\section{Introduction}

Deep, near-infrared $K$-band ($2.2\mu$m) galaxy surveys are a powerful tool for
studying galaxy evolution, eg. Gardner, Cowie \& Wainscoat \shortcite{gcw93},
Cowie et al.\ \shortcite{cghshw94}, Glazebrook et al.\ \shortcite{gpmc95}.
Compared to blue-optical light, near-infrared light is a better tracer
of mass in evolved stars and the correction for redshift
dimming (the ``$k$-correction'') is approximately independent of
morphological type.
The rapid evolution in galaxy number counts apparent in the $b_J$ band is not
seen in the $K$ band, eg. Koo \& Kron \shortcite{kk92}.
However, it is vital to have a reliable determination of the $K$-band
luminosity function (LF) for {\em nearby} galaxies in order to interpret
faint galaxy counts and to calculate the clustering of $K$-selected galaxy
samples.
Local $K$-band luminosity functions have been measured from optically-selected
samples by Mobasher, Sharples \& Ellis \shortcite{mse93} and
Szokoly et al. \shortcite{sscm98}
and from $K$-selected samples by Glazebrook et al.\ \shortcite{gpmc95} and
Gardner et al. \shortcite{gsfc97}.
Since all of these surveys are flux-limited, the majority of galaxies have
$K$-band luminosities close to $L_K^*$.
Even in the largest sample (Gardner et al., 510 galaxies), it was
feasible to measure the LF over a range
of only 5 magnitudes, ${M^*}^{+3}_{-2}$, and so the the faint-end slope of
the $K$-band luminosity function, so important for predicting galaxy number
counts, is still rather poorly constrained.

In this paper we present a new estimate of the $K$-band LF over a range
of ten magnitudes, ${M^*}^{+8}_{-2}$, based on a subsample of galaxies
selected from the Stromlo-APM galaxy survey \cite{lpme96}.
This survey is an ideal source
for estimating the $K$-band LF since redshifts have already
been measured for 1797 galaxies with $b_J < 17.15$
over a very large volume of space.
The solid angle of the survey is 1.3 sr and the median redshift is
about 15,300 km/s.
The key to measuring the faint-end of the $K$-band LF without a huge
investment of telescope time is to observe galaxies selected by their
intrinsic luminosity rather than their apparent flux.
One can make use of the fact that near-infrared and optical
luminosities are correlated \cite{mes86,sci96},
in order to preferentially select galaxies of high and low luminosity and
thus sample the luminosity range more uniformly than a flux-limited sample.
One is thus able to measure the luminosity function to fainter
luminosities than from a flux-limited sample of similar size.
This sampling strategy is described in \S\ref{sec:sample} and
the observations and data reduction are discussed in \S\ref{sec:obs}.
Our method of estimating $\phi(M_K)$ from a $b_J$-selected sample is given
in \S\ref{sec:method} and we test this method in \S\ref{sec:test}.
Our results are presented in \S\ref{sec:results} and we conclude in
\S\ref{sec:concs}.
For notational convenience, we will denote absolute magnitudes
in the $K$ and $b_J$ bands by $M_K$ and $M_B$ respectively.
Throughout, we assume a Hubble constant of $H_0 = 100$ km/s/Mpc.

\section{Sample Selection} \label{sec:sample}

The aim in selecting a subset of Stromlo-APM galaxies for which to obtain
$K$-band photometry was to sample the magnitude range
$-22 \le M_B \le -13$ (the full range of $b_J$ absolute magnitudes
in the Stromlo-APM survey) as uniformly as possible.
An added complication in defining the sample arose because we wished
to obtain optical CCD images for the same sample of galaxies.
One planned use of this optical imaging is to measure morphological
parameters for a representative sample of galaxies at low redshift
in order to compare with HST observations of galaxies
at high redshift, $z > 0.4$, \cite{brinch98}.
To obtain comparable linear resolution to the HST data required observing
galaxies at $z < 0.04$, assuming ground-based seeing of 1.3 arcsecond.
Our primary sample thus consists of galaxies at redshifts $z < 0.04$.
We divided the magnitude interval $-22 \le M_B \le -13$ into 90 bins
each of width 0.1 mag.
We then randomly selected up to six galaxies from the Stromlo-APM survey with
$z < 0.04$ in each bin.
Due to the redshift limit of $z < 0.04$, this primary sample
contains rather few galaxies brighter than $M_B = -20$.
We therefore formed a supplementary sample, consisting of galaxies at
$z > 0.04$ to top up each magnitude bin, where possible,
to six galaxies.
This supplementary sample consists entirely of galaxies with $M_B < -20$.
The primary sample contains 283 galaxies, and the supplementary sample
contains 80 galaxies, giving a total sample size of 363 galaxies.

\section{Observations and data reduction} \label{sec:obs}

Imaging of the above sample of galaxies was carried out at the
Cerro Tololo Interamerican Observatory (CTIO) 1.5m telescope using
the CIRIM infrared array through the standard $K$ filter
over the nine nights 1996 August 31 -- September 4 and 1997 October 19--22.
The pixel size at $f/7.5$ is $1.16''$, allowing most galaxies
to be observed at 9 non-overlapping positions on the $256 \times 256$ array.
Two frames were taken with the galaxy at the central position,
thus yielding 10 frames per galaxy.
Total integration time for each galaxy was 300 seconds.
For 11 galaxies with angular size more than 100 arcsec,
we obtained four on-source and four offset-sky integrations of 75 seconds each.
In the following, these will be referred to as ``biggrid'' observations.
Standard stars were observed from the list of Elias et al.\ \shortcite{efmn82}.
Dark frames and dome flats were taken at the start of each night.

