\section{Introduction}

Significant numbers of massive galaxies have been identified
in the near-infrared with colors so red that they are not
found in surveys that select galaxies at visual wavelengths.  Such
objects have been named {\em extremely red objects} (EROs).  With
surface densities similar to that of bright quasars (Hu \& Ridgeway,
1994), they represent an important component of the population of
high-redshift galaxies.

The identification of EROs coincided with the development and
implementation of infrared arrays for astronomical research, the first
EROs being noted in the K-band surveys of Elston, Rieke, \&
Rieke (1988, 1989).  Since then, many groups have identified objects
with extremely red colors (McCarthy, Persson, \& West 1992; Eisenhardt
\& Dickinson 1992; Persson et al. 1993, Hu \& Ridgeway 1994, Soifer et al.
1994, Dey, Spinrad \& Dickinson 1995, Djorgovski, et al. 1995, among
others), though often with widely different definitions of {\em extremely
red}.  These objects were primarily identified in the fields around
high-redshift active galactic nuclei, i.e. around radio galaxies and
quasars.  In this paper, we define the selection criteria for bright EROs
to be $K^\prime$\ $\leq$ 19\fm0 and with ($R$-$K^\prime$) $\geq$ 6$^{\rm m}$.

The broad-band spectral energy distributions (SEDs) of these galaxies
were initially thought to be fitted best by an old stellar population,
such as found in present-day elliptical galaxies.  In this scenario, a
redshifted ($z > 0.85$) strong 4000\AA-break falling between the $R$\ and
$K$ filter bandpasses and the lack of any appreciable restframe UV light
from a young population of stars are responsible for the extremely
red colors.  An alternative interpretation is that these galaxies are
starbursts or active galactic nuclei (AGN), perhaps triggered by a
merging event.  In this scenario, the presence of significant quantities
of interstellar dust hides the star-forming regions or AGN, considerably
reddening the observed SEDs.  These would also lie at redshifts $1 < z
< 2$, which shifts the restframe UV into the optical bandpass while
still sampling the restframe optical/near-IR wavelengths in the K filter.

Archetypal examples of each of these interpretations have been
identified.  LBDS 53W091 (Dunlop et al. 1996) is a
faint radio galaxy at a redshift of 1.55, derived from an absorption
line spectrum obtained at restframe UV wavelengths.  Spectral
modeling suggests an elliptical galaxy spectrum with an age of up to
3.5 Gyr (Dunlop et al. 1996, Spinrad et al. 1997, but see Bruzual \&
Magris 1997 for an alternative viewpoint), implying a much higher
redshift of formation for this system.  In contrast, HR10 (Hu \&
Ridgeway 1994) has been identified as a dusty starforming galaxy,
perhaps with an active nucleus, at $z=1.44$ (Graham \& Dey 1996).
Detections of this galaxy at submillimeter and millimeter wavelengths
support the interpretation that the SED is dominated by the presence
of dust (Cimatti et al. 1998).

In both cases, the EROs are most likely to lie in the redshift range $1
< z < 2$.  It is difficult to produce such red colors at lower
redshifts, while higher redshift objects become exceedingly luminous.
If EROs are dominated by old stellar populations, then massive galaxy
formation was well underway at $z > 3$.  Dust-dominated EROs, however,
imply that much of the massive galaxy formation could actually occur at
late times, supporting hierarchical galaxy formation models.  The
contrasting implications these two scenarios have for the history of
galaxy formation is the primary driver behind continuing work on these
objects.

We have therefore started a wide-area field survey to select
an unbiased and statistically significant sample of EROs for further
study.  This paper presents initial results from the first of seven
fields, obtained as part of the Calar Alto Deep Imaging Survey (CADIS,
Meisenheimer et al. 1997), each of which will cover $\sim$160 square
arcminutes.  Section 2 describes the observations and reductions.
Section 3 includes a discussion of the sample selection, morphology,
and surface and volume density of the EROs identified here.  Section
4 gives a brief summary of this paper.

\section{Observations and Reductions}

The field observed in this paper is centered on $\alpha_{1950}$ $=$
16$^h$ 23$^m$ 29$^s$ $\delta_{1950}$ $=$ $+55^\circ$ 50$^\prime$
47$^{\prime\prime}$ (hereafter the 16\,h field).  It is one of 7 fields
in the CADIS survey (Meisenheimer et al. 1997), chosen to be free of
bright stars or galaxies, known galaxy clusters, and in regions of
minimal extinction due to Galactic cirrus.  All magnitudes quoted
in this paper are referred to $\alpha$ Lyra.

\subsection{K$^\prime$ imaging}

The $K^\prime$\ imaging was obtained with the Omega-Prime camera (Bizenberger
et al. 1998) at the prime focus of the Calar Alto 3.5m telescope on UT
1996 June 7-8 and UT 1996 August 22-29.  The majority of these nights
were photometric.  Omega-Prime is a direct imaging camera equipped with
a HAWAII 1024$^2$ HgCdTe array.  The image scale is 0.396"/pixel (6.75
arcminute square field of view).  Data were taken in a 2$\times$2 mosaic
to cover as much of the CADIS survey field (16$^\prime$ diameter) as
possible.  The $K^\prime$\ filter ($\lambda_{\rm cent} = 2.12$\,$\mu$m, $\Delta
\lambda = 0.35$\,$\mu$m, see Wainscoat and Cowie 1991) was used because
its bluer bandpass relative to a standard K filter significantly reduces
the thermal background seen by the array.

The data were reduced with a standard infrared reduction algorithm in
IRAF\footnote{IRAF is distributed by the National Optical Astronomy
Observatories, which are operated by the Association of Universities
for Research in Astronomy, Inc., under cooperative agreement with the
National Science Foundation.}.  A sky frame was constructed from 6 to
12 temporally adjacent images, scaled to have the same median counts
before averaging all but the 1 or 2 most extreme values in each pixel,
effectively removing stars from the final sky image.  The sky frame was
then scaled to and subtracted from the data frame and the result divided
by a normalized dome flat to remove pixel-to-pixel fluctuations in
quantum efficiency.

On photometric nights, standard stars from the UKIRT faint standards
list (Casali \& Hawarden 1992) were observed every 2-3 hours at several
different airmasses to determine both a photometric zero point and an
extinction correction.  Individual images were then corrected for
extinction and put onto a common photometric scale.  Night to night
variations of the zero points were under 0.1 magnitude, giving a
measure of the systematic uncertainty in the $K^\prime$\ calibrations.

As IRAF was not able to handle either the large number of frames (1209)
or the large size of the data set (4.8\,GB), a procedure for IDL was
developed to stack the individual images into the final mosaic.  Known
bad pixels were ignored and only good pixels were included in the
mosaic.  The software used a multi-pass sigma-clipping algorithm to
reject cosmic rays and other transient phenomena (e.g meteors and
satellites).

The final mosaic is 2419 $\times$ 2546 pixels (15\farcm96 $\times$
16\farcm80), including regions of higher noise around the edge where
coverage was incomplete.  The deepest portion of the mosaic, with per
pixel exposure times greater than 6000 seconds, covers the central
$\sim$12.5 arcmin square with a seeing of 1.1 arcsec full-width at
half-maximum (FWHM).  A grayscale image of the exposure time map, along
with the relative fields of view for the optical and infrared data, is
shown in Figure~\ref{xtmap}.  The central portion of the $K^\prime$\ mosaic
has a 5$\sigma$ point source detection limit of $K^\prime$\ $=$ 20\fm5, within
an aperture diameter of 2\farcs2, twice the seeing FWHM.  Limits for
extended sources are up to 0\fm50--0\fm75 brighter, depending on the
object morphology.

\placefigure{xtmap}

\subsection{$R$-band imaging}

The $R$-band imaging was obtained with the CAFOS reimaging camera at the
Calar Alto 2.2\,m telescope on UT 1996 May 15 (5$\times$500\,s,
1\farcs3 seeing) and UT 1997 February 5 (4$\times$700\,s 1\farcs6
seeing).  CAFOS was equipped with a LORAL 2048$^2$ CCD
(0.33$^{\prime\prime}$/15\,$\mu$m pixel) for the May observations, and
with a SITe 2048$^2$ CCD (0.53$^{\prime\prime}$/24\,$\mu$m pixel
covering a circular field of view $\sim$16 arcminutes in diameter) for
the remaining observations.  The $R$\ filter used in CADIS ($\lambda_{\rm
cent} = 648$\,nm, $\Delta \lambda = 171$\,nm) is slightly narrower and
bluer than a standard Johnson $R$\ filter.  This difference has no
significant impact on the goals of this survey.

The data were reduced with the MIDAS software package.  Following bias
subtraction and flatfielding, multiple exposures obtained on a single
night were shifted into coincidence, bad pixels and cosmic rays masked
out, then summed to form a deep image for that night.

Flux calibration of the $R$-band data was established relative to two
stars in the 16\,h field for which we have good spectrophotometry.  These
stars were calibrated against an HST UV spectrophotometric standard
(AGK$+81^\circ266$) as well as 3 stars from Oke (1990).  The magnitudes
of the two 16\,h field reference stars in the $R$\ filter were determined by
integrating their calibrated spectra over the $R$\ filter bandpass and
instrumental transmission curve.  The 2$\sigma$ limits were established
within an aperture diameter twice the seeing FWHM, scaling from the
per-pixel noise in the stacked data.  A full description of the
calibration of the CADIS optical data set can be found in Fockenbrock
(1998).

\subsection{Object Detection}

