format: text

text:
# Mean Field Calculation of Thermal Properties of Simple Nucleon Matter on a Lattice

## I Introduction

Nuclear excitations are complicated dynamical phenomena, depending on the detailed structure of the individual nucleus, and must be examined based on the specific structure of the nucleus, such as whether it is closed or open shell. As the excitation energy gets higher, however, the excitations depend less on specific nuclear structure and start to exhibit more common features among (heavy) nuclei. These features are expected to be reasonably well represented by excitations of nuclear matter. Furthermore, dynamics of supernovae and neutron stars, which have been of much astronomical interest, are expected to be better understood through the study of excitations of neutron matter dean ; arpns . The gross features of the thermal properties of nucleon matter have been examined by means of statistical models stat  and lattice gas models latgas . More realistic descriptions of the thermal properties have been provided through applications of various approaches in the nuclear many-body theories TTT ; panda-t ; kyoto-t ; other-t ; jkuo ; gapref . Applications of the traditional nuclear many-body theories regarding nucleon matter at zero temperature have been extensive arpns ; panda ; kyoto ; fantoni , and provide the most reliable information on the properties of nucleon matter at low temperatures.

Previously, one of us (R. S.) has collaborated on a Monte Carlo calculation of nuclear matter on a lattice MKSV , which provides a new framework for studying the thermal properties of nucleon matter. Though the computational space was small and the nucleon-nucleon interaction was simple, the calculation has proven to be of much promise, demonstrating the occurrence of a phase transition around 15 MeV using the parameters adjusted to reproduce the saturation properties. A similar calculation of nucleon matter of a single species has also been initiated in the same work, using the same form of the Hamiltonian. Though a phase transition appeared to take place at a few MeV, the evidence for it was not quite solid owing to statistical fluctuations, which are enhanced at low temperature (a sign problem). The phase transitions may correspond to those expected through paired nucleons (Cooper pairs) in nucleon matter dean ; arpns .

In order to gain a better understanding of the possible, latter phase transition, we apply in this work the analytic means of a mean field approach to the problem in the same lattice formulation. From this work, we do not expect to be able to draw precise quantitative conclusions, but rather we will try to learn the nature of the phase transition at a semi-quantitative level. For this purpose, we take the thermodynamical (infinite volume) limit for numerical results, so as to obtain a clear signal of the phase transition. The mean field results in this work will also serve as a reference for the more extensive Monte Carlo calculation that we are currently carrying out.

Our Hamiltonian for single-species nucleon matter turns out to be an extended, attractive Hubbard model, which has been studied as a simple model of high-temperature superconductivity RMP:MRR . The mean field calculation shows that the low-temperature, low-density state is a superfluid state and undergoes a continuous (second-order) phase transition to the normal state as the temperature and/or density increases. As the Hamiltonian is not yet fully realistic and the values of the interaction parameters are uncertain, our results are not quite comparable to those for neutron matter, except perhaps at a very low density. But we demonstrate that the approach is promising for the study of low-density neutron matter.

Furthermore, we find that a density-wave state coexists with the superfluid state, suggesting that the state of neutron matter may be more complicated than the simple description of a superfluid state as it is often characterized.

The outline of this work is as follows: After the introduction in Sec. I, the Hamiltonian and its discretized form in the coordinate space are presented and are identified as an extended Hubbard model in Sec. II. The Hartree-Fock Bogoliubov approximation is applied and the Hamiltonian is diagonalized in Sec. III. Thermodynamical properties numerically calculated are shown in Sec. IV, an attempt to apply our calculation to the problem of low-density neutron matter is discussed in Sec. V, and the possibility of a density wave phase is examined in Sec. VI. Discussions and conclusion are presented in Sec. VII.

## II Simplified Hamiltonian and Extended Hubbard Model

The Hamiltonian consists of the kinetic and potential terms $`\widehat{K}`$ and $`\widehat{V}`$, respectively:

$`\widehat{H}`$ $`=`$ $`\widehat{K}+\widehat{V}`$ (1)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m_N}}{\displaystyle \underset{\sigma \tau }{}}{\displaystyle d𝐫\widehat{\psi }_{\sigma \tau }^{}(𝐫)^2\widehat{\psi }_{\sigma \tau }(𝐫)}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{\sigma \sigma ^{}}{\tau \tau ^{}}}{}}{\displaystyle d𝐫d𝐫^{}\widehat{\psi }_{\sigma \tau }^{}(𝐫)\widehat{\psi }_{\sigma ^{}\tau ^{}}^{}(𝐫^{})V(𝐫𝐫^{})\widehat{\psi }_{\sigma ^{}\tau ^{}}(𝐫^{})\widehat{\psi }_{\sigma \tau }(𝐫)},`$
where $`m_N`$ is the nucleon mass, and $`\sigma =\pm 1/2`$ and $`\tau =\pm 1/2`$ are the spin ($``$ or $``$) and the isospin ($`p`$ or $`n`$), respectively. $`\widehat{\psi }_{\sigma \tau }^{}(𝐫)`$ and $`\widehat{\psi }_{\sigma \tau }(𝐫)`$ are the creation and annihilation operators of the nucleon, with the spin $`\sigma `$ and isospin $`\tau `$ at the position $`𝐫`$. As in the previous Monte Carlo lattice calculation MKSV , we include only the central and spin-exchange interactions, $`V_c`$ and $`V_\sigma `$, respectively:

$$V(𝐫𝐫^{})=V_c(𝐫𝐫^{})+V_\sigma (𝐫𝐫^{})𝝈𝝈^{}.$$
(2)
$`V_c`$ and $`V_\sigma `$ are taken to consist of on-site and next-neighbor interactions,

$$\begin{array}{c}V_c\left(𝐫𝐫^{}\right)=V_c^{(0)}\delta \left(𝐫𝐫^{}\right)+V_c^{(2)}\left[_𝐫^2\delta \left(𝐫𝐫^{}\right)\right]\hfill \\ V_\sigma \left(𝐫𝐫^{}\right)=V_\sigma ^{(0)}\delta \left(𝐫𝐫^{}\right)+V_\sigma ^{(2)}\left[_𝐫^2\delta \left(𝐫𝐫^{}\right)\right],\hfill \end{array}$$
(3)
where $`V^{(2)}`$ terms can be written explicitly exhibiting their hermiticity.

As physics of the lattice description is more apparent in the coordinate space, we consider the discretized coordinate with an inter-nucleon spacing $`a`$ in the cubic lattice with the torus boundary conditions. We thus focus our interest on the physics of the spatial separation greater than $`a`$ in each direction, or of the momentum component roughly between $`\pi /a`$ and $`\pi /a`$, by eliminating (or integrating out) the physics of the shorter distance. The discretization corresponds to

$`𝐫`$ $``$ $`a𝐧_i,`$
$`{\displaystyle d𝐫}`$ $``$ $`a^3{\displaystyle \underset{i}{}},`$
$`{\displaystyle \underset{\sigma \tau }{}}{\displaystyle 𝑑𝐫\widehat{\psi }_{\sigma \tau }^{}(𝐫)\widehat{\psi }_{\sigma \tau }(𝐫)}`$ $``$ $`{\displaystyle \underset{i\sigma \tau }{}}\widehat{c}_{i\sigma \tau }^{}\widehat{c}_{i\sigma \tau },`$
where $`i`$ denotes a lattice site specified by $`𝐧_i`$ with its component ranging $`[aN/2,aN/2]`$. Here, $`N`$ is the number of sites in each spatial direction. Note that the creation and annihilation operators, $`\widehat{c}_{i\sigma \tau }^{}`$ and $`\widehat{c}_{i\sigma \tau }`$, have no dimension as defined. We also apply the identity

$$\underset{i}{}\sigma _{\alpha \beta }^{(i)}\sigma _{\gamma \delta }^{(i)}=2\delta _{\alpha \delta }\delta _{\beta \gamma }\delta _{\alpha \beta }\delta _{\gamma \delta },$$

where $`i`$ denotes the spatial components $`(x,y,z)`$, and the Greek indices denote the components of the Pauli spin matrix, 1 or 2.

In this work, we study the simplified case of nucleon matter consisting of a single nucleon species, such as neutron matter. The Hamiltonian Eq. (1) is then expressed in a spatially discretized form:

$`\widehat{H}`$ $`=`$ $`t{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+6t{\displaystyle \underset{i\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }+\stackrel{~}{U}{\displaystyle \underset{i}{}}\widehat{c}_i^{}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i`$ (4)
$`+V_1{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }^{}\widehat{c}_{j\sigma }\widehat{c}_{i\sigma }+V_2{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }^{}\widehat{c}_{j\sigma }\widehat{c}_{i\sigma }+V_3{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }^{}\widehat{c}_{j\sigma }\widehat{c}_{i\sigma },`$
where $`i,j`$ denotes the pairs of next-neighbor sites, and $`t`$ is the hopping (kinetic energy) parameter defined as

$$t=\frac{\mathrm{}^2}{2m_Na^2}.$$
(5)
Here, the potential parameters $`\stackrel{~}{U}`$ and $`V`$’s are expressed in terms of linear combinations of $`V^{(0)}`$’s and $`V^{(2)}`$’s. The Hamiltonian Eq. (4) is now in the form of an extended Hubbard model, which is the Hubbard model with the on-site spin-pairing interaction of the $`\stackrel{~}{U}`$ term, modified by the next-neighbor interaction of the $`V`$ terms. As the $`\stackrel{~}{U}`$ value will be taken to be negative, our model is an extended, attractive Hubbard model. The repulsive Hubbard model has been well studied in condensed matter physics as a model of strongly correlated electron systems Fradkin , but the attractive model is generally less studied. In recent years, however, the extended, attractive Hubbard model has drawn much attention as the model describing the essential features of high-temperature superconductivity RMP:MRR . Note that the extended, attractive (negative-$`U`$) Hubbard model used in condensed matter physics, however, is usually of the simpler form shown below and has no $`6t`$ term as a part of the kinetic energy tara .

When the spin-dependent next-neighbor interaction is taken to be small and negligible,

$$V_\sigma ^{(2)}=0,$$
(6)
the Hamiltonian is simplified:

$$\widehat{H}=t\underset{i,j\sigma }{}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+6t\underset{i\sigma }{}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }+U\underset{i}{}\widehat{c}_i^{}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i+V\underset{i,j\sigma \sigma ^{}}{}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }\widehat{c}_{j\sigma ^{}}^{}\widehat{c}_{j\sigma ^{}},$$
(7)
where

$`U`$ $`=`$ $`{\displaystyle \frac{1}{a^3}}\left(V_c^{(0)}6{\displaystyle \frac{V_c^{(2)}}{a^2}}3V_\sigma ^{(0)}\right),`$
$`V`$ $`=`$ $`{\displaystyle \frac{1}{2a^5}}V_c^{(2)}.`$ (8)
Apart from the lattice spacing $`a`$, the Hamiltonian of Eq. (7) now describes dynamics with two parameters. The Hamiltonian Eq. (4) \[and thus also Eq. (7)\] possesses an underlying particle-hole symmetry, which affects thermodynamical properties as discussed in Sec. VI. The symmetry is not an explicit property in our original Hamiltonian, Eqs. (1)–(3). We elaborate on the symmetry in Appendix A.

## III Hartree-Fock Bogoliubov Approximation and Gap Equations

We now apply the mean field method in the Hartree-Fock Bogoliubov approximation HFB . Here, we expect effects of the $`U`$ term to dominate the thermal properties of the single species matter as in the standard BCS description schrieffer ; Leg ; Noz , but we also wish to treat their single-particle aspects in the Hartree-Fock approximation on the same footing. Since the method is well known, we limit the description of the formalism to the key steps that are specifically relevant to our calculation.

The nature of the mean field approximation is apparent in the spatial representation. Through the application of the Wick theorem, the decoupling scheme for the $`U`$ term is

$$\widehat{c}_i^{}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i\mathrm{\Delta }_i^{}\widehat{c}_i\widehat{c}_i\mathrm{\Delta }_i\widehat{c}_i^{}\widehat{c}_i^{}\left|\mathrm{\Delta }_i\right|^2+n_i\widehat{c}_i^{}\widehat{c}_i+n_i\widehat{c}_i^{}\widehat{c}_in_in_i,$$

where

$`\mathrm{\Delta }_i\widehat{c}_i\widehat{c}_i=\widehat{c}_i\widehat{c}_i,`$
$`n_i\widehat{c}_i^{}\widehat{c}_i\mathrm{and}n_i\widehat{c}_i^{}\widehat{c}_i.`$
Here, $`\mathrm{}`$ denotes the expectation value in the BCS-like ground state, which is to be determined in a self-consistent way. The order parameters $`\mathrm{\Delta }_i`$’s are related to the local density of the condensate of nucleon pairs, while $`n_i`$’s are the average number of nucleons. Note that the inclusion of $`n_i`$ as a variational parameter distinguishes the present treatment from the standard BCS HFB . For the $`V`$ terms, we have

$$\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }\widehat{c}_{j\sigma ^{}}^{}\widehat{c}_{j\sigma ^{}}n_{i\sigma }\widehat{c}_{j\sigma ^{}}^{}\widehat{c}_{j\sigma ^{}}+n_{j\sigma ^{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }n_{i\sigma }n_{j\sigma ^{}}.$$

In the following, $`\mathrm{\Delta }_i`$ and $`n_i`$’s will be assumed to be independent of the site $`i`$, or global, and will be formally treated as the variational parameters:

$`\mathrm{\Delta }2U\mathrm{\Delta }_i=2U\mathrm{\Delta }_i^{},`$
$`n2n_i=2n_i.`$
Our $`\mathrm{\Delta }`$ is defined to be the gap energy itself, with the dimension of energy in the unit of $`2U`$, and carries an extra factor of 2 in comparison to the often-used $`\mathrm{\Delta }`$. We also note that we ignore the $`V`$-term contribution to the hopping term, as they merely change somewhat the strength of the hopping term and of the constant part of the energy, without affecting the physics of the phase transition. The mean field approximation then yields in a cubic lattice with six neighboring sites

$`\widehat{H}`$ $`=`$ $`t{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+\left(6t+{\displaystyle \frac{n}{2}}\overline{U}\right){\displaystyle \underset{i\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$ (9)
$`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{i}{}}\left(\widehat{c}_i^{}\widehat{c}_i^{}+\widehat{c}_i\widehat{c}_i\right){\displaystyle \frac{N^3}{4}}\left({\displaystyle \frac{\mathrm{\Delta }^2}{U}}+n^2\overline{U}\right),`$
where $`\overline{U}U+24V`$.

For convenience, we will carry out the mean field calculations in the momentum space. The momentum representation is introduced by discretizing the momentum as $`𝐩2\pi 𝐤/(Na)`$, with each component of $`𝐤`$ being an integer, ranging $`[N/2,N/2]`$. Note that we now have $`𝐫𝐩2\pi 𝐧_i𝐤/N`$. The coordinate and momentum representations of the operator $`\widehat{c}`$ are related through the Fourier transformations

$$\begin{array}{c}\widehat{c}_{j\sigma }=\frac{1}{\sqrt{N^3}}\underset{𝐤}{}e^{2i\pi 𝐤𝐧_j/N}\widehat{c}_{𝐤\sigma },\hfill \\ \widehat{c}_{𝐤\sigma }=\frac{1}{\sqrt{N^3}}\underset{j}{}e^{2i\pi 𝐤𝐧_j/N}\widehat{c}_{j\sigma },\hfill \end{array}$$

and similarly for the nucleon creation operator, $`\widehat{c}^{}`$. Note that the discretized orthonormality relation is

$$\underset{j}{}\mathrm{exp}\left(2i\pi 𝐤𝐧_j/N\right)=N^3\delta _{𝐤,\mathrm{𝟎}}.$$
(10)
With the chemical potential $`\mu `$, $`\widehat{H}\mu N^3\widehat{n}`$ is in the momentum space

$`\widehat{H}\mu N^3\widehat{n}`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\left(\widehat{c}_𝐤^{}\widehat{c}_𝐤\right)\left(\begin{array}{cc}ϵ_𝐤\overline{\mu }& \frac{1}{2}\mathrm{\Delta }\\ \frac{1}{2}\mathrm{\Delta }& (ϵ_𝐤\overline{\mu })\end{array}\right)\left(\begin{array}{c}\widehat{c}_𝐤\\ \widehat{c}_𝐤^{}\end{array}\right)`$ (16)
$`{\displaystyle \frac{N^3}{4}}\left({\displaystyle \frac{\mathrm{\Delta }^2}{U}}+n^2\overline{U}\right),`$
where

$$\overline{\mu }=\mu \left(6t+n\overline{U}/2\right),$$
(17)
and

$$ϵ_𝐤=t\underset{𝐞}{}\mathrm{exp}\left(2i\pi 𝐤𝐞/N\right)=2t\underset{j=x,y,z}{}\mathrm{cos}\left(2\pi k_j/N\right)$$
(18)
is a part of the kinetic energy of a quasi-particle expressed in terms of the unit vector $`𝐞`$ showing a next-neighbor site. Note that $`ϵ_𝐤=ϵ_𝐤`$ and $`_𝐤ϵ_𝐤=0`$. We see that in the Hamiltonian, $`V`$ appears only as $`\overline{U}=U+24V`$, and merely shifts the chemical potential and the total energy: It does not actively participate in the generation of the phase transition.

As is well known, the spin density $`s`$,

$$s=\widehat{s}\frac{1}{2N^3}\underset{i}{}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i^{}\widehat{c}_i=\frac{1}{2N^3}\underset{𝐤}{}\widehat{c}_𝐤^{}\widehat{c}_𝐤\widehat{c}_𝐤^{}\widehat{c}_𝐤,$$
(19)
is conserved in the mean field approach, while $`n`$,

$$n=\widehat{n}\frac{1}{N^3}\underset{i\sigma }{}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }=\frac{1}{N^3}\underset{𝐤\sigma }{}\widehat{c}_{𝐤\sigma }^{}\widehat{c}_{𝐤\sigma }$$
(20)
(and thus the total nucleon number), is not, as the spin number operator $`\widehat{s}`$ commutes with $`\widehat{H}`$ of Eq. (9), but the number density operator $`\widehat{n}`$ does not. We remedy this problem by the standard method of introducing a Lagrange multiplier corresponding to the chemical potential $`\mu `$, by subtracting a term $`\mu N^3\widehat{n}=\mu _{i\sigma }\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$ from the $`\widehat{H}`$. $`\mu `$ will be adjusted to achieve the desired value of the conjugate parameter, $`n`$. The formalism is thus essentially the canonical ensemble method.

We now diagonalize $`\widehat{H}\mu N^3\widehat{n}`$ by the use of the Bogoliubov-Valatin transformation,

$$\begin{array}{c}\widehat{\beta }_{𝐤+}^{}=u_𝐤\widehat{c}_𝐤^{}v_𝐤\widehat{c}_𝐤,\hfill \\ \widehat{\beta }_𝐤^{}=u_𝐤\widehat{c}_𝐤^{}+v_𝐤\widehat{c}_𝐤,\hfill \end{array}$$
(21)
where $`u_𝐤`$ and $`v_𝐤`$ are taken to be real and are given by

$$u_𝐤^2=\frac{1}{2}\left(1+\frac{ϵ_𝐤\overline{\mu }}{E_𝐤}\right)\mathrm{and}v_𝐤^2=\frac{1}{2}\left(1\frac{ϵ_𝐤\overline{\mu }}{E_𝐤}\right),$$

satisfying $`u_𝐤^2+v_𝐤^2=1`$. $`\widehat{H}\mu N^3\widehat{n}`$ is now expressed as that of a system of free quasi-particles:

$$\widehat{H}\mu N^3\widehat{n}=\underset{𝐤\lambda =\pm }{}E_𝐤\widehat{\beta }_{𝐤\lambda }^{}\widehat{\beta }_{𝐤\lambda }+N^3(E_{\mathrm{GS}}\mu n),$$
(22)
with the energy of a quasi-particle

$$E_𝐤=\sqrt{(ϵ_𝐤\overline{\mu })^2+\mathrm{\Delta }^2/4}$$
(23)
and the ground-state energy of the system

$$E_{\mathrm{GS}}=\frac{1}{4}\left(\frac{\mathrm{\Delta }^2}{U}+n^2\overline{U}\right)\frac{1}{N^3}\underset{𝐤}{}E_𝐤\overline{\mu }+\mu n.$$
(24)
Equation (21) shows that $`\widehat{\beta }_{𝐤\lambda }`$’s obey anti-commutation relations and that the quasi-particles are fermions. Furthermore, Eq. (22) implies that they form a system of free fermions. As a consequence of the thermal average, the internal energy is then given by

$$E\widehat{H}=\underset{𝐤\lambda }{}E_𝐤n_{𝐤\lambda }+N^3E_{\mathrm{GS}},$$
(25)
with the energy per lattice site being $`E/N^3`$. Here, $`n_{𝐤\lambda }`$ is the momentum distribution of the quasi-particles, $`n_{𝐤\lambda }\widehat{\beta }_{𝐤\lambda }^{}\widehat{\beta }_{𝐤\lambda }`$, and is determined by the requirement that the free energy $`F`$ introduced below is minimized by a variation of $`n_{𝐤\lambda }`$, $`\delta F/\delta n_{𝐤\lambda }=0`$ for $`\lambda =\pm `$. We obtain

$$n_𝐤n_{𝐤+}=n_𝐤=\left[\mathrm{exp}\left(E_𝐤/T\right)+1\right]^1,$$
(26)
which have the limiting values, $`n_𝐤0`$ and $`1/2`$ as $`T0^+`$ and $`\mathrm{}`$, respectively. Note that throughout this work, we denote the temperature $`T`$ in the unit of the Boltzmann constant.

By combining with the entropy $`S`$, the (Helmholtz) free energy is expressed as

$`F(T,a;\mathrm{\Delta },n)`$ $``$ $`ETS`$ (27)
$`=`$ $`{\displaystyle \underset{𝐤\lambda }{}}E_𝐤n_{𝐤\lambda }+N^3E_{\mathrm{GS}}`$
$`+T{\displaystyle \underset{𝐤\lambda }{}}\left[n_{𝐤\lambda }\mathrm{ln}n_{𝐤\lambda }+\left(1n_{𝐤\lambda }\right)\mathrm{ln}\left(1n_{𝐤\lambda }\right)\right],`$
depending on $`a`$ through the $`𝐤`$ sum because the spatial volume $`(aN)^3`$ depends on $`a`$ with $`N`$ fixed. $`F`$ is a function of $`\mathrm{\Delta }`$ and $`n`$. $`\mathrm{\Delta }`$ and $`\mu `$ are determined so as to minimize $`F`$ for variations of $`\mathrm{\Delta }`$ and $`n`$, while $`T`$ and $`a`$ are fixed. The conditions

$$\frac{1}{N^3}\frac{F}{n}=\mu ,\frac{F}{\mathrm{\Delta }}=0$$

provide the gap equations

$$\begin{array}{c}n1=\frac{1}{N^3}\underset{𝐤}{}\frac{ϵ_𝐤\overline{\mu }}{E_𝐤}\left(2n_𝐤1\right),\hfill \\ \mathrm{\Delta }\left(1\frac{U}{2N^3}\underset{𝐤}{}\frac{2n_𝐤1}{E_𝐤}\right)=0\hfill \end{array}$$
(28)
from which $`\mu `$ and $`\mathrm{\Delta }`$ are determined.

## IV Thermodynamical properties

We now apply the formalism so far described, to compute various thermodynamical quantities. For clarifying our presentation, we place some expressions of the thermodynamical variables in Appendix B. All numerical results are calculated in the thermodynamical limit $`N\mathrm{}`$. In the limit, the summation over the discretized momentum space of each component of $`𝐤`$ ranging in $`[N/2,N/2]`$ is replaced by the integral over the first Brillioun zone with each component of the momentum $`𝐩`$ ranging $`[\pi /a,\pi /a]`$:

$$\frac{1}{N^3}\underset{𝐤}{}\left(\frac{a}{2\pi }\right)^3_{\pi /a}^{\pi /a}_{\pi /a}^{\pi /a}_{\pi /a}^{\pi /a}d^3𝐩.$$

### IV.1 Potential parameters

We apply the parameter values used in the previous Monte Carlo lattice calculation for nuclear matter MKSV :

$`\begin{array}{c}V_c^{(0)}=181.5\mathrm{MeV}\mathrm{fm}^3,\hfill \\ V_c^{(2)}=37.8\mathrm{MeV}\mathrm{fm}^5,\hfill \\ V_\sigma ^{(0)}=31.25\mathrm{MeV}\mathrm{fm}^3,\hfill \\ V_\sigma ^{(2)}=0,\hfill \end{array}`$ (33)
with the lattice spacing $`a=1.842\mathrm{fm}`$. These parameter values give

$$\begin{array}{c}t=6.11\mathrm{MeV},\hfill \\ U=24.74\mathrm{MeV},\hfill \\ V=0.89\mathrm{MeV}.\hfill \end{array}$$

The parameter values of Eq. (33) were chosen in Ref. MKSV  so as to reproduce the saturation density and energy of nuclear matter on a finite $`4\times 4\times 4`$ lattice for the same Hamiltonian as ours, Eqs. (1)–(3). Our Hamiltonian has no explicit $`\tau `$-dependent term, and the parameter values effectively include the strong neutron-proton interactions for nuclear matter. The use of the parameters is thus not quite adequate as a realistic description of the nucleon matter of single species, such as the neutron matter. Furthermore, finite lattice volume effects make the thermodynamical limit ($`N\mathrm{}`$) calculation different from the finite volume calculation. For comparison purposes with the previous and future Monte Carlo calculations, however, we use the above parameter values except when the $`U`$ dependence of $`\mathrm{\Delta }`$ is examined. Any conclusion that we could draw from the numerical results in this section is then qualitative.

### IV.2 Gap parameter $`\mathrm{\Delta }`$

Equations (28) determine $`\mathrm{\Delta }`$ and $`\mu `$. Figure 2 illustrates $`\mathrm{\Delta }`$ as a function of the temperature $`T`$ for $`n=`$ 0.5 (one-quarter filling), 1.0 (one-half filling), and 1.5 (three-quarter filling). In the figure, we see that $`\mathrm{\Delta }`$ vanishes at $`T=T_c`$,

$$T_c=\{\begin{array}{cc}\hfill 0.66t\mathrm{or}\mathrm{\hspace{0.33em}\hspace{0.33em}4.0}\mathrm{MeV}& \mathrm{for}n=1.0,\hfill \\ \hfill 0.55t\mathrm{or}\mathrm{\hspace{0.33em}\hspace{0.33em}3.3}\mathrm{MeV}& \mathrm{for}n=0.5\mathrm{and}\mathrm{\hspace{0.33em}\hspace{0.33em}1.5}.\hfill \end{array}$$

The temperature dependence of $`\mathrm{\Delta }`$ is the same for $`n=0.5`$ and 1.5. This is a consequence of the symmetry with respect to $`n=1`$ and is discussed further in Sec. IV D and Appendix A. Figure 2 also shows

$$\mathrm{\Delta }(T=0)(22.5)t1215\mathrm{MeV}.$$

The explicit values of $`T_c`$ and $`\mathrm{\Delta }(T=0)`$ depend sensitively on the parameter values as discussed below, but $`T_c`$ and $`\mathrm{\Delta }(T=0)`$ satisfy

$$\mathrm{\Delta }(T=0)3.6T_c,$$

which is the well-known relation at the weak-coupling limit except for the latter to have a slightly smaller coefficient 3.54 schrieffer . In comparison to the weak limit, our calculation thus somewhat underestimates $`T_c`$ in relation to $`\mathrm{\Delta }(T=0)`$. Note that different mean field calculations have been reported to yield the coefficient smaller than the weak-limit value gapref  and also even much larger jkuo  than ours.

Near $`T_c`$, the gap equations Eqs. (28) yield

$$\mathrm{\Delta }\{\begin{array}{cc}(TT_c)^\beta & \mathrm{for}T<T_c,\hfill \\ 0& \mathrm{for}T>T_c,\hfill \end{array}$$

with $`\beta 0.45`$. Note that the well-known mean field value of the critical exponent $`\beta `$ in the simple BCS theory is 1/2 renorm .

The physics of the phase transition depends on the strength of the potential parameter $`U`$. Figure 2 illustrates how sensitively the value of $`\mathrm{\Delta }`$ depends on $`U`$. We see that $`\mathrm{\Delta }0`$ as $`U0`$. As is well known, $`\mathrm{\Delta }`$ does not vanish for a finite $`U`$. In fact, Eq. (28) yields the well-known dependence of $`\mathrm{\Delta }`$ on $`|U|`$ for $`U0`$:

$$\mathrm{\Delta }Ae^{B/|U|},$$

where $`A`$ and $`B`$ are constant and independent of $`U`$.

### IV.3 Second-order (continuous) phase transition

The temperature dependence of $`\mathrm{\Delta }`$ in Fig. 2 is a well-known dependence of the order parameter for a second-order phase transition. The variation of $`\mathrm{\Delta }`$ as shown in Fig. 2 implies that the phase transition takes place from a superfluid state generated by spin pairing to the normal state, as the temperature increases. The features of the second-order phase transition are clearly seen in the temperature dependence of the thermodynamic quantities expressed in terms of the temperature derivatives of the free energy $`F`$ in a successive order. We consider the internal energy $`E`$, the entropy $`S`$, and the heat capacity $`C_v`$,

$$\begin{array}{c}E=T^2\left[\frac{(F/T)}{T}\right]_{a,n},\hfill \\ S=\left[\frac{F}{T}\right]_{a,n},\hfill \\ C_v=T\left[\frac{S}{T}\right]_{a,n}=T\left[\frac{^2F}{T^2}\right]_{a,n}.\hfill \end{array}$$
(34)
The temperature dependence of the quantities is calculated using Eqs. (25), (27), (53), and (54), and is illustrated in Figs. 4 and 4. $`E`$ and $`S`$ are continuous at the critical temperature $`T_c`$ as seen in Fig. 4, while $`C_v`$ has a jump at $`T_c`$ as in Fig. 4. These behaviors demonstrate the generic features of the second-order phase transition. The amount of the discontinuity in $`C_v`$ at $`T_c`$, $`\mathrm{\Delta }C_v`$, is relative to $`C_v`$ of the normal phase,

$$\mathrm{\Delta }C_v/C_v(\mathrm{normal})\{\begin{array}{cc}\hfill 1.74& \mathrm{for}n=1.0,\hfill \\ \hfill 1.43& \mathrm{for}n=0.5\mathrm{and}\mathrm{\hspace{0.33em}\hspace{0.33em}1.5},\hfill \end{array}$$

while the BCS mean field value is $`12/7\zeta (3)1.43`$, independent of $`n`$ schrieffer . ($`\zeta `$ is the zeta function.)

The thermal quantities involving volume derivatives form a set of quantities similar to the temperature derivatives. We consider the pressure $`P`$ and the isothermal compressibility $`\kappa _T`$. Here, following the common practice in nuclear physics, we examine the incompressibility $`K9/(\kappa _T\rho )`$, defined in terms of the density $`\rho =𝒱/(aN)^3=n/a^3`$ with the spatial volume $`𝒱=(aN)^3`$. A volume derivative is then a derivative with respect to the lattice spacing $`a`$. A derivative with respect to $`a`$ requires, however, the knowledge of $`a`$ dependence of $`U`$ and $`V`$, that is, their renormalization flow when $`a`$ is varied. In this work, for simplicity, we assume that their $`a`$ dependence is small, at least around the value of $`a`$ we use. $`P`$ and $`K`$ are written as

$`P`$ $`=`$ $`\left[{\displaystyle \frac{F}{𝒱}}\right]_{T,a},`$ (35)
$`K`$ $`=`$ $`9{\displaystyle \frac{𝒱}{\rho }}\left[{\displaystyle \frac{P}{𝒱}}\right]_{T,a}=9{\displaystyle \frac{𝒱}{\rho }}\left[{\displaystyle \frac{^2F}{𝒱^2}}\right]_{T,a}.`$ (36)
$`P`$ is calculated using Eq. (56) and $`K`$ is obtained numerically from the temperature dependence of $`P`$. The temperature dependence of $`P`$ and $`K`$ confirms the second-order phase transition, as seen in Fig. 4.

There are many other quantities that describe thermodynamic properties of the system described by our Hamiltonian, but they are either related to the quantities already shown, or their features depend strongly on the explicit form of the Hamiltonian. We thus do not show them in this exploratory work. For example, the temperature dependence of double occupancy per site,

$$D\frac{1}{N^3}\underset{i}{}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i^{}\widehat{c}_i,$$
(37)
provides the amount of the spin pairing that participates in the phase transition. But $`D`$, as well as the kinetic energy per site, is related to $`\mathrm{\Delta }`$ and $`n`$ in the mean field approximation armen  and provides no new information, as shown in Appendices A and B.

### IV.4 Particle-hole symmetry and phase diagram

For illustrative purposes, however, we show the density dependence of $`\mu `$, $`D`$, $`E`$, and $`KE`$ at $`T=0`$, in Fig. 6. As noted above, Fig. 2 suggests a symmetric dependence of $`\mathrm{\Delta }`$ on the density $`n`$, which is more clearly illustrated in Fig. 6. The symmetry is generated as a consequence of the particle-hole symmetry, as seen from the fact that Eq. (28) is invariant under the particle-hole conjugation though the symmetry is implicit in our Hamiltonian, Eqs. (4) and (7). There is a group of Hamiltonians, in which the particle-hole symmetry is implicit, yet yielding (in the mean field results) the energy spectrum of the system with the explicit symmetry. The Hamiltonians do yield different behaviors of thermodynamic variables. Figure 6 is an example. We elaborate on the issue of this symmetry in Appendix A.

Some thermodynamic variables are made to exhibit this symmetry explicitly by modifying the Hamiltonian to possess the explicit symmetry. Whether it is explicit or not, however, the thermodynamic properties obtained from the Hamiltonian are affected by the symmetry. The symmetry is a consequence of our computational method using the lattice configuration, and it is an artifact. In order to extract physically realistic results, we should therefore stay away from the region of the symmetry and should confine ourselves to a small value of $`n`$ by appropriately adjusting the value of the lattice spacing $`a`$, so as to simulate the desired density of the nucleon matter. We discuss this point again in the following section, where we attempt to apply our calculation to a case of low-density neutron matter.

Combining the variations of the thermodynamic quantities, some of which have been presented so far, we obtain the phase diagrams of the present system described in the mean field theory. Figure 8 shows the phase diagram in the region of small densities where we expect the above-mentioned symmetry to be generating less distortion.

## V Low-density Neutron matter

In the previous section, we have used the parameter values most appropriate as a description of nuclear matter and have examined the nature of the single species nucleon matter described by our model. In this section, we discuss whether our model could be made a realistic description of neutron matter.

First, we have the question of whether our lattice would meet the basic momentum requirement imposed by the lattice spacing. A lattice description can be made realistic when the lattice spacing is less than the momentum scale of the system. By taking the Fermi momentum of the neutron matter as an estimate of the momentum scale, we have

$$\pi /a>p_F.$$
(38)
Our lattice spacing, $`a=1.842`$ fm, yields the density of $`\rho =1/a^30.160\mathrm{fm}^3`$ for $`n=1.0`$ (the lattice space being half full). The Fermi momentum corresponding to this density is

$$p_F1.68\mathrm{fm}^1$$

using $`\rho =p_{F}^{}{}_{}{}^{3}/3\pi ^2`$. The value of $`p_F`$ is practically the same as

$$\pi /a1.71\mathrm{fm}^1,$$

and thus the lattice with the above lattice spacing is applicable to a density much smaller than $`\rho 0.160\mathrm{fm}^3`$. Note that the preceding discussion yields that the condition

$$\pi /3>n$$

meets Eq. (38) independently of $`a`$.

Second, there is the question of whether our Hamiltonian is appropriate for a realistic description of low-density neutron matter. The nucleon-nucleon $`{}_{}{}^{1}S_{0}^{}`$ phase shift is much greater than the nucleon-nucleon phase shifts of other states below the laboratory energy $`E_{lab}`$ 100 MeV. The nucleon momentum in the center-of-mass coordinate system $`p_{cm}`$ corresponding to this $`E_{lab}`$ is about 1.2 fm <sup>-1</sup> through $`E_{lab}=4(p_{cm}^2/2m_N)`$, and is smaller than the above $`\pi /ap_{cutoff}`$ 1.7 $`\mathrm{fm}^1`$. We thus infer that our Hamiltonian form of the $`S`$-wave should be reasonable for neutron matter of a density less than 0.17 $`\mathrm{fm}^3`$, which corresponds to the Fermi momentum of 1.7 $`\mathrm{fm}^1`$.

As to the parameter values in the Hamiltonian, it would be best to determine them for our lattice size from experimental $`{}_{}{}^{1}S_{0}^{}`$ phase shifts by applying the method of effective field theory MS ; lstv . Instead, as an exploratory study, we simply adjust the $`U`$-parameter value so as to see whether our approach could come close to other mean field calculations of low-density neutron matter in the literature. Figure 8 illustrates that we could obtain a somewhat reasonable density dependence of $`\mathrm{\Delta }`$ by increasing the magnitude of $`U`$. We leave a more serious determination of the parameters for our future work.

## VI Density wave phase

Our analysis has so far been strictly based on the mean field approximation applied to a spin-pairing phase at the same site. In the approximation, the $`V`$ term merely shifts the effective chemical potential and is inactive in generating the phase transition. As the $`V`$ term represents the pairing of the spin densities at the adjacent sites, such a role may be a reasonable one in this phase transistion. Would the $`V`$ term ever play an active role in generating a different phase transition? In this section we briefly examine this possibility.

The most likely phase in which the $`V`$ term would play the major role would be a density wave phase generated by a coupling of the densities of the opposite spins at the adjacent sites. We examine how the $`V`$ term could generate such a phase transition, again in the mean field approximation, and see whether the phase transition would occur with our parameter values.

We follow the Hartree-Fock Bogoliubov approximation with the same decoupling scheme as before. For simplicity, however, we ignore the superfluid phase. The density wave phase is introduced by making a replacement,

$`\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$ $``$ $`n_{i\sigma }+(\delta /2)\mathrm{cos}(2\pi 𝐧_i𝐪/N)`$ (39)
$``$ $`n/2+(\delta /2)\mathrm{cos}(2\pi 𝐧_i𝐪/N).`$
Here, the inhomogeneous order parameter for the density wave depends on the amplitude $`\delta `$ and the wave number vector $`𝐪`$. In the following, we examine the wave modes in which the adjacent sites are the maximum and minimum of the amplitude in a cubic lattice. That is, at least one component of $`2\pi \widehat{q}/N`$ is $`\pm \pi `$. The number of the non-zero components, or the dimension of the density wave $`d`$, provides a convenient parameter $`\eta `$,

$$\eta =32d.$$
(40)
As in the previous spin-pairing case, $`n`$ and $`\delta `$ are treated as independent parameters, and $`\mu `$, the parameter conjugate to $`n`$, is adjusted to the desired value of $`n`$.

The steps to the gap equations are similar to the case of the superfluid phase discussed in Sec. III, and are shown in Appendix C. By minimizing $`F\mu N^3n`$, we obtain

$$\begin{array}{c}n=\frac{1}{N^3}\underset{𝐤\lambda }{}n_{𝐤\lambda },\hfill \\ \delta \left(1\frac{U+8\eta V}{N^3}\underset{𝐤}{}\frac{n_{𝐤+}n_𝐤}{\xi _𝐤}\right)=0\hfill \end{array}$$
(41)
by varying $`n`$ and $`\delta `$, respectively. Equation (41) determines $`\delta `$ and $`\mu `$. Furthermore, as $`n_{𝐤+}<n_𝐤`$ and

$$U+8\eta V<0$$
(42)
for our parameter values, Eq. (41) shows that the density waves of all dimensions are expected to occur owing to the strong, attractive $`U`$.

## VII Discussion and Conclusion

The form of the nucleon-nucleon interaction, Eqs. (1)–(3), shows that the nucleon-nucleon interaction used in this work is of the S states. As is well known, the realistic nucleon-nucleon interaction is highly state-dependent. The relevant nucleon energy of interest to us here is a few hundred MeV, corresponding to the Fermi energy region of the nuclear matter density. In this energy region, the attractive neutron-neutron interaction is known to be dominated by the $`{}_{}{}^{3}P_{2}^{}`$ interaction driven by the spin-orbit force, coupled with the $`{}_{}{}^{3}F_{2}^{}`$ interaction associated with the tensor force NN ; TT . Our Hamiltonian accommodates none of these features of the interaction. Our objective as noted in Sec. I is to understand the essential physics associated with the thermal properties of nucleon matter, but our finding in this work is limited in this sense and is perhaps most applicable to low-density neutron matter.

There is a serious question of how good the mean field calculation is in our case. A mean field approximation ignores most features of particle correlations. As the particle correlations are the vital ingredient of critical phenomena, the mean field approximation is generally believed to be only of qualitative use, and some times not even qualitative, as fluctuations could alter the nature of the phase transition. The situation, however, depends on the nature of the problem renorm , as the prominent success of the Ginzburg-Landau/BCS theory shows superc , especially at temperatures not too close to the critical one. Our problem is in three dimensions, close to the usual upper critical dimension of four under the Ginsburg criterion renorm . We are hoping that our calculation, being similar to the BCS theory, is not far off, but this remains to be seen. This issue is under further investigation by incorporating a renormalization approach, as has recently been done at zero temperature schwenk .

Our Hamiltonian has a form similar to that of the Skyrme interaction skyrme  (though ours is a truncated form.) The Skyrme interaction is one of the effective interactions that are phenomenologically introduced to achieve quantitative agreement with experiments, usually by the use of a mean field approximation such as the Hartree-Fock calculations. Though it is still not quite realistic, our Hamiltonian has a justification in this sense. A lattice calculation such as the previous Monte Carlo lattice calculation MKSV , however, accounts for all the complexity of the many-body interaction with no approximation other than the numerical, in the lattice framework. The nucleon-nucleon interaction used in it should not be then an effective interaction like the Skyrme interaction, but an interaction in free space. The parameters $`U`$ and $`V`$ are expected to be determined from scattering data through the use of effective field theory by extending Lüscher’s formula luescher ; wash  for the large $`{}_{}{}^{1}S_{0}^{}`$ scattering length. This issue is presently under investigation lstv . Note that the values of $`U`$ and $`V`$ will depend on the lattice spacing $`a`$.

In conclusion, our mean field calculation on a lattice with a simple nucleon-nucleon interaction suggests a second-order phase transition taking place at a low temperature in single species nucleon matter described by a simple Hamiltonian. Thermodynamic variables show their dependence on the temperature and density variations as expected under the phase transition in the mean field approach. The transition changes the phase of the matter from a superfluid state due to nucleon pairing, to a normal state as the temperature increases and the density decrease. This dependence is in qualitative agreement with the findings that have been reported in the literature dean .

ACKNOWLEDGMENTS

This work is supported by the U.S. Department of Energy under Grant No. DE-FG03-87ER40347 at CSUN, and the U.S. National Science Foundation under grants PHY0071856 and PHY0244899 at Caltech. T. A. is supported under the 21st Century COE Program at the Tokyo Institute of Technology, “Nanometer-Scale Quantum Physics,” by the Ministry of Education, Culture, Sports, Science and Technology in Japan. R. S. acknowledges Satoshi Okamoto for patiently explaining various aspects of condensed matter BCS theory. For his warm hospitality, R. S. thanks Masahiko Iwasaki of RIKEN, where a part of this work was done. T. A. and R. S. acknowledge for its hospitality the Institute for Nuclear Theory at the University of Washington, where the last stage of this work was carried out.

## Appendix A The particle-hole symmetry

The Hamiltonian of Eq. (7) does not explicitly exhibit particle-hole symmetry, but it can be modified to do so by adding a single-particle Hamiltonian,

$$\mathrm{\Delta }\widehat{H}=N^3\overline{U}/4\left(6t+\overline{U}/2\right)\underset{i\sigma }{}\widehat{n}_{i\sigma },$$
(43)
where $`\widehat{n}_{i\sigma }\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$. The new Hamiltonian,

$`\widehat{H}^{}`$ $``$ $`\widehat{H}+\mathrm{\Delta }\widehat{H}`$ (44)
$`=`$ $`t{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+U{\displaystyle \underset{i}{}}\left({\displaystyle \frac{1}{2}}\widehat{n}_i\right)\left({\displaystyle \frac{1}{2}}\widehat{n}_i\right)`$
$`+V{\displaystyle \underset{i,j\sigma \sigma ^{}}{}}\left({\displaystyle \frac{1}{2}}\widehat{n}_{i\sigma }\right)\left({\displaystyle \frac{1}{2}}\widehat{n}_{j\sigma ^{}}\right),`$
is symmetric under the particle-hole conjugation (or with respect to half filling, $`n=1`$):

$`\widehat{c}_{i\sigma }^{}`$ $``$ $`\widehat{c}_{i\sigma }`$
$`\widehat{c}_{i\sigma }`$ $``$ $`\widehat{c}_{i\sigma }^{},`$
or $`n_{i\sigma }1n_{i\sigma }`$, which is $`n2n`$ with $`n_{i\sigma }n_{i\sigma }`$. Under the conjugation, the Hamiltonian becomes that of holes, with $`tt`$.

When we repeat the BCS formulation in Sec. III, after the Bogoliubov transformation, we obtain the diagonalized form of the new Hamiltonian, Eq. (44), as

$$\widehat{H}^{}=\underset{𝐤\lambda }{}E_𝐤^{}\widehat{n}_{𝐤\lambda }^{}+N^3E_{\mathrm{GS}}^{}$$
(45)
with the quasi-particle energy

$$E_𝐤^{}=\sqrt{(ϵ_𝐤\overline{\mu }^{})^2+\mathrm{\Delta }^2/4}$$
(46)
and

$$n_𝐤^{}n_{𝐤+}^{}=n_𝐤^{}=\left[\mathrm{exp}\left(E_𝐤^{}/T\right)+1\right]^1.$$
(47)
Here, the ground-state energy is

$$E_{\mathrm{GS}}^{}=\frac{1}{4}\left(\frac{\mathrm{\Delta }^2}{U}+\frac{n^21}{4}\overline{U}\right)\frac{1}{N^3}\underset{𝐤}{}E_𝐤^{}\overline{\mu }^{}+\mu n,$$
(48)
with

$$\overline{\mu }^{}=\mu (n1)\overline{U}/2.$$
(49)
Equations (17) and (49) show that the $`\mu `$ value is shifted, and Eqs. (24) and (48) tell us that the expression of the ground-state energy is altered. The kinetic energy is also changed:

$$KE^{}\frac{t}{N^3}\underset{i,j,\sigma }{}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }=\frac{1}{N^3}\underset{𝐤}{}\frac{\left(ϵ_𝐤\overline{\mu }^{}\right)ϵ_𝐤}{E_𝐤^{}},$$
(50)
in comparison to Eq. (59). We also see that $`E_{\mathrm{GS}}^{}`$ and $`KE^{}`$ for the new Hamiltonian explicitly exhibit the particle-hole symmetry about half filling ($`n=1`$).

The double occupancy per site for the new Hamiltonian $`D^{}`$,

$$D^{}\frac{1}{N^3}\underset{i}{}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i^{}\widehat{c}_i=\frac{\mathrm{\Delta }^2}{4U^2}+\frac{n^2}{4},$$
(51)
differs from that for the original Hamiltonian, Eq. (58). Note that neither $`D`$ nor $`D^{}`$ is symmetric about $`n=1`$.

The gap equations for the new Hamiltonian become

$$\begin{array}{c}n1=\frac{1}{N^3}\underset{𝐤}{}\frac{ϵ_𝐤\overline{\mu }^{}}{E_𝐤^{}}\left(2n_𝐤^{}1\right),\hfill \\ \mathrm{\Delta }\left(1\frac{U}{2N^3}\underset{𝐤}{}\frac{2n_𝐤^{}1}{E_𝐤^{}}\right)=0.\hfill \end{array}$$
(52)
These gap equations are the same as Eqs. (28) except $`\overline{\mu }`$ is replaced by $`\overline{\mu }^{}`$. As we solve the gap equations for $`\mathrm{\Delta }`$ and $`\mu `$ (or $`\overline{\mu }`$) for a fixed $`n`$, the gap equations for the two Hamiltonians yield the same set of $`\mathrm{\Delta }`$ and $`\overline{\mu }`$, thus the same $`T_c`$. Both Hamiltonians thus provide the same excitation energy spectrum.

## Appendix B Thermodynamic variables

We list here the expressions of the thermodynamic variables, which are used in Sec. IV. The entropy is given in terms of $`n_{𝐤\lambda }`$ in Eq. (26) by

$$S=\underset{𝐤\lambda }{}\left[n_{𝐤\lambda }\mathrm{ln}n_{𝐤\lambda }+\left(1n_{𝐤\lambda }\right)\mathrm{ln}\left(1n_{𝐤\lambda }\right)\right].$$
(53)
By the use of Eqs. (23), (24), and (27), the heat capacity $`C_v`$ is expressed as

$`C_v`$ $`=`$ $`T\left[{\displaystyle \frac{^2F}{T^2}}\right]_{a,n}`$ (54)
$`=`$ $`{\displaystyle \frac{1}{N^3}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Delta }_T}{4E_𝐤}}\left(2n_𝐤1\right){\displaystyle \frac{2}{N^3}}{\displaystyle \underset{𝐤}{}}n_𝐤^2\left({\displaystyle \frac{\mathrm{\Delta }\mathrm{\Delta }_T}{4}}{\displaystyle \frac{E_𝐤^2}{T^2}}\right)e^{E_𝐤/T}{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Delta }_T}{2U}}.`$
Here, $`\mathrm{\Delta }_T`$ is

$$\mathrm{\Delta }_T\frac{\mathrm{\Delta }}{T},$$
(55)
and is numerically calculated from the solution of the gap equations, Eq. (28).

The pressure $`P`$ is written in terms of the space volume $`𝒱=(aN)^3`$ as

$`P`$ $`=`$ $`\left[{\displaystyle \frac{F}{𝒱}}\right]_{T,n}`$ (56)
$`=`$ $`{\displaystyle \frac{1}{N^3}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{E_𝐤}{a^3}}(2n_𝐤1)+{\displaystyle \frac{n}{2}}\left({\displaystyle \frac{n}{2}}1\right)\left({\displaystyle \frac{U}{a^3}}+24{\displaystyle \frac{V}{a^3}}\right),`$
where

$`{\displaystyle \frac{E_𝐤}{a^3}}`$ $`=`$ $`{\displaystyle \frac{ϵ_𝐤\overline{\mu }}{E_𝐤}}\left[{\displaystyle \frac{ϵ_𝐤}{a^3}}+{\displaystyle \frac{n}{2}}\left({\displaystyle \frac{U}{a^3}}+24{\displaystyle \frac{V}{a^3}}\right)\right],`$
$`{\displaystyle \frac{ϵ_𝐤}{a^3}}`$ $`=`$ $`{\displaystyle \frac{2t}{3a^2}}{\displaystyle \underset{i=x,y,z}{}}k_i\mathrm{cos}\left({\displaystyle \frac{2\pi }{N}}k_i\right),`$
$`{\displaystyle \frac{U}{a^3}}`$ $`=`$ $`{\displaystyle \frac{1}{a^6}}\left[\left(V_c^{(0)}5{\displaystyle \frac{V_c^{(2)}}{a^2}}\right)3\left(V_\sigma ^{(0)}5{\displaystyle \frac{V_\sigma ^{(2)}}{a^2}}\right)\right],`$
$`{\displaystyle \frac{V}{a^3}}`$ $`=`$ $`{\displaystyle \frac{5}{6a^8}}V_c^{(2)}.`$ (57)
The double occupancy per site $`D`$ is expressed as

$`D`$ $``$ $`{\displaystyle \frac{1}{N^3}}{\displaystyle \underset{i}{}}\widehat{c}_i^{}\widehat{c}_i\widehat{c}_i^{}\widehat{c}_i={\displaystyle \frac{1}{N^3}}{\displaystyle \frac{E}{U}}`$ (58)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }^2}{4U^2}}+{\displaystyle \frac{n^2}{4}}{\displaystyle \frac{n}{N^3}}{\displaystyle \underset{𝐤}{}}\left(ϵ_𝐤\overline{\mu }\right){\displaystyle \frac{n_𝐤^2}{T}}e^{E_𝐤/T}.`$
The kinetic energy per site $`KE`$ is written as

$`KE`$ $``$ $`{\displaystyle \frac{t}{N^3}}{\displaystyle \underset{i,j,\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+6{\displaystyle \frac{t}{N^3}}{\displaystyle \underset{i\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }={\displaystyle \frac{t}{N^3}}{\displaystyle \frac{E}{t}}`$ (59)
$`=`$ $`{\displaystyle \frac{1}{N^3}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\left(ϵ_𝐤\overline{\mu }\right)\left(ϵ_𝐤+6t\right)}{E_𝐤}}\left[2n_𝐤\left(1{\displaystyle \frac{n_𝐤}{T}}E_𝐤e^{E_𝐤/T}\right)1\right]+6t.`$
## Appendix C Density wave phase

The Hamiltonian in the coordinate space is reduced to

$`\widehat{H}`$ $`=`$ $`t{\displaystyle \underset{i,j\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+(6t+n\overline{U}/2){\displaystyle \underset{i\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$ (60)
$`+(U/2+4\eta V)\delta {\displaystyle \underset{i\sigma }{}}\mathrm{cos}(2\pi 𝐧_i𝐪/N)\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$
$`\left[(U/4)(n^2+\delta ^2/2)+V(6n^2+\eta \delta ^2)\right]N^3,`$
where the parameter $`\eta `$ depends on the dimension of the density wave $`d`$, $`\eta =32d`$. Here, we have used

$`{\displaystyle \underset{i,j\sigma }{}}(n_j+n_j)\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$ $``$ $`{\displaystyle \underset{i,j\sigma }{}}(n_i+n_i)\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }+\delta {\displaystyle \underset{i𝐞\sigma }{}}\mathrm{cos}(2\pi (𝐧_i+𝐞)𝐪/N)\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$ (61)
$``$ $`6n{\displaystyle \underset{i\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }+2\eta \delta {\displaystyle \underset{i\sigma }{}}\mathrm{cos}(2\pi 𝐧_i𝐪/N)\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }`$
and

$`{\displaystyle \underset{i,j}{}}(n_i+n_i)(n_j+n_j)`$ $``$ $`{\displaystyle \underset{i,j}{}}(n_i+n_i)(n_j+n_j)+2\delta {\displaystyle \underset{ij\sigma }{}}\mathrm{cos}(2\pi 𝐧_i𝐪/N)(n_j+n_j)`$ (62)
$`+\delta ^2{\displaystyle \underset{i𝐞}{}}\mathrm{cos}(2\pi 𝐧_i𝐪/N)\mathrm{cos}(2\pi (𝐧_i+𝐞)𝐪/N)`$
$``$ $`(6n^2+\eta \delta ^2)N^3,`$
where the nucleon densities $`n`$’s are defined in the same way as in Sec. III. In the last step of Eq. (62), we have used the identities

$$\begin{array}{c}\underset{i}{}\mathrm{cos}(2\pi 𝐧_i𝐪/N)=0\hfill \\ \underset{i}{}\mathrm{cos}^2(2\pi 𝐧_i𝐪/N)=\frac{N^3}{2}\hfill \end{array}$$
(63)
by applying the discretized orthonormality relation Eq. (10).

Using

$$\underset{i\sigma }{}\mathrm{cos}\left(\frac{2\pi }{N}𝐧_i𝐪\right)\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }=\frac{1}{2}\underset{𝐤\sigma }{}\left(\widehat{c}_{𝐤𝐪/2\sigma }^{}{}_{}{}^{}\widehat{c}_{𝐤+𝐪/2\sigma }+\widehat{c}_{𝐤+𝐪/2\sigma }^{}{}_{}{}^{}\widehat{c}_{𝐤𝐪/2\sigma }\right),$$
(64)
we obtain the momentum space representation

$`\widehat{H}\mu N^3\widehat{n}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐤\sigma }{}}\left(ϵ_+\overline{\mu }\right)\widehat{c}_{+\sigma }^{}\widehat{c}_{+\sigma }+{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐤\sigma }{}}\left(ϵ_{}\overline{\mu }\right)\widehat{c}_\sigma ^{}\widehat{c}_\sigma `$ (70)
$`+{\displaystyle \frac{1}{8}}\overline{\delta }{\displaystyle \underset{𝐤\sigma }{}}\left(\widehat{c}_\sigma ^{}\widehat{c}_{+\sigma }+\widehat{c}_{+\sigma }^{}\widehat{c}_\sigma \right){\displaystyle \frac{1}{2}}\delta \overline{\delta }N^3+E_0N^3`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐤\sigma }{}}\left(\widehat{c}_{+\sigma }^{}\widehat{c}_\sigma ^{}\right)\left(\begin{array}{cc}ϵ_+\overline{\mu }& \frac{1}{4}\overline{\delta }\\ \frac{1}{4}\overline{\delta }& ϵ_{}\overline{\mu }\end{array}\right)\left(\begin{array}{c}\widehat{c}_{+\sigma }\\ \widehat{c}_\sigma \end{array}\right)`$
$`{\displaystyle \frac{1}{16}}\delta \overline{\delta }N^3+E_0N^3,`$
where the subscripts $`+`$ and $``$ denote $`𝐤+𝐪/2`$ and $`𝐤𝐪/2`$, respectively; that is, $`\widehat{c}_{\pm \sigma }=\widehat{c}_{𝐤\pm 𝐪/2\sigma }`$. $`ϵ_\pm `$ is defined as

$$ϵ_\pm =t\underset{𝐞}{}\mathrm{exp}\left(i2\pi (𝐤\pm 𝐪/2)𝐞/N\right)=2t\underset{j=x,y,z}{}\mathrm{cos}\left(2\pi (k_j\pm q_j/2)/N\right).$$

We have also defined

$$\overline{\delta }=\delta \left(2U+16\eta V\right).$$
(71)
$`\widehat{H}\mu N^3\widehat{n}`$ is diagonalized through the transformation

$$\begin{array}{c}\widehat{\alpha }_{𝐤+}=u_𝐤\widehat{c}_{+\sigma }+v_𝐤\widehat{c}_\sigma ,\hfill \\ \widehat{\alpha }_𝐤=u_𝐤\widehat{c}_\sigma v_𝐤\widehat{c}_{+\sigma },\hfill \end{array}$$
(72)
where $`u_𝐤`$ and $`v_𝐤`$ are taken to be real and are given by

$$u_𝐤^2=\frac{1}{2}\left(1+\frac{ϵ_+ϵ_{}}{\xi _𝐤}\right)\mathrm{and}v_𝐤^2=\frac{1}{2}\left(1\frac{ϵ_+ϵ_{}}{\xi _𝐤}\right),$$

satisfying $`u_𝐤^2+v_𝐤^2=1`$. Here,

$$\xi _𝐤=\sqrt{(ϵ_+ϵ_{})^2+\overline{\delta }^2/4}.$$
(73)
$`\widehat{H}\mu N^3\widehat{n}`$ is now expressed as that of a system of free quasi-particles:

$$\widehat{H}\mu N^3\widehat{n}=\frac{1}{2}\underset{𝐤\lambda }{}E_{𝐤\lambda }\widehat{\alpha }_{𝐤\lambda }^{}\widehat{\alpha }_{𝐤\lambda }+N^3(E_{\mathrm{GS}}\mu n),$$
(74)
where $`\lambda =\pm `$ and

$$E_{𝐤\pm }=\frac{1}{2}\left[(ϵ_++ϵ_{}2\overline{\mu })\pm \xi _𝐤\right]$$
(75)
is the energy of the quasi-particles, and

$$E_{\mathrm{GS}}=\delta \overline{\delta }/16n^2\overline{U}/4+\mu n$$
(76)
is the ground-state energy of the system.


text:
# 1 Problem

## 1 Problem

Use the Weizsäcker-Williams method to deduce the radiated power, and its angular distribution, emitted by an electron of charge $`e`$ that undergoes oscillatory motion, $`z=z_0\mathrm{cos}\omega t`$.

## 2 Solution

This problem is a continuation of an earlier note on the Weizsäcker-Williams method .

The basic idea of the method is that if a rapidly moving electron is sufficiently perturbed, it will emit $`\alpha =e^2/\mathrm{}c`$ photons per unit bandwidth during per formation time $`t_0`$, where the latter is the time it takes for the radiation to pull one wavelength ahead of the charge. Thus, the photon number spectrum of the radiation is

$$\frac{dn}{d\omega dt_0}\alpha .$$
(1)
Underlying the method is the fact that the electromagnetic fields of a relativistic charge are largely transverse to the direction of motion of the charge, and hence much like a free electromagnetic wave (pulse).

Here, we seek to apply this method to nonrelativistic radiation by an oscillating charge. Because the charge has a maximum velocity (relative to the speed of light $`c`$)

$$\beta =\frac{v_{\mathrm{max}}}{c}=\frac{z_0\omega }{c}$$
(2)
much less than one, we anticipate that the radiation will not have the full strength attainable in the relativistic limit. Rather, the amplitude of the radiation will be scaled by a factor $`\beta `$, and so the radiated power will be scaled down by a factor $`\beta ^2`$.111A related argument was given in Sec. IIID of  for the case of radiation by an electron passing through a weak undulator.

Because the motion of the charge is periodic, the radiation spectrum is a line at angular frequency $`\omega `$. We estimate that the rate of photon emission at frequency $`\omega `$ is roughly the same as that in one unit of bandwidth in the continuum version of the model,

Because the motion of the electron is bounded (with amplitude small compared to the wavelength $`\lambda =2\pi c/\omega `$), the formation length is simply a wavelength, and the formation time $`t_0`$ is simply one period, $`2\pi /\omega `$.

Applying, the basic principle (1) of the Weizsäcker-Williams method, with the modifications described above for nonrelativistic, period motion of the charge, the number of photons radiated per unit time is

$$\frac{dN}{dt}\alpha \beta ^2\frac{\omega }{2\pi }=\frac{e^2}{\mathrm{}c}\frac{z_0^2\omega ^2}{c^2}\frac{\omega }{2\pi }.$$
(3)
The radiated power is $`\mathrm{}\omega `$ times this:

$$P=\mathrm{}\omega \frac{dN}{dt}\frac{(ez_0)^2\omega ^4}{2\pi c^3}=\frac{p_0^2\omega ^4}{2\pi c^3}.$$
(4)
where $`p_0=ez_0`$ is the (maximum) dipole moment of the electron. The result (4) is about 1/2 that of the usual expression for the time-averaged power from an oscillation dipole $`p=p_0\mathrm{cos}\omega t`$,

$$P=\frac{p_0^2\omega ^4}{3c^3}.$$
(5)
The factor $`2\pi /3`$ difference between forms (4) and (5) indicates that the effective bandwidth of the continuum that is “squeezed” into the line spectrum is only about 1/2 unit, rather than one as assumed above.

Some care is required to determine the angular distribution of the radiation from the Weizsäcker-Williams perspective. One might be tempted to argue that the radiation will be emitted preferentially in the direction of velocity v of the charge, because wavelike fields that surround the charge (its cloud of “virtual photons”) have $`𝐄\times 𝐁`$ largely parallel to v. And indeed, in the relativistic limit the radiation is emitted primarily in the v direction. However, the present problem is in the nonrelativistic limit, in which the average velocity of the charge is zero, so we must look elsewhere for a description of the angular distribution.

A useful argument is based on the radiation reaction. Power is emitted by the oscillating charge, so we expect a back reaction on the charge. This reaction must come from the interaction of the charge with the radiation fields. The power of the reaction force is

$$P_{\mathrm{react}}=𝐅_{\mathrm{react}}𝐯=e𝐄_{\mathrm{react}}𝐯=e𝐄_{\mathrm{rad}}𝐯.$$
(6)
To provide the required radiation reaction, the electric field $`𝐄_{\mathrm{rad}}`$ of the radiation must have a large component along the direction of v, i.e., along the $`z`$ axis in the present problem. This indicates that the electric field of radiation emitted at angle $`\theta `$ measured with respect to the $`z`$ axis varies as $`\mathrm{sin}\theta `$ (rather than $`\mathrm{cos}\theta `$ as would be too naively inferred from the previous paragraph).

The radiated power goes as the square of the electric field, so the angular distribution corresponding to eq. (4) is

$$\frac{dP}{d\mathrm{\Omega }}=\frac{3p_0^2\omega ^4}{8\pi ^2c^3}\mathrm{sin}^2\theta .$$
(7)
Or, applying the factor $`2\pi /3`$, we recover the usual Hertzian expression,

$$\frac{dP}{d\mathrm{\Omega }}=\frac{p_0^2\omega ^4}{4\pi c^3}\mathrm{sin}^2\theta .$$
(8)
K.T.M. wishes to thank Dan Handelsman, N2DT, and David Jeffries for the e-discussions that encouraged a different way of thinking about radiation from antennas.


text:
# Star-Forming Galaxies in the Sloan Digital Sky Survey – The View from Pittsburgh

## 1 Indroduction

The CMU-Pitt Value Added Catalog (VAC) database (http://astrophysics.phys. cmu.edu/vac/) allows users to conduct research on SDSS spectroscopic galaxy data. The current public VAC contains 133,440 SDSS DR1 galaxies, of which 124,411 are unique. This enables the pursuit of projects which require very large galaxy samples. We here describe two such projects; first, the search for rare, X-ray emitting normal SFGs, and second, an investigation of the L-Z relation for galaxies of different sizes – compact to large – and its evolution with redshift.

## 2 Sample Definition

The sample was drawn from all objects spectrally classified as galaxies. We required that the H$`\beta `$ rest equivalent width of an object be larger than 5Å . The threshold of 5Å was set because the SDSS pipeline does not account for underlying Balmer absorption when calculating equivalent widths, and was motivated by the work of Kobulnicky et al. (1999), who recommend a statistical correction of 3$`\pm `$2Å for underlying H$`\beta `$ absorption to integrated spectra of late-type galaxies. We adopted a conservative criterion for separating SFGs and AGN. We required that the galaxies fell into the SFG regime in all three diagnostic diagrams of Kewley et al. (2001). The SDSS SFGs/AGN separation is also described in Miller et al. (2003). In addition, we checked that all galaxies classified as SFGs conform to the same relation of H$`\beta `$ versus \[OIII\]$`\lambda \lambda `$5007 line width. The SDSS wavelength range and the SFGs/AGN separation together, limit the redshift range of the sample from about 0.02 to 0.3. We use a concordance cosmology with $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`H_0`$=70 km s$`^1`$Mpc-1.

## 3 SDSS/ROSAT SFGs

We searched for positional matches to our objects with RASS Position Sensitive Proportional Counter (PSPC) data, with catalogs of pointed PSPC observations, and with catalogs of pointed High Resolution Imager (HRI) observations, maintained at MPE. The search radius was 60” for PSPC, and 20” for HRI catalogs. This resulted in a list of 146 unique SDSS DR1 SFGs which have positional coincidence with a ROSAT source (or source list). The ROSAT error circle may include a cluster of galaxies, an AGN or an X-ray bright star, as well as the SFG, but we are interested only in secure X-ray/SFG detections. AGN are X-ray bright and fairly numerous, and a match of our list against NED and SDSS QSO catalogs (Richards et al. 2002) immediately revealed that 1/4 could be QSOs, while another 1/10 could be clusters. We used a soft-band luminosity of 1041 erg s-1 to separate SFGs from luminous, AGN and clusters in the remaining 2/3 of sample. In order to transform the ROSAT count rates into an unabsorbed flux, we used the Galactic HI column densities of Dickey & Lockman (1990), and assumed an intrisic powerlaw spectrum with a photon index of 2.0. We converted redshifts to luminosity distances and used these to convert fluxes to luminosities. There are twelve objects with soft X-ray luminosities $``$ 1041 erg s-1. In addition, the data suggest two possible matches based on small fX-over-fR ratios. While these have LX $`>`$ 1041 erg s-1, they were added to our candidate sample, bringing it up to 14 objects. We investigated the properties of all 14 cases in great detail (but did not pursue the nature of the rest of the sources).

Four SFGs were previously known to be soft X-ray emitters, namely the edge-on SB galaxy IC 2233, the Blue Compact Dwarf (BCD) I Zw 18, the SA galaxy NGC 4030, and the edge-on SA galaxy NGC 5907. There are four newly identified X-ray emitting SFGs: SDSS J014143.17+134033.16, a face-on Spiral galaxy also known as NPM1G +13.0066, SDSS J081313.18+455940.7, the BCD Mrk 86, SDSS J130039.24+023002.81, the face-on Virgo-cluster SB galaxy NGC 4900, and SDSS J165817.7+644220.11, an inclined Spiral galaxy. The fluxes in the 0.1–2.4 keV band range from about 0.5 to 1.4$`\times `$10-14 erg s-1 cm-2; fairly deep because of pointed ROSAT observations. The logarithmic luminosities range from about 37.8 to 41.5. Matches coinciding with the known clusters Abell 119, Abell 318 and RDCS J0910+5422 were not attributed to SFGs. We suggest that there may be galaxy groups or clusters at J0233+0041, with a redshift of about 0.022, and at J1717+6420, with a redshift of about 0.035.

Table 1 lists the basic data for the X-ray emitting SFGs identified in this work. Table 2 gives fluxes and luminosities which can be calibrated to a SFR. Our objects follow the soft-X-ray—FIR and soft-X-ray—radio relations derived by Ranalli et al. (2003). They propose a soft-X-ray—SFR calibration based on 23 galaxies, which is nicely confirmed by our added sample of 8 SDSS/ROSAT/SFGs, as shown by Fig. 1. The last column of Table 2 gives the X-ray SFRs of our sample based on equation (14) of Ranalli et al. (2003).

$``$ Our results underscore the utility of soft X-ray fluxes to derive SFRs.

## 4 SDSS L-Z Relation

Local galaxies are observed to conform to a luminosity-metallicity relation: luminous galaxies exhibit higher metallicites than do dwarfs (Melbourne & Salzer 2002, and references therein). The L-Z relation is usually shown as total absolute B magnitude versus O/H ratio. The L-Z relation is increasingly being used to study galaxy evolution at intermediate (e.g. Lilly et al. 2003) and high (e.g. Pettini et al. 2001) redshifts, when galaxies were younger and smaller.

In order to study the dependence of the L-Z relation on galaxy size as well as on redshift, we divided our SDSS DR1 SFG sample into three groups depending on $`r`$-band size: compact galaxies with 0$`<`$Re\[kpc\]$``$2, small galaxies with 2$`<`$Re\[kpc\]$``$5; this is the largest group and includes Milky-Way–sized objects, and large galaxies having 5$`<`$Re\[kpc\]$``$100. Photometry in $`g`$ and $`r`$ was transformed to the B band following Smith et al. (2002); galactic foreground reddening was removed and K-corrections were applied (Blanton et al. 2003). In many cases, the SDSS pipeline magnitudes referred to HII regions within galaxies, rather than to entire galaxies. This phenomenon is also known as “shredding”. Of course, the compact galaxy sample was most severely affected by shredding, and we used in our analysis only “parents” which had not been deblended into “children”. This has eliminated a number of real compact galaxies, e.g., the two HII regions of I Zw 18 are listed as separated galaxies in the DR1. We used the N2 parameter (Denicoló et al. 2002) in order to derive O/H ratios. The N2 parameter has two advantages compared with the R23 index (Pagel et al. 1979): it is not affected by reddening, and its calibration to O abundance is single-valued. The SDSS fibers have a size of 3” and the O/H ratios usually refer to the inner regions of the galaxies. The L-Z relation for about 5,000 SDSS SFGs with high signal-to-noise data is shown in Fig. 2.

The results of the KISS survey for $``$500 local SFGs (Melbourne & Salzer 2002) indicate a flattening of the L-Z relation for dwarfs. Fig. 2 suggests that small SDSS SFGs have a flatter L-Z relations than large ones. We also see a large amount of scatter among compact, dwarf galaxies. We interpret this to indicate that these galaxies are affected more by feedback than luminous/large galaxies.

The top left quadrant of each panel in Fig. 2 shows where SDSS SFGs (of all sizes) in four redshift slices lie relative to the KISS L-Z relation (Melbourne & Salzer 2002). We see a strong evolution, in the sense that the SDSS SFGs systematically shift across the KISS L-Z relation as redshift increases. Specifically, the data suggest that at a given metallicity (say, 12+log(O/H) = 8.7), galaxies are less luminous now than they were just 2-3 Gyrs ago; or, that at a given L (say, MB = -19.5), galaxy metallicity declines with increasing redshift.

The data in each of the redshift slices, do not extend over the luminosity range of the entire sample: in the lowest-z bin the volume is too small to include luminous galaxies, while in the highest-z bin, the limited depth of the survey precludes the detection of dwarfs. Therefore, the main concern is that luminosity introduces some kind of bias. In order to start to address this question, we used a magnitude limited sample with $`r`$ $``$17.77, as in the SDSS main galaxy sample (Strauss et al. 2002). This has the effect of removing galaxies at low O/H and decreasing the scatter among compact galaxies, but also, of cutting galaxies from the two highest-z bins. In principle, the L-Z releation could vary in both, slope and zero point, but because our galaxies cover only a small fraction of the entire L range in each redshift slice, we assumed a constant slope is retained and tested only for variation in zero point with redshift. The difference in zero point of the L-Z relation when comparing the lowest with the now highest, 0.1$`<`$z$``$0.2 bin, is significant, and amounts to over 1.5 magnitudes.

$``$ SDSS data indicate a very rapid evolution in the L-Z relation for late-type galaxies since z$`<0`$.3.

Acknowledgments Drs. M. Bernardi and D. Vanden Berk helped us understand SDSS parameters. RS-L thanks the MPE for the invitation to Munich in fall of 2002. AMH gratefully acknowledges support provided by NASA through Hubble Fellowship grant HST-HF-01140.01-A awarded by the Space Telescope Science Institute (contract NAS 5-26555). Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. We made extensive use of NED.


text:
# The world-wide air transportation network: Anomalous centrality, community structure, and cities’ global roles

## The data

Many measures—including total number of passengers, total number of flights, or total amount of cargo—quantifying the importance of the world airports are compiled and publicized . We study here the OAG MAX database111For more information on the OAG MAX database, visit http://oagdata.com/solutions/max.aspx., which comprises flight schedule data of more than 800 of the world’s airlines for the period November 1, 2000 to October 31, 2001. This database is compiled by OAG, a division of Reed Business Information Group, and includes all scheduled flights and scheduled charter flights, both for big aircrafts—air carriers—and small aircrafts—air taxis.

We focus our analysis on a network of cities, not of airports—for example, the Newark, JFK and La Guardia airports are all assigned to New York City. We further restrict our analysis to passenger flights operating in the time period November 1, 2000 to November 7, 2000. Even though this data is more than four years old, the resulting world-wide airport network is virtually indistinguishable from the network one would obtain if using data collected today. The reason is that air traffic patterns are strongly correlated with: (i) socio-economic factors, such as population density and economic development; and (ii) geo-political factors such as the distribution of the continents over the surface of the Earth and the locations of borders between states . Clearly, the time scales associated to changes in these factors are much longer than the lag in the data we analyze here.

During the period considered, there are 531,574 unique non-stop passenger flights, or flight segments, operating between 3883 distinct cities. We identify 27,051 distinct city pairs having non-stop connections. The fact that the database is highly redundant—that is, that most connections between pairs of cities are represented by more than one flight—adds reliability to our analysis. Specifically, the fact that unscheduled flights are not considered does not mean, in general, that the corresponding link between a certain pair of cities is missing in the network, since analogous scheduled flights may still operate between them. Similarly, even if some airlines have canceled their flights between a pair of cities since November 2000, it is highly unlikely that all of them have.

We create the corresponding adjacency matrix for this network, which turns out to be almost symmetrical. The very minor asymmetry stems from the fact that a small number of flights follow a “circular” pattern, i.e, a flight might go from A to B to C, and back to A. To simplify the analysis, we symmetrize the adjacency matrix.

Further, we build regional networks for different geographic regions (Table 1). Specifically, we generate twenty-one regional networks at different aggregation levels. At the highest-aggregation level, we generate six networks; one each for Africa, Asia and Middle East, Europe, Latin America, North America, and Oceania. For each of these regions, except for North America and Oceania, we generate between two and five sub-networks. For instance, the Asia and Middle East network is further subdivided into South Asia, Central Asia, Southeast Asia, Northeast Asia, and Middle East.

## Large-scale structure of the air transportation network

A ubiquitous characteristic of complex networks is the so-called “small-world” property . In a small-world network, pairs of nodes are connected by short paths as one expects for a random graph . Crucially, nodes in small-world networks also have a high-degree of cliquishness, as one finds in low-dimensional lattices but not in random graphs.

In the air transportation network, the average shortest path length $`d`$ is the average minimum number of flights that one needs to take to get from any city to any other city in the world. We find that for the $`719`$ cities in the Asia and Middle East network $`d=3.5`$, and that the average shortest path length between the $`3663`$ cities in the giant component of the world-wide network is only about one step greater, $`d=4.4`$. Actually, most pairs of cities (56%) are connected by four steps or less. More generally, we find that $`d`$ grows logarithmically with the number $`S`$ of cities in the network, $`d\mathrm{log}S`$. This behavior is consistent with both random graphs and small-world networks, but not with low-dimensional networks, for which $`d`$ grows more rapidly with $`S`$.

Still, some pairs of cities are considerably further away from each other than the average. The farthest cities in the network are Mount Pleasant, in the Falkland Islands, and Wasu, in Papua New Guinea: To get from one city to the other, one needs to take fifteen different flights. From Mount Pleasant, one can fly to Punta Arenas, in Chile, and from there to some hubs in Latin America. At the other end of the path, from Wasu one needs to fly to Port Moresby, which requires a unique sequence of eight flights. In the center of the path, between Punta Arenas and Port Moresby, six different flights are needed. In contrast with what happens the ends of the path, in the central region of the path there are hundreds of different flight combinations, all of them connecting Punta Arenas and Port Moresby in six steps.

The clustering coefficient $`C`$, which quantifies the local cliquishness of a network, is defined as the probability that two cities that are directly connected to a third city are also directly connected to each other. We find that $`C`$ is typically larger for the air transportation network than for a random graph and that it increases with size. These results are consistent with the expectations for a small-world network but not those for a random graph. For the world-wide network, we find $`C=0.62`$ while its randomization yields $`C=0.049`$. Therefore, we conclude that the air transportation network is, as expected, a small-world network .

Another fundamental aspect in which real-world networks often deviate from the random graphs typically considered in mathematical analysis  is the degree distribution, that is, the distribution of the number of links of the nodes . In binomial random graphs, all nodes have similar degrees, while many real world networks have some nodes that are significantly more connected than others. Specifically, many complex networks, termed scale-free, have degree distributions that decay as a power law. A plausible mechanism for such a phenomenon is preferential attachment , that is, the tendency to connect preferentially to nodes with already high degrees.

To gain greater insight into the structure and evolution of the air transportation network, we calculate the degree distribution of the cities. The degree of a city is the number of other cities to which it is connected by a non-stop flight. In Fig. 1a, we show the cumulative degree distribution 222The cumulative degree distribution $`P(>k)`$ gives the probability that a city has $`k`$ or more connections to other cities, and is defined as $`P(>k)=_{k^{}=k}^{\mathrm{}}p(k^{})`$, where $`p(k)`$ is the probability density function.  for the world-wide air transportation network. The data suggest that $`P(>k)`$ has a truncated power-law decaying tail

$$P(>k)k^\alpha f(k/k_\times )$$
(1)
where $`\alpha =1.0\pm 0.1`$ is the power law exponent, $`f(u)`$ is a truncation function, and $`k_\times `$ is a crossover value that depends on the size of the network. The measured value of the exponent $`\alpha `$ would imply that, as one increases the size of the network, the average degree of the cities are also expected to increase .

The degree of a node is a source of information on its importance. However, the degree does not provide complete information on the role of the node in the network. To start to address this issue, we consider the “betweenness centrality” of the cities comprising the world-wide air transportation network. The betweenness $`B_i`$ of city $`i`$ is defined as the number of shortest paths connecting any two cities that involve a transfer at city $`i`$ . We define the normalized betweenness as $`b_i=B_i/B`$, where $`B`$ represents the average betweenness for the network. We plot, in Fig. 1b, the cumulative distribution $`P(>b)`$ of the normalized betweenness for the world-wide air transportation network. Our results suggest that the distribution of betweennesses for the air transportation network obeys the functional form

$$P(>b)b^\nu g(b/b_\times )$$
(2)
where $`\nu =0.9\pm 0.1`$ is the power law exponent, $`g(u)`$ is a truncation function, and $`b_\times `$ is a crossover value that depends on the size of the network.

A question prompted by the previous results regarding the degree and the centrality of cities is: “Are the most connected cities also the most central?” To answer this question, we analyze first the network obtained by randomizing the world-wide air transportation network (Fig. 1b). We find that the distribution of betweennesses still decays as a power law but, in this case, with a much larger exponent value $`\nu =1.5\pm 0.1`$. This finding indicates the existence of anomalously large betweenness centralities in the air transportation network.

For the randomized network, the degree of a node and its betweenness centrality are strongly correlated—i.e., highly connected nodes are also the most central (Fig. 2a). In contrast, for the world-wide air transportation network it turns out that there are cities that are not hubs—i.e., have small degrees—but that nonetheless have very large betweennesses (Fig. 2a).

To better illustrate this finding, we plot the 25 most connected cities and contrast such a plot with another of the 25 most “central” cities according to their betweenness (Figs. 2b and c). While the most connected cities are located mostly in Western Europe and North America, the most central cities are distributed uniformly across all the continents. Significantly, each continent has at least one central city, which is typically highly-connected when compared to other cities in the continent—Johannesburg in Africa or Buenos Aires and São Paulo in South America. Interestingly, besides these cities with relatively large degree, there are others—such as Anchorage (Alaska, U.S.) and Port Moresby (Papua New Guinea)—that, despite having small degrees are among the most central in the network (Table 1).

## Degree-betweenness anomalies and multi-community networks

Nodes with small degree and large centrality can be regarded as anomalies. Other complex networks that have been described in the literature, like the Internet , do not display such a behavior, and nodes with the highest degree are also those with the highest betweenness . It is, in principle, easy to construct a network in which a node has small degree and large centrality—think, for example, of a network formed by two communities that are connected to one another through a single node with only two links. The relevant question is, however, “what general and plausible mechanism would give rise to scale-free networks with the obtained anomalous distribution of betweenness centralities?”

To answer this question it is useful to consider a region such as Alaska. Alaska is a sparsely populated, isolated region with a disproportionately large—for its population size—number of airports. Most Alaskan airports have connections only to other Alaskan airports. This fact makes sense geographically. However, distance-wise it would also make sense for some Alaskan airports to be connected to airports in Canada’s Northern Territories. These connections are, however, absent. Instead, a few Alaskan airports, singularly Anchorage, are connected to the continental US. The reason is clear, the Alaskan population needs to be connected to the political centers, which are located in the continental US, while there are political constraints making it difficult to have connections to cities in Canada, even to ones that are close geographically . It is now obvious why Anchorage’s centrality is so large. Indeed, the existence of nodes with anomalous centrality is related to the existence of regions with a high density of airports but few connections to the outside. The degree-betweenness anomaly is therefore ultimately related to the existence of “communities” in the network.

The unexpected finding of central nodes with low degree is a very important one because central nodes play a key role in phenomena such as diffusion and congestion , and in the cohesiveness of complex networks . Therefore, our finding of anomalous centralities points to the need to (i) identify the communities in the air transportation network and (ii) establish new ways to characterize the role of each city based on its pattern of intra- and inter-community connections and not merely on its degree.

Community structure— To identify communities in the air transportation network, we use the definition of modularity introduced in Refs. . The modularity of a given partition of the nodes into groups is maximum when nodes that are densely connected among them are grouped together and separated from the other nodes in the network. To find the partition that maximizes the modularity, we use simulated annealing . We display in Fig. 3 the communities identified by our algorithm in the world-wide air transportation network.333We do not know the geographical coordinates of approximately 10% of the 3663 cities in the giant component of the world-wide air transportation network. Those cities are not plotted in the map. Also, some small cities may be misplaced due to duplications in the three-letter code of the corresponding airport.

As we surmised, both Alaska and Papua New Guinea form separate communities.444Alaska and Papua New Guinea are small communities compared to most of the others, which confirms the idea that these are very well defined communities—otherwise, they would be incorporated in a larger community. This fact is particularly important taking into consideration that the community identification algorithm does not take into account the betweenness of the nodes at all. This fact explains the large betweenness centrality of Anchorage and Port Moresby, as they provide the main links to the outside world for the other cities in their communities.

Another significant result is that even though geographical distance plays a clear role in the definition of the communities, the composition of some of the communities cannot be explained by purely geographical considerations. For example, the community that contains most cities in Europe also contains most airports in Asian Russia. Similarly, Chinese and Japanese cities are mostly grouped with cities in the other countries in Southeast Asia, but India is mostly grouped with the Arabic Peninsula countries and with countries in Northeastern Africa. These facts are consistent with the important role of political factors in determining community structure .

Global role of cities— We characterize the role of each city in the air transportation network based on its pattern of intra- and inter-community connections. We first distinguish nodes that play the role of hubs in their communities from those that are non-hubs. Note that cities like Anchorage are hubs in their communities but they are not hubs if one considers all the nodes in the network. Thus, we define the within-community degree of a node. If $`\kappa _i`$ is the number of links of node $`i`$ to other nodes in its community $`s_i`$, $`\overline{\kappa }_{s_i}`$ is the average of $`\kappa `$ over all the nodes in $`s_i`$, and $`\sigma _{\kappa _{s_i}}`$ is the standard deviation of $`\kappa `$ in $`s_i`$, then

$$z_i=\frac{\kappa _i\overline{\kappa }_{s_i}}{\sigma _{\kappa _{s_i}}}$$
(3)
is the so-called $`z`$-score. The within-community degree $`z`$-score measures how “well-connected” node $`i`$ is to other nodes in the community.

We then distinguish nodes based on their connections to nodes in communities other than their own. For example, two nodes with the same $`z`$-score will play different roles if one of them is connected to several nodes in other communities while the other is not. We define the participation coefficient $`P_i`$ of node $`i`$ as

$$P_i=1\underset{s=1}{\overset{N_M}{}}\left(\frac{\kappa _{is}}{k_i}\right)^2$$
(4)
where $`\kappa _{is}`$ is the number of links of node $`i`$ to nodes in community $`s`$, and $`k_i`$ is the total degree of node $`i`$. The participation coefficient of a node is therefore close to one if its links are uniformly distributed among all the communities, and zero if all its links are within its own community.

We hypothesize that the role of a node can be determined, to a great extent, by its within-module degree and its participation coefficient . We define heuristically seven different “universal roles,” each one corresponding to a different region in the $`zP`$ phase-space. According to the within-module degree, we classify nodes with $`z2.5`$ as module hubs and nodes $`z<2.5`$ as non-hubs. Both hub and non-hub nodes are then more finely characterized by using the values of the participation coefficient .

We divide non-hub nodes into four different roles: (R1) ultra-peripheral nodes, i.e., nodes with all their links within their module ($`P0.05`$); (R2) peripheral nodes, i.e., nodes with most links within their module ($`0.05<P0.62`$); (R3) non-hub connector nodes, i.e., nodes with many links to other modules ($`0.62<P0.80`$); and (R4) non-hub kinless nodes, i.e., nodes with links homogeneously distributed among all modules ($`P>0.80`$).

We divide hub nodes into three different roles: (R5) provincial hubs. i.e., hub nodes with the vast majority of links within their module ($`P0.30`$); (R6) connector hubs, i.e., hubs with many links to most of the other modules ($`0.30<P0.75`$); and (R7) kinless hubs, i.e., hubs with links homogeneously distributed among all modules ($`P>0.75`$).

For each city in the world-wide air transportation network, we calculate its within-community degree $`z_i`$ and its participation coefficient $`P_i`$. Then, we assign each city a role according to the definitions above (Figs. 4a and c). Significantly, 95.4% of the cities in the world-wide air transportation network are classified as either peripheral or ultra-peripheral. Additionally, there is a small fraction of non-hub connectors (0.5%). This result suggests that cities which are not hubs in their respective communities rarely have links to many other communities in the air transportation network. This situation is in stark contrast to what happens in some biological networks, in which non-hub connectors seem to be relatively frequent and to play an important role .

The remaining 4.1% of the nodes are hubs. We find approximately equal fractions of provincial and connector hubs. The former include cities that, for historical, political, or geographical reasons, are comparatively not well-connected to other communities. Examples are Denver, Philadelphia, and Detroit, in North America; Stuttgart, Copenhagen, Istanbul, and Barcelona, in the community formed by Europe, North Africa and the former Soviet Union; Adelaide and Christchurch in Oceania; Brasilia in South America; Fairbanks and Juneau in Alaska; and the already discussed case of Port Moresby. Connector hubs include the most recognizable airport hubs in the word: Chicago, New York, Los Angeles, and Mexico City in North America; Frankfurt, London, Paris, and Rome in Europe; Beijing, Tokyo, and Seoul in the South-Eastern Asian community; Delhi, Abu Dhabi, and Kuwait in the community comprising India, the Arabic Peninsula and North-Eastern Africa; Buenos Aires, Santiago, and São Paulo in South America; Melbourne, Auckland, and Sydney in Oceania; and Anchorage in Alaska.

The fractions of cities with each role in the world-wide air transportation network contrast with the corresponding fractions in a randomization of the network (Fig. 4b). In this case, the community identification algorithm still yields certain communities, but the network lacks “real” community structure. The identification of roles enables one to realize that these communities are somehow artificial. Indeed, many cities are either kinless hubs or kinless non-hubs due to the absence of a real community structure, and the network contains essentially no provincial or connector hubs.

## Discussion

We carried out a “systems” analysis of the structure of the world-wide air transportation network. The study enables us to unveil a number of significant results. The world-wide air transportation network is a small-word network in which: (i) the number of non-stop connections from a given city, and (ii) the number of shortest paths going through a given city have distributions that are scale-free. Surprisingly, the nodes with more connections are not always the most central in the network. We hypothesize that the origin of such a behavior is the multi-community structure of the network. We find the communities in the network and demonstrate that their structure can only be understood in terms of both geographical and political considerations.

Our analysis of the community structure of the air transportation network is important for two additional reasons. First, it allow us to identify the most efficient ways in which to engineer the structure of the network. Specifically, having identified the communities, one can identify which ones are poorly connected and the ways to minimize that problem. Second, cities that connect different communities play a disproportionate role in important dynamic processes such as the propagation of infections such as SARS. As, we show, finding the communities is the first step toward identifying these cities.

The existence of communities and the understanding that different cities may have very different impacts on the global behavior of the air transportation system, call for the definition of the role of each city. We address this issue by classifying cities into seven roles, according to their patterns of inter- and intra-community connections. We find that most of the nodes (95%) are peripheral, that is, the vast majority of their connections are within their own communities. We also find that nodes that connect different communities are typically hubs within their own community, although not necessarily global hubs. This finding is in stark contrast with the behavior observed in certain biological networks, in which non-hub connectors are more frequent .

The fact that different networks seem to be formed by nodes with network-specific roles points to the more general question of what evolutionary constraints and pressures determine the topology of complex networks, and how the presence or absence of specific roles affects the performance of these networks.

Acknowledgments

We thank A. Arenas, A. Barrat, M. Barthélémy, A. Díaz-Guilera, A. A. Moreira, R. Pastor-Satorras, M. Sales-Pardo, D. Stouffer, and A. Vespignani for stimulating discussions and helpful suggestions. We also thank OAG for making their electronic database of airline flights available to us, and Landings.com for providing us with the geographical coordinates of the world airports.

Tables

Figures