The infrafred frames were reduced using {\sc iraf},
mosaiced with the {\sc dimsum}\footnote{{\sc dimsum} is the Deep Infrared
Mosaicing Software package developed by
Peter Eisenhardt, Mark Dickinson, Adam Stanford, and John Ward, and is
available via ftp from
ftp://iraf.noao.edu/iraf/contrib/dimsumV2/dimsum.tar.Z} package and
image detection and photometry was performed using SExtractor \cite{ba96}.

\subsection{Basic reduction}

The basic reduction process consisted of the following steps:
\begin{enumerate}
\item Non-linearity correction using the {\sc irlincor} task with parameters
$c1 = 0.9997$, $c2 = 0.0257$, $c3 = 0.0158$ (Mike Keane, private communication),
\item Subtraction of dark frame,
\item Flatfielding by dome flat,
\item Masking of bad pixels.
\end{enumerate}

\subsection{Mosaicing}

At this stage of the reduction we had ten frames per galaxy, with the galaxy
at a different position on each frame.
Mosaicing of the frames and subtraction of sky background was performed
using the {\sc dimsum} package.
As described by Stanford, Eisenhardt \& Dickinson \shortcite{sed95},
{\sc dimsum} employs a two-pass procedure in order to mask out faint as
well as bright images when constructing the sky background,
resulting in a much flatter sky than is obtained from median filtering of
the galaxy frames with outlier rejection.
Individual galaxy frames were block replicated by a factor
of 4 in each dimension, allowing alignment to be performed by integer offsets
without interpolation.
Co-added images and corresponding exposure weight maps were made by
summing the aligned images and were finally block averaged $2 \times 2$
in order to economize on disc space.

For the 11 ``biggrid'' observations, the four offset sky frames were median
filtered with outlier rejection in order to estimate the sky background.
The sky background was subtracted from each of the four on-source frames,
which were then aligned and coadded.

\subsection{Image detection and photometry}

We used SExtractor 2.0.15 to detect and measure images in the mosaiced frames.
For both standard stars and galaxies we used the {\sc mag\_best} estimate
of magnitude.
This yields a pseudo-total magnitude \cite{kron80} except in crowded fields,
when a corrected isophotal magnitude is measured instead.
Magnitude errors were estimated by combining in quadrature SExtractor's
estimate of the error from photon statistics and the difference between
magnitudes measured using local and global estimates of the sky background.
Of the selected sample of 363 galaxies, 351 were observed under photometric
conditions and 345 yielded a
$K$-band magnitude with an estimated error of less than 0.1 mag
(rms mag error = 0.04 mag).

\subsection{Calibration}

Since we observed only single band infrared imaging, we used the following
simple relation to convert observed magnitudes $k$ to standard CIT
\cite{efmn82} magnitudes $K$:
\begin{equation}
K = k + k_0 + k_X X, \label{eqn:trans}
\end{equation}
where $k_0$ is a zero-point offset, $k_X$ is the extinction coefficient in
the $K$-band and $X$ is the airmass of the observation.
We made a total of 74 standard star observations during the runs
(an average of 8 per night) and we initially
fitted the parameters $k_0$ and $k_X$ for all standard star observations
combined, obtaining $k_0 = 7.588 \pm 0.020$ and $k_X = 0.081 \pm 0.018$
with rms residual magnitude error of 0.018.
Six standard star observations with large residuals were omitted from the
fitting procedure, three of these were from the night of 1996 August 31,
which was partially non-photometric.
Holding the zero-point term fixed at $k_0 = 7.588$ we then fitted the
extinction coefficient separately for each night, with results shown
in Table~\ref{tab:calib}.
Galaxy magnitudes were converted to the CIT system using (\ref{eqn:trans})
with $k_0 = 7.588$ and $k_X$ coefficients from the Table.
Galaxies observed during non-photometric conditions on the nights of
1996 August 31 and 1997 October 20 were rejected: these cases were obvious
from the rapidly varying sky background.

\begin{table}
 \caption{Standard star calibration data.}
 \label{tab:calib}
 \begin{tabular}{lrrcc}
 \hline
 \hline
 Night & $N_{\rm std}^1$ & $N_{\rm del}^2$ & $k_X^3$ & rms\\
 \hline
1996 Aug 31 & 11 & 3 & $0.105 \pm 0.006$ & 0.017\\
1996 Sep 01 &  4 & 1 & $0.101 \pm 0.002$ & 0.003\\
1996 Sep 02 &  9 & 0 & $0.081 \pm 0.004$ & 0.014\\
1996 Sep 03 &  9 & 0 & $0.082 \pm 0.004$ & 0.015\\
1996 Sep 04 &  7 & 0 & $0.092 \pm 0.004$ & 0.011\\
1997 Oct 19 &  9 & 0 & $0.083 \pm 0.003$ & 0.010\\
1997 Oct 20 &  7 & 1 & $0.067 \pm 0.004$ & 0.010\\
1997 Oct 21 & 11 & 1 & $0.075 \pm 0.002$ & 0.007\\
1997 Oct 22 &  7 & 0 & $0.066 \pm 0.005$ & 0.013\\
  \hline
  \hline
 \end{tabular}

 \medskip
 $^1$ Number of standard stars observed.\\
 $^2$ Number of outliers deleted.\\
 $^3$ Extinction coefficient in (\ref{eqn:trans}).
\end{table}

\subsection{Photometric repeatability}

\begin{figure}
\centerline{\psfig{figure=repeats.eps,width=\linewidth}}
\caption{Repeatability of galaxy photometry.
The top panel plots the rms $K$ magnitude between repeated observations
against $K$ magnitude from the coadded data.
The symbols indicate the rms magnitude error and the length of
the line shows the estimated magnitude error.
The lower panel plots the same information but as a function of mean
$K$-band surface brightness $\mu_K$ within the measurement aperture.
Open symbols denote sources deblended by SExtractor, stars indicate a
``biggrid'' observation.
\label{fig:repeats}}
\end{figure}

Thirty-one galaxies, mostly of low $K$-band surface brightness, were observed
on more than one occasion, which allows us to assess the repeatability of
our photometry.
In these cases, we ran SExtractor on the mosaiced galaxy frames both before
and after coadding observations.
We obtain final photometry from the coadded frame and use
the two individual frames and the coadded frame to estimate the rms error.
In four cases, one observation was made using ``biggrid'': in these cases
we do not coadd the observations and instead use the {\sc dimsum}-reduced
observation for the final magnitude and estimated error.

In Figure~\ref{fig:repeats} we plot both our estimated errors and the rms
errors between repeated observations as a function of $K$ magnitude and mean
surface brightness $\mu_K$.
The symbols indicate the rms magnitude error and the length of
the line shows the estimated magnitude error from the coadded image.
Thus if our estimated errors are a good estimate of the true rms then the lower
ends of the error bars should reach zero.
In most cases, the rms error is within 2--3 times the estimated error
and in a few cases the error is slightly overestimated.
The five points with rms $> 0.4$ all correspond to images which were broken
up by SExtractor's deblender (open symbols).
We attempted to sum the flux from the components, but clearly the
$K$ magnitude for these objects is accurate to only $\sim 0.5$ mag.
We see that most of the discrepant points are of extremely low surface
brightness, $\mu_K \ga 21.5$ mag arcsec$^{-2}$.
For unbroken images, all rms errors are less than 0.2 mag and the estimated
errors provide a reasonable estimate of the rms.
The stars denote objects in which one of the observations used separate
sky exposures (``biggrid'' observations).
For two out of four of these objects, we see an rms magnitude error larger
than 0.1 but with a negligible estimated error (from the non-biggrid
observation).
The reason for this is the poorer sky subtraction of the biggrid observations
compared with the {\sc dimsum} reduction of the majority of our observations.

\subsection{Matching to APM}

Since the mosaiced frames cover an area of sky around 8 arcmin on a side,
albeit not to uniform depth, they provide $K$-band photometry for
many objects in each field in addition to the target galaxy.
We therefore matched the images detected by SExtractor with images in the
APM scans.
Using the matched objects in each frame, we calculated a 6-parameter transform
from CIRIM pixel coordinates to APM plate coordinates and thence to RA \& Dec.

The $K$-band photometry for the target galaxies is presented in the Appendix
to this paper.

\subsection{Sample properties}

\begin{figure}
\centerline{\psfig{figure=ksel.eps,width=\linewidth}}
\caption{Illustration of observational selection effects in our sample.
The upper panel plots observed $b_j - K$ colour versus apparent $K$ magnitude.
The line shows the apparent magnitude limit for the Stromlo-APM survey,
$b_J = 17.15$.
The lower panel plots redshift against apparent $K$ magnitude.
The line in this case separates the primary and supplementary samples
at $z = 0.04$.}
\label{fig:ksel}
\end{figure}

Since our source survey is limited by $b_J$ flux, our $K$-band sample is
not a complete one, as illustrated in Figure~\ref{fig:ksel}.
At faint $K$ magnitudes, only blue objects will be in the Stromlo-APM sample.
We are roughly complete to $K \approx 12$, the bluest galaxies can be seen
as faint as $K \approx 15$.
The luminosity function estimator described in \S\ref{sec:method} corrects
for this $K$-band incompleteness and allows us to use all galaxies with
$K$-band photometry.
In the lower panel of this Figure we plot the redshift distribution as
a function of apparent $K$ magnitude.
The requirement that $z < 0.04$ for the primary sample has a noticeable
effect on the overall redshift distribution: this effect is also accounted
for by our luminosity function estimator.

\begin{figure}
\centerline{\psfig{figure=colmag.eps,width=\linewidth,angle=-90}}
\caption{Rest-frame $M_B - M_K$ versus $M_k$ colour-magnitude plot.
	Plus signs represent early-type galaxies, asterisks late-type
	galaxies and dots represent unclassified galaxies.
	The solid line shows a least-squares fit to all galaxies,
	the dashed line a fit to early types and the dot-dashed line a fit
	to late type galaxies.}
\label{fig:colmag}
\end{figure}

To calculate absolute magnitudes we use galactocentric recession velocities
and assume k-corrections in the $b_J$ band as given in Table~2 of
Efstathiou, Ellis \& Peterson \shortcite{eep88}.
We assume a $K$-band k-correction of $-2.5 z$ for all galaxy types
(ie. the near-infrared k-correction {\em brightens} galaxies with redshift).
This is a very good approximation to the $K$-band k-correction of
Glazebrook et al.\ \shortcite{gpmc95} for redshifts $z < 0.15$.

Figure~\ref{fig:colmag} shows the rest-frame $(M_B - M_K)$ versus $M_K$
colour-magnitude relation for our data.
The fit to all galaxies is given by
\begin{equation}
(M_B - M_K) = -0.212 \times M_K - 1.20,\ \sigma = 0.77. \label{eqn:colmag}
\end{equation}
The slopes in this relation for early and late-type galaxies
($-0.135 \pm 0.036$ and $-0.209 \pm 0.024$ respectively) are consistent
with those found for E-S0 ($-0.095 \pm 0.013$) and Sa-Sdm ($-0.24 \pm 0.03$)
galaxies by Mobasher, Ellis \& Sharples \shortcite{mes86}, although we find a
significantly larger scatter about the relation.
At least part of this scatter is due to saturated APM $b_J$ magnitudes
for galaxies $b_J \la 15$ and a few instances of the APM scans breaking up a
very large, bright galaxy into fragments.
An improved colour-magnitude relation will be available when we incorporate
optical CCD magnitudes for these galaxies (Loveday \& Lilly, in preparation).

\section{Estimating the $K$-Band Luminosity Function}
\label{sec:method}

When one has a sample selected on optical $b_J$ magnitude and wishes to
estimate the $K$-band luminosity function,
the best way to proceed is to calculate a bivariate luminosity function
(BLF, $\phi(M_K, M_B)$) allowing for known selection in $b_J$ flux and $M_B$ absolute
magnitude and then to integrate over $M_B$ to obtain $\phi(M_K)$.
One can estimate the shape of $\phi(M_K, M_B)$,
independently of inhomogeneities in the galaxy distribution, using
the maximum likelihood method of Sandage, Tammann \& Yahil \shortcite{sty79}.
The probability of seeing a galaxy with $K$-band luminosity $L_K^i$
and $B$ band luminosity $L_B^i$ at redshift $z_i$ in our sample is given by
\begin{equation}
p_i = \frac{\phi(L_K^i, L_B^i) S(L_B^i)}
      {\int_{{L_K}_{\rm min}(z_i)}^{{L_K}_{\rm max}(z_i)}
      \int_{{L_B}_{\rm min}(z_i)}^{{L_B}_{\rm max}(z_i)}
      \phi(L_K, L_B) S(L_B) dL_K dL_B}.
\end{equation}
The function $S(L_B)$ accounts for the known selection in absolute $B$
magnitude and the luminosity limits ${L_B}_{\rm min}(z_i)$ and
${L_B}_{\rm max}(z_i)$
are the minimum and maximum $B$-band luminosities observable at
redshift $z_i$ in a sample limited by apparent $b$ magnitude.
For the sample analysed here, there are no flux limits in the $k$-band,
and so the integral over $K$-band luminosity runs from
0 to $+\infty$.
The maximum-likelihood shape of the BLF $\phi(M_K, M_B)$ is estimated by maximizing
the likelihood ${\cal L} = \prod_{i=1}^{N_g} p_i$ (the product of the
individual probabilities $p_i$ for the $N_g$ galaxies in the sample)
with respect to the parameters describing the BLF.

In practice, we do not have a good {\em a priori} parametric model for
$\phi(M_K, M_B)$, and so instead we measure $\phi(M_K, M_B)$\
in a non-parametric way using an extension of the Efstathiou, Ellis \& Peterson
\shortcite{eep88} stepwise maximum likelihood (SWML) method.
Sodr\'{e} \& Lahav \shortcite{sl93} have extended the SWML method to estimate
the bivariate diameter-luminosity function and to allow for sample
incompleteness.
We adopt their extension of the SWML estimator here, including the sampling
function $S(M_B)$ separately for the primary and supplementary galaxy samples.
These sampling functions are illustrated in Figure~\ref{fig:samp}, which shows
histograms of $M_B$ for the entire Stromlo-APM sample, for the selected
galaxies with good $K$-band photometry and the sampling functions $S(M_B)$.

\begin{figure}
\centerline{\psfig{figure=samp.eps,width=\linewidth}}
\caption{Sample selection as a function of $M_B$ absolute magnitude for the
primary (thick lines) and supplementary (thin lines) samples.
The upper panel shows histograms of $M_B$ for all 1797 galaxies in the
Stromlo-APM survey.
The middle panel shows the same for the 345 selected galaxies with $K$-band
photometry.
The lower panel shows the sampling functions $S(M_B)$, the ratio of selected
to total galaxies in each bin.}
\label{fig:samp}
\end{figure}

We measure the bivariate luminosity function in bins of absolute
$K$ and $B$ magnitude,
\begin{equation}
\phi(M_K, M_B) = \phi_{jk},\ \ j = 1, \ldots, N_{M_K},\ \ k = 1, \ldots, N_{M_B},
\end{equation}
where $M_K^j - \Delta M_K/2 < M_K < M_K^j + \Delta M_K/2$ and
$M_B^k - \Delta M_B/2 < M_B < M_B^k + \Delta M_B/2$.

The log-likelihood is given by
\begin{eqnarray}
\ln{\cal L} & = & \sum_{i=1}^{N_g} \sum_{j=1}^{N_{M_K}} \sum_{k=1}^{N_{M_B}}
	      W_{ijk} \ln[\phi_{jk} S(M_B^k)] - \nonumber \\
	      & & \sum_{i=1}^{N_g} \ln\left[
	      \sum_{j=1}^{N_{M_K}} \sum_{k=1}^{N_{M_B}}
	      H_{ijk} \phi_{jk} \right] + {\rm const}.
\label{eqn:logL}
\end{eqnarray}
Here $W_{ijk} \equiv W(M_K^i - M_K^j,\ M_B^i - M_B^k)$, with,
\begin{equation}
W(x,y) = \left\{ \begin{array}{ll}
		1 & \mbox{if} -\Delta M_K/2 \le x \le \Delta M_K/2\\
		  & \mbox{and} -\Delta M_B/2 \le y \le \Delta M_B/2,\\
			0 & \mbox{otherwise.}
		    \end{array}
	    \right.
\end{equation}
The sampling function $S(M_B)$ (Fig.~\ref{fig:samp})
is incorporated into the ramp function $H$.
Writing $M_K^- = M_K^j - \Delta M_K/2$,\ $M_K^+ = M_K^j + \Delta M_K/2$,\
$M_B^- = M_B^k - \Delta M_B/2$,\ and $M_B^+ = M_B^k + \Delta M_B/2$, then
\begin{equation}
H_{ijk} = \frac{1}{\Delta M_K \Delta M_B} \int_{M_K^-}^{M_K'} dM_K
	  \int_{M_B^-}^{M_B'} S(M_B) dM_B,
\end{equation}
where $M_K' = \max[M_K^-, \min(M_K^+, {M_K}_{\rm lim}^i)]$ and
$M_B' = \max[M_B^-, \min(M_B^+, {M_B}_{\rm lim}^i)]$.

To fix the otherwise arbitrary normalisation constant we apply the constraint
\begin{equation}
g = \sum_j \sum_k \phi_{jk} \left(\frac{L_K}{L_f}\frac{L_B}{L_f}\right)^\beta
    \Delta M_K \Delta M_B - 1 = 0
\end{equation}
with $\beta = 1.5$ and the fiducial luminosity $L_f$ corresponding to
$M = -20$ using a Lagrangian multiplier $\lambda$.
The new likelihood $\ln {\cal L'} = \ln {\cal L} + \lambda g$ is maximised
with respect to the $\phi_{jk}$ and $\lambda$, requiring that $\lambda = 0$.

Setting $\partial\ln{\cal L'}/\partial\phi_{jk} = 0$, one arrives
at a maximum likelihood estimate for $\phi_{jk}$,
\begin{equation}
\phi_{jk} = n_{jk}\left/\sum_i^{N_g}\left[
	    \frac{H_{ijk} \Delta M_K \Delta M_B}
	    {\left(\sum_{l=1}^{N_{M_K}} \sum_{m=1}^{N_{M_B}} \phi_{lm} H_{ilm}
	    \Delta M_K \Delta M_B \right)} \right] \right.,
\label{eqn:phiKB}
\end{equation}
where $n_{jk} = \sum_{i=1}^{N_g} W_{ijk}$ is the number of galaxies
in the $(j,k)th$ bin.
Errors in $\phi_{jk}$ are estimated via the inverse of the information matrix
(Efstathiou et al.\ 1988).

As with all density-independent estimators, information about the overall
normalisation is lost.
We therefore normalise our BLF to the mean density of galaxies with
$-22 \le M_B \le -13$ in the full Stromlo-APM sample,
$\bar{n} = 0.071 h^3{\rm Mpc}^{-3}$, calculated as described by Loveday et al.\
\shortcite{lpem92}.
(Note that the density $\bar{n} = 0.047 h^3{\rm Mpc}^{-3}$ quoted by
Loveday et al.\ (1992) is for the restricted magnitude range
$-22 \le M_B \le -15$.)

Finally, we obtain the $K$-band luminosity function by summing over
$B$ magnitude bins,
\begin{equation}
\phi(M_K^j) = \sum_{k=1}^{N_{M_B}} \phi_{jk}.
\label{eqn:phiK}
\end{equation}
One can then fit a given functional form, eg. a Schechter \shortcite{schec76}
function, to the stepwise estimate by least-squares.

\section{Test of the Method} \label{sec:test}

We have tested the above procedure by using it to
estimate the $K$-band luminosity function from a set of Monte Carlo
simulations.
We generated nine mock Stromlo surveys by a Soneira \& Peebles \shortcite{sp78}
hierarchical clustering simulation with similar clustering properties to that
measured from the Stromlo-APM Survey \cite{lmep95}.
Each galaxy in the simulation was assigned a $K$-band luminosity
drawn at random from
a Schechter function with $\alpha = -1.20$ and $M_K^* = -23.6$
and then assigned a $b_J$ magnitude according to our observed
colour-luminosity relation (\ref{eqn:colmag}).
Galaxies were selected on their apparent $b_J$ magnitude, $b_J < 17.15$.
This process was repeated until each simulation contained 2000 galaxies.
We then sampled each simulation by absolute $M_B$ magnitude as described
in \S\ref{sec:sample}, finally yielding an average of 359 galaxies
per simulation.
We calculated the K-band luminosity function $\phi(M_K)$ for each
simulation as described in \S\ref{sec:method} and fit a Schechter function
to each by least squares.
Averaging over the nine simulations, and estimating the BLF in bins of width
0.5 mag, we measure mean and rms Schechter
function parameters $\alpha = -1.17 \pm 0.07$, $M^* = -23.6 \pm 0.2$.
The errors on the mean values are $\sqrt 9$ times smaller than the quoted
rms scatter between the simulations.
Our estimates of $\alpha$ and $M^*$ are in excellent agreement with the
input luminosity function.

\begin{figure}
\centerline{\psfig{figure=test.eps,width=\linewidth}}
\caption{Test of our method to estimate the $K$-band luminosity function
from $b_J$-limited samples.
The upper panels show the $K$-band LF estimated from each of nine
Soneira-Peebles simulations and the lower panel shows the mean and rms
(error bars offset to right).
In each case the curve shows the {\em input} LF.
\label{fig:test}}
\end{figure}

The SWML estimates of the $K$-band LF for the simulations
are shown in Figure~\ref{fig:test}.
The points show the SWML estimate of $\phi(M_K)$ from
the nine simulations and the curve shows the input Schechter function with
shape $\alpha = -1.20$, $M^*_K = -23.6$.
The error bars going through the data points show the errors determined
from the covariance matrix.
The lower panel shows the average over the 9 simulations.
The rms scatter between realisations for each data point (shown by
the error bars offset slightly to the right) are in satisfactory agreement
with the predicted errors.

Overall, we find that our procedure for estimating $\phi(M_K)$ from a
sample limited by apparent $b$ magnitude and further selected by absolute
$B$ magnitude provides a robust and unbiased estimate of the
$K$-band LF over a wide range of absolute magnitudes.

\section{Results} \label{sec:results}

\begin{figure*}
\centerline{\psfig{figure=klf.eps,width=12cm,angle=-90}}
\caption{The $K$-band luminosity function estimated from our sample
(symbols) together with the best-fit Schechter function.
The inset shows the 1 and 2 $\sigma$ likelihood contours for the shape
parameters $\alpha$ and $M_K^*$.}
\label{fig:klf}
\end{figure*}

Our estimated $K$-band luminosity function is shown in
Figure~\ref{fig:klf}.
The curve shows a Schechter function fit to the SWML estimate using least
squares.
We allow for finite bin width by calculating the mean predicted
$\langle \phi_j \rangle$ at the absolute magnitude of each galaxy in each bin,
rather than simply calculating $\phi(M)$ at the bin centre.
The inset shows 1 and 2 $\sigma$ likelihood contours for the shape parameters
$\alpha$ and $M_K^*$, where the normalisation $\phi^*$ is adjusted to maximize
the likelihood at each grid point.
The best-fit Schechter parameters are $\alpha = -1.16 \pm 0.19$,
$M_K^* = -23.58 \pm 0.42$ and
$\phi^* = 0.0121  \pm 0.0082 h^3 {\rm Mpc}^{-3}$.
The quoted errors on $\alpha$ and $M_K^*$ come from the bounding-box of the
1$\sigma$ likelihood contour.

Clearly the characteristic magnitude $M_K^*$ is rather poorly constrained
from these data, and so in order to measure the faint-end slope $\alpha$
independent of $M_K^*$, we have also fit a straight line to our estimated LF
over the restricted luminosity range $M > -22$.
We measure a slope $-1.08 \pm 0.15$, consistent with that measured
from the Schechter fit.
Furthermore, this slope does not change significantly varying the bright
luminosity cut from $-20.5$ to $-23$.

Since we normalised the LF to the same number density as the full Stromlo-APM
sample, the error in $\phi^*$
is dominated by the uncertainty in shape of the LF.
The Schechter function provides a reasonable fit to the SWML estimate over the
full range of ten magnitudes; there is no indication of any faint-end turnup.

\begin{figure*}
\centerline{\psfig{figure=comp.eps,width=15cm,angle=-90}}
\caption{Comparison of our new estimate of the $K$-band luminosity function
(dotted line) with estimates from (a) Mobasher et al., (b) Glazebrook et al.,
(c) Gardner et al.\ and (d) Szokoly et al.}
\label{fig:comp}
\end{figure*}

\begin{table*}
\begin{minipage}{12cm}
 \begin{center}
 \caption{Schechter function fits to $K$-band LFs.}
 \label{tab:comp}
 \begin{math}
 \begin{array}{lccc}
 \hline
 \hline
 {\rm Sample} & \alpha & M^*_K & \phi^* / h^{-3}{\rm Mpc}^3\\
 \hline
\mbox{Mobasher et al.\ 1993} & -1.0 \pm 0.3 & -23.4 \pm 0.3 & 0.0112 \pm 0.0016\\
\mbox{Glazebrook et al.\ 1995 ($z < 0.2$)} & -1.04 \pm 0.31 & -23.02 \pm 0.23 & 0.029 \pm 0.007\\
\mbox{Gardner et al.\ 1997} & -0.91 \pm 0.24 & -23.12 \pm 0.17 & 0.0166\\
\mbox{Szokoly et al.\ 1998} & -1.3 \pm 0.2 & -23.6 \pm 0.3 & 0.012 \pm 0.004\\
\mbox{This work (all galaxies)} & -1.16 \pm 0.19 & -23.58 \pm 0.42 & 0.012 \pm 0.008\\
\mbox{This work (non-ELG)} & -1.3 \pm 0.4 & -23.5 \pm 0.7 & 0.001 \pm 0.001\\
\mbox{This work (ELG)} & -1.2 \pm 0.4 & -22.5 \pm 0.8 & 0.011 \pm 0.013\\
  \hline
  \hline
 \end{array}
 \end{math}
\end{center}
\end{minipage}
\end{table*}

In Figure~\ref{fig:comp} and Table~\ref{tab:comp}, we compare our new estimate
of the $K$-band luminosity function with those of other workers.
Following Glazebrook et al.\ (1995) and Gardner et al.\ (1997), we have added
0.22 mag to the points from Figure~2 of Mobasher et al.\ (1993) to account for
their method of calculating K-corrections.
We have also applied an aperture correction of $-0.30$ mag to the
Glazebrook et al.\ data (see Gardner et al.\ 1997).
Ours is the first estimate of $\phi(M_K)$ fainter than $M_K = -20$.
Brighter than this all estimates are in reasonable agreement,
particularly given the uncertainties in the normalisation of the LF from
small samples of galaxies.
The normalisation of the Glazebrook et al.\ (1995) sample is about twice
that of all the other samples.
This is likely to be due to sampling fluctuations, since their survey covers
only 552 arcmin$^2$ and contains a total of only 124 galaxies.

Note that the bright-end points of our estimated LF have dropped significantly
since the preliminary analysis of Loveday \shortcite{love98car}.
This is largely due to the improved K-band photometry obtained by mosaiccing
the galaxy frames with {\sc dimsum} compared with a cruder algorithm used
earlier.

\subsection{ELG versus non-ELG}

As an aid to understanding the processes of galaxy formation and
evolution, it is of great help
to have estimates of the luminosity function for different types of galaxies.
The most objective way of separating Stromlo-APM galaxies is via the equivalent
width of the \mbox{H$\alpha$}\ emission line \cite[hereafter LTM]{ltm99}.
Due to the small size of our $K$-band sample, we subdivide into just two
subsamples, those with $\mbox{EW\,(H$\alpha$)} < 10$\AA\ (non-ELG, 138 galaxies) and those with
$\mbox{EW\,(H$\alpha$)} \ge 10$\AA\ (ELG, 134 galaxies).
For the remaining 71 galaxies, no \mbox{EW\,(H$\alpha$)}\ measurement is available (LTM).
For both subsamples we recalculated the sampling function $S(M_B)$ appropriate
for galaxies with \mbox{EW\,(H$\alpha$)}\ less or greater than 10\AA\ as appropriate.
The LFs were normalised to mean densities of $\bar{n} = 0.018\ h^{3}{\rm Mpc}^{-3}$
(non-ELG) and $\bar{n} = 0.098\ h^{3}{\rm Mpc}^{-3}$ (ELG) as determined for galaxies
with these emission line properties from the full Stromlo-APM sample.
The fact that the estimated mean density of ELGs is larger than that for
galaxies of any \mbox{EW\,(H$\alpha$)}\ ($\bar{n} = 0.071 h^{3}{\rm Mpc}^{-3}$, cf. \S\ref{sec:method})
is a consequence of the mean density of subsamples of Stromlo-APM galaxies
being determined to only $\sim 30\%$ accuracy (LTM).

\begin{figure}
\centerline{\psfig{figure=ew.eps,width=\linewidth,angle=-90}}
\caption{Comparison of the $K$-band luminosity function estimated from our
non-ELG sample (filled symbols, continuous line) and ELG sample (open symbols,
dashed line).}
\label{fig:ew}
\end{figure}
\begin{figure}
\centerline{\psfig{figure=ew_like.eps,width=\linewidth,angle=-90}}
\caption{$1\sigma$ likelihood contours for the best fit Schechter parameters
for the non-ELG and ELG samples.}
\label{fig:ew_like}
\end{figure}

The estimated LFs for each subsample are shown in Figure~\ref{fig:ew}
and the $1\sigma$ likelihood contours for the best fit Schechter parameters
are shown in Figure~\ref{fig:ew_like}.
The Schechter function parameters themselves are listed in Table~\ref{tab:comp}.
Compared with the optical luminosity functions for \mbox{EW\,(H$\alpha$)}-selected samples (LTM),
the $K$-band LFs for ELGs and non-ELGs are very similar in shape,
with just a $\sim 1\sigma$ indication that ELGs have characteristic magnitudes
$M^*$ about one magnitude fainter than non-ELGs.
Note that the faint-end slopes for the two samples are completely consistent.
These results are in accord with the LFs estimated for E/S0 and spiral
galaxies by Mobasher et al.\ (1993).

\section{Conclusions} \label{sec:concs}

We have presented a new estimate of the $K$-band luminosity
function obtained from infrared array imaging of a subsample of galaxies
from the Stromlo-APM survey.
We measure $\phi(M_K)$ over a range of 10 magnitudes, ${M^*}^{+8}_{-2}$,
a significantly greater range of luminosities than has been measured before now.
Our LF is consistent with earlier estimates, and shows that a
Schechter function of moderate faint-end slope ($\alpha = -1.16$) provides
a good fit to the LF as faint as $M^* + 8$ with no hint of any upturn
at the faint-end.
This observation appears to rule out the suggestion of Glazebrook et al.\
(1995) that there may be an extra dwarf component to the $K$-band LF, or
at least that if there is such a component, it must have characteristic
magnitude $M_{\rm dwarf} \ga M^*_K + 8$.
In contrast, the faint-end upturn seen in the field $b_J$ LF \cite{love97}
is apparent by $M^*_B + 6$, and in the CfA Redshift Survey \cite{mhg94} by
$M^*_{Zw} + 3$.
In the Coma cluster, an upturn in the $H$-band
($1.6\mu m$) LF is also seen at $M^*_H + 3$ \cite{desd98}.
Note that we do see a slight upturn in $\phi(M_K)$ near $M^*_K + 3$,
faintward of the three low points around $M_K = -22$, but the LF flattens
out again by $M_K = -20$.

The $K$-band LFs for non-ELGs and ELGs are very similar in shape, with no
difference in faint-end slopes, contrary to what is seen in the optical LF
(LTM).
The steep optical LF for ELGs might be consistent with a $\sim$ flat mass
function in a fading burst model \cite{hp97},
although the particular model that Hogg \& Phinney consider
would steepen the $K$-band LF even more than the optical LF.

The main limitation with our sample is that the $b_J$ selection of the
source survey introduces a bias against red galaxies compared to a
$K$-selected sample.
However, the good agreement between the $K$-band LFs determined from $K$
and $b_J$-limited samples in Fig.~\ref{fig:comp} suggests that our
$b_J$ selection does not significantly affect the estimated $K$-band LF.
In particular, it is unlikely that we are underestimating the faint-end slope
of the $K$-band LF since red galaxies tend to be luminous in $K$
\cite{mes86,sci96}.
We will further address the issue of colour selection in a future paper
in which we plan to parameterize the bivariate LF
$\phi(M_K,M_B)$ using CCD $b_J$ magnitudes in addition to the $K$-band data
presented here.

\section*{Acknowledgments}

It is a pleasure to thank George Efstathiou for suggesting this project,
the CTIO staff for their excellent support, and Bahram Mobasher, Jon Gardner,
Karl Glazebrook and Gyula Szokoly for sending me their $K$-band LF data points.


\section{INTRODUCTION}

The applicability of perturbative QCD to exclusive processes at large
momenta is an
interesting research problem. The Brodsky-LePage \cite{BL} pQCD based
factorization
has been only partially successful.
In this case the process is factorized into a perturbatively calculable
hard scattering piece and the soft distribution amplitude.
The pion electromagnetic
form factor \cite{BL,FJ,ER} at momentum transfer $q^2 = -Q^2$, for example,
can be written as
\begin{equation}
F_\pi(Q^2) = \int dx_1 dx_2 \phi(x_2,Q) H(x_1,x_2,Q) \phi(x_1,Q)
\end{equation}
where $\phi(x,Q) $ are the distribution amplitudes which can
be expressed in terms of the pion wave function $\psi(x,\vec k_T)$ as
\begin{equation}
 \phi(x,Q) = \int^Q d^2k_T\psi(x,\vec k_T) .
\end{equation}
Here $x$ is the longitudinal momentum fraction and $\vec k_T$ the
transverse momentum carried by the quark. The factorization is
possible provided the external photon momentum $Q^2$ is much larger
than the intrinsic quark transverse momentum $k_T^2$, in which case
the $k_T$ dependence of the hard scattering $H$ can be neglected.

The formalism predicts that at large momenta the cross section for
exclusive processes
$d\sigma/dt$, where $t$ is the momentum transfer squared, scales
like $1/t^{n-2}$ up to logs, where $n$ is the total number of elementary
partons participating in the process. The underlying reason
for the power law is
scale invariance of the fundamental
theory. The extra logarithmic
dependence is given by QCD scaling violations.
The dominant contribution to this scattering arises from the valence
quark, since every additional parton leads to an additional suppression
factor of $1/t$. Physically the scattering probes the short distance
part of the hadron wave function. Dominance by the short distance
wave functions leads to several predictions such as helicity conservation,
color transparency  \cite{BM,JPR} etc.

The successes and failures of this scheme
are well known. The predicted momentum dependence of exclusive
processes, in particular the hadronic electromagnetic form factors,
 have generally been found to be in good agreement with data.
However more detailed dynamical predictions such as helicity conservation
in hadron-hadron collisions fail to agree. Calculation
of electromagnetic form factors using this factorization
scheme has been criticised by several authors. The basic problem
is that the momentum scales of the exchanged gluons tend to become
rather small, and the applicability of pQCD becomes doubtful. The
normalization of form factors is largely unknown; use of asymptotic
distribution
amplitudes tends to give small normalizations compared to data. Form
factor magnitudes can
be enhanced by use of model distribution amplitudes which peak closer to the
endpoints, namely $x\rightarrow 0,1$, which then exacerbates the problem
of small internal momentum transfers.

\section{THE SUDAKOV FORM FACTOR}

In order to investigate this problem
in more detail, Botts and Sterman \cite{BS}
and Li and Sterman \cite{LS} developed an alternate factorization which
does not neglect the $k_T$ dependence of the hard scattering. This
formalism also
includes use of a Sudakov form factor. For the case of pion form factor
\cite{LS}
the starting point is,
\begin{equation}
F_\pi(Q^2) = \int dx_1 dx_2 d\vec k_{T1} d\vec k_{T2} \psi^*(x_2,\vec
k_{T2},P_2)
H(x_1,x_2,Q^2,\vec k_{T1},\vec k_{T2}) \psi(x_1,\vec k_{T1},P_1),
\end{equation}
where it is again assumed that the process is factorizable into
hard scattering and soft hadronic wave functions $\psi(x,\vec k_T,P)$. The
calculation is simplified by dropping the $k_T$ dependence in the
quark propagators in hard scattering kernel $H$,
in which case only the combination $\vec k_{T1} + \vec k_{T2}$
appears in the calculation.
The authors \cite{LS} work in configuration space where this can be written
as
\begin{equation}
F_\pi(Q^2) = \int dx_1 dx_2 {d^2\vec b\over (2\pi)^2}  {\cal P}(x_2,b,P_2,\mu)
\tilde H(x_1,x_2,Q^2,\vec b,\mu) {\cal P}(x_1,b,P_1,\mu),
\end{equation}
where ${\cal P}(x,b,P,\mu)$ and $\tilde H(x_1,x_2,Q^2,\vec b,\mu)$ are the
Fourier transforms of the wave function and hard scattering
respectively;
$\vec b$ is conjugate to $\vec k_{T1} + \vec k_{T2}$, $\mu$ is the
renormalization
scale and $P_1$, $P_2$ are the initial and final momenta of the pion.

Sudakov form factors are obtained
by summing the leading and next
to leading logarithms using renormalization group (RG) techniques.
The wave function at small $b$ with a transverse momentum $k_T$ cutoff
equal to $\omega=1/b$ can be approximated by the distribution amplitude
$\phi(x,1/b)$.
Large $k_T$ corrections can be evaluated perturbatively, which result in the
Sudakov form factor. The final result is given by:
\begin{eqnarray}
{\cal P}(x,b,P,\mu) &=& \exp\left[-s(x,\omega,Q)  - s(1-x,\omega,Q) -
2\int^\mu_\omega
{d\bar \mu\over \bar \mu} \gamma_q(\alpha_s(\bar\mu))\right] \nonumber\\
&\times &
\phi(x,1/b) + O(\alpha_s(\omega)).
\end{eqnarray}
where $\gamma_q(\alpha_s)$ is the quark anomalous dimension.
The explicit formula for the function $s(x,\omega,\mu)$ is given
in Li and Sterman \cite{LS}.
Here $\omega=1/b$ plays the role of the factorization scale,
above which QCD corrections give the perturbative evolution of the wave
function $P$, and below which QCD corrections are
absorbed into the nonperturbative distribution amplitude $\phi$. For the
case of the pion,
$1/b$ is the natural choice for this scale. However as discussed below,
for the proton the relevant scale is not obvious and several possibilities
exist in the literature.

\begin{figure} [t,b]
\hbox{\hspace{6em}
 \hbox{\psfig{figure=pion2.ps,height=6cm}}}
  \caption {\em The Sudakov form factor $\exp(-S)$ with $Q^2 = 4 $ GeV$^2$.
For this calculation the QCD scale parameter $\Lambda$ was taken to
be $0.1$ GeV.
    }
  \label{pion2}
\end{figure}

The final formula for the form factor, after incorporating the
renormalization group evolution of the hard scattering from the
renormalization scale
$\mu$ to $t$, $t={\rm max}(\sqrt{x_1x_2}Q,1/b)$, is given by\cite{LS},

\begin{eqnarray}
F_\pi(Q^2) &=& 16\pi C_F \int^1_0 dx_1 dx_2 \phi(x_1)\phi(x_2)\int_0^\infty
b db\alpha_s(t) K_0(\sqrt{x_1x_2}Qb)\nonumber \\
&\times & \exp[-S(x_1,x_2,\omega,Q)],
\end{eqnarray}

where $$S(x_1,x_2,b,Q) = \sum_{i=1}^2\left[s(x_i,b,Q) + s(1-x_i,b,Q)\right]
-4\int_\omega^t{d\bar \mu\over\bar\mu}\gamma_q(\alpha_s(\bar\mu)).$$
The function $e^{-S}$ is plotted in fig. \ref{pion2}. It cuts off large
$b$ regions of the integral and hence the calculation is infrared finite,
without
needing any arbitrary infrared cutoff such as a gluon mass. At
small $b$ the function has been set equal to one.

\begin{figure} [t,b]
\hbox{\hspace{6em}
 \hbox{\psfig{figure=pion1.ps,height=6cm}}}
  \caption {\em The pion form factor $F_\pi (Q^2)$
using the asymptotic (dotted line) and the CZ
    (solid line) distribution amplitudes.
    The experimental data, taken from Ref. 9 with errorbars are
also shown.}
  \label{pion1}
\end{figure}

The resulting form factor using asymptotic as well as  CZ \cite{CZ}
distribution amplitudes is shown in fig. \ref{pion1}. A remarkable
fact is that the correct asymptotic $Q^2$ behavior is seen beyond the scale
of about $Q=1$ GeV,
irrespective of the choice of wave function. In contrast to the Brodsky-LePage
factorization, the $k_T$ dependence of the hard scattering is not
neglected, and
hence this $Q^2$ dependence
does not follow trivially.  It is rather
a detailed dynamical prediction of the theory and depends on the relative
size of intrinsic $k_T^2$ and $x_1 x_2 Q^2$. The prediction is robust, since
it is independent of the details of the distribution amplitude used.
This simple yet important point
justifies the basic idea of Brodsky-LePage factorization, namely that $k_T$
can be treated as negligible in the hard scattering.

\begin{figure} [t,b]
\hbox{\hspace{6em}
 \hbox{\psfig{figure=pion3.ps,height=6cm}}}
  \caption {\em Dependence of $Q^{2}F_{\pi}(Q^2)$ on the transverse
distance cutoff $b_{c}$
    with the asymptotic (dotted line),
    and CZ (solid line) distribution amplitudes for $Q^2=4$ GeV$^2$.
	The QCD scale parameter $\Lambda$ was chosen to be $0.1$ GeV for this
calculation.}
  \label{pion3}
\end{figure}

We note that the normalization of the theoretical result falls below the
experimental data for both choices of distribution amplitude.
However, the large difference between theory and experiment at high momenta
should be interpreted with caution, since as emphasized by Sterman
and Stoler \cite{SS}, there may be large systematic errors in the experimental
extraction of the form factor which are not shown in the figure.
Further theoretical issues in this extraction have been raised in
Ref. \cite{CM}.

In any event, the theoretical normalization of the form factor is comparatively
murky, because it
depends on the details of the distribution amplitude.
Furthermore, the leading order pQCD amplitudes that are practical to
calculate may not give a very
reliable estimate of the
normalization. One can investigate this further by considering the
transverse separation
cutoff ($b_c$) dependence of the form factor.  This can give an idea about the
integration regions important for the calculation. We show the $b_c$ dependence
in fig. \ref{pion3} as originally discussed in Ref. \cite{LS}. Based on
this plot Li and
Sterman argue that roughly 50\% of the contribution can regarded
as perturbative, since it is obtained
from the region where $\alpha_s/\pi < 0.7$.
The observation also implies that higher order contributions in $\alpha_s$
are not
negligible, and the leading order predictions for the normalization of the
form factor cannot be regarded as accurate. The next to leading order
calculation \cite{passek} of the pion form factor also leads to the same
conclusion.

\medskip

We are left with the following interesting situation: Although the basic
Brodsky-LePage
factorization is correct, one may need to go to higher orders in $\alpha_s$
in order to obtain an accurate prediction for the form factor normalizations.
However the predicted $Q^2$ dependence appears
to be quite robust, and independent of the theoretical uncertainties
such as the choice of distribution amplitude.

\section{THE PROTON}

The improved factorization has also been applied to the proton Dirac form
factor $F^p_1(Q^2)$\cite{L}.
The calculation is considerably more complicated compared to the pion.
Here also it is necessary to use distribution amplitudes which peak close to
the end points in order to obtain the experimental normalization
of the form factor. In contrast to pion, there is no
natural choice for the infrared cutoff $\omega$ in the Sudakov
exponent, due to the presence of three quarks and resulting three
distances.

The Sudakov resummation of large logarithms in $\cal P$ leads to
\begin{equation}
{\cal P}(x_i,{\bf b}_i,P,\mu)
=\exp\left[-\sum_{l=1}^3 s(x_l,w,Q)-3\int_{w}^{\mu}
\frac{d\bar{\mu}}{\bar{\mu}}\gamma_{q}\left(\alpha_s(\bar{\mu})
\right)\right]
\times\phi(x_{i},w)\; ,
\label{8}
\end{equation}
where the quark anomalous dimension $\gamma_{q}(\alpha_s)=-\alpha_{s}/\pi$
in axial gauge governs the RG evolution of $\cal P$.
The function $\phi$ is the standard proton distribution amplitude.
The exponent $s$ is given in Ref. \cite{L}.

In equation \ref{8} we use the same infrared parameter $\omega$ for all
the three $s(x_l,\omega,Q)$ for $l=1,2,3$ as well as in the integrals
over the anomalous dimension. Earlier
Li \cite{L} chose to use different infrared cutoffs
 $b_l$
for each exponent $s$ and for each integral involving $\gamma_q$.
As pointed out in \cite{JKB} this choice does not
does not suppress the soft divergences
from $b_l\to 1/\Lambda$ completely, where $\Lambda$ is the QCD
scale parameter. For example, the divergences from
$b_1\to 1/\Lambda$, which appear in $\phi(x_i,w)$ as $w\to \Lambda$,
survive as $x_1\to 0$, since
$s(x_1,b_1,Q)$ vanishes and $s(x_2,b_2,Q)$ and $s(x_3,b_3,Q)$ remain finite.
On the other hand, $w$ should play the role of the factorization scale,
above which QCD corrections give the perturbative evolution of the wave
function ${\cal P}$ in Eq.~(\ref{8}), and below which QCD corrections are
absorbed into the initial condition $\phi$. It is then not
reasonable to choose the cutoffs $b_l$ for the Sudakov resummation
different from $w$.

A modified choice of the cutoffs,
 $w=1/b_{max}$, $b_{max}
=\max(b_l)$, $l=1,2,3$,
was proposed in Ref. \cite{JKB}.
This choice was found to suppress
the soft enhancements, and the form factor was found
to saturate as $b_c \rightarrow 1/\Lambda$.
The authors also included a model non-perturbative
soft wave function in the calculation.
Unfortunately, it turned out that the normalization of the
resulting $Q^4F_1$ was
found to be
less than half of that of the data for all the distribution amplitudes
explored \cite{JKB}. Bolz et al \cite{JKB} then concluded that pQCD
is unable to fit the experimental form factor.

Kundu et al \cite{KLSJ} reexamined the situation.  The group argued that the
 form factor normalization is sensitive to the precise
choice of the infrared cutoff $w$. They used $w=c/b_{max}$, $b_{max}
=\max(b_l)$, $l=1,2,3$, as the infrared cutoff in the Sudakov
exponent,
instead of the one used by Bolz et al \cite{JKB}.
The introduction of  parameter $c$ is
natural from the viewpoint of the resummation, since the scale $cw$,
with $c$ of order unity, is as appropriate as $w$ \cite{BS}.
Kundu et al \cite{KLSJ} find that the calculation is in good agreement
with data
using the King-Sachrajda (KS)
\cite{KS} distribution amplitude and setting $c=1.14$.

The choice $c=1.14$ can also be motivated physically
by considering the proton as a quark-diquark type configuration.
The diquark constituents are the two quarks closest to
each other in the transverse plane. Let $d_{\rm typ}$ be the distance
between the
center of mass of the diquark and the remaining third quark. Then
the infrared cutoff scale $\omega$ can be taken to
be $1/d_{\rm typ}$.
We choose $c$ such that
for a large number of randomly chosen triangles,
 we get for the average
$\langle d_{typ}/b_{max} \rangle=1/c$.
Defining $c$ in such a way, gives
$c\approx1.14$.

\begin{figure} [t,b]
\hbox{\hspace{6em}
 \hbox{\psfig{figure=fig2.ps,height=6cm}}}
  \caption {\em Dependence of $Q^4F_1^p$ on $Q^2$ using the
    KS distribution amplitude ($c=1.14$, solid line; $c=1$, dense-dot line)
     and the CZ distribution amplitude ($c=1.14$, dashed line;
               $c=1$, dotted line).
    The experimental data with error bars, taken from Ref. 17,
 are also shown.}
  \label{proton1}
\end{figure}

\begin{figure} [t,b]
\hbox{\hspace{6em}
 \hbox{\psfig{figure=fig3.ps,height=6cm}}}
  \caption {\em Dependence of $Q^{4}F_{1}^{p}$ on the cutoff $b_{c}$
    with the KS distribution amplitude for $Q^2=12$ GeV$^2$ (dotted line),
    $Q^2=16$ GeV$^2$ (dashed line), $Q^2=25$
    GeV$^2$ (dense-dot line), and $Q^2=36$ GeV$^2$ (solid line).
	The QCD scale parameter $\Lambda$ was taken to be $0.2$
GeV for this calculation.
}
  \label{proton2}
\end{figure}

\begin{figure} [t,b]
\hbox{\hspace{6em}
 \hbox{\psfig{figure=soft.ps,height=6cm}}}
  \caption {\em $Q^2$ dependence of $Q^4F_1^p$ including a soft model
	wave function with $<k_T^2> = 0.271$ GeV$^2$.
The four upper curves at $Q^2 = 35$ GeV$^2$ use the KS
distribution amplitude with the infrared parameter
$c=1.4$ (solid curve), $c=1.2$ (dashed
curve) and $c=1$ (dashed-dotted curve). The dotted curve shows the result
without including model soft wave function with $c=1$. The two lower curves
at $Q^2=35$ GeV$^2$ use CZ distribution amplitude
with $c=1$. The lowest curve includes the
soft model wave function, whereas the upper curve does not. The
    experimental data  with error bars  [17] are also shown.}
  \label{proton3}
\end{figure}

The results of the calculation using KS \cite{KS}
and CZ \cite{CZ} distribution amplitudes and $c=1$ and 1.14
are shown in fig. \ref{proton1}.
It is found that with $c=1.14$ the KS distribution amplitude gives good
agreement with data. The $b_c$ dependence of the form factor is shown in
fig. \ref{proton2}, which shows saturation at about
$b_c = 0.8/\Lambda$.
The result after including
a model nonperturbative soft wave function
are displayed in fig. \ref{proton3}. Again we find that choosing
$c$ of order unity gives pQCD calculations in agreement with data.
For all choices of the distribution
amplitude and parameter $c$, independent of whether the model soft
wave function is included or not, the $Q^2$ dependence of the form factor is
in good agreement with data. An analogous situation was found for pion
form factor.

\medskip

The natural agreement of $Q^{2}$ dependence of the pQCD calculations
is in contrast to data fits
obtained using soft model \cite{Rady,Diehl}. In such models the
$Q^2$ dependence depends on the details of the
model wave function. Soft model predictions at high momentum have a
tendency to fall
more strongly than experimental data. As also noted for the pion,
the approximate power law behaviour in  $Q^2$ is
not directly implied by the factorization and is a detailed dynamical
prediction of the calculation. This could have been achieved only
if the intrinsic $k_T$ were negligible in the hard scattering kernel.
Thus the observation of power-law dependence in the data
lends considerable support to the basic
factorization of Brodsky-LePage. Nevertheless, the relatively small
magnitude of
internal momentum
scales and the sensitive dependence of final result on the
infrared scale $w$ indicates that calculation of normalizations using
leading order diagrams is not reliable. While higher order
contributions are propably non-negligible, there is every reason to
believe that the power-law dependence of the calculations is robust.

We have reviewed the current status of pQCD
calculations of hadronic electromagnetic
form factors. We argue that the normalization of the form factor
cannot be predicted reliably by a leading order calculation in $\alpha_s$.
Detailed calculations including the soft $k_T$ dependence, however, support
the basic factorization scheme. One finds the correct asymptotic $Q^2$
evolution
of the form factors for $Q$ as small as 2-3 GeV, independent of the choice
of distribution amplitude and other theoretical uncertainties.
Hence we conclude that agreement of quark counting scaling predictions is
not accidental but
is well supported by detailed dynamical calculations.

PJ and BK thank the staff of ICTP, Trieste, for hospitality during a
visit where this paper was written. This work was supported by
BRNS grant No. DAE/PHY/96152,
the National Science Council of the
Republic of China Grant No. NSC-87-2112-M-006-018,
 the Crafoord Foundation,
the Helge Ax:son Johnson Foundation, DOE Grant number
DE-FGO2-98ER41079,
the KU General Research Fund and NSF-K*STAR
Program under the Kansas Institute for Theoretical and Computational
Science.


\section{Introduction}

Millisecond oscillations in the X-ray brightness during thermonuclear bursts,
``burst oscillations", have been observed from six low
mass X-ray binaries (LMXB) with the Rossi X-ray Timing Explorer (RXTE) (see
\markcite{SSZ}Strohmayer, Swank \& Zhang {\it et al.} 1998 for a recent review).
Considerable evidence points to rotational modulation as the source of these
pulsations (see for example, \markcite{SZS}Strohmayer, Zhang \& Swank 1997;
\markcite{SM99}Strohmayer \& Markwardt 1999). Anisotropic X-ray emission caused by
either localized or
inhomogeneous nuclear burning produces either one or a pair of hot spots on the
surface which are then modulated by rotation of the neutron star. A remarkable
property of these oscillations is the frequency evolution which occurs in the
cooling tail of some bursts. Recently, \markcite{SM99}Strohmayer \& Markwardt (1999)
have shown that the frequency in the cooling tail of bursts from 4U 1728-34 and
4U 1702-429 is well described by an exponential chirp model whose frequency increases
asymptotically toward a limiting value. \markcite{SJGL}Strohmayer et. al (1997)
have argued this evolution results from angular momentum conservation of the
thermonuclear shell, which cools, shrinks and spins up as the surface radiates
away the thermonuclear energy. To date, only frequency increases have been reported
in the cooling tails of bursts, consistent with settling of the shell as its energy is
radiated away.

In this Letter we report observations of a {\it decreasing} burst oscillation
frequency in the tail of an X-ray burst. We find that an episode of spin down in
the cooling tail of a burst observed on December 31, 1996 at 17:36:52 UTC
(hereafter, burst A, or the ``spin down'' burst)
from 4U 1636-53 is correlated with the presence of an extended. In \S 1 we present
an analysis of the frequency evolution in this burst, with emphasis on the spin down
episode. In \S 2 we present time resolved energy spectra of the spin down burst, and we
investigate the energetics of the extended tail. Throughout, we compare the temporal
and spectral behavior of the spin down burst with a different burst observed on
December 29, 1996 at 23:26:46 UTC (hereafter, burst B) which does not show either a
spin down episode nor an extended tail of emission, but which is similar to the spin
down burst in most other respects. We conclude in \S 3 with a summary and discussion of
the spin down episode and extended emission in the context of an additional, delayed
thermonuclear energy release which might re-expand the thermonuclear shell and perhaps
account for both the spin down and the extended tail of thermal emission.

\section{Evidence for Spin Down}

Oscillations at 580 Hz were discovered in thermonuclear bursts from
4U 1636-53 by Zhang et al. (1996). More recently, Miller (1999a) has reported evidence
during the rising phase of bursts of a significant modulation at half the 580 Hz
frequency  suggesting that 580 Hz is twice the neutron star spin frequency
and that a pair of antipodal spots produce the oscillations. Here we focus on a
burst from 4U 1636-53 which shows a unique decrease in the $\approx 580$ Hz oscillation
frequency. To study the evolution in frequency of burst oscillations we employ the
$Z_n^2$ statistic (\markcite{B83}Buccheri et al. 1983). We have described this
method previously, and details can be found in \markcite{SM99}Strohmayer \&
Markwardt (1999). We first
constructed for both bursts A and B a dynamic ``variability'' spectrum by computing
$Z_1^2$ as a function of time on a grid of frequency values in the vicinity of 580 Hz.
We used 2 second intervals to compute $Z_1^2$ and started a new interval every 0.25
seconds. This variability spectrum is very similar to a standard dynamic power spectrum,
however, the $Z_1^2$ statistic allows for a more densely sampled frequency grid than a
standard Fast Fourier Transform power spectrum. The results are shown in Figure 1
(bursts A and B are in the top and bottom panels, respectively) as contour maps of
constant $Z_1^2$ through each burst. The contour map for the spin down burst (top
panel) suggests that the oscillation began with a frequency near 579.6 Hz at burst
onset, reappeared later in the burst after ``touchdown'' of the photosphere at an
increased frequency, $\approx 580.7$ Hz, but then beginning near 11 seconds dropped to
$\approx 579.6$ Hz over several seconds. For comparison, we also show in Figure 1
(bottom panel) a similar variability spectrum for burst B which also shows strong
oscillations near 580 Hz, but shows no evidence of a similar spin down episode.

To investigate the evolution of the oscillation frequency more quantitatively we fit
a model for the temporal evolution of the frequency, $\nu(t)$, to the 4.5 second
interval during which the oscillation is evident in the dynamic variability spectrum
(Figure 1, top panel). Our model is composed of two linear segments, each with its
own slope, joined continuously at a break time $t_b$. This is similar to the model
employed by \markcite{M99b}Miller (1999b), and has four free parameters, the initial
frequency, $\nu_0$, the two slopes, $d_{\nu}^1$ and $d_{\nu}^2$, and the break time,
$t_b$. We used this frequency model to compute phases $\phi_{t_j}$ for each X-ray
event, {\it viz.} $\phi_{t_j} = \int_0^{t_j} \nu(t') dt'$, where $t_j$ are the
photon arrival times, and then varied the model parameters to maximize the $Z_1^2$
statistic. We used a downhill simplex method for the maximization (see Press et al.
1989). Figure 2 compares $Z_1^2$ vs. parameter $\nu_0$ for the best fitting two
segment model (solid histogram) and a simple constant frequency model ($\nu(t) =
\nu_0$, dashed histogram). The two segment model produces a significant increase in the
maximum $Z_1^2$ of about 40 compared with no frequency evolution, and it also
yields a single, narrower peak in the $Z_1^2$ distribution.  The
increase of 40 in $Z_1^2$, which for a purely random process is distributed as $\chi^2$
with 2 degrees of freedom, argues convincingly that the frequency drop is significant.
We note that \markcite{M99b}Miller (1999b) has also identified the same spin down
episode during this burst using a different, but related method. The best fitting two
segment model is shown graphically as the solid curve in Figure 1 (top panel).

\section{Time history, spectral evolution and burst energetics}

A comparison of the 2 - 20 keV time history of the spin down burst with other bursts
from the same observations reveals that this burst is also unique in having an
extended tail of thermal emission. This is well illustrated in Figure 3,
which compares the 2 - 20 keV time histories of the spin down burst and burst B.
To further investigate the energetics of the
thermal burst emission we performed a spectral evolution analysis.
We accumulated energy spectra for varying time intervals through both bursts. Using
XSPEC we fit blackbody spectra to each interval by first subtracting a pre-burst interval
as background, and then investigated the temporal evolution of the blackbody temperature,
$kT$, inferred radius, $R_{BB}$ and bolometric flux, $F$. In most intervals we obtained
acceptable fits with the blackbody model. The results for both bursts are summarized
in Figure 4, We have aligned the burst profiles in time
for direct comparison. Both bursts show evidence for radius expansion shortly after
onset in that $kT$ drops initially and then recovers. Their peak fluxes are also
similar, consistent with being Eddington limited. Out to about 7
seconds post-onset both bursts show the same qualitative behavior, after this,
however, the spin down burst (solid curve) shows a much more gradual decrease
in both the blackbody temperature $kT$ and the bolometric flux $F$ than is evident in
burst B. We integrated the flux versus time profile for each burst in order to estimate
fluences and establish the energy budget in the extended tail. We find fluences of $1.4
\times 10^{-6}$ and $5.1 \times 10^{-7}$ ergs cm$^{-2}$ for bursts A and B
respectively. That is, the spin down burst has about 2.75 times more energy than burst
B. Put another way, most of the energy in the spin down burst is in the extended tail.
In figure 4 we also indicate with a vertical dotted line the time $t_b$ associated with
the beginning of the spin down episode based on our modelling of the 580 Hz
oscillations. The spectral evolution analysis indicates that at about the same
time the spin down episode began there was also a change in its spectral evolution
as compared with that of burst B. This behavior is evident in figure 5 which
shows the evolution of $kT$ (dashed curve) and the inferred blackbody radius $R_{BB}$
(solid curve) for the spin down burst. Notice the secondary increase in $R_{BB}$ and an
associated dip in $kT$ near time $t_b$ (vertical dotted line). This behavior
is similar to the signature of radius expansion seen earlier in both bursts, but at
a weaker level, it suggests that at this time there may have been an additional
thermonuclear energy input in the accreted layers, perhaps at greater depth, which
then diffused out on a longer timescale, producing the extended tail. That this
spectral signature ocurred near the same time as the onset of the spin down
episode suggests that the two events may be causally related.

\section{Discussion and Summary}

The observation of thermonuclear bursts with extended tails is not a new phenomenon.
\markcite{CCG}Czerny, Czerny, \& Grindlay (1987) reported on a burst from the soft X-ray
transient Aql X-1 which showed a long, relatively flat X-ray tail. Bursts following
this one were found to have much shorter, weaker tails. \markcite{F92}Fushiki et al.
(1992) argued that such long tails were caused by an extended phase of hydrogen burning
due to electron captures at high density ($\rho \approx 10^7$ g
cm$^{-3}$) in the accreted envelope. Such behavior is made possible because of the long
time required to accumulate an unstable pile of hydrogen-rich thermonuclear fuel when
the neutron star is relatively cool ($\approx 10^7$ K) prior to the onset of accretion.
This, they argued, could occur in transients such as Aql X-1, which have long quiescent
periods during which the neutron star envelope can cool, and thus the first
thermonuclear burst after the onset of an accretion driven outburst should show the
longest extended tail. Other researchers have shown that the thermal state of the
neutron star envelope, the abundance of CNO materials in the accreted matter, and
variations in the mass accretion rate all can have profound effects on the character of
bursts produced during an accretion driven outburst (see \markcite{T93}Taam et al.
1993; \markcite{WW85}Woosley \& Weaver 1985; and \markcite{AJ82}Ayasli \& Joss 1982,
and references therein). For example, \markcite{T93}Taam et al. (1993) showed that for
low CNO abundances and cool neutron star envelopes the subsequent
bursting behavior can be extremely erratic, with burst recurrence times varying by as
much as two orders of magnitude. They also showed that such conditions produce dwarf
bursts, with short recurrence times and peak fluxes less than a tenth Eddington, and
that many bursts do not burn all the fuel accumulated since the last X-ray burst. Thus
residual fuel, in particular hydrogen, can survive and provide energy for subsequent
bursts. These effects lead to a great diversity in the properties of X-ray bursts
observed from a single accreting neutron star. Some of these effects were likely
at work during the December, 1996 observations of 4U 1636-53 discussed here, as
both a burst with a long extended tail, as well as a dwarf burst were observed (see
\markcite{M99b}Miller 1999b).

The spin up of burst oscillations in the cooling tails of thermonuclear bursts from
4U 1728-34 and 4U 1702-429 has been discussed in terms of angular momentum
conservation of the thermonuclear shell (see \markcite{SJGL}Strohmayer et al 1997;
\markcite{SM99}Strohmayer \& Markwardt 1999). Expanded at burst onset by the initial
thermonuclear energy release, the shell spins down due to its larger rotational
moment of inertia compared to its pre-burst value. As the accreted layer subsequently
cools its scale height decreases and it comes back into co-rotation with the neutron
star over $\approx$ 10 seconds. To date the putative initial spin down at burst onset
has not been observationally confirmed, perhaps due to the radiative diffusion delay,
on the order of a second, which can hide the oscillations until after the shell has
expanded and spun down (see, for example, \markcite{B98}Bildsten 1998). We continue to
search for such a signature, however. Although the initial spin down at burst onset has
not been seen, the observation of a spin down episode in the {\it tail} of a burst begs
the question; can it be understood in a similar context, that is, by invoking a second,
thermal expansion of the burning layers?  The supporting evidence
is the presence of spin down associated with the spectral evidence for an
additional energy source, the extended tail, as well as the spectral variation observed
at the time the spin down commenced (see Figure 5). Based on these
observations we suggest that the spin down began with a second episode of thermonuclear
energy release, perhaps in a hydrogen-rich layer underlying that responsible for the
initial instability, and built up over several preceeding bursts. Such a scenario is not
so unlikely based on previous theoretical work (see \markcite{T93}Taam et al. 1993,
\markcite{F92}Fushiki et al. 1992). The observed rate of spin down, $d_{\nu}^2 =
-1.01\times 10^{-3}$ s$^{-1}$, interpreted as an increase in the height of the angular
momentum conserving shell gives $dr/dt \approx (\Delta\nu / 2\Delta T \nu_0)R \approx
5.25$ m s$^{-1}$, for a neutron star radius of $R = 10$ km. Calculations predict
increases in the scale height of the bursting layer on the order of 20-30 m during
thermonuclear flashes (see \markcite{J77}Joss 1977; and \markcite{B95}Bildsten 1995).
Based on this and the energy evident in the extended tail, the additional expansion of
about 12 m does not appear overly excessive. If correct this scenario would require
that the oscillation frequency eventually increase again later in the tail.
Unfortunately the oscillation dies away before another increase is seen.

It is interesting to note that bursts from 4U 1636-53 do not appear to show the same
systematic evolution of the oscillation frequency as is evident in bursts from
4U 1728-34 and 4U 1702-429 (see for example, \markcite{M99b}Miller 1999b; and
\markcite{SM99}Strohmayer \& Markwardt 1999). In particular, there is no strong
evidence for an exponential-like recovery that is often
seen in 4U 1728-34 and 4U 1702-429. Rather, in 4U 1636-53, when the burst oscillation
frequency reappears after photospheric touchdown in many bursts it appears almost
immediately at the higher frequency. In the context of a spinning shell this might
suggest that the shell recouples to the underlying star more quickly than in 4U 1728-34
or 4U 1702-429. Interestingly, 4U 1636-53 is also the only source to show significant
pulsations at the sub-harmonic of the strongest oscillation frequency, and this has
been interpreted as revealing the presence of a pair of antipodal hot spots (see
\markcite{M99a}Miller 1999a). These properties may be related and could, for example,
indicate the presence of a stronger magnetic field in 4U 1636-53 than the other sources.

Another physical process which could alter the observed frequency is related to
general relativistic (GR) time stretching. If the burst evolution modulates the
location in radius of the photosphere,
then the rotation at that radius, as seen by a distant observer, is affected by a
redshift such that, $\Delta r /R = (R/r_g)(1-r_g/R)(1-1/(\nu_{h}/\nu_{l})^2)$, where
$\Delta r$ is the change in height of the photosphere, $r_g = 2GM/c^2$ is the
schwarzchild radius, and $\nu_h/\nu_l$ is the ratio of the highest and the lowest
observed frequencies.  If this were the sole cause of the frequency changes, then
it would imply a height change for the photosphere of
$\approx 120$ m, which is much larger than the increases predicted theoretically
for bursting shells.  Note that this effect works counter to angular momentum
conservation of the shell, since increasing the height makes the frequency
higher compared to deeper layers. Since the thicknesses of pre- and post-burst shells
are on the order of 20 - 50 m, we estimate from the above that the GR correction
amounts to about 10 - 20\% of the observed frequency change, and, if the angular
momentum conservation effect is at work, requires a modest {\it increase} in the
height of the shell over that estimated non-relativistically.

We have reported in detail on the first observations of a spin down in the frequency
of X-ray brightness oscillations in an X-ray burst. We have shown that this event
is coincident with the ocurrence of an extended tail as well as a spectral signature
which both suggest a secondary release of thermonuclear energy in the accreted layer.
It is always difficult to draw conclusions based on a single event, however,
if the association of spin down episodes with an extended X-ray tail can be
confirmed in additional bursts this will provide strong evidence in support of the
hypothesis that angular momentum conservation of the thermonuclear shell is
responsible for the observed frequency variations during bursts. The combination
of spectral and timing studies during bursts with oscillations can then give us a
unique new probe of the physics of thermonuclear burning on neutron stars.

\acknowledgements

We thank Craig Markwardt and Jean Swank for many helpful discussions.

\vfill\eject