Objects were separately identified on the deep $R$\ and $K^\prime$\ images using
the Source Extractor software (Bertin \& Arnouts 1996) and the
resulting lists merged into a master catalog.   This catalog was then
used to extract photometric magnitudes for all objects in each
individual exposure (mosaics from each night for the $K^\prime$\ data) using
photometry software (Meisenheimer \& R\"oser 1986) incorporating a
Gaussian weighting function, designed to compensate for variations in
seeing between individual exposures.  The procedure recovers accurate
{\em colors} for all objects, but underestimates the total magnitude
for extended objects.  For the faint galaxy images from which the
sample of EROs was drawn, typically subtending only
1-2$^{\prime\prime}$, the underestimate is of order 0.1 magnitude.

A final list of objects was selected from the master catalog,
containing objects with $K^\prime$\ magnitudes above the 5$\sigma$ limit of
20\fm5 that also appear on at least one of the two $R$\ images (from the
LORAL and SITe CCDs).  Figure~\ref{xtmap} shows a contour map and image
of the exposure time, overlaid with the outline of the survey area.
The full area imaged in both $K^\prime$\ and $R$\ covers 154 square arcminutes.

\section{Discussion}

\subsection{ERO Selection}

We select objects with extremely red colors from the color-magnitude
diagram for the full survey area, shown in Figure~\ref{fullcmd}.  We
define the selection criteria for bright EROs as objects with $K^\prime$\
$\leq$ 19\fm0 and ($R$-$K^\prime$) $\geq$ 6$^{\rm m}$.  A more detailed plot of
this region of the color-magnitude diagram is shown in Figure~\ref{zoomcmd}.
Eight objects in our survey field meet these selection
criteria.  Their positions, magnitudes and colors are summarized in
Table~\ref{ERO_data}, while finding charts from the $K^\prime$\ data are
given in Figure~\ref{FCs}.

\placefigure{fullcmd}

\placefigure{zoomcmd}

\placetable{ERO_data}

\placefigure{FCs}

The magnitude limit defining the bright sample of EROs is partially driven
by the current data: most EROs brighter than $K^\prime$\ $\leq$ 19\fm0 are
detected above the 2$\sigma$ limit in $R$\, while most of those fainter
than this are not.  This also represents the practical limit for obtaining
followup near-infrared spectroscopy at 8m-class telescopes (Moorwood \&
Spyromilio 1997) with exposure times of order 1 hour.

In the absence of significant dust extinction, the reddest galaxy colors
will come from the oldest stellar populations at any redshift.  A
present-day elliptical galaxy SED under a no-evolution model peaks in
color at $1 \lesssim z \lesssim 3$ (ignoring the Lyman break at higher
redshift), with an $R$-$K^\prime$ color of about 6.  For an old stellar population,
this color is produced by a combination of the cosmological k-corrections
and a strong 4000\AA-break lying between the two bandpasses.  Since all
galaxies should be bluer than this limit, selecting objects with $R$-$K^\prime$
$\geq$ 6$^{\rm m}$ will isolate a sample of galaxies with extremely red
colors.  The specific choice of $R$-$K^\prime$ $\geq$ 6\fm0 is, however, also
motivated by historical precedence, approximating the limits used to
identify other EROs in the literature.  Such objects were identified as
lying away from the locus of other field objects in the color-magnitude
diagram.  However, the areas surveyed, and thus the number of objects
detected, were generally too small to address this point statistically.

Figure~\ref{colorhist} displays histograms of the $R$-$K^\prime$ colors for objects
with 16$^{\rm m}$\ $\leq$ $K^\prime$\ $\leq$ 20$^{\rm m}$\ in 1$^{\rm m}$\ bins.
Table~\ref{colorhist_data} lists the mean color and FWHM for the
Gaussian least-square fits which are also plotted in Figure~\ref{colorhist}.
For example, in the range 18\fm5 $\leq$ $K^\prime$\ $\leq$ 19\fm5, the bin containing
most of the EROs, the mean color is 3\fm92 and the standard deviation is
1\fm11.  An ($R$-$K^\prime$) $\geq$ 6$^{\rm m}$ selection criterion therefore corresponds
to colors that deviate by more than about two standard deviations from the
mean color.  The extremely red galaxies selected here are thus the reddest
2\% of the population and should consist of the oldest and/or the dustiest
high redshift galaxies.

\placefigure{colorhist}

\placetable{colorhist_data}

There is a distinct trend to redder colors at fainter magnitudes, with
the mean ($R$-$K^\prime$) color becoming 0\fm7 redder between $K^\prime$\ $=$
16\fm5 and $K^\prime$\ $=$ 19\fm5.  This continues the trend noted at
brighter magnitudes in (I-K) colors by Huang et al. (1997).  Over
the range sampled by our data, this trend probably reflects an increasing
fraction of high redshift galaxies.

There are an additional 43 objects with ($R$-$K^\prime$) $\geq$ 6$^{\rm m}$ in the
range 19\fm0 $<$ $K^\prime$\ $\leq$ 20\fm0.  These fainter EROs lie only in
the central 120 square arcminutes of the survey field, where the
overlap of the two sets of $R$-band data reach deep enough to select
EROs with $K^\prime$\ $\leq$ 20\fm0.  We note that since both the color and
color dispersion for the objects in the 16$^{\rm h}$ field increase
with increasing magnitudes, the EROs in this fainter magnitude range
may no longer have truly extreme colors.  The current optical data do
not have sufficient depth to address this more quantitatively.

We chose to set the detection limit in the $R$\ data at 2$\sigma$ to
ensure that all of the EROs are selected, at the expense of some possible
contamination of the sample.  In practice, contamination is likely to be
low.  Consider a Malmquist bias vector defined on the color-magnitude
diagram (Figure~\ref{zoomcmd}).  For objects near the ERO selection
region, it would be dominated by the uncertainties in the $R$-band
magnitudes, since most of the objects there are well-detected at $K^\prime$\
but near the 2$\sigma$ detection limit in $R$.  The Malmquist bias vector
would thus point to {\em bluer} colors and slightly brighter $K^\prime$\
magnitudes.  Inspection of Figure~\ref{zoomcmd} shows that there are
a comparable number of objects which could scatter into the ERO region
as are already there, so contamination of the sample is likely to be
low and would be coming from objects with redder colors.

\subsection{Morphology and Environment}

The radial brightness profiles of the EROs were compared to those of
several stars near each object to determine whether the EROs are
resolved.  Five of the 8 bright EROs are easily resolved, with FWHMs
ranging from 1\farcs4 to 2\farcs3 in 1\farcs1 seeing, and are therefore
likely to be galaxies.  Subtracting the intrinsic seeing in quadrature
gives deconvolved FWHMs ranging from 0\farcs9 to 2\farcs0.  Such sizes
are similar to both HR10 (0\farcs7, Hu \& Ridgeway 1994) and 53W091
(1\farcs3, Dunlop et al. 1996).  A sixth object, C16-ERO3, is formally
resolved, but at a low statistical significance, with a deconvolved
FWHM $\leq$ 0\farcs4.

The remaining 2 bright EROs in the 16\,h field are unresolved.  Their
broadband colors indicate that they are low mass Galactic stars.  One
has subsequently been spectroscopically confirmed as a low mass star of
spectral type L1 (Wolf et al. 1998).  The second is well-detected in
the CADIS filters ($K^\prime$ $=$ 16\fm32), with colors fully consistent with
a star of spectral type M9 (Wolf et al. 1998), though without spectroscopic
confirmation we cannot completely rule out that the object may be a compact
galaxy or AGN at high redshift.  Finding charts for these two objects
can be found in Wolf et al. (1998), and they are not considered further
in this paper.

Only one of the six resolved EROs is significantly non-circular (C16-ERO1).
The deconvolved FWHMs along the major and minor axes are 1\farcs9 and
1\farcs1, respectively, with the major axis oriented at a position angle
(PA) of 143$^\circ$, east of north.  There is also a companion ERO with
$K^\prime$\ $>$ 19\fm0, 2\farcs8 away at PA $=$ 126$^\circ$.  The proximity of
the second ERO and the apparent alignment between the PA of C16-ERO1 and
the direction to the fainter galaxy suggest that this might be an
interacting system.  C16-ERO1 and C16-ERO2 also make the closest pair
of EROs in the current sample, separated by 39\farcs2 at PA $=$ 125\fdg5,
possibly indicating the presence of a cluster or group.  One other object
(C16-ERO4) has a fainter ERO companion nearby (4\farcs9 away at PA $=$
127$^\circ$), though there is no evidence for interaction in the current
data set.

Of the fainter 43 EROs, 37 are either clearly or marginally resolved.
Below K$^\prime \sim$ 19\fm5, it becomes difficult to distinguish
between unresolved and slightly extended sources on the survey data.
Nevertheless, the fainter EROs in this field are clearly dominated by
galaxies.

\subsection{Surface Density}

The 6 bright, extragalactic EROs identified here were found in an area
of 154\,arcmin$^2$, yielding a surface density of 0.039 $\pm$ 0.016
arcmin$^{-2}$ (poisson uncertainties).  This density is a factor of
four higher than found by Hu \& Ridgeway (1994), estimated from
100\,arcmin$^2$ of I and $K^\prime$\ data from the Hawaii Deep K-band Surveys
(Cowie et al. 1994).  In their survey, Hu \& Ridgeway found no galaxies
as red as HR10, implying a surface density of $\lesssim$\,0.01\,arcmin$^{-2}$.
The difference is not statistically significant due to the currently large
uncertainties from the small number statistics as well as a possible
contribution from field-to-field variations.  Combining these data with
those from future fields obtained as part of this survey will reduce this
uncertainty to $\sim$15\%.  Still, the surface density of EROs is
comparable to or higher than that of bright ($B < 22$$^{\rm m}$) quasars, making
the EROs a significant population.

Dey, Spinrad, and Dickinson (1995) suggest that EROs are more common
along the line of sight to or clustered with high redshift radio
galaxies and quasars.  Estimates of the surface density range up to
10-100 times that of the field EROs.  If the EROs are at the same
redshifts as the high-redshift AGN, then many of them would have
luminosities significantly in excess of $L^*$.  With the EROs in the
foreground there is the possibility of a magnification bias, where the
ERO (or a group or cluster containing it) is acting as a gravitational
lens for the background AGN.  This higher surface density in the fields
of high-redshift AGN may also reflect what we term an ``observational
bias.''  This is a combination of ({\em i}) AGN were a common target for
observations in the near-IR, ({\em ii}) unusual objects of any sort were
noted (here objects with extremely red colors) when the data were published,
and ({\em iii}) earlier generations of near-IR imagers had fields of
$\sim$\,1\,arcmin$^2$.  A systematic measurement of the surface density
in AGN fields still remains to be done.

Guiderdoni et al. (1997, their model C) predict the surface density
($\sim 0.05$\,arcmin$^{-2})$ of high-redshift dusty star-forming
galaxies observed at far-infrared wavelengths that is remarkably
similar to the surface density of bright EROs found here, suggesting
that EROs may be dusty starbursts.  The predictions were based on local
ultra-luminous infrared galaxies, actively forming stars but completely
shrouded by dust.  The detection of significant quantities of dust in
HR 10 (Cimatti et al. 1998) supports this picture.

Including the additional 43 EROs with $K^\prime$\ $\leq$ 20$^{\rm m}$, the
surface density is 0.33$\pm$0.05 arcmin$^{-2}$.  As noted above,
the fainter sample is not well defined over the survey area, and
this is only an estimate to the true surface density of faint EROs.
The depth of the $R$\ band data is currently the limiting factor, and
deeper imaging will allow us to make a stronger statement about
these fainter EROs.

\subsection{Volume Density}

We can derive a volume density of EROs if we can estimate the
redshift range over which they occur.  This volume density would
then be directly comparable to known populations of galaxies, perhaps
providing some insight into their nature.

If the ERO population is dominated by old (elliptical) galaxies, the
red color would be coming from a strong 4000\AA\ break.  This spectral
break comes primarily from Balmer continuum limit and Ca H$+$K
absorption lines and can significantly redden an ($R$-$K^\prime$) color when
the break lies between the two filter bandpasses.  For the $R$\ filter
used here, this implies a redshift beyond $z=0.85$.

If, on the other hand, the ERO population is dominated by dusty
starburst galaxies, the red color would be due to dust extinction
reddening a flat SED.  A color of ($R$-$K^\prime$) $\sim$ 6$^{\rm m}$\ can be produced
at $z = 0.85$ by a V-band extinction of $A_V \sim 5^{\rm m}$, assuming
the Seaton (1979) Milky Way extinction law.  In the grayer extinction
law of Calzetti et al. (1994), derived empirically from a sample of
local starburst galaxies, a factor of two more extinction would be
required to give similar results.  The steeply rising extinction curves
at shorter wavelengths, however, imply that less dust extinction would
be required to produce the same reddening at higher redshifts.  In
practice, a moderate continuum decrement across the 4000\AA\ break is
also seen in relatively young stellar populations, contributing to the
red color and relaxing the extinction requirements.  We therefore adopt
a redshift of $z=0.85$ as the lower redshift limit for both types of
galaxies.

For the high redshift limit, consider an L$^*$ galaxy at $z = 2$.  It
would have an apparent magnitude of $K^\prime$\ $\simeq$ 20\fm3 (H$_0 =
70$\,km\,s$^{-1}$\,Mpc$^{-1}$ and q$_0 = 0.1$), assuming only passive
evolution (evolution and k-corrections from Poggianti 1997 and Fioc \&
Rocca-Volmerange 1997 give similar results).  Our selection
of bright EROs at $K^\prime$\ $\leq$ 19$^{\rm m}$\ is thus already sampling
the luminosity function well above L$^*$ at $z = 2$.  Beyond this
redshift, the galaxies would have to be intrinsically as luminous as
bright quasars to be selected as EROs.  We therefore adopt the range
$z \in [0.85,2.0]$ where EROs are most likely to be found with the
selection criteria used here.

The comoving volume covered in this redshift range over our survey
field in the assumed cosmology is 3.55$\times10^5$\,Mpc$^3$.  For
the six bright EROs in our sample, this yields a volume density of
$1.7 \pm 0.7 \times 10^{-5}$\,Mpc$^{-3}$.  Extending the redshift
range down to $z = 0.5$ would only increase the comoving volume by
15\%.  Increasing the upper redshift limit to $z = 3$ would
approximately double the sampled volume.  The volume density derived
here is, therefore, not very sensitive to the assumed redshift range.

We can compare this estimate of the volume density of EROs with that
of other known types of objects.  For quasars (M$_{\rm B \leq -23^m}$),
the volume density in a similar redshift range (Boyle 1991) is only $2
\times 10^{-6}$\,Mpc$^{-3}$, a factor of 8 lower than the bright ERO
density.  A simple extrapolation of the Boyle luminosity function to
M$_{\rm B} \sim -19^m$ (it is strictly only valid for M$_{\rm B} \leq
-21^m$) is required to reach space densities similar to those of the
EROs.  Thus the {\em bright} ERO galaxies appear to be considerably
more numerous than quasars, and perhaps as common as Seyfert Is in
the local universe.  The EROs, however, have been compared to both
old ellipticals and young dusty starbursts, and it is perhaps more
interesting to compare them to these groups of objects.

The red galaxies from the Canada-France Redshift Survey (CFRS, Lilly
et al. 1995) represent a population of objects which to not undergo
significant evolution from $z \sim 1$ to the present, similar to what
would be expected of the EROs if they are dominated by old elliptical
galaxies.  The integrated luminosity function for the CFRS red galaxies
reaches a volume density comparable to the EROs at a luminosity of
$\sim$4\,L$^*$.  Thus, if the ERO population is dominated by old
elliptical galaxies, then the EROs are very massive and suggest that
the most massive systems collapsed early in the history of the universe.
Such galaxies at $z \sim 1.5$ could also mark the location of
high redshift (forming?) clusters or groups of galaxies.

Luminous infrared galaxies are often powered by strong starbursts or
AGN, but emit most of their bolometric luminosity at restframe far-infrared
wavelengths due to the presence of significant quantities of dust.  The
ERO population could be high-redshift analogs to these luminous infrared
galaxies if the EROs are dominated by objects like HR10.  The integrated
luminosity function for the luminous infrared galaxies (Soifer \&
Neugebauer 1991) reaches a volume density comparable to the EROs at a
bolometric luminosity of a few\,$\times\,10^{11}L_\odot$, corresponding
to $\sim$10\,L$^*$.  Thus, if the ERO population is dominated by young,
dusty starbursts similar to the luminous infrared galaxies, then many
massive galaxies were undergoing formative mergers at relatively late
times.

\section{Summary}

We have surveyed an area of 154 sq. arcmin to a 5$\sigma$
point-source limit of $K^\prime$\ $=$ 20\fm5, representing the largest
contiguous field imaged in the near infrared at this depth to date.
From these data, we have selected 8 EROs with $K^\prime$\ $\leq$ 19\fm0 and
with ($R$-$K^\prime$) $\geq$ 6$^{\rm m}$.

Six of the 8 bright EROs are resolved at $K^\prime$\, and are, therefore, most
likely to be galaxies at high redshift ($z \geq 0.85$).  These EROs
have a surface density of 0.039 $\pm$ 0.016 arcmin$^{-2}$, or, with
simple assumptions about their range of redshifts, a volume density
comparable to that of Seyfert galaxies.  Interpretations as to the
nature of these galaxies range from old, evolved elliptical galaxies to
young dusty starbursts or AGN.

The remaining 2 EROs are low-mass Galactic stars, either main sequence
M stars or brown dwarfs (Wolf et al. 1998).  Such low-mass stars can be
seen to distances well beyond 100 pc, depending on their intrinsic
luminosity.  Morphological considerations alone rule out the
possibility that the majority of the EROs with $K^\prime$\ $\leq$ 20$^{\rm m}$\ are
low-mass stars.  The true count of such stars in high-latitude fields
will provide interesting constraints on Galactic structure and the
faint end of the local stellar luminosity function.

EROs appear to be a significant population of galaxies absent from
samples derived from observations at visual wavelengths.  Future
work on these objects will distinguish between the old elliptical
and young dusty starburst interpretations.  Observations at submillimeter
(ground based) and mid- to far-infrared (SIRTF, NGST) wavelengths will
be able to quantitatively determine how much dust, if any, is present.
In either case, the EROs are likely drawn from a population of massive
galaxies at high redshift.   Whichever interpretation prevails, the EROs
will provide important clues to the formation and evolution of massive
galaxies.

\acknowledgments

We are grateful to the team that constructed the OMEGA camera, particularly
M. McCaughrean and P. Bizenberger, without whom the wide field, deep
K$^\prime$ survey would not have been possible.  We also thank the staff
at the Calar Alto observatory for support during the observing runs, and
MPIA for support of the CADIS key project.  This research was supported
by the Max-Planck-Society.

\clearpage


\section{Introduction}

The aim of this paper is to show that random partitions have
a very natural and direct connection to various structures which are
well known in integrable systems.  This
connection is arguably even more natural than, for example, in the case
of random matrices. In a sense,
we show that solitaire (which is related to increasing subsequences in
random permutations and thus to the Plancherel measure on partitions \cite{AD,BDJ})
and soliton have much more in common than the general notion of solitude.

The other character in the title, the infinite wedge space,
is our main technical tool. Starting from the fundamental work
of the Kyoto school, this object
 plays a prominent role in integrable systems.

\subsection{}

More concretely, we consider two types of measures on partitions. The first
one is the Schur measure $\mathfrak{M}$ for which $\mathfrak{M}(\lambda)$ is proportional to $s_\lambda(x)\, s_\lambda(y)$,
where $s_\lambda$ is the Schur function and $x$ and $y$ are two independent sets
of parameters. This measure has a natural representation-theoretic interpretation,
see Section \ref{ha}. The Plancherel measure is a very particular specialization
of $\mathfrak{M}$.

The other measure is the uniform measure which gives equal weight to all
partitions of a given number. Our main findings are the following.

\subsection{Schur measure}
We consider the correlation functions of $\mathfrak{M}$, that is, the probability that
the set $\mathfrak{S}(\lambda)=\{\lambda_i - i+\frac12\}$ contains a given set $X\subset \mathbb{Z}+\frac12$.
Here the $\frac12$'s are just to make certain formulas more symmetric. We prove,
see Theorem \ref{t1},  that
\begin{equation}\label{e01}
\mathfrak{M}\big(\lambda, X\subset\mathfrak{S}(\lambda)\big) = \det \big[ K(x_i,x_j) \big]_{x_i,x_j\in K} \,,
\end{equation}
where the kernel $K$ has a nice generating function, and thus a contour integral
representation, in terms of the parameters of the Schur measure $\mathfrak{M}$, see Theorem \ref{t2}.
Note the
similarity to situation in random matrices \cite{Me}, but also observe that $K$
does not involve any objects even remotely as complicated as polynomials
orthogonal with respect to an arbitrary measure.

We also prove that, as functions of the parameters of the measure $\mathfrak{M}$, the
correlation functions satisfy an infinite hierarchy of PDE's, namely the
Toda lattice hierarchy of Ueno and Takasaki \cite{UT}, see Theorem \ref{TDL}.
Again, this is a well
known phenomenon in mathematical physics that various correlation functions tend to
be $\tau$-functions of integrable hierarchies.

Both of these results are quite straightforward once one interprets the
correlation functions as certain matrix elements in the infinite wedge
space.

\subsection{Uniform measure}

We give a new, more simple, and conceptual proof of the formula from \cite{BO},
which is reproduced in Theorem \ref{t3} below,
for the following averages, called the $n$-point functions,
$$
F(t_1,\dots,t_n) = \sum_{\lambda} q^{|\lambda|} \,
\prod_{k=1}^n \sum_{i=1}^\infty t_k^{\lambda_i-i+\frac12}\,.
$$
These $n$-point functions are sums of
determinants involving genus $1$ theta functions
and their derivatives. The key to these $n$-point functions is
a system of $q$-difference equations \cite{BO} the representation-theoretic
derivation of which is given in Theorem \ref{t5}.

We also investigate the local structure of a typical large partition
and find that, in contrast to the Plancherel measure case \cite{BOO},
it is quite trivial. Locally, a typical partition is a trajectory
of a random walk, with probability to make a vertical (resp.\ horizontal)
step depending on the global position on the limit shape.

\subsection{Connections and applications}

\subsubsection{Asymptotic problems}
Our formula for the correlation functions of the Schur measure
generalizes the exact formulas
for the correlations functions for the Plancherel measure used in \cite{BOO,J}
to prove the conjecture of Baik, Deift, and Johansson \cite{BDJ} about increasing
subsequences in a random permutation and also used in \cite{BOO} to analyze
the local structure of a Plancherel typical partition in the ``bulk'' of the
limit shape. It also generalizes
the more general formula of Borodin and Olshanski \cite{BO2} for the correlation
functions of the so called $z$-measures, see Section \ref{P&z}. The
asymptotics of $z$-measures is important for the harmonic analysis on
the infinite symmetric group $S(\infty)$, see \cite{BO1,KOV}.  Since our
kernel $K$ has a simple integral representation, the formula \eqref{e01} is
particularly suitable for asymptotic investigations. We also point out
that our proofs are considerably simpler and arguably much more
conceptual than the ones given in \cite{BOO,BO2}.

\subsubsection{Toeplitz determinants and Fredholm determinants}
The results of the present paper were used in \cite{BorO}
to solve a problem proposed, independently, by A.~Its and P.~Deift.
The problem was to find a general identity of the form
\begin{equation}\label{e02}
\det(\phi_{i-j})_{1\le i,j\le n} = \det (1-K_\phi)_{\ell_2(\{n,n+1,\dots\})}\,,
\end{equation}
where the kernel $K_\phi$ in the Fredholm determinant admits
an integral representation in terms of the generating function $\phi(z)= \sum \phi_n z^n$
for the entries of the Toeplitz determinant.  See also \cite{BW} for alternative proofs
of the result of \cite{BorO}. Again, \eqref{e02} seems to
be particularly suitable for asymptotic analysis.

\subsubsection{Ramified coverings, moduli spaces, and ergodic theory}
It is known, see e.g.\ \cite{D}, that the Schur and uniform measures
on partitions are related to the enumeration of ramified coverings
of the sphere and the torus, respectively.

In the case of the torus, the exact formula \eqref{e22} was used in \cite{EO}
to compute the asymptotics of the number of ramified coverings with
given ramification type as the degree goes to infinity. These numbers,
which are certain rather complicated polynomials in Bernoulli numbers,
can be identified with volumes of certain moduli spaces and are
important in ergodic theory \cite{EM}. For the coverings of the sphere, see \cite{O2}.

\subsection{Acknowledgments}

I very much benefited from the  discussions with S.~Bloch, A.~Borodin,
P.~Deift, A.~Eskin, M.~Kashiwara, S.~Kerov, G.~Olshanski,
Ya.~Pugaj, T.~Spencer, A.~Vershik and others. Much of this paper has grown
out of my attempt to better understand the work Borodin, Olshanski, and Kerov
on the $z$-measures. I would also
like to thank NSF for financial support
under grant DMS-9801466.

\section{Correlation functions of Schur measure}\label{s2}

\subsection{The Schur measure}

\subsubsection{Definition}

Consider the set of all partitions $\lambda$ and introduce
the following
function of $\lambda$
$$
\mathfrak{M}(\lambda)=\frac1Z\, s_\lambda(x)\, s_\lambda(y) \,
$$
where $s_\lambda$ are the Schur functions \cite{M}
in auxiliary variables $x_1,x_2,\dots$ and $y_1,y_2,\dots$,
and $Z$ is the sum in the Cauchy identity for the Schur
functions
$$
Z= \sum_\lambda s_\lambda(x)\, s_\lambda(y) = \prod_{i,j} (1-x_i y_j)^{-1} \,.
$$
It is clear that if, for exmaple, $\{y_i\}=\overline{\{x_i\}}\subset\mathbb{C}$ and $Z<\infty$ then
$\mathfrak{M}$ is a probability measure on the set of all partitions
which we shall call the \emph{Schur measure}. This measure depends
on countably many parameters. Since our
treatment will be purely algebraic, we shall not require
positivity from the measure $\mathfrak{M}$.

\subsubsection{Harmonic analysis interpretation}\label{ha}

It is well known that the Schur functions $s_\lambda$ are characters of
irreducible representations of the general linear group $GL$ and also
encode characters of the symmetric group $S(n)$ as in \eqref{e007} below.
Hence one can think of $\mathfrak{M}$ as of a representation-valued measure such that
$$
\mathfrak{M}(\lambda)\propto\bla\boxtimes\bla \,.
$$
where $\bla$ is a representation of $GL$ or $S(n)$ corresponding to the
partition $\lambda$ and $\boxtimes$ denotes the outer tensor product.

Recall that the classical spaces of noncommutative harmonic analysis
\begin{equation}\label{dual}
\mathbb{C}[\textup{Mat}(n,m)]\,,\quad  \left(\mathbb{C}^n\right)^{\otimes m}\!, \quad
\mathbb{C}[S(n)]
\end{equation}
all decompose into a multiplicity-free direct sum of representations
of the form $\bla\boxtimes\bla$ with respect to the natural actions of $GL(n)\times GL(m)$,
$GL(n)\times S(m)$, and $S(n)\times S(n)$, respectively. In either of these cases,
$\mathfrak{M}(\lambda)$ can be interpreted as the portion of the space which transforms according
to the symmetry type $\lambda$.

\subsubsection{Power-sum parameters}

It is convenient to introduce another
parameters (also known as the Miwa variables) for the Schur measure
\begin{equation}\label{miwa}
t_k = \frac1k\sum_i x_i^k\,, \quad t'_k=\frac1k\sum_i y_i^k\,,
\quad k=1,2,\dots \,.
\end{equation}
Since the power-sum symmetric functions are free
commutative generators of the algebra of the symmetric
functions, the Schur function $s_\lambda(x)$ is a homogeneous
polynomial in the $t_k$'s of degree $|\lambda|$, where
$$
\deg t_k = k\,.
$$
More precisely, one has \cite{M}
\begin{equation}\label{e007}
s_\lambda(x)= \sum_{\rho=1^{r_1} 2^{r_2} 3^{r_3} \dots} \chi^\lambda_\rho\,
\prod_k \frac{t_k^{r_k}}{r_k!}\,,
\end{equation}
where $\rho=1^{r_1} 2^{r_2} 3^{r_3} \dots$ means that the
partition $\rho$ has $r_k$ parts of length $k$ and $\chi^\lambda_\rho$
is the character of a permutation with cycles $\rho$ in the
representation labeled by $\lambda$.
Also, we have
$$
Z=\exp\left(\sum_{k>0} k\, t_k\, t'_k\right) \,.
$$

\subsubsection{Plancherel measure and $z$-measures}\label{P&z}

In particular, if we set
$$
t=t'=(\sqrt\xi,0,0,\dots)
$$
then $\mathfrak{M}$ specializes to the Poissonized Plancherel measure
\cite{BDJ,BOO,J}
$$
\mathfrak{M}(\lambda)= e^{-\xi}\, \xi^{|\lambda|}
\left(\frac{\dim\lambda}{|\lambda|!}\right)^2\,,
$$
where $\dim\lambda$ is the dimension of the irreducible
representation labeled  by $\lambda$ and $\xi$ is the
parameter of the Poissonization.

The Plancherel measure
corresponds to taking just the dimensions of irreducible pieces
in the last space in \eqref{dual}. Similarly, taking dimensions
in the first space in \eqref{dual} and allowing $z=n$ and $z'=m$
be arbitrary complex one obtains the $z$-measure.

More precisely, if one sets
$$
t_k = \frac{\xi^{k/2} z}k \,, \quad t'_k = \frac{\xi^{k/2} z'}k\,,
\quad k=1,2,\dots\,,
$$
where $z$, $z'$, and $\xi$ are parameters,
 then one obtains the $z$-measure \cite{KOV,BO1,BO2}
\begin{equation}\label{mz}
\mathfrak{M}(\lambda) = (1-\xi)^{zz'} \xi^{|\lambda|} \,
s_\lambda(\underbrace{1,\dots,1}_{\textup{$z$ times}})\,
 s_\lambda(\underbrace{1,\dots,1}_{\textup{$z'$ times}})\,.
\end{equation}
Here $s_\lambda(\underbrace{1,\dots,1}_{\textup{$z$ times}})$ stands for
the polynomial in $z$ which for positive integer values of $z$
specializes to $s_\lambda(x)$ evaluated at $x_1=\dots=x_z=1$ and $x_i=0$,
$i>z$. This polynomial is well known to be \cite{M}
$$
s_\lambda(\underbrace{1,\dots,1}_{\textup{$z$ times}}) =
\prod_{\square\in\lambda} \frac{z+c(\square)}
{h(\square)}\,,
$$
where $c(\square)$ is the content of a square $\square\in\lambda$ and
$h(\square)$ is its hook-length. More precisely, \eqref{mz}
is the mixture of the $z$-measures on partitions of a
fixed number $n=|\lambda|$ by the negative binomial distribution
on $n=0,1,2,\dots$ with parameter $\xi$.
The Plancherel measure is a limit case of the $z$-measures.

\subsection{Formula for correlation functions}\label{s23}

\subsubsection{Definition of the correlation functions}\label{s214}

It is convenient to introduce the following coordinates
on the set of partitions. To a partition $\lambda$ we
associate a subset
$$
\mathfrak{S}(\lambda)=\{\lambda_i-i+1/2\}\subset \mathbb{Z}+{\textstyle \frac12} \,.
$$
For example, $\mathfrak{S}(\emptyset)=\{-\frac 12,-\frac 32,-\frac 52\}$. The geometric
meaning of the set $\mathfrak{S}(\lambda)$ can be seen from the following picture where
the diagram of $\lambda$ is rotated $135^\circ$ and the elements of $\mathfrak{S}(\lambda)$ are
marked by $\bullet$
\begin{center}
\scalebox{0.5}{\includegraphics{part3.EPS}}
\end{center}
Subsets $S\subset\mathbb{Z}+\frac12$ of the form $S=\mathfrak{S}(\lambda)$ can be characterized
by
$$
|S_+|=|S_-|<\infty
$$
where
$$
S_+ = S \setminus \left(\mathbb{Z}_{\le 0} - {\textstyle \frac12}\right) \,, \quad
S_- = \left(\mathbb{Z}_{\le 0} - {\textstyle \frac12}\right) \setminus S \,.
$$
The number $|\mathfrak{S}_+(\lambda)|=|\mathfrak{S}_-(\lambda)|$ is the number of
squares in the diagonal of the diagram of $\lambda$ and the
finite set $\mathfrak{S}_+(\lambda)\cup \mathfrak{S}_-(\lambda) \subset \mathbb{Z}+\frac 12$
is known as the \emph{modified Frobenius coordinates of
$\lambda$}, see \cite{KO}.

Given a finite subset $X\in\mathbb{Z}+{\textstyle \frac12}$, define the correlation
function
$$
\rho(X)=\mathfrak{M}\big(\{\lambda, X\subset \mathfrak{S}(\lambda)\}\big) \,.
$$
Our goal in the present section  is to prove a determinantal
formula for the correlation function
$$
\rho(X)=\det\Big[K(x_i,x_j)\Big]_{x_i,x_j\in X}
$$
where $K$ is a certain kernel which we shall
obtain momentarily.

\subsubsection{Correlation functions as matrix coefficients}

In what follows, we assume that the reader is familiar with
the basics of the infinite wedge space and our notational conventions
as summarized for the reader's convenience in the Appendix.

By the definition of the correlation functions $\rho(X)$ and
the equations \eqref{schur} and \eqref{e14} we have
\begin{equation}\label{e110}
\rho(X)=\frac1Z\sum_{\mathfrak{S}(\lambda)\supset X} s_\lambda(x)\, s_\lambda(y) =
\frac 1Z \left(\Gamma_+(t) \left( \prod_{x\in X} \psi_x \psi^*_x \right) \Gamma_-(t') \, v_\emptyset,v_\emptyset\right).
\end{equation}
We also have, see \eqref{Gamma_vac} and \eqref{Gamma_Gamma},
$$
\Gamma_+(t) \, v_\emptyset = \Gamma_-(t')^*\,  v_\emptyset = v_\emptyset\,, \quad  \Gamma_+(t)\, \Gamma_-(t')  =
Z\, \Gamma_-(t')  \, \Gamma_+(t) \,.
$$
{F}rom this and \eqref{e110} it follows that
\begin{equation}\label{e11}
\rho(X)=
\left( \prod_{x\in X} \Psi_x \Psi^*_x \, v_\emptyset,v_\emptyset\right)\,,
\end{equation}
where
$$
\Psi_k = G\, \psi_k \, G^{-1} \,, \quad
\Psi^*_k = G\, \psi^*_k \, G^{-1} \,, \quad G= \Gamma_+(t)\, \Gamma_-(t')^{-1}\,.
$$
The presentation \eqref{e11} results in the following

\begin{Theorem}\label{t1} We have
\begin{equation}\label{e12}
\rho(X)=\det\Big[K(x_i,x_j)\Big]_{x_i,x_j\in X} \,,
\end{equation}
where $K(x,y)= \left( \Psi_x \Psi^*_y \, v_\emptyset,v_\emptyset\right)$.
\end{Theorem}

 \noindent\emph{Proof.}
The formula \eqref{e12} can be seen as a particular case of Wick's
theorem, see e.g.\  \cite{BS}. An equivalent direct combinatorial
argument goes as follows.  Suppose $|X|=n$ and $X=\{x_1,\dots,x_n\}$. Since
$\psi_{x_i} \psi^*_{x_j} = -  \psi^*_{x_j} \psi_{x_i}$ for $i\ne j$,
the operators $\Psi_{x_i}$ and $\Psi^*_{x_j}$ also anti-commute.
Therefore
$$
\rho(X)=\left(\Psi_{x_1} \cdots \Psi_{x_n}
\Psi^*_{x_n} \cdots \Psi^*_{x_1}   \, v_\emptyset,v_\emptyset\right)\,.
$$
Applied to the vacuum $v_\emptyset$,
the operator
$
\Psi^*_{x_n} \cdots \Psi^*_{x_1}
$
removes an $n$-tuple of the vectors $\underline{s_i}$ and then the operator
$
\Psi_{x_1} \cdots \Psi_{x_n}
$
has to put all of them back, in $n!$ possible orders. Depending on
parity of corresponding permutation, these $n!$ terms
appear with a $\pm$ sign and this produces the determinant.

\subsubsection{Generating function for the kernel $K$}

Form the following generating function
\begin{align*}
\mathsf{K}(z,w)&=\sum_{i,j\in\mathbb{Z}+\frac12} z^i w^{-j}\, K(i,j) \\
&=\left(G\, \psi(z)\, \psi^*(w) \, G^{-1} v_\emptyset,v_\emptyset\right)\,.
\end{align*}
Using the formulas \eqref{e111} we compute
$$
G\, \psi(z) \, \psi^*(w) \, G^{-1} =
\frac{J(z)}{J(w)} \,
 \psi(z)\, \psi^*(w) \,,
$$
where
\begin{equation}\label{J(z)}
J(z)=\frac{\gamma(z,t)}{\gamma(z^{-1},t')}=
\exp\left(\sum_{k\ge 1} t_k\, z^k - \sum_{k\ge 1} t'_k\, z^{-k}\right)=
\prod_i \frac{1-y_i/z}{1-x_i\,z } \,.
\end{equation}
{F}rom this and the formula
$$
\left(\psi(z) \,\psi^*(w)\, v_\emptyset,v_\emptyset \right)= \sum_{j\in\mathbb{Z}_{\ge 0}+\frac12}
\left(\frac wz\right)^j = \frac{\sqrt{zw}}{z-w} \,, \quad |w|<|z| \,,
$$
we obtain the following

\begin{Theorem}\label{t2}  We have for $|w|<|z|$
$$
\mathsf{K}(z,w) =  \frac{\sqrt{zw}}{z-w} \,  \frac{J(z)}{J(w)}\,,
$$
where $J(z)$ is defined by \eqref{J(z)}.
\end{Theorem}

\noindent
The following results generalize Propositions 2.9 and
2.8 of \cite{BOO}, respectively.

\begin{Corollary}\label{c1} Define $J_n(t,t')$ by
$J(z)=\sum_{n\in\mathbb{Z}} J_n \, z^n$. We have
$$
K(i,j)=\sum_{k\in\mathbb{Z}_{\ge 0}+\frac12} J_{i+k} (t,t') \, J_{-j-k} (-t,-t')  \,.
$$
\end{Corollary}

\begin{Corollary}\label{c2} We have
$\displaystyle
\frac{\partial}{\partial t_1} \mathsf{K}(z,w) = \frac{1}{\sqrt{zw}}  \frac{J(z)}{J(w)}\,.
$
\end{Corollary}

In other words, the generating function $\mathsf{K}(z,w)$ factors
after taking the partial with respect to
 $t_1$.  Similar formulas hold for partials
with respect to other parameters.

\subsubsection{$K$ as an integrable kernel}

Recall that a kernel $K$ is called integrable
\cite{IIKS,De} if the kernel $(x-y)\, K(x,y)$
has finite rank. In more invariant terms, $K$ is integrable, if its
commutator with the operator of multiplication by the independent
coordinate has finite rank.

In our case, we have
\begin{multline}
\sum_{i,j\in\mathbb{Z}+\frac12} z^i w^{-j}\, (i-j) K(i,j) =
\left(z\frac{\partial}{\partial z} + w\frac{\partial}{\partial w}\right) \mathsf{K}(z,w)=\\
\frac{\sqrt{zw}}{z-w} \,\frac{J(z)}{J(w)} \,
\left(z\frac{\partial}{\partial z} -w\frac{\partial}{\partial w}\right)
\log (J(z) J(w))  \,.
\end{multline}
Therefore, if the Schur measure is specialized in such a way that
$$
\frac1{z-w} \left(z\frac{\partial}{\partial z} -w\frac{\partial}{\partial w}\right)
\log (J(z) J(w)) = \sum_{i=1}^N f_i(z)\, g_i(w)\,,
$$
for some $N$ and some functions $f_1,\dots,f_N$ and $g_1,\dots,g_N$, then
$K$ is integrable. This happens  if
$$
\frac{\partial}{\partial z} \log J(z) \in \mathbb{C}(z)\,,
$$
where $\mathbb{C}(z)$ stands for the field of rational functions.

\subsection{Remarks}

\subsubsection{Multivariate Bessel functions}
Remark that the functions $J_n(t,t')$  are the most natural multivariate generalization
of the classical Bessel functions and specialize
to the classical Bessel functions in the Plancherel specialization of
the Schur measure \cite{BOO,J}.  The following picture, where the $z^4$ coefficient in
$\exp\left[\frac{x}{2}(z-z^{-1})+\frac{y}{2}(z^2-z^{-2})\right]$ is plotted
as a function of $x$ and $y$, illustrates how naturally yet nontrivially
this bivariate Bessel function interpolates between the classical Bessel
functions of order 2 and 4.
\begin{center}
\scalebox{0.4}{\includegraphics{j4.eps}}
\end{center}
These Bessel functions should not be confused with the (much more complicated)
multivariate Bessel functions which arise as suitable limits of the multivariate
Jacobi polynomials.

\subsubsection{Toda lattice hierarchy}

The operator $\psi_k \psi_k^*$ satisfies \eqref{commOm} and so do the operators
$$
A_X = \prod_{x\in X} \psi_x \psi_x^* \,.
$$
Denote by $X+n$ the translation of the set $X$ by an integer number $n$. We have
\begin{multline*}
(\Gamma_+(t)\, A_X  \, \Gamma_-(t) \, v_n,v_n) =
  (\Gamma_+(t)\, R^{-n} \, A_X \, R^n \, \Gamma_-(t) \, v_\emptyset, v_\emptyset) \\ =
 (\Gamma_+(t)\, A_{X-n} \, \Gamma_-(t) \, v_\emptyset, v_\emptyset) = Z\, \rho(X-n)
\end{multline*}
Thus, we have the following
\begin{Theorem}\label{TDL} For any set $X\subset \mathbb{Z}+\frac12$ the sequence
$\left\{Z\,\rho(X-n)\right\}_{n\in\mathbb{Z}}$ of functions of $t$ and $t'$
is a sequence of $\tau$-functions for the Toda lattice hierarchy.
\end{Theorem}

In the particular case when
$$
X=\{-\frac12,-\frac32,-\frac52,\dots\}\,,
$$
these $\tau$-functions are Fredholm determinants with the kernel $K$.
The Toda lattice equations for these Fredholm determinants are
related to the results of \cite{TW1,TW2} and of the recent paper \cite{AM}.
See \cite{BorO} for
a discussion of the relationship between these Fredholm determinants and
Toeplitz determinants.

The times $t$ and $t'$ in the Toda lattice hierarchy admit the
following combinatorial interpretation. We have
$$
\al_{-n} \, v_\mu = \sum_{\textup{$\lambda/\mu$ is a rim hook of size $n$}}
 (-1)^{\textup{height of $\lambda/\mu$ + 1}} \, v_\lambda \,,
$$
where by a rim hook we mean a skew diagram $\lambda/\mu$ which is connected and lies
on the rim of $\lambda$ and the height of $\lambda/\mu$ is the number of rows it
occupies. The action of the operator $\Gamma_-(t')= \prod_n \exp(t'_n  \al_{-n})$
can be interpreted as holding our initial partition for the time $t_n$ under
a Poisson random stream of $n$-hooks, where each hook attaches to $\al$ with
a $\pm$ sign according to its height.

Also remark that in the Plancherel case the correlations depend
only on the product $t_1 t'_1$ and so \eqref{Toda} becomes an ODE
equivalent to the usual Toda lattice \cite{T}.

\subsubsection{Other algebras and Kerov's construction}\label{s241}
Instead of the operators $\al_n$, which generate an action of
the Heisenberg algebra, one can consider representation of
other algebras on $\Lambda^{\frac\infty2}V$. Sandwiching products of $\psi_k \psi^*_k$
between some raising and lowering operators will again
produce formulas similar to ours. The analog of the kernel $K$
will involve matrix elements of this representation.
In particular, an $\mathfrak{sl}_2$-action on $\Lambda^{\frac\infty2}V$ leads again to $z$-measures
and provides another proof of the Borodin-Olshanski formula, see
\cite{O3}. The action of $\mathfrak{sl}_2$ in the basis $\{v_\lambda\}$ gives
certain operators on partitions considered by S.~Kerov
in his analysis of the $z$-measure \cite{Ke}.

\section{The uniform measure}\label{s3}

In this section we shall consider the uniform measure
on partitions on $n$ and the related measure
$$
\mathcal{M}_q(\lambda)= (q;q)_\infty \, q^{|\lambda|} \,, \quad q\in [0,1)\,,
$$
on the set of all partitions. The normalization
factor
$$
(q;q)_\infty=(1-q)(1-q^2)\dots
$$
makes $\mathcal{M}_q(\lambda)$ a probability measure;  this is clear
from the classical identity
$$
\sum_\lambda q^{|\lambda|} = \prod_{n\ge 1} \frac1{1-q^n} \,.
$$

\subsection{The $n$-point functions}

\subsubsection{Definition}

Again, since our treatment will be purely algebraic, we
shall not require $\mathcal{M}_q$ to be positive and we shall allow
$q$ to be any number from the unit disk $|q|<1$.

Let $f(\lambda)$ be a function on the set of partitions. We set
$$
\langle f\rangle_q\ =(q ;q)_\infty \sum_\lambda f(\lambda)\, q^{|\lambda|}\,,
$$
assuming this expectation converges.

Fix an integer $n\ge 1$ and variables $t_1,\ldots,t_n$. These
are not in any way related to the variables $t$ of Section \ref{s2}.
We consider the following $n$-point function
\begin{equation}\label{e21}
F(t_1,\ldots,t_n) := \left\langle \prod_{k=1}^n\left(\sum_{i=1}^\infty
t_k^{\lambda_i-i+\frac12}\right)\right\rangle_q.
\end{equation}

This series converges if $1<|t_i|$ and $|t_1\cdots t_n|<|q|^{-1}$ and admits a
meromorphic analytic continuation onto $\mathbb{C}^n$, see \cite{BO}
and below. The Laurent coefficients of \eqref{e21} near the point
$t_1=\dots=t_n=1$ are the averages of the form $\langle f \rangle_q$ for
polynomial functions $f$ and all averages of polynomials can
be reconstructed in this fashion.

\subsubsection{Formula for the $n$-point function}

The following formula for the $n$-point functions
was obtained in \cite{BO}. Consider the following
genus $1$ theta function
\begin{align*}
\Theta(x) = \Theta_{11}(x;q)&=\sum_{n\in{\mathbb{Z}}}
(-1)^n q^{\frac{(n+1/2)^2}2} x^{n+1/2} \\
&=q^{1/4}  (x^{1/2}-x^{-1/2}) (q;q)_\infty\, (q x;q)_\infty (q/x;q)_\infty \,.
\end{align*}
This is the only odd genus $1$ theta function and its precise
normalization is not important because the main formula
will be homogeneous in $\Theta$.  Set
\begin{equation*}
\Theta^{(k)}(x):= \left(x\,
\frac{d}{dx}\right)^k \Theta(x;q)\,.
\end{equation*}
We have the following result from \cite{BO}

\begin{Theorem}[\cite{BO}]\label{t3}
\begin{equation}\label{e22}
\vspace{-3 \jot}
F(t_1,\dots,t_n)=
\sum_{\sigma\in S(n)}\, \frac
{\displaystyle \det\left( \frac{\displaystyle \Theta^{(j-i+1)}(t_{\sigma(1)}\cdots
    t_{\sigma(n-j)})}{\displaystyle (j-i+1)!}
\right)_{i,j=1}^n}
{\displaystyle \Theta(t_{\sigma(1)})\, \Theta(t_{\sigma(1)} t_{\sigma(2)}) \dots
\Theta(t_{\sigma(1)}\cdots t_{\sigma(n)})}
\end{equation}
Here $\sigma$ runs through all permutations $S(n)$ of $\{1,\dotsc,n\}$, the
matrices in the numerator have size $n\times n$, and we define
$1/(-n)!=0$ if $n\ge 1$.
\end{Theorem}

\subsubsection{Strategy of proof}

The proof of \eqref{e22} given in \cite{BO} is based on the following
two properties of the $n$-point functions. First, these functions
satisfy a $q$-difference equation:
\begin{Proposition}[\cite{BO}]\label{p1} We have
\begin{multline}\label{e23}
F(qt_1,t_2,\dots,t_n)=-q^{1/2} t_1 \dots t_n \times \\ \sum_{s=0}^{n-1}
(-1)^s
\sum_{1<i_1<\dots < i_s \le n} F(t_1 t_{i_1} t_{i_2} \cdots t_{i_s},
\dots,\widehat{\,t_{i_1}}, \dots,\widehat{\,t_{i_s}}, \dots) \,,
\end{multline}
where hats mean that the corresponding term should be
omitted.
\end{Proposition}
Observe that since the $n$-point functions are obviously
symmetric,
we have parallel formulas for $F(t_1,\dots,qt_k,\dots,t_n)$.

Another property of the $n$-point functions
needed in the proof of \eqref{e22} is the description of
their singularities. All singularities
lie on the divisors
$$
q^m t_{i_1} t_{i_2} \dots t_{i_k} = 1 \,, \quad m\in\mathbb{Z}\,,
$$
where $\{i_1,\dots,i_k\}\subset \{1,\dots,n\}$ is an
arbitrary subset. By symmetry, it suffices to consider
the divisor $q^m t_1 \dots t_k=1$ on which we have
\begin{equation}\label{e24}
F(t_1,\dots,t_n)=(-1)^{m}
 \frac{q^{m^2/2} m^{k-1}}{q^m t_1\dots t_k-1}
\frac{F(t_{k+1},\dots,t_n)}{(t_{k+1}\dots t_n)^m}+ \dots \,,
\end{equation}
where dots stand for terms regular on the divisor
 $q^m t_1 \dots t_k=1$  and we
assume that
$$
F(t_{k+1},\dots,t_n)=1\,,\quad k=n\,.
$$

Since we have a $q$-difference equation for the
$n$-point functions, it suffices to prove the following
\begin{Proposition}[\cite{BO}]\label{p2} We have
$$
F(t_1,\dots,t_n)=
 \frac{1}{t_1-1}\,
{F(t_{2},\dots,t_n)}+ \dots \,,
$$
and $F(t_1,\dots,t_n)$ is regular on the
divisors $t_1 t_2 \cdots t_k=1$ for $k>0$.
\end{Proposition}

In this section  we shall give a representation-theoretic
proof of Propositions \ref{p1} and \ref{p2}.

\subsection{Singularities of the $n$-point functions}

\subsubsection{Convergence of traces}

 Introduce the following operators on $\Lambda^{\frac\infty2}V$
\begin{align*}
\mathsf{T}(t)&=\sum_{k\in\mathbb{Z}+\frac12} t^k\, \psi_k\, \psi^*_k\,, \\
\nr{\mathsf{T}(t)}&=\sum_{k\in\mathbb{Z}+\frac12} t^k\,
\nr{\psi_k\, \psi^*_k}\,,
\end{align*}
where the normal ordering was defined in \eqref{e112}.
The action of the operator $\nr{\mathsf{T}(t)}$ in the infinite wedge
module is well defined for any $t$. The action of $\mathsf{T}(t)$
makes sense only if $|t|>1$ in which case we have
\begin{equation}\label{e25}
\nr{\mathsf{T}(t)} = \mathsf{T}(t) - \sum_{k=-\infty}^{-1/2} t^k  = \mathsf{T}(t) -
\frac1{t^{1/2}-t^{-1/2}} \,.
\end{equation}

Set $\sigma_{\!S}= \sum_{k\in S_+} t^k - \sum_{k\in S_-} t^k$. It is clear that
$$
\nr{\mathsf{T}(t)}\, f_S = \sigma_{\!S} \, f_S
\,.
$$
We have the following estimates
\begin{alignat*}{2}
\left|\sigma_{\!S}\right|& \le |t|^{\|S_+\|}+|S_-|\,,&\quad &|t|\ge 1\,,\\
\left|\sigma_{\!S}\right|& \le |t|^{-\|S_-\|}+|S_+|\,,&\quad &|t|\le 1\,,
\end{alignat*}
where
$$
\|S\|=\sum_{k\in S} |k| \,.
$$
Since
$$
q^H \, f_S = q^{\|S_+\| + \|S_-\|} \, f_S
$$
it follows that the following trace
\begin{equation}\label{e26}
\nr{\FF}(t_1,\dots,t_n)=
\tr_{\Lambda^{\frac\infty2}V} \,
\left(q^H\, \nr{\mathsf{T}(t_1)}\, \nr{\mathsf{T}(t_2)} \cdots  \nr{\mathsf{T}(t_n)}
\right)
\end{equation}
converges absolutely provided
$$
|q| < |t_{i_1} t_{i_2} \dots t_{i_k}| < |q|^{-1}
$$
for any subset $\{i_1,\dots,i_k\}\subset \{1,\dots,n\}$.
Let us denote by $\Delta$ this domain of the convergence
of \eqref{e26}.

\subsubsection{Poles of the $n$-point functions}

Consider the trace
\begin{equation}\label{e27}
\mathsf{F}(t_1,\dots,t_n)=
\tr_{\Lambda^{\frac\infty2}V}
\left(q^H\, \mathsf{T}(t_1)\, \mathsf{T}(t_2)\cdots \mathsf{T}(t_n)
\right)
\end{equation}
which converges absolutely provided $|t_i|>1$ and
$|t_1\cdots t_n|<|q|^{-1}$. It is clear that using \eqref{e25}
one obtains a meromorphic continuation of $\mathsf{F}$ onto the
domain $\Delta$.

In particular, consider the zero charge subspace
$\Lambda_0 \subset \Lambda^{\frac\infty2}V$, that is, the kernel of
the charge operator $C$. This subspace is preserved
by operators $\nr{\mathsf{T}(t)}$ and $H$. It is spanned by the vectors
$v_\lambda$, where $\lambda$ is a partition, and therefore
$$
F(t_1,\dots,t_n)=(q;q)_\infty\, \tr_{\Lambda_0}
\left(q^H\, \mathsf{T}(t_1)\, \mathsf{T}(t_2)\cdots \mathsf{T}(t_n)
\right) \,.
$$
In fact, this is the original definition of the $n$-point
function in \cite{BO}. From the above discussion, we
have the following immediate conclusion

\begin{Theorem} The $n$-point function $F(t_1,\dots,t_n)$
admits a meromorphic continuation onto the domain $\Delta$
with simple poles along the divisors $t_i=1$ and no
other singularities. We have
$$
F(t_1,\dots,t_n)=
 \frac{1}{t^{1/2}_1-t^{-1/2}_1}\,
{F(t_{2},\dots,t_n)}+ \dots \,,
$$
where dots stand for terms regular on the divisor $t_1=1$.
\end{Theorem}

\subsection{The $q$-difference equation}

\subsubsection{The effect of the translation operator}

Recall that $R$ denotes the translation operator \eqref{defR}.
It is clear that
\begin{equation}\label{e210}
R^{-1} \, \mathsf{T}(t) \,R = t\, \mathsf{T}(t) \,.
\end{equation}
Therefore, from \eqref{e29} and  \eqref{decC}   we obtain
\begin{align}
\mathsf{F}(t_1,\dots,t_n)&=\sum_k \tr_{R^k \, \Lambda_0}
\left(q^H\, \mathsf{T}(t_1)\, \mathsf{T}(t_2)\cdots \mathsf{T}(t_n)
\right)\notag\\
&=\sum_{k} q^{\frac{k^2}2} \, (t_1 t_2 \cdots t_n)^k
\tr_{\Lambda_0}
\left(q^H\, \mathsf{T}(t_1)\, \mathsf{T}(t_2)\cdots \mathsf{T}(t_n)
\right)\notag\\
&=\Theta_3(t_1 t_2 \cdots t_n;q)\, (q;q)_\infty \,F(t_1,\dots,t_n) \label{e211}\,,
\end{align}
where
\begin{equation}\label{e212}
\Theta_3(t;q)=\sum_{n\in\mathbb{Z}} q^{\frac{n^2}2} \, t^n \,.
\end{equation}
The function $\Theta_3$ satisfies the following $q$-difference
equation
$$
\Theta_3(qt;q)=q^{-1/2}t^{-1} \, \Theta_3(t;q) \,.
$$
We conclude that the $q$-difference equation \eqref{e23} for
the $n$-point function is equivalent to the following
equation
\begin{multline*}
\mathsf{F}(qt_1,t_2,\dots,t_n)=\\
\sum_{s=0}^{n-1}
(-1)^{s+1}
\sum_{1<i_1<\dots < i_s \le n} \mathsf{F}(t_1 t_{i_1} t_{i_2} \cdots t_{i_s},
\dots,\widehat{\,t_{i_1}}, \dots,\widehat{\,t_{i_s}}, \dots) \,,
\end{multline*}
which will now be established.

\subsubsection{Proof of the $q$-difference equation}

Note that
$$
[\mathsf{T}(t),\psi^*_k]=-t^k\, \psi^*_k \,.
$$
It follows that
\begin{equation}\label{e214}
\mathsf{T}(t_1) \cdots \mathsf{T}(t_n) \, \psi^*_k =
\sum_{P\subset\{1,\dots,n\}} (-1)^{|P|} \left(\prod_{i\in P}
t^k_i\right) \,
\psi_k^* \, \prod_{i\notin P} \mathsf{T}(t_i) \,.
\end{equation}

Let us assume that $|t_i|>1$ and $|t_1\cdots t_n|<q^{-1}$ so that
the trace \eqref{e27} converges absolutely. Let us write simply
``$\tr$'' instead of $\tr_{\Lambda^{\frac\infty2}V}$. From \eqref{e214} we have
\begin{multline}\label{e215}
\tr q^H \, \psi_k \mathsf{T}(t_2) \cdots \mathsf{T}(t_n)  \, \psi^*_k  =
\\
\sum_{P\subset\{2,\dots,n\}} (-1)^{|P|} \left(\prod_{i\in P}
t^k_i\right)\,
\tr q^H \, \psi_k\, \psi_k^* \, \prod_{i\notin P} \mathsf{T}(t_i) \,.
\end{multline}
On the other hand, because of the absolute convergence, we have
\begin{align*}
\tr q^H \, \psi_k\,  \mathsf{T}(t_2) \cdots \mathsf{T}(t_n) \, \psi^*_k
&= \tr \psi^*_k\, q^H \, \psi_k\,  \mathsf{T}(t_2) \cdots \mathsf{T}(t_n) \\
&= q^k  \tr q^H \,\psi^*_k\,  \psi_k\,  \mathsf{T}(t_2) \cdots \mathsf{T}(t_n) \,,
\end{align*}
where we used the relation
$
q^H\, \psi^*_k\, q^{-H} = q^{-k} \, \psi^*_k
$\,.
Summing the equations \eqref{e215} we obtain
\begin{multline*}
\tr q^H \, \widetilde{\TT}(qt_1)\,  \mathsf{T}(t_2) \cdots \mathsf{T}(t_n) =\\
\sum_{s=0}^{n-1}
(-1)^{s}
\sum_{1<i_1<\dots < i_s \le n} \mathsf{F}(t_1 t_{i_1} t_{i_2} \cdots t_{i_s},
\dots,\widehat{\,t_{i_1}}, \dots,\widehat{\,t_{i_s}}, \dots) \,,
\end{multline*}
where
$$
\widetilde{\TT}(t)=\sum_k t^k \,\psi^*_k\,  \psi_k\,.
$$
From the relations
\begin{alignat*}{2}
\mathsf{T}(t) &= \nr{\mathsf{T}(t)}+\frac1{t^{1/2}-t^{-1/2}}\,,& \quad &|t|>1 \,, \\
\widetilde{\TT}(t) &= -\nr{\mathsf{T}(t)}-\frac1{t^{1/2}-t^{-1/2}}\,,& \quad &|t|<1\,,
\end{alignat*}
it follows that
$$
\tr q^H \, \widetilde{\TT}(qt_1)\,  \mathsf{T}(t_2) \cdots \mathsf{T}(t_n) =
- \mathsf{F}(qt_1,t_2,\dots,t_n)\,,
$$
where the right-hand side is computed by analytic continuation.
This establishes the following

\begin{Theorem}\label{t5}
 The meromorphic function $\mathsf{F}(t_1,\dots,t_n)$
satisfies the following $q$-difference equation
\begin{multline*}
\mathsf{F}(qt_1,t_2,\dots,t_n)=\\
\sum_{s=0}^{n-1}
(-1)^{s+1}
\sum_{1<i_1<\dots < i_s \le n} \mathsf{F}(t_1 t_{i_1} t_{i_2} \cdots t_{i_s},
\dots,\widehat{\,t_{i_1}}, \dots,\widehat{\,t_{i_s}}, \dots) \,.
\end{multline*}
\end{Theorem}

In view of \eqref{e211}, this establishes the $q$-difference equation
\eqref{e23} and this concludes the proof of Propositions \ref{p1} and \ref{p2}

\subsection{An $\mathcal{M}_q$-typical  partition is locally a random walk}\label{s4}

The $n$-point functions encode global properties of a measure $\mathcal{M}_q$.
In this subsection we turn to local properties of an $\mathcal{M}_q$-typical
partitions and find that their asymptotics is quite trivial. More precisely,
as $q\to 1$, the local structure of an  $\mathcal{M}_q$-typical partition
becomes a random walk with probability to go horizontally or vertically
 depending on the global position on the rim
of the limit shape.

\subsubsection{Correlation functions for Frobenius coordinates}

Recall that the modified Frobenius coordinates of a partition $\lambda$
are defined as follows
$$
\Fr\lambda = \mathfrak{S}(\lambda)_+ \cup \mathfrak{S}(\lambda)_- \subset \mathbb{Z}+{\textstyle \frac12}\,,
$$
where the set $\mathfrak{S}(\lambda)=\{\lambda_1-1/2,\lambda_2-3/2,\dots\}$ was defined in \eqref{e15}.
Given a set $X\subset \mathbb{Z}+{\textstyle \frac12}$ define the corresponding
\emph{correlation function} by
\begin{align*}
\varrho(X,q)&=\mathcal{M}_q\left(\left\{\lambda, X\subset\Fr\lambda\right\}\right)\\
&=(q;q)_\infty\, \sum_{\Fr\lambda\supset X}  q^{|\lambda|} \,,
\end{align*}
here the summation is over all $\lambda$ such that $X\subset\Fr\lambda$.

Recall the following triple product formula for the theta function
$\Theta_3(z,q)$ defined in \eqref{e212}
$$
\Theta_3(z,q)=(q;q)_\infty \, (q^{1/2} z; q)_\infty \, (q^{1/2} /z; q)_\infty \,.
$$
Also recall that one way to prove this triple product formula uses
the identity
\begin{equation}\label{e31}
[z^0] \,  (q^{1/2} z; q)_\infty \, (q^{1/2} /z; q) = \sum_{\lambda} q^{|\lambda|}\,,
\end{equation}
where the symbol $[z^0]$ in the left-hand side means that we take
the constant term in $z$. The combinatorial meaning of \eqref{e31} is the
following: on the right, we have all partitions, on the left, we have
all possible modified Frobenius coordinates.  The same argument establishes
the following formula. If  $X=\{x_1,x_2,\dots\}$ then
$$
 \varrho(X,q) =[z^0] \left( \prod_{x_i>0} \frac{q^{x_i}\, z}{1+q^{x_i}\, z} \,
\prod_{x_i<0} \frac{q^{-x_i}\, z^{-1}}{1+q^{-x_i}\, z^{-1}} \,  \Theta_3(z,q) \right) \,.
$$
More compactly, this can be written as follows
$$
\varrho(X,q) =[z^0] \left( \prod_{x\in X} \frac{(q^x\, z)^{\sgn(x)/2}}{q^{x/2}\, z^{1/2}+q^{-x/2}\, z^{-1/2}} \,
 \,  \Theta_3(z,q) \right) \,.
$$
Replacing the constant term extraction by an integral we obtain
\begin{equation}\label{e32}
\varrho(X,e^{-2\pi r}) = \int_{-1/2}^{1/2} \,\,\prod_{x\in X} \frac{e^{\sgn(x)\, \pi(is-rx)}}{2\cosh(\pi(is-rx))}
\, \Theta_3(e^{2\pi i s}, e^{-2\pi r}) \, ds\,,
\end{equation}
where $r>0$.

\subsubsection{Asymptotics of the correlation functions}

Let $N=N(q)$ denote the expected size of a partition with
respect to the measure $\mathcal{M}_q$
\begin{align*}
N(q)&=(q;q)_\infty  \sum_\lambda |\lambda|\, q^{|\lambda|}\\
&= - q \, \frac{d}{dq} \log (q;q)_\infty = \frac 1{24} + G_2(q) \,,
\end{align*}
where $G_k(q)=\frac{\zeta(1-k)}2+\sum_n q^n \sum_{d|n} d^{k-1}$ is the
Eisenstein series. This is, in fact, the simplest example of a polynomial
average with respect to $\mathcal{M}_q$ which, in the general case, are computed
by the formula \eqref{e22}. From the quasi-modular property of $G_2$ it
follows that
$$
N(e^{-2\pi r})  \sim \frac{\zeta(2)}{(2\pi r)^2}=\frac1{24 r^2} \,, \quad r\to +0 \,.
$$

Let us assume that the set $X$ changes with $q\to 1$ in such a way that
all limits
$$
\frac{x_i}{\sqrt N} \to  \xi_i \,, \quad x_i\in X\,, \quad q\to 1 \,,
$$
exist. That is, we consider the scaling of our typical partition by
the factor of $\sqrt N$ in both directions, where $N$ is the area of the
typical diagram. We are interested in the asymptotics of the
correlation functions \eqref{e32}. Observe that
by our assumption about the growth of $X$ we have
$$
x_i \, r \to  \frac{\xi_i}{2\sqrt 6} \,, \quad q=e^{-2\pi r} \,, \quad r\to 0 \,.
$$

The Jacobi imaginary transformation \cite{HTF} gives
$$
\Theta_3(e^{2\pi i s}, e^{-2\pi r}) = r^{-1/2}
\sum_{n\in\mathbb{Z}} \exp\left(-\frac{\pi\,(s-n)^2}r\right) \,.
$$
Clearly, as $r\to +0$, only the $n=0$ summand in the above formula is
relevant and the asymptotics of the integral \eqref{e32} is determined
by the value of the integrand at $s=0$. Therefore, we obtain the following

\begin{Theorem}\label{t6} Suppose that, as $q\to1$, a finite set $X\subset\mathbb{Z}+{\textstyle \frac12}$
varies in such a way that all limits
$$
\frac{x_i}{\sqrt N} \to  \xi_i \,, \quad x_i\in X\,, \quad q\to 1 \,,
$$
exist, where $N=\frac 1{24} + G_2(q)$ is the expectation of $|\lambda|$. Then
\begin{equation}\label{e33}
\lim_{q\to 1} \varrho(X,q)  = \prod_{x_i\in X}
\left(1+\exp\left(\frac{\pi|\xi_i|}{\sqrt{6}}\right)\right)^{-1} \,.
\end{equation}
\end{Theorem}

The error term in this formula can be also found by the Laplace
method.
This theorem says that locally a typical random
partition is just a trajectory of a  random walk.

Theorem \ref{t6} is in agreement with Vershik's theorem \cite{V} about
the limit shape of a typical partition with respect to the uniform
measure. Namely, Vershik's theorem asserts that after the
scaling by the square root of the area in both direction, a typical
partition converges to the following limit shape
$$
\exp\left(-\frac{\pi x}{\sqrt{6}}\right)+\exp\left(-\frac{\pi y}{\sqrt{6}}\right)=1
$$
In new coordinates, $u=x-y$, $v=x+y$ this limit shape becomes
$$
v=\Upsilon(u) \,,  \quad \Upsilon(u)=
\frac{2\sqrt{6}}{\pi} \log \left(2 \cosh \left(\frac{\pi u}{2\sqrt 6}\right)\right) \,.
$$
From this we obtain the formula
$$
\frac{1-|\Upsilon'(u)|}{2} = \left(1+\exp\left(\frac{\pi|u|}{\sqrt{6}}\right)\right)^{-1}\,,
$$
which is precisely what we see in the right-hand side of \eqref{e33}.
See the discussion in in Section 1.3 of  \cite{BOO} for why this should be the case.


\section{Abstract}
We consider the thermal properties of cold, dense clouds of molecular
hydrogen and atomic helium.  For cloud masses below
$10^{-1.7} \mathrm{M_{\sun}}$, the internal pressure is sufficient to permit
the existence of particles of solid or liquid hydrogen at temperatures
above the microwave background temperature.  Optically thin thermal
continuum emission by these particles can balance cosmic-ray heating
of the cloud, leading to equilibria which are thermally {\it stable\/}
even though the heating rate is independent of cloud temperature.  For
the Galaxy, the known heating rate in the disk sets a minimum mass of
order $10^{-6}\mathrm{M_{\sun}}$ necessary for survival.  Clouds of this
type may in principle comprise most of the dark matter in the Galactic
halo.  However, we caution that the equilibria do not exist at
redshifts $z\ga1$ when the temperature of the microwave background was
substantially larger than its current value; the formation and
survival of such clouds to the present epoch therefore remain open
questions.
\end{abstract}

\keywords{ISM: clouds --- galaxies: halos --- dark matter}

\section{Introduction}

Walker \& Wardle (1998) showed that a population of neutral, AU-sized
clouds in the Galactic halo could be responsible for the ``Extreme
Scattering Events'' (ESEs) observed in the radio flux towards several
quasars (Fiedler et al.  1987, 1994).  In this model the cloud
surfaces are exposed to UV radiation from hot stars in the Galactic
disk, producing a photo-ionised wind.  When one of these clouds
crosses the line of sight to a compact radio source, the flux varies
as a result of refraction by the ionised gas (cf.  Henriksen \& Widrow
1995).  This model explains the observed flux variations quite
naturally; but if the clouds are self-gravitating, then the ESE event
rate implies that the cloud population comprises a significant
fraction of the Galaxy's mass.

This halo cloud population cannot contain much dust mixed with the gas
as this would lead to optical extinction events of distant
stars: either the clouds have extremely low metallicity, or any dust
grains have sedimented to the cloud centre.  Given this, several
factors make the clouds difficult to detect (Pfenniger, Combes
\& Martinet 1994): cold molecular hydrogen is, by and large,
invisible; the clouds are small; they are transparent in most regions
of the electromagnetic spectrum; and they cover a small fraction of the
sky. The clouds are not sufficiently compact to cause gravitational
lensing towards the LMC, although Draine (1998) has shown that there
is substantial optical refraction by the neutral gas, so that
microlensing experiments (Paczy\'nski 1996) already place useful
constraints on the properties of low-mass halo clouds.

Given that this hypothesised cloud population does not violate
observational constraints, the primary issues that need to be
addressed are theoretical: (i) how and when did these clouds form?
and (ii) how do they resist gravitational collapse?  The second of
these is addressed in this \emph{Letter}.

We begin by writing down equations describing a simple ``one-zone''
model of a cloud, characterised by a single temperature and
pressure (\S\ref{sec:model}), and show that particles of solid H$_2$ may
exist in the clouds (cf. Pfenniger \& Combes 1994).  At temperatures
above the microwave background temperature, these particles cool the
cloud by thermal continuum radiation, admitting equilibria in which this
cooling balances heating by cosmic rays.  In \S\ref{sec:stability}, we
show that (for optically thin emission) these equilibria are thermally
stable: if the cloud contracts the coolant is destroyed by the increase
in temperature, and the power deposited by cosmic rays causes the cloud
to expand and the temperature to return to its original value. We conclude
that, within the context of our one-zone model, the viable mass range for
Galactic clouds is $10^{-6}$--$10^{-1.7}\mathrm{M_{\sun}}$.

\section{Cloud model}
\label{sec:model}

Virial equilibrium implies that for a self-gravitating cloud
characterised by mass $M$, temperature $T$, and radius $R$,
\begin{equation}
	R \approx GM\mu/kT
	\label{eq:R}
\end{equation}
where $\mu$ is the mean molecular weight.
The pressure in the cloud can be related to the temperature upon
noting that $P\sim GM^2/R^4$, yielding
\begin{equation}
	P = \frac{q}{G^3M^2} \, \left( \frac{kT}{\mu} \right)^4 \,,
	\label{eq:P}
\end{equation}
where $q$ depends on the cloud's structure.  For polytropes $q$ rises
monotonically from 9.2 to 40 as the polytropic index runs from 3/2
to 9/2, so we adopt $q=20$.

\begin{figure*}
\centerline{\epsfxsize=8cm \epsfbox{coolants.eps}}
\figcaption{\emph{Top:} cooling by rotational transitions of H$_2$, HD
and LiH in the gas phase, and thermal continuum emission from
particles of solid H$_2$ for a $10^{-3}\mathrm{M_{\sun}}$ cloud in virial
equilibrium.  \emph{Bottom:} solid H$_2$ cooling rates for cloud
masses $10^{-3}$, $10^{-5}$ and $10^{-7}\mathrm{M_{\sun}}$.  Dotted portion of the
$10^{-7}\mathrm{M_{\sun}}$ curve indicates where the coolant is liquid H$_2$.  The
dashed line in both plots indicates the cooling rate that would give a
Kelvin-Helmholtz time scale of 10 Gyr.
\label{fig:cooling}}
\end{figure*}

At sufficiently high pressures a fraction $x$ of the molecular
hydrogen assumes solid (or liquid) form.
Then $\mu=(1-x+2y)m/(1-x+y)$, where $m$
is the mass of an H$_2$ molecule, and $y\approx1/6$ is the
abundance ratio He:H$_2$ by number.
Neglecting the temperature difference between the
phases, in equilibrium the partial pressure of H$_2$ equals the
saturated vapour pressure, i.e.
	\begin{equation}
		\frac{1-x}{1-x+y}P = \left(\frac{2\pi m}{h^2}\right)^{3/2}
		(kT)^{5/2} e^{-T_v/T} \,,
		\label{eq:Psat}
	\end{equation}
(valid for $0<x<1$), where $k T_v$ is the heat of vapourisation for
H$_2$ (Phinney 1985).  With $T_v = 91.5 \u K $, the vapour
pressure given by the RHS of eq. (\ref{eq:Psat}) is within 20 \% of
the available experimental data (Souers 1986).

Hydrogen grains can cool the gas in a manner similar to dust
grains in molecular clouds: the gas cools via collisions with slightly
colder solid particles, which in turn cool by thermal continuum
emission.  To calculate the cooling by solid H$_2$, first consider the
net power radiated by a single particle (we employ an `escape probability'
formulation of radiative transfer):
\begin{equation}
	L_s = 4\sigma\left(\frac{C(T)T^4}{1+\tau} -
   \frac{C(T_b)T_b^4}{1+\tau_b}\right) \,,
	\label{eq:L_s}
\end{equation}
where $C(T)$ is the Planck-mean absorption cross-section, $T_b$ is the
cosmic microwave background temperature, and $\tau$ and $\tau_b = \tau
C(T_b)/C(T)$ are the Planck-mean optical depth of the cloud to thermal
radiation characterized by $T$ and $T_b$ respectively.  Assuming that
the particle size is $\ll \lambda $ ($ \sim 0.1 \u cm $ at the
temperatures of relevance here), we may write $C(T) = C_m(T) m_s$,
$m_s$ being the particle mass, and for spherical grains we have
(Draine \& Lee 1984)
\begin{equation}
	C_m(T) \approx \frac{15 (4\pi)^3}{28\rho_s}
	\frac{\lambda_2}{(\varepsilon_1+2)^2}
	\left(\frac{kT}{hc}\right)^2\, ,
	\label{eq:Cm}
\end{equation}
where $\rho_s = 0.087 \u g \ut cm -3 $ is the density of
solid H$_2$  (Souers 1986),the complex dielectric function of the
solid is $\varepsilon_1 + i\varepsilon_2$, and we have assumed that
$\varepsilon_2 =
\lambda_2/\lambda$ as expected at low frequencies. The net cooling rate
per unit mass of cloud material (gas and solid) is then
\begin{equation}
	\Lambda = \frac{4\pi R^2\sigma}{M} \left(\frac{\tau T^4}{1+\tau} -
    \frac{\tau_b T_b^4}{1+\tau_b}\right) \,,
	\label{eq:Lambda}
\end{equation}
where $\tau_b = \tau\,(T_b/T)^2$, and
\begin{equation}
	\tau = C_m \frac{x}{1+2y}\, \frac{M} {\pi R^2} \,.
	\label{eq:tau}
\end{equation}

To evaluate $\Lambda$, we require optical constants for solid $\hh$ in
the microwave.  The particles are expected to be almost pure
para-hydrogen as an ortho-para mixture of the solid relaxes to para
($J=0$) form in a few days (Souers 1986).  The low frequency value of
$\varepsilon_1$ for para-hydrogen has been measured (Souers 1986) as
$\varepsilon_1\simeq1.25$; the low frequency limit of $\varepsilon_2$
is less certain.  Jochemsen et al.\ (1978) measured the extinction
coefficient of a single crystal of solid para-hydrogen in the region
of interest ($\lambda\sim0.1\;{\rm cm}$), but could not determine
whether this continuum extinction was due to absorption or scattering
within the crystal.  Because these measurements do not conform to the
anticipated low-frequency behaviour ($\propto1/\lambda$), and absorption
bands are not expected below the S(0) line, it is likely that the
absorption of pure crystalline $\hh$ is much smaller than the measured
extinction, and we can only infer a limit: $\varepsilon_2\la 1.8\times10^{-3}$.
However, the low-frequency absorption of solid H$_2$ grains could be
strongly enhanced by impurity species and lattice defects. For the
purposes of this paper we adopt $\varepsilon_2 = \lambda_2/\lambda$
and $\lambda_2 = 10^{-4} \u cm $.  Within the confines of the model,
this assumption represents one of our main areas of  uncertainty.

\begin{figure*}
	\centerline{\epsfxsize=8cm \epsfbox{tau_contours.eps}}
	\figcaption{Solid $\hh$ cooling rate versus cloud mass for
	different optical depths.  The curves, labelled by $\tau$, are
	truncated at high masses when the solutions become thermally
	unstable.  The optically thick limit, corresponding to black-body
	emission, is labelled ``$\infty$''.  The dashed curve
	delineates where half of the hydrogen is in solid form, an
	approximate upper limit for thermal stability.
	\label{fig:tau}}
\end{figure*}

For a given cloud mass, the fraction of H$_2$ in the solid phase can
be determined from $T$ using eqs (\ref{eq:P}) and (\ref{eq:Psat}).
This allows the cooling by solid particles to be calculated as a
function of $T$.  The upper panel of Fig \ref{fig:cooling} illustrates
this for a cloud of mass $10^{-3}\mathrm{M_{\sun}}$.  For comparison, the cooling
contributed by rotational lines of gas-phase H$_2$, HD and LiH is also
plotted; energies and A-values from Turner, Kirby-Docken \& Dalgarno
(1977), Abgrall, Roueff \& Viala (1982), and Gianturco et al.\ (1996).
The adopted deuterium and LiH abundances are $3\times10^{-5}$, and
$1.2 \times 10^{-10}$ respectively (Schramm \& Turner 1998).  We
employ an escape-probability formulation of radiative transfer, in
which the optically-thin cooling rates are divided by $(1+\tau)$ where
$\tau$ is the optical depth at line centre (rotational transitions) or
the Planck-mean (continuum emission).

There is a critical temperature $T_c$ at which H$_2$ in the cloud lies on
the border between the solid and gaseous phases, i.e.  the partial
pressure of H$_2$ is equal to its saturated vapour pressure.  From
this point $x$ and $\Lambda$ increase precipitously as $T$ is reduced,
until $\tau \sim 1$ and the cooling is then roughly black-body.
$\Lambda$ then increases more slowly, peaking when $x\sim 1/2$ and
subsequently dropping to zero as the temperature approaches that of
the microwave background.  At lower temperatures the cloud is heated
by the background radiation.

The solid H$_2$ cooling curves for cloud masses of $10^{-3}$,
$10^{-5}$ and $10^{-7} \mathrm{M_{\sun}}$ are compared in the lower panel of
figure \ref{fig:cooling}.  Decreasing the cloud mass has two
consequences: the critical temperature $T_c$ decreases, whereas the maximum
value of $\Lambda$ decreases (because in this circumstance the
emission is optically thick, and cloud surface area is proportional to
$M^2$ at a given temperature).  For cloud masses $\la 10^{-7} \mathrm{M_{\sun}}$,
$T_c$ is above the H$_2$ triple point (13.8\,K) and liquid droplets of
H$_2$ form instead.  This does not qualitatively affect the
calculations as the density and saturated vapour pressure of the
liquid are within 50\% of those of the solid for $T\la 20 \,K $
(Souers 1986), and the optical properties are similar (in the sense
that $\varepsilon_2$ is small and uncertain), thus in Fig.
\ref{fig:cooling} we continue the cooling curve for a $10^{-7}\mathrm{M_{\sun}}$
cloud above the triple point as a dotted curve.  Thermally stable
solutions do not exist for masses $\la 10^{-7.5}\mathrm{M_{\sun}} $, as the
partial pressure of H$_2$ exceeds the saturated vapour pressure unless
$x > 0.5$ (see \S\ref{sec:stability}).  On the other hand, for $M\ga
10^{-1.7} \mathrm{M_{\sun}}$ $T_c$ is below the CMB temperature and the solid phase
warms the cloud rather than cooling it.

\section{Thermal Equilibrium}
\label{sec:stability}

In thermal equilibrium, the cloud temperature is set by the balance
between cooling and heating.  We assume that clouds in the Galactic
halo are heated primarily by cosmic rays.  The local interstellar
cosmic-ray ionisation rate in the Galactic disc, $\sim 3\ee -17 \ut s
-1 \ut H -1 $ (Webber 1998), implies a heating rate $\Gamma \sim
3\times10^{-4} \mathrm{erg\,g^{-1}\,s^{-1}}$ (Cravens \& Dalgarno
1978).  The cosmic-ray heating in the halo is uncertain but should be
somewhat lower, say $\sim 10^{-5}\mathrm{erg\,g^{-1}\,s^{-1}}$.  In
Fig.  \ref{fig:tau}, we show the cooling rates for clouds of masses
$10^{-7}$--$10^{-1.7}\mathrm{M_{\sun}}$.  The solid curves are contours of constant
optical depth; the dashed curve shows the optically-thick limit and
represents the maximum cooling rate for each cloud mass.  It appears
that solid hydrogen can provide the necessary cooling for
planetary-mass gas clouds at the cosmic-ray heating rates expected in
the Galactic disk and halo.

The upper panel of Fig.  \ref{fig:cooling} shows that there are
typically three equilibrium temperatures available for cloud masses
between $10^{-7.5}$ and $10^{-1.7}\mathrm{M_{\sun}}$ at the expected heating rates:
solid H$_2$ provides one barely above $T_b$ and one a few degrees
higher; the gas-phase coolants provide an equilibrium above 30 K. We
now show that thermal stability requires $\Lambda$ to be a decreasing
function of $T$, and therefore only the second of these three
equilibria is stable.

In virial equilibrium the total energy per unit mass is approximately
$-\frac{3}{2}kT/\mu$ (the internal excitation of the gas is
negligible at the low temperatures of interest here), so the thermal
evolution of the cloud is determined by
\begin{equation}
	\frac{3k}{2\mu}\frac{{\rm d}T}{{\rm d}t} =\Lambda - \Gamma.
	\label{eq:energy}
\end{equation}
In the absence of heating the cloud contracts on the Kelvin-Helmholtz
time-scale $t_{\mathrm{KH}} = \frac{3}{2} kT / (\mu\Lambda) $.  Note that this
time-scale can be a substantial fraction of the Hubble time for temperatures
of a few Kelvin: the dashed curves in figure \ref{fig:cooling}
show the cooling rate that yields $t_{\mathrm{KH}} = 10 \u Gyr $.  In thermal
equilibrium cosmic-ray heating replaces the energy radiated away by
the cloud implying, for example, $t_{\mathrm{KH}} \sim2\times10^6$~yr for $\Gamma
\sim 10^{-5} \u erg \ut g -1 \ut s -1 $ (at $T\sim10$~K).
This is much greater than
the sound crossing time ($\sim10^2$~yr), so the response of a cloud
to dynamical perturbations is adiabatic, to a good approximation, and
dynamical stability is assured.  However, eq.  (\ref{eq:energy}) shows
that perturbations to the
cloud temperature grow or decay as $e^{\alpha t}$ where
\begin{equation}
	\alpha= t_{\mathrm{KH}}^{-1} \frac{T}{\Lambda}\,\frac{\rm d\;\;}{{\rm d}T}
	(\Lambda - \Gamma),
	\label{eq:growth_rate}
\end{equation}
and the right-hand side of this equation is evaluated at the
equilibrium temperature.
Thus a cloud is thermally stable only if a decrease (increase) in
cloud temperature leads to cooling outstripping (lagging) heating.
For cosmic ray heating, $\Gamma$ is independent of $T$ if the column
through the cloud is insufficient to cause significant attenuation of
cosmic rays (changes in temperature affect the cloud's column density
through the virial relationship $R\propto 1/T$).  Thermal stability
then requires that $\Lambda$ be a decreasing function of $T$, and
we conclude that {\it only the equilibrium on the high-temperature shoulder
of the solid hydrogen cooling curve is stable.\/}
In fact the column density of each cloud ($\sim10^2\;{\rm g\,cm^{-2}}$:
Walker 1999) is sufficient to stop sub-GeV cosmic-ray protons (and
all electrons), leading to a dependence of $\Gamma$ on $T$;
this dependence is too weak to affect our conclusions concerning
stability.

\section{Discussion}
\label{sec:discussion}

The suggestion that cold gas could comprise a significant fraction of
the Galaxy's dark matter is not new, previous proposals include: a
fractal medium in the the outer reaches of the Galactic disk
(Pfenniger et al.\ 1994); isolated halo clouds (Gerhard \& Silk 1996);
and mini clusters of clouds in the halo (de Paolis et al.\ 1995; Gerhard
\& Silk 1996).  However, to date there has been no compelling reason
to believe that isolated, cold gas clouds -- as inferred by Walker \&
Wardle (1998) -- could support themselves for long periods against
gravitational collapse.  We have shown that such clouds can be
stabilised by the precipitation/sublimation of particles of solid
hydrogen (or by the condensation/evaporation of droplets of liquid
hydrogen) if these particles dominate the radiative cooling of the
cloud.  The key feature which confers thermal stability is that these
particles are destroyed, hence cooling becomes less efficient, as the
cloud temperature increases.  This feature will be present in any
model where condensed hydrogen is the principal coolant, and consequently
we expect that more sophisticated structural treatments will also
admit stable solutions.

The masses of thermally stable clouds lie in the approximate range
$10^{-7.5}$--$10^{-1.7} \mathrm{M_{\sun}}$.  The lower limit is increased to
$10^{-6}\mathrm{M_{\sun}}$ if subject to cosmic ray heating similar to that in the
Galactic disc for an interval $\ga kT/\Gamma\sim10^5 \u yr $.  As halo
clouds take much longer than this to pass through the cosmic-ray disc,
this limit is appropriate even for a halo cloud population.  For cloud
masses in the range $10^{-6}$--$10^{-1.7}\;\mathrm{M_{\sun}}$ the radiative cooling
simply readjusts, on the timescale $t_{\mathrm{KH}}$, to maintain equilibrium as
$\Gamma$ varies through the orbit ($\sim10^8$~yr) of a cloud around
the Galaxy.

The typical particle radius $a$ is constrained by the requirement that
the clouds not produce significant extinction at optical wavelengths:
the geometrical optical depth of the particles, $\tau_g \approx
\tau/C_m a \rho_s $ should be less than 1.  Adopting $T=5$\,K and $\tau
= 0.01$ this translates to $a(\u mm ) \ga (10^{-4}\u cm )/\lambda_2$.
Millimeter-size particles settle to the centre of the cloud in $\sim
10^4 \u yr $, this time is shortened if $\lambda_2$ is
significantly less than $10^{-4}\u cm $ and the particles are required to be
larger.  Settling may be counteracted by convective motions or
sublimation resulting from the higher temperatures deeper within the
cloud --- issues that must await a more sophisticated treatment of the
cloud structure.

If the hypothesised population of clouds exists, their thermal microwave
emission may be detectable as a Galactic continuum background at temperatures
just above the Cosmic Microwave Background; a Galactic component of this kind
has in fact been isolated in the COBE FIRAS data (Reach et al.\ 1985).  One would
like to compare these data with the theory presented here, but it is
difficult to predict the total microwave intensity for our model because
the distribution of cosmic-ray density away from the Galactic plane is
only loosely constrained (see Webber et al.\ 1994). A similar uncertainty afflicts
the modelling of $\gamma$-ray production from baryonic material in the
Galactic halo (cf. de Paolis et al.\ 1995; Salati et al.\ 1996; de Paolis
et al.\ 1999; Kalberla, Shchekinov \& Dettmar 1999).
Nevertheless, microwave and $\gamma$-ray emissivities are each
proportional to the local cosmic-ray flux (assuming the cosmic-ray spectrum
does not vary greatly), so we can write $I_\mu\simeq I_\gamma j_\mu/
j_\gamma$, for emissivities $j$ and intensities $I$. At high latititudes
the Galactic $\gamma$-ray background is $I_\gamma\sim10^{-6}\;{\rm ph\,
cm^{-2}\,s^{-1}\,sr^{-1}}$, above 1~GeV (Dixon et al.\ 1998); local to the Sun
the corresponding (optically thin) emissivity is $1.1\times10^{-3}\;{\rm
ph\;s^{-1}\,g^{-1}\,sr^{-1}}$ (Bertsch et al.\ 1993). Thus for a cosmic-ray
heating rate (again, local to the Sun) of $\Gamma=4\pi j_\mu\sim3\times
10^{-4}\;{\rm erg\,s^{-1}\,g^{-1}}$ (Cravens \& Dalgarno 1978), we expect
$I_\mu\sim2\times10^{-8}\;{\rm erg\,cm^{-2}\,s^{-1}\,sr^{-1}}$.
This is roughly 1\% of that observed in the FIRAS cold component
at high latitude (Reach et al.\ 1985), so the microwave data do not
exclude the possibility of a cold cloud population heated by cosmic
rays.  Sciama (1999) proposed that all of the cold excess may be
accounted for by cosmic-ray heating of cold clouds, but this appears
to be based on an overestimate of the gamma-ray flux, and an
underestimate of the high-latitude FIRAS flux.

\section{Conclusions}
We have demonstrated that, by virtue of the solid/gas phase transition
of hydrogen,  cold, planetary-mass Galactic gas clouds can be thermally
stable even when they are heated at a temperature-independent rate.
Our analysis applies to the present epoch, with the microwave background
temperature at $T_b < 3$~K; for background temperatures $T_b\ga 6$~K,
our model admits no stable mass range.  Consequently the longevity of
the clouds at redshifts $z\ga 1$ is problematic. We cannot, however,
hope to address this issue until a firm theoretical basis for the
formation of such clouds has been established.

\acknowledgements

We thank Sterl Phinney for making available a copy of an unpublished
preprint and Bruce Draine for thoughtful comments on the manuscript.
The Special Research Centre for Theoretical Astrophysics is funded by
the Australian Research Council under its Special Research Centres
programme.