text:
# Statistical mechanics of granular gases in compartmentalized systems

## I Introduction

In recent years considerable progress has been achieved in the understanding of granular gases, i.e. assemblies of moving inelastic solid grains mutually colliding and loosing a little energy in each collision gases ; jaeger . They exhibit a fascinating and rich phenomenology which comprehends clustering goldhirsch ; gold2 ; mcnamara  spontaneous vortex formation baldassa , and breakdown of kinetic energy equipartition in granular mixtures Dufty ; Ioepuglio . More fundamentally granular gases represent one of the prototypical non equilibrium systems.

Major progress have been achieved by considering the simplest situations, i.e. spatially homogeneous and time-independent systems. In spite of that, inhomogeneities do not always represent a nuisance or an undesired complication, but, on the contrary, they can be a source of new insight. Within such a perspective, external fields have been introduced on purpose in order to probe the inhomogeneous behavior of granular gases. The seminal work of Schlichting and Nordmeier experiment  has stimulated a vivid interest in the behavior of the so called compartmentalized systems. A box, whose base moves periodically up and down, is separated by a vertical barrier into two compartments which communicate through an orifice. For strong shaking the two halves are equally populated, as it would occur in the case of a standard molecular fluid, whereas for weaker driving the symmetry is spontaneously broken. The resulting scenario resembles that of an equilibrium thermodynamical phase transition, with the population difference between the two compartments playing the role of the order parameter and the driving intensity that of temperature. Our understanding of the problem has increased since then due to a series of studies. These range from new experiments Lohse ; Lohse1 , phenomenological flux models Eggers ; Droz , hydrodynamic equations and Molecular Dynamics simulation Brey .

The scope of the present contribution is to show that by choosing a simplified yet significant model of compartmentalized granular gas it is possible not only to derive a set of equations describing the evolution of the macro-state of the system, but also to obtain information about its microscopic fluctuations. The present approach represents a bridge between the phenomenological level of refs. Lohse ; Lohse1 ; Droz  and the statistical mechanical level dererumpallettarum .

The present paper is organized as follows: in the first section we define the model and introduce the statistical description of the system, based on a Boltzmann equation for the distribution functions, modified to take into account the stochastic driving. At this stage we follow the strategy of integrating out microscopic fluctuations going from a microscopic description based on the Boltzmann equation to a macroscopic level, where only the occupation numbers and the granular temperatures of the two compartments, i.e. the first two moments of the distribution function are taken into account. This is equivalent to neglect inhomogeneities of the system at scales smaller than the linear size of the compartments. To obtain qualitative insight, we first try a Gaussian solution. In other words, we assume that the velocity PDF’s are Gaussians, whose normalizations and variances (related to the occupation numbers of the compartments and to their granular temperatures, respectively), can be determined by means of a set of self-consistent differential equations. The analysis shows the existence of a transition from a symmetric phase, where the compartments are equally populated and are at the same granular temperature, to an asymmetric phase, where the two compartments have different properties. Various predictions are made: we locate the bifurcation point, obtain relations between the asymptotic values of the temperatures and densities in the compartments, estimate the characteristic times to observe the symmetry breaking, etc. In the second section, we relax the Gaussian hypothesis about the nature of the velocity fluctuations and let the populations in the two compartments to fluctuate as well. This is done by solving numerically the full Boltzmann equation by means of the Direct Simulation Monte Carlo (DSMC) puglisi . We notice deviations with respect to the treatment of section I, in particular near the critical point. Finally in the last section we present our conclusions.

## II The model

Let us consider a system, composed of two compartments A and B of the same volume, $`V_A=V/2`$, and containing $`N_A`$ and $`N_B`$ particles, respectively. Particle pairs belonging to the same compartment may collide inelastically loosing a fraction of their kinetic energy, but conserving their total momentum. The post-collisional velocities ($`𝐯_1^{},𝐯_2^{}`$) are determined by the transformation:

$`𝐯_1^{}=𝐯_1\frac{1}{2}(1+\alpha )(𝐯_{12}\widehat{𝝈})\widehat{𝝈}`$ (1)
where $`𝐯_{12}=𝐯_1𝐯_2`$ and $`\alpha `$ is the restitution coefficient.

In addition, we assume that all the grains are subject to the action of an external stochastic driving force and their motion between two successive collisions is described by the following Ornstein-Uhlenbeck process:

$$\frac{\mathrm{d}𝐯_i}{\mathrm{d}t}=\frac{1}{\tau _b}𝐯_i+\widehat{\xi _𝐢},$$
(2)
where $`\tau _b^1𝐯_i`$ is a friction term and $`\xi _𝐢`$ a Gaussian random acceleration, of zero average and variance given by:

$$\xi _{i\alpha }(t)\xi _{j\beta }(t^{})=2\frac{T_bm}{\tau _b}\delta _{ij}\delta _{\alpha \beta }\delta (tt^{}),$$
(3)
where $`T_b`$ is a measure of the intensity of the driving. Notice that in the elastic case ($`\alpha =1`$), the average kinetic energy per particle moving in a d-dimensional space is $`K=dT_b/2`$, which is just the ideal gas value. The simultaneous presence of frictional dissipation and random “kicks” renders the kinetic energy of the system stationary even in the absence of collisional dissipation.

To complete the model, we allow the particles contained in A (B) to move into B (A) with a probability per unit time, $`\tau _s^1`$, provided their kinetic energy exceeds a given threshold, let say $`T_s=\frac{1}{2}mu_s^2`$. Such a mechanism schematizes the jump process by which particles may pass from one compartment to the other by overcoming a vertical barrier of height $`h`$ only if $`\frac{1}{2}mv^2>mgh`$.

We shall make the key assumption, dictated by mathematical convenience, of ignoring the spatial gradients which characterize experimental situations. To be specific we stipulate that within the spatial domain, that represents a compartment, the system is homogeneous. The effect of the dividing vertical wall is schematized by a selection rule of probabilistic nature which allows only the more energetic particles to cross the barrier.

To make analytic progress we shall transform the stochastic evolution equations for the velocities of the particles into a system of deterministic equations for the distribution functions. In order to achieve this description, we shall treat the collisions at the level of the Boltzmann molecular chaos approximation. This approximation, widely employed even in the case of granular systems, is equivalent to neglect correlations among the colliding particles.

In order to study the statistical evolution of the system we assume that the single particle phase-space distribution function $`f(𝐫,𝐯,t)`$ satisfies the following properties:

$$N_{A(B)}(t)=_{V_{A(B)}}d𝐫d𝐯f(𝐫,𝐯,t)$$
(4)
where $`N_{A(B)}(t)`$ is the average number of particles in compartment A (B) at instant $`t`$. The average kinetic energy per particle, the granular temperature $`T_{A(B)}`$, is defined as:

$$T_{A(B)}(t)=\frac{1}{N_{A(B)}(t)d}_{V_{A(B)}}d𝐫d𝐯m𝐯^2f(𝐫,𝐯,t)$$
(5)
In the following we shall make the assumption that to describe the essential properties of the system, we do not need the detailed information contained in $`f(𝐫,𝐯,t)`$, but these can can be captured by the following coarse grained distributions obtained by eliminating the $`𝐫`$ dependence of the original distributions

$$f_{A(B)}(𝐯,t)=\frac{1}{V_{A(B)}}_{V_{A(B)}}d𝐫f(𝐫,𝐯,t)$$
(6)
The change of $`f_{A(B)}(𝐯,t)`$ is given by:

$$_tf_A(𝐯_1,t)=I(f_A,f_A)+\frac{T_b}{\tau _b}\left(\frac{}{𝐯_1}\right)^2f_A(𝐯_1,t)+\frac{1}{\tau _b}\frac{}{𝐯_1}f_A(𝐯_1,t)\frac{1}{\tau _s}\theta (|𝐯_1|u_s)[f_A(𝐯_1,t)f_B(𝐯_1,t)]$$
(7)
where $`\theta (x)`$ is the Heaviside function. The first term represents the change of the distribution due to collisions, the second and the third terms are due to the interaction with the heat-bath and the last describes the population change due to the particles migrating from one compartment to the other puglisi .

In order to obtain an explicit expression for the collision integral $`I(f,f)`$ in the case of inelastic hard spheres, we shall follow closely the derivation of Ernst and Van Noije vannoije  and set:

$$I(f_A,f_A)=\sigma ^{d1}d𝐯_2^{}d\widehat{𝝈}(𝐯_{12}\widehat{𝝈})\left\{\frac{1}{\alpha ^2}f_A(𝐯_1^{},t)f_A(𝐯_2^{},t)f_A(𝐯_1,t)f_A(𝐯_2,t)\right\}$$
(8)
The prime on the $`\widehat{𝝈}`$ integration enforces the condition $`𝐯_{12}\widehat{𝝈}>0`$, where $`\widehat{𝝈}`$ is a unit vector along the line of centers of the colliding spheres at contact, whereas $`𝐯_i^{}`$ represent the precollisional velocities, which can be computed by inversion of eq. (2).

By integrating over the spatial coordinate and over the velocity it is straightforward to obtain the equation for the rate of change of the occupation numbers $`N_A`$ and $`N_B`$. It reads:

$$\frac{\mathrm{d}N_A(t)}{\mathrm{d}t}=\frac{V_A}{\tau _s}d𝐯_1[f_A(𝐯_1,t)f_B(𝐯_1,t)]\theta (|𝐯_1|u_s).$$
(9)
We see that when the typical time $`\tau _s`$ diverges, there is no particle exchange between the two compartments, thus $`N_A,N_B`$ are constant, and the temperatures reach a constant value, which depends on the heat bath properties, the physical properties of the particles and on their number.

In order to obtain the energy equation, instead, we multiply both sides by $`m𝐯_1^2`$ and integrate over coordinate and velocity space:

$$\frac{\mathrm{d}[N_A(t)T_A(t)]}{\mathrm{d}t}=\frac{2}{\tau _b}N_A[T_bT_A(t)]+\frac{2}{d}\frac{N_A}{V_A}\sigma ^{d1}v_A\mu _2N_A(t)T_A(t)$$
(10)
$$\frac{V_A}{\tau _sd}d𝐯_1m𝐯_1^2[f_A(𝐯_1,t)f_B(𝐯_1,t)]\theta (|𝐯_1|u_s)$$
(11)
where we have introduced, following ref. vannoije , the non dimensional quantity $`\mu _2`$ :

$$\mu _2=\frac{1}{v_A^3}\frac{V^2}{N^2}d𝐯_1𝐯_1^2I(f_A,f_A)=$$
(12)
$$\frac{1}{v_A^3}\frac{V^2}{N^2}d𝐯_1𝐯_1^2d𝐯_2^{}d\widehat{𝝈}(𝐯_{12}\widehat{𝝈})\left\{\frac{1}{\alpha ^2}f_A(𝐯_1^{},t)f_A(𝐯_2^{},t)f_A(𝐯_1,t)f_A(𝐯_2,t)\right\}$$

with $`v_{A(B)}^2=2T_{A(B)}/m`$. In the following we shall set $`m=1`$. A more compact notation is achieved if one introduces the collision frequency $`\omega _A(t)`$ and the non dimensional spontaneous cooling rate $`\gamma `$:

$$\omega _A=\frac{\mathrm{\Omega }_d}{\sqrt{2\pi }}\frac{N_A}{V_A}\sigma ^{d1}v_A$$
(13)
$$\gamma =\frac{\sqrt{2\pi }}{d\mathrm{\Omega }_d}\mu _2$$
(14)
with $`\mathrm{\Omega }_d=\frac{2\pi ^{d/2}}{\mathrm{\Gamma }(d/2)}`$. the surface area of a d-dimensional unit sphere. In the remainder of the paper we shall measure the temperature in units of $`T_s`$, time in units of $`\tau _b`$ and length in units of $`\sigma `$.

To proceed further analytically we choose $`d=2`$ and assume an approximate Gaussian form sonine  of the distribution function:

$$f_{A(B)}(𝐯,t)=\frac{N_{A(B)}}{V_{A(B)}}\frac{1}{\pi 2T_{A(B)}}\mathrm{exp}(\frac{𝐯^2}{2T_{A(B)}})$$
(15)
In this case one can evaluate explicitly the collision term vannoije :

$$\mu _2=(1\alpha ^2)\frac{\pi ^{d/21}}{\sqrt{2}\mathrm{\Gamma }(d/2)}$$
(16)
$$\gamma \omega _A=\sigma (1\alpha ^2)\frac{N_A}{2V_A}\sqrt{\frac{T_A}{m}}$$
(17)
We find:

$$\frac{\mathrm{d}N_A(t)}{\mathrm{d}t}=\frac{1}{\tau _s}[N_Be^{T_s/T_B}N_Ae^{T_s/T_A}]$$
(18)
where the temperature $`T_s`$ is given by $`T_s=\frac{1}{2}mu_s^2`$. Notice that the right hand side of the equation above represents the difference between the incoming flux and the outgoing flux of compartment A. A similar expression has been chosen by Droz and Lipowski Droz  on phenomenological grounds. The model we study allows us to obtain self-consistently the temperatures $`T_A`$ and $`T_B`$, a feature which was not present in most of the studies dedicated to compartmentalized systemsextension . In fact, we can write:

$$\frac{\mathrm{d}[N_A(t)T_A(t)]}{\mathrm{d}t}=\frac{2}{\tau _s}[N_A(T_A+T_S)e^{T_s/T_A}N_B(T_B+T_S)e^{T_s/T_B}]2\gamma \omega _AN_AT_A+\frac{2}{\tau _b}N_A(T_bT_A)$$
(19)
$$\frac{\mathrm{d}[N_B(t)T_B(t)]}{\mathrm{d}t}=\frac{2}{\tau _s}[N_B(T_B+T_S)e^{T_s/T_B}N_A(T_A+T_S)e^{T_s/T_A}]2\gamma \omega _BN_BT_B+\frac{2}{\tau _b}N_B(T_bT_B)$$
(20)
Notice that the Arrhenius type of behavior of the transition rate $`\tau _se^{T_s/T_A}`$ results from the combination of two assumptions: a) the fact that only particles whose energy exceeds the threshold $`T_s`$ can overcome the barrier; b) the Gaussian ansatz for the velocity distribution functions. The latter ingredient, although very convenient for numerical work, could result too crude in some physically relevant situations, because the dynamics represent a severe probe of the extreme value statistics of the system. It is well known that the velocity distribution of driven and undriven granular assemblies might be characterized by fat velocity tails puglisi .

## III Mean field analysis

It is convenient to rewrite the coupled equations as:

$`{\displaystyle \frac{\mathrm{d}N_A(t)}{\mathrm{d}t}}`$ $`={\displaystyle \frac{1}{\tau _s}}[N_Be^{T_s/T_B}N_Ae^{T_s/T_A}]`$ (21a)
$`\begin{array}{cc}\hfill N_A{\displaystyle \frac{\mathrm{d}T_A(t)}{\mathrm{d}t}}& ={\displaystyle \frac{1}{\tau _s}}[2(N_AT_Ae^{T_s/T_A}N_BT_Be^{T_s/T_B})+(N_Ae^{T_s/T_A}N_Be^{T_s/T_B})(2T_sT_A)]\hfill \\ & 2\gamma \omega _AN_AT_A+{\displaystyle \frac{2}{\tau _b}}N_A(T_bT_A)\hfill \end{array}`$ (21b)
$`\begin{array}{cc}\hfill N_B{\displaystyle \frac{\mathrm{d}T_B(t)}{\mathrm{d}t}}& ={\displaystyle \frac{1}{\tau _s}}[2(N_BT_Be^{T_s/T_B}N_AT_Ae^{T_s/T_A})+(N_Be^{T_s/T_B}N_Ae^{T_s/T_A})(2T_sT_B)]\hfill \\ & 2\gamma \omega _BN_BT_B+{\displaystyle \frac{2}{\tau _b}}N_B(T_bT_B)\hfill \end{array}`$ (21c)
Let us observe that in the temperature equations there appear two types of fluxes: the first due to the unbalance of kinetic energies, and the second due to and the population difference. We identify the first with a heat conduction process and the second with a particle diffusion process.

One sees by inspection that the choice $`N_A=N_B=N^{}=N/2`$ and $`T_A=T_B=T^{}`$ represents a symmetric solution for all values of the control parameters. The granular temperature of the symmetric state, $`T^{}`$, is given by the nonlinear equation:

$$T^{}\left[1+\tau _b\sigma (1\alpha ^2)\frac{N^{}}{2V_A}\sqrt{\frac{T^{}}{m}}\right]=T_b$$
(22)
In order to ascertain the stability of such a symmetric solution we assume $`T_A=T^{}+\delta T_A`$, $`T_B=T^{}+\delta T_B`$ and $`N_A=N^{}+\delta N_A`$ ($`\delta N_B=\delta N_A`$) and expand the equations to linear order about the fixed point $`T^{},N^{}`$ with the result:

$$\delta \dot{N}_A=\frac{1}{\tau _s}\text{e}^{T_s/T^{}}\left[2\delta N_A+\frac{N^{}T_s}{(T^{})^2}(\delta T_A\delta T_B)\right]$$
(23)
$$\delta \dot{T}_A=\frac{1}{\tau _s}\text{e}^{T_s/T^{}}\left[2+\frac{T_s}{T^{}}+\left(\frac{T_s}{T^{}}\right)^2\right](\delta T_A\delta T_B)\left(3\gamma \omega ^{}+\frac{2}{\tau _b}\right)\delta T_A\frac{2}{N^{}}\left[\frac{1}{\tau _s}\text{e}^{T_s/T^{}}(T^{}+2T_s)+\gamma \omega ^{}T^{}\right]\delta N_A$$
(24)
$$\delta \dot{T}_B=\frac{1}{\tau _s}\text{e}^{T_s/T^{}}\left[2+\frac{T_s}{T^{}}+\left(\frac{T_s}{T^{}}\right)^2\right](\delta T_A\delta T_B)\left(3\gamma \omega ^{}+\frac{2}{\tau _b}\right)\delta T_B+\frac{2}{N^{}}\left[\frac{1}{\tau _s}\text{e}^{T_s/T^{}}(T^{}+2T_s)+\gamma \omega ^{}T^{}\right]\delta N_A$$
(25)
The eigenvalues of the associated $`3\times 3`$ matrix of coefficients give the three relaxation modes. By adding and subtracting the last two equations the linear system factorizes into a decoupled equation for the average value of the two temperatures ($`(\delta T_A+\delta T_B)/2`$) with eigenvalue $`\lambda _3=(3\gamma \omega ^{}+\frac{2}{\tau _b})`$ and a system of rank 2 involving the temperature difference and the occupation number. The more negative of the remaining eigenvalues (we call it $`\lambda _2`$) is of the same order as $`\lambda _3`$, while the second, say $`\lambda _1`$, is smaller in absolute value and vanishes at the “special” temperature, $`T_{cr}`$. This is obtained by solving the transcendental equation:

$$\tau _s\gamma \omega ^{}\left(\frac{T_s}{T_{cr}}\frac{3}{2}\right)\frac{\tau _s}{\tau _b}=2e^{T_s/T_{cr}}$$
(26)
Fig. 1 shows that above the temperature $`T_{cr}`$, $`\lambda _1`$ is negative, thus making the symmetric solution stable, whereas for temperatures below it $`\lambda _1`$ is positive. In such a case the symmetric fixed point is unstable and the solution flows away, to an asymmetric fixed point. Hereafter, we shall use the word “critical temperature” or “critical line” instead of “special” facciamoipedanti .

Let us also notice that when the typical collision time is much shorter than the characteristic times $`\tau _c`$ and $`\tau _b`$ the temperature $`T_{cr}2T_s/3`$, which represents an upper bound for such a quantity.

At this stage few comments are in order. The transition results from the competition between two effects: the diffusion due to the “thermal” agitation of the particles, i.e. the tendency to fill all the available space, and the dissipation of energy during collisions which favors clusterization.

When the external drive is sufficiently weak there appears a second fixed point. This is found by imposing the detailed balance:

$$\frac{N_B}{N_A}=\frac{e^{T_s/T_A}}{e^{T_s/T_B}}$$
(27)
and the energy balance

$$\frac{2}{\tau _s}(N_AT_Ae^{T_s/T_A}N_BT_Be^{T_s/T_B})2\gamma \omega _AN_AT_A+\frac{2}{\tau _b}N_A(T_bT_A)=0$$
(28)
$$\frac{2}{\tau _s}(N_AT_Ae^{T_s/T_A}N_BT_Be^{T_s/T_B})2\gamma \omega _BN_BT_B+\frac{2}{\tau _b}N_B(T_bT_B)=0$$
(29)
A good approximation to the highest eigenvalue in such a case reads:

$$\lambda _1=\frac{2}{3}e^{T_s/T^{}}(\omega ^{}\gamma (T_s/T_{}3/2)1/\tau _b2e^{T_s/T^{}}).$$

One observes two regions where the solutions are apparently metastable. The first is the low temperature region, where a detailed inspection shows that the symmetric state, in reality, is very weakly unstable. From fig. 1 we see that at low temperatures ($`T^{}0`$) the positive eigenvalue $`\lambda _1`$ vanishes exponentially as $`e^{T_s/T^{}}`$, thus indicating that the symmetric solution might appear stable if the observation time is finite. In fact, the inter-well diffusion for $`T_b0`$ is almost completely suppressed by the Arrhenius factor.

The second region where the system displays genuine (meta)stability occurs just above $`T_{cr}`$. A branch with $`N_AN_B`$ and $`T_AT_B`$ is observed by integrating nmerically the system (21) by means of an Euler scheme. These asymmetric stationary states are obtained by preparing the system in a subcritical configuration ($`T^{ast}<T_{cr}`$), and successively increasing the temperature above $`T_{cr}`$. The width of the region is about a ten per cent of $`T_{cr}`$. The observed hysteresis agrees qualitatively with the results relative to the model of ref. Droz .

The numerical solutions of the coupled equations is displayed in figs. 2 and 3. We observe that a perturbation about the symmetric solution is re-adsorbed for $`T^{ast}>T_{cr}`$, whereas for temperatures below $`T_{cr}`$ the perturbation grows initially at an exponential rate, before saturating about a finite value. Let us notice that at low temperatures, due to the presence of the Arrhenius factor, the saturation process occurs extremely slowly, a phenomenon which is not related to the slowing down which occurs only in the vicinity of $`T_{cr}`$.

In fig. 4 we display the line of critical points $`T_{cr}`$ as a function of the inelasticity $`(1\alpha )`$ for two different values of the total number of particles. Above the line the solutions are symmetric ($`N_A=N_B`$), below are asymmetric. Let us observe that for $`\alpha =1`$ one finds $`T_{cr}=0`$, because there is no clustering instability in systems of elastic particles. The smaller the value of $`\alpha `$ the higher the value of the critical temperature. We also notice that the slope of the critical line is non negative. When the total number of particles decreases, also the critical temperature decreases, because of the lower collision rate.

Figure 5 illustrates the behavior of the order parameter $`ϵ=|N_AN_B|/N`$ versus the heat bath temperature for two different values of the total number of particles.

Let us notice that when particles are added to the system the critical point shifts up-wards. This explains why the steady number of particles in the less populated compartment decreases when more particles are added Brey .

Finally, we comment that, within our perspective, the approach of Lohse et al. Lohse  and of Lipowski and Droz Droz  is equivalent to an adiabatic approximation for the temperature variables. In their cases, the two granular temperatures are assumed to be slaved by the occupation number variables. In other words, one postulates that the temperatures depend on the instantaneous values of the occupation numbers. We implemented this idea within our approach, but we found a large discrepancy with the previous approximation. In reality, the temperature and the occupation variable vary on the same time scale and no slaving principle seems to occur.

## IV Monte Carlo simulations

In deriving the evolution equations of the previous section we have implicitly stipulated that the values of the occupation number and of the kinetic energy per particle in each box have narrow distributions around their mean values and that these quantities are sufficient to characterize the state of the system. Do statistical fluctuations play any role or even modify the picture presented above? In the present section, which is not intended to be exhaustive, we solve the model of section I by means of the DSMC method. The motivation is twofold: we want to validate the picture and the approximations introduced above and unveil some phenomena which are not accounted for by a mean field description.

We simulated two ensembles of $`N_A`$ and $`N_B`$ particles, respectively, subject to a Gaussian forcing, viscous friction and inelastic collisions. In addition, the particles of high energy can change compartment with probability per unit time $`\tau _s^1`$.

The scheme consists of the following ingredients:

* time is discretized, i.e. $`t=ndt`$
* update all the velocities to simulate the random forcing and the viscous damping:

$$v_i^\alpha (t+dt)=v_i^\alpha (t)e^{\frac{dt}{\tau _b}}+\sqrt{T_b\left(1e^{\frac{2dt}{\tau _b}}\right)}W(t)$$
(30)
where $`W(t)`$ is a normally distributed deviate with zero mean and unit variance.
* At every time step and for both compartments, a collision step is performed, i.e. an adequate number of randomly chosen pairs of velocities are updated with the collision rule (1): the pairs are chosen with a probability proportional to their relative velocity and normalized in order to have a mean collision frequency as calculated in eq. (13).
* Place in the other compartment with probability $`dt/\tau _s`$ particles with kinetic energy greater than $`T_s`$.
* Change the time counter $`n`$ and restart.

To summarize, at every step each particle experiences a gaussian kick and receives energy from the bath, but dissipates a fraction of its kinetic energy by collision and by damping. With respect to previous DSMC simulations of granular gases particles can migrate to another compartment, whenever their energy exceeds a fixed threshold. No packing effects are included.

The parameters chosen in the simulation are $`T_s=1`$, $`V=100`$, $`\sigma =1`$, $`\tau _b=1`$, $`\tau _s=0.5`$, $`\alpha =0.7`$, whereas the total number of particles $`N`$ and $`T_b`$ have been varied.

First, we analyze the evolution of the order parameter $`|N_AN_B|/N`$, when the system is prepared in a symmetric configuration. The order parameter displays a behavior similar to that obtained by means of the mean field theory. In addition, it displays fluctuations around its asymptotic value. In fig. 6, an average over $`200`$ realizations of the evolution of the order parameter $`|N_AN_B|/N`$ is shown. The numerical simulations are in agreement with the results of the mean field theory. For $`T_b>0.45`$ the system remains substantially homogeneous. For $`T_b<0.45`$ the homogeneous state becomes unstable and the stable configuration that is reached is strongly asymmetric.

Interestingly, we observe a much slower growth (see inset) near the critical temperature, $`T_{cr}0.45`$. The same kind of phenomenon appears at low temperatures ($`T<T_s`$), because the transitions from one compartment to the other represent very rare events.

In figure 7 the temperature evolution $`T_A(t)`$ and $`T_B(t)`$ (after averaging over $`200`$ realizations and relabeling the compartments in such a way that $`A`$ is always the most populated) are displayed. The plateau values of $`T_A`$ and $`T_B`$ verify (with small deviations, not larger than $`10\%`$) the equation (27).

The velocity fluctuation can be appreciated by considering the velocity PDF’s of the two compartments fig. 8. Appreciable deviations from a Gaussian have not been detected observed in all cases considered (neither above, neither below, nor in the proximity of the critical temperature).

The order parameter, instead, shows less trivial fluctuations. A typical trajectory is shown in fig. 9. The associated distribution is shown in fig. 10. We observe that at high temperature the population distribution $`P(N_A)`$ is well approximated by a gaussian distribution, but displays two symmetric peaks about the central value $`N_A=500`$ when the system approaches the critical temperature. Finally, at low temperature we find two very narrow peaks, well separated and centered around $`N_A=950`$ and $`N_B=50`$.

Let us consider more closely the behavior of the order parameter near the critical point. According to mean field description of section II, at $`T_{cr}`$ the less negative eigenvalue of the matrix of evolution vanishes. The picture looks somehow different when microscopic fluctuations are accounted for. One observes a remarkable behavior in the order parameter evolution, namely in the earlier regime it grows in a power law fashion as $`(N_AN_B)^2t`$, and not exponentially. To explain this phenomenology we remark that the evolution of the order parameter near $`T_{cr}`$ can be assimilated to that of a particle undergoing a random walk. We assume that the early evolution of $`\mathrm{\Delta }N(t)=N_A(t)N/2`$ can be effectively described by the following equation:

$$\frac{\mathrm{\Delta }N(t)}{t}=\gamma (T_b)\mathrm{\Delta }N(t)+A(T_b)\eta $$
(31)
with $`\gamma `$ lesser than zero above the critical temperature and greater than zero under the critical temperature and $`\eta `$ a white gaussian noise. In this case

$$C(\tau ,t)=\mathrm{\Delta }N(t+\tau )\mathrm{\Delta }N(t)=\left[\mathrm{\Delta }N(0)\mathrm{\Delta }N(0)+\frac{A^2}{2\gamma }\right]e^{\gamma (2t+\tau )}\frac{A^2}{2\gamma }e^{\gamma \tau }$$
(32)
which means that if $`\gamma <0`$ (asymptotically in a time greater than $`1/|\gamma |`$) the equal time correlation function becomes:

$$\mathrm{\Delta }N(t)\mathrm{\Delta }N(t)\frac{A^2}{2\gamma }$$
(33)
or equivalently $`A=\sqrt{2\gamma \mathrm{\Delta }N^2}`$, where $`\mathrm{\Delta }N^2=lim_t\mathrm{}C(0,t)`$. On the other hand, when $`\gamma =0`$ (i.e. in correspondence of the critical point), the equal time correlation function displays a diffusive behavior:

$$\mathrm{\Delta }N(t)\mathrm{\Delta }N(t)A^2t$$
(34)
that means that the diffusion coefficient for the variable $`\mathrm{\Delta }N`$ is given by $`D=A^2/2`$.

Eventually in the late regime the growth saturates due to the onsert of non linear effects.

In order to verify the plausibility of eq. (31) we determine $`A`$ and $`\gamma `$ by extracting them from numerical measures of $`lim_t\mathrm{}C(\tau ,t)`$ and imposing the asymptotic behavior $`\frac{A^2}{2\gamma }e^{\gamma \tau }`$. In fig. 11 we display the values of the constants obtained from simulations above the critical temperature. We also added to the plot the value of the diffusion coefficient, $`D=A^2/2`$, as measured at the critical temperature.

Interestingly, when the characteristic time $`\tau _s`$ decreases (namely $`\tau _s=0.1`$) a new phenomenon appears: the temperature of the less populated compartment exhibits an initial regime, during which it increases, due to the decrement of the number of collisions the particles experience, followed by a second regime where the temperature instead decreases and eventually reaches a temperature lower than the temperature of the more populated container. This kind of anomaly is due to the fact that the fastest and more energetic particles are removed causing a large negative fluctuation of the associated average kinetic energy. The positive energy flux provided by the heat bath is not sufficient to compensate the energy loss due to the removal of the fast particles. This phenomenon is not observed in the mean field model, because it results from the large non-gaussian fluctuations of the velocity PDF. We checked the fourth cumulant and observed that whereas the cumulant relative to the populated compartment corresponds to a gaussian, the cumulant of the small population strongly fluctuates. Experimentally it might be hard to observe such a phenomenon, which is probably an artifact of the model, because the limit of small $`\tau _s`$ is not very realistic.

Finally, our schematization of the compartmentalized granular gas recalls the Gibbs ensemble method of equilibrium statistical mechanics, whereby it is possible to study first order phase coexistence without interfaces Frenkel .

## V Conclusions

To summarize, we introduced a model for a compartmentalized granular gas which allows to bridge between the microscopic level to the hydrodynamic level. In the first part of the paper we have derived a Fokker-Planck-Boltzmann description starting from the stochastic evolution of the particles coordinates. Next, by employing a gaussian ansatz for the velocity distribution function, we have obtained a closed set of equations for the slowly varying fields, namely the granular temperatures and occupation numbers of each compartment. Let us comment that with respect to existing flux models, our approach treats analytically the granular temperature on equal footing as the occupation variables. The solution of the resulting equations shows the existence of two different regimes: at strong shakings the populations in the two compartments are symmetric and for weak shakings the two populations are asymmetric together with their granular temperatures. A critical point separates these two behaviors. The dynamics has been characterized in the various regimes. In the last part of the paper we have solved the full model by means of DSMC. The presence of stochastic fluctuations leads to properties not observed in a mean field description. These are the presence of critical fluctuations, of anomalies in the dynamics of the populations and in population fluctuations.

## VI Acknowledgments

U.M.B.M. acknowledges the support of the Cofin MIUR “Fisica Statistica di Sistemi Classici e Quantistici”. A.P. thanks INFM Center for Statistical Mechanics and Complexity.


format: text

text:
# The Munich Near-Infrared Cluster Survey (MUNICS) – II. The 𝐾-band luminosity function of field galaxies to 𝑧∼1.2

## 1 Introduction

The luminosity function (LF) is the most basic statistic used to study galaxy populations and their evolution. Its dependence on wavelength and look-back time provides important constraints on the evolution of the global properties of the galaxy population.

Recently, progress in measuring the optical galaxy LF was made both locally and at higher redshift. Locally, results from the 2dF survey (Folkes et al., 1999; Madgwick et al., 2002) and 2MASS (Kochanek et al., 2001) provide much improved measurements of the LF as a function of galaxy morphology, environment, and wavelength. In spite of some controversy regarding the normalization and the very faint end slope of the local LF, it is a well established result that the LF depends on galaxy type, although the correlations are not always very tight. Generally, the faint end is populated by galaxies of smaller mass, later morphology, bluer colors and later spectral type, and stronger line emission. The bright end is dominated by early-type spirals and ellipticals, and the very bright end by giant ellipticals (e.g. Marzke et al., 1994; Bromley et al., 1998; Madgwick et al., 2002).

A number of studies of the luminosity function of galaxies at $`0z1`$ using optically selected redshift surveys have been published in recent years (e.g. Efstathiou et al., 1988; Lilly et al., 1995b; Heyl et al., 1997; Lin et al., 1997; Liu et al., 1998; Ratcliffe et al., 1998; Lin et al., 1999). These surveys (with samples of typically hundreds of objects) have consistently found similar trends in the evolution of the rest-frame $`B`$-band field galaxy LF. The main result was the contrast between the rapid evolution of the blue, star-forming sub-population and the mild change in the redder, early-type population. Simultaneously, first results at $`z3`$ became available through the study of Lyman-break galaxies (Shapley et al., 2001).

The ground-breaking CFRS (Lilly et al., 1995a) found that the luminosity function in the rest-frame $`B`$-band of the blue field population brightens by roughly 1 mag to $`z1`$ and also becomes steeper at the faint end at $`z0.5`$ (see also Heyl et al., 1997). In contrast, the red population was found to show very little change in either number density or luminosity (Lilly et al., 1995b). Lin et al. (1999) attempted to discriminate explicitly between number-density and luminosity evolution in the CNOC2 sample consisting of roughly 2000 galaxies. They found that the early-type population shows positive luminosity evolution (1.6 mag) which is nearly compensated by negative density evolution (factor of 0.5), so that there is little net change in their overall LF. The intermediate-type population shows positive luminosity evolution (0.9 mag) plus weak positive density evolution (factor of 1.7), resulting in mild positive evolution in its luminosity density. The amount of luminosity evolution in the early-type and intermediate-type populations was found to be consistent with expectations from models of passive evolution of their stellar populations. In contrast, the late-type population is best fit by a strong increase in number density at high redshift (factor of 4.1), accompanied by only little positive luminosity evolution (0.2 mag). The overall $`B`$-band luminosity density of late-type objects was found to increase rapidly in both the CFRS and CNOC2 samples, while the luminosity density of early-type objects was found to be nearly constant.

In contrast to the blue and optical wave-bands, the near infrared $`K`$-band light of galaxies is much less affected by dust, much less sensitive to ongoing star-formation and much less dependent on galaxy type. Furthermore, the near-IR luminosity is as a reasonable tracer of the stellar mass of a galaxy, since the near-IR light is dominated by the light of old and evolved stars. Therefore, by investigating the $`K`$-band properties of galaxies as a function of redshift, we may hope to be able to move from a picture dominated by the evolution of star formation to one which focuses on the assembly history of mass in these systems – one of the most fundamental predictions of cold dark matter based structure formation theories.

Most near-IR selected surveys so far have been either small in size or shallow (e.g. Glazebrook et al., 1995; Gardner et al., 1997; Cowie et al., 1996; Szokoly et al., 1998). Furthermore, optical followup spectroscopy of near-IR selected samples to $`z1`$ is difficult because of the large range of optical to near-IR (e.g. $`RK`$) colors of galaxies, although spectroscopic samples of a few hundred objects are now available (Feulner et al., 2003; Pozzetti et al., 2003; Cimatti et al., 2002). If any change in the LF was seen at all, only mild luminosity evolution at $`z>0.5`$ has been significantly detected.

In this work, we aim at measuring the rest-frame $`K`$-band luminosity function of galaxies and its evolution with redshift in the range $`0.4<z<1.2`$ using the $`K`$-band selected, BVRIJK multicolor data set of the MUNICS project (Drory et al., 2001b). The sample contains more than 5000 objects and covers a large solid angle. We use spectroscopically calibrated photometric redshifts (using more than 500 spectroscopic redshifts) to estimate distances and spectral energy distribution (SED) fitting techniques to extrapolate the SED to rest-frame $`K`$.

This paper is organized as follows. First, we briefly introduce the galaxy sample we use in this work in Sect. 2. Next, in Sect. 3, we describe the photometric redshift estimation technique along with the construction of a set of semi-empirical SEDs best matching the SEDs present in the sample. In Sect. 4 we discuss the construction of the rest-frame $`K`$-band luminosity function from photometric redshifts and assess the influence of errors in $`z`$ by means of Monte-Carlo simulations. The resulting luminosity function is presented in Sect. 5, and its evolution is investigated in Sect. 6. Sect. 7 discusses the results and Sect. 8 summarizes this work.

We assume $`\mathrm{\Omega }_\mathrm{M}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ throughout this paper. We write Hubble’s Constant as $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, using $`h=0.65`$ unless the quantities in question can be written in a form explicitly depending on $`h`$.

## 2 The galaxy sample

MUNICS is a wide-area, medium-deep, photometric and spectroscopic survey selected in the $`K`$ band, reaching $`K19.5`$. It covers an area of roughly one square degree in the $`K`$ and $`J`$ bands with optical follow-up imaging in the $`I`$, $`R`$, $`V`$, and $`B`$ bands in 0.4 square degrees. Drory et al. (2001b, hereafter MUNICS I) discusses the field selection, object extraction and photometry in the $`K,J,I,R`$, and $`V`$ bands. The $`B`$ band imaging was completed more recently and the data were processed in an analogous way as described in MUNICS I. Detection biases, completeness, and photometric biases of the MUNICS data are analyzed in detail in Snigula et al. (2002, hereafter MUNICS IV).

The MUNICS photometric survey is complemented by spectroscopic follow-up observations of all galaxies down to $`K17.5`$ in 0.25 square degrees, and a sparsely selected deeper sample down to $`K19`$. It contains more than 550 secured redshifts. The spectra cover a wide wavelength range of $`40008500`$Å at $`13.2`$Å resolution, and sample galaxies at $`0<z<1`$. These observations are described in detail in Feulner et al. (2003, hereafter MUNICS V).

The galaxy sample used in this work is a subsample of the MUNICS survey Mosaic Fields (see MUNICS I), selected for best photometric homogeneity, good seeing, and similar depth. Furthermore, in each of the survey patches, areas close to the image borders in any passband, areas around bright stars, and regions suffering from blooming are excluded. The subsample covers 0.28 square degrees in $`B`$, $`V`$, $`R`$, $`I`$, $`J`$, and $`K`$. Table 1 lists the Mosaic Fields used in the subsequent analysis. It is almost identical to the sample used in Drory et al. (2001a, hereafter MUNICS III), except for the additional $`B`$-band imaging data (and hence refined photometric redshifts).

Stars are identified following the procedure described in MUNICS I, adding a color criterion in the $`JK`$ vs. $`VI`$ plane to exclude faint stellar sources which cannot be separated from galaxies morphologically (see Fig. 10 in MUNICS I). This color criterion may also exclude $`z0.2`$ compact blue galaxies. Such galaxies are very unlikely to be present in the $`K`$-selected sample, given our magnitude limit. We will restrict our analysis of the luminosity function to $`z>0.4`$ such that this is anyway not a problem.

To make sure that our star-galaxy separation does not cause any systematic biases in our analysis, we double check the method by including stellar SEDs in the template library used to estimate photometric redshifts (see below). We find that the objects identified as stars are indeed better fit by a pure stellar SED than by any galactic SED (see right panel in Fig. 5). Most exceptions (a few percent of the stars) are most probably binaries. Using only the SED fit to discriminate stars from galaxies does not change any of the results of this work.

## 3 Photometric redshifts

Photometric redshifts are derived using the method presented in Bender et al. (2001) and Bender (2003). We only briefly outline the procedure in this section. The method is a template matching algorithm rooted in Bayesian statistics closely resembling the method presented by Benítez (2000). However, instead of relying on a predetermined set of template SEDs, semi-empirical templates matching the photometric properties of the sample are used. The templates are derived by fitting stellar population models of Maraston (1998) of different age and dust extinction and Kinney-Calzetti (Kinney et al., 1996) spectra to combined broad-band energy distributions of MUNICS and FDF (Heidt et al., 2003) galaxies having spectroscopic redshifts. In this way, representative galaxy templates of mixed stellar populations (variable age, metallicity, and dust extinction) optimized for the MUNICS dataset are obtained.

The total spectroscopic sample is divided into 2 groups of objects. The first subsample (all objects with redshifts in the field S2F1) is used for constructing SED templates, the second subsample (all other objects with spectroscopic redshifts) is used for comparing spectroscopic and photometric redshifts and thereby calibrating the SED library.

The observed-frame apparent magnitudes of objects in the first subsample are transformed to rest-frame redshift zero, and fitted by an initial set of stellar population synthesis models. Objects best fitting the same model are grouped together. Since each group will contain objects from a variety of redshifts, a densly sampled SED from the broad band photometry of these objects is obtained. This procedure is illustrated in Fig. 1.

This initial set of SEDs is used to determine photometric redshifts for the total sample of objects having spectroscopic redshifts. The photometric redshifts are compared to the spectroscopic ones, and, additionally, the same de-redshifting procedure is applied to the spectroscopic sub-sample not used for the initial construction of the SEDs, only that now we group the objects by the SED that gave the best fit during the determination of the photometric redshift. Using this comparison, deficiencies in the set of SEDs can be identified as those become apparent through systematic offsets between the de-redshifted magnitudes and the SED templates. This is the case since such deficiencies lead to a wrongly determined photometric redshift and therefore to the assignment of a wrong SED.

This procedure is repeated with a refined set of SEDs, by changing SEDs, abolishing some and adding others, until a satisfactory library of template SEDs is found. Fig. 2 shows the final template SED library used to derive photometric redshifts in what follows.

In Fig. 3 we plot the difference between spectroscopic and photometric redshift vs. spectroscopic redshift for the subsample used to construct the SED templates and for the rest of the sample which was used to test the procedure. There is no apparent difference between these two error distributions.

Finally, Fig. 4 compares photometric and spectroscopic redshifts for all $`500`$ objects within five MUNICS Mosaic Fields and shows the distribution of redshift errors. The typical scatter in the relative redshift error $`\mathrm{\Delta }z/(1+z)`$ is 0.055. The mean redshift bias is negligible. The distribution of the errors is roughly Gaussian. There is no visible difference between the distributions among the survey fields. Although this performance is encouraging, it is important to say that the spectroscopic data become sparse at $`z0.6`$ and there are only two spectroscopic redshifts at $`z>1`$.

To illustrate the use of photometric redshifts, we show two instructive examples of photometric redshift determinations of galaxies and the identification of an M-star in Fig. 5. These examples help to understand how the technique works and the uncertainties involved.

Firstly, a spiral-like system at redshift around unity. Here the Balmer break is redshifted beyond the $`R`$-band filter, and only one SED contributes significantly to the global peak in the redshift probability distribution. The rather broad probability distribution in redshift is due to the fact that the Balmer break is in the rather large sampling gap between the $`R`$ and $`I`$ bands, and therefore its position is not very well determined. Note that this object, although very bright in the optical, is very faint in the $`K`$ band, and therefore probably a rather low-mass system.

Secondly, an early-type object at redshift $`z0.82`$. The redshift determination can be regarded as quite secure although there are competing SEDs around the global peak of the redshift probability function. In this case, the redshift is determined by the 4000Å break and the steep decline in flux bluer thereof. Because the object is undetected in $`B`$ and only barely detected in $`V`$, the rest-frame UV and blue slopes of the spectrum are not firmly determined and hence slightly differing effective ages are giving reasonable fits. This uncertainty in age or type translates into an uncertainty in redshift. Therefore there are competing SEDs at similar but not identical redshifts contributing to the total redshift probability distribution.

Common to both above examples is that the position in redshift of the global maximum of the total redshift probability distribution is always compatible with the redshift one would derive by looking at the probability distribution of the most likely SED.

The third panel in Fig. 5 shows an object best fit by the SED of an M2 star. We have used a criterion based on a comparison of the best $`\chi ^2`$ for redshifted galaxy template SEDs and stellar SEDs to discriminate between star and galaxies, to test the robustness of the morphological and color based star-galaxy separation presented in MUNICS I, hoping to be able to improve the procedure at faint magnitudes. The main problem hereby is the lack of reliable SEDs for cool stars with the necessary broad wavelength coverage. Using only this SED based star-galaxy separation instead of the morphological method does not change the results of this work. Note that we present an independent measure of the reliability of our star-galaxy separation in MUNICS V, where we test it against blind spectroscopy.

It is interesting to compare the photometric redshifts obtained using the present six-color photometry to the ones used in MUNICS III, where we did not have B-band imaging yet and the spectroscopic calibration sample was smaller (310 objects as opposed to 550 now). This comparison is shown in Fig.6. In spite of overall reasonable agreement, there are some important differences. The most notable one is that objects tend to scatter to lower redshifts if $`B`$-band photometry is missing. Objects are scattered preferably to the region in redshift below which the 4000Å break enters the observed filter set (because this peak in the $`z`$-probability function cannot be excluded by the missing blue photometry). Another set of objects is redistributed from $`z1.2`$ (six-colors) to lower redshifts (five-colors). These objects all have blue colors and rather late type SEDs. Both effects are small in numbers, affecting $`3`$ per cent of the sample. Therefore they do not show up in the spectroscopic calibration and are difficult to eliminate systematically. The influence of such an effect on the luminosity function can be significant although the number of affected objects is small if rare high-luminosity objects are preferably affected. This is discussed below in Sect. 5.

In Fig. 7 we plot absolute $`K`$-band magnitude vs. photometric redshift and the redshift distribution of the total sample discussed here, containing 5132 galaxies. The redshift distribution peaks around $`z0.5`$ and has a tail extending to $`z3`$. An analytical fit of the form

$$\frac{dN}{dz}=\frac{\beta z^2}{\mathrm{\Gamma }(3/\beta )z_0^3}e^{(z/z_0)^\beta }$$
(1)
is also shown. The best-fitting values are $`z_0=0.101,\beta =0.746`$.

Finally, we plot the cumulative redshift distribution of bright galaxies in the $`K`$-band in Fig. 8. We select galaxies with apparent $`K`$-band magnitudes in the interval $`16<m_K<18`$, where the present sample is definitively complete. The figure also compares the data in the single survey fields to the sample average. We find that cosmic variance is significant on scales of $`100`$ square arc minutes, the typical size of previous infrared selected surveys.

## 4 Computing the luminosity function

### 4.1 The $`V_{\mathrm{max}}`$ method

We estimate the luminosity function, $`\mathrm{\Phi }(M)dM`$, the comoving number density of galaxies with absolute magnitude in the range $`[M,M+dM)`$, using the $`V_{\mathrm{max}}`$ formalism (Schmidt, 1968) to account for the fact that some fainter galaxies are not visible in the whole survey volume.

Each galaxy in a given redshift bin $`[z_l,z_h)`$ contributes to the number density an amount inversely proportional to the volume in which the galaxy is detectable in the redshift bin given all relevant observational constraints:

$$V_i=d\mathrm{\Omega }_{z_l}^{\mathrm{min}(z_h,z_{max})}\frac{dV}{dz}𝑑z,$$
(2)
where $`dV/dz`$ is the comoving volume element, $`d\mathrm{\Omega }`$ is the survey area, $`z_{max}`$ is the maximum redshift at which galaxy $`i`$ having absolute magnitude $`M_{K,i}`$ is still detectable given the limiting apparent magnitude of the survey and the galaxy’s SED (the best-fit SED from the photometric redshift determination in our case).

Additionally, the contribution of each galaxy $`i`$ is weighted by the inverse of the detection probability, $`P(m_{K,i})`$, where we assume that the detection probability is independent of the galaxy type and can be approximated by that of point-like sources. We only include objects with $`P(m_{K,i})>0.75`$, such that this correction is always small. We have checked that this correction does not bias our results by comparing to what we get for higher completeness limits. Also, the results of completeness simulations discussed in MUNICS IV demonstrate that although there exist profile-dependent surface-brightness selection biases, these are under control for redshifts up to $`z1`$, and that the onset of incompleteness is roughly independent of profile type.

The comoving number density of objects in a given absolute magnitude interval and redshift bin is finally calculated as

$$\mathrm{\Phi }(M)dM=\underset{i}{}\frac{1}{V_{\mathrm{max}}^i}\frac{1}{P(m_{K,i})}dM,$$
(3)
where the sum is to be taken over all objects $`i`$ in the bin.

The advantage of the $`V_{\mathrm{max}}`$ method is that it is non-parametric, i.e. no assumption on the form of the LF is made and the estimate of the LF in each redshift and magnitude bin is independent. There is no need to compute the normalization of the LF separately. It is quick and easy to compute. The major disadvantage is that this method is sensitive to clustering in the sample. In our case we might hope that this is not a major concern since we average over a large volume and probe multiple independent lines of sight.

Takeuchi et al. (2000) performed a systematic comparison of different estimators for the luminosity function, finding that the $`V_{\mathrm{max}}`$ estimator yields a completely unbiased result if there is no inhomogeneity. Earlier claims that the $`V_{\mathrm{max}}`$ estimator is biased even without clustering were not confirmed. The $`V_{\mathrm{max}}`$ estimator was also found to give consistent results with other statistical estimators analyzed, despite of its sensitivity on large scale structure.

The last missing ingredient is the absolute $`K`$-band magnitude. It is computed by extrapolation from the observed-frame colors, using the best-fit SED from the photometric redshift code. The near-IR slopes of the SEDs are fairly uniform across galaxy types (the $`K`$-band $`k`$-corrections are small and almost type-independent) and so the uncertainty introduced by this extrapolation is rather small, of the order of $`\mathrm{\Delta }M_K0.1`$ mag in the mean, and thus small compared with the uncertainty of the total rest-frame $`K`$-band magnitude coming from the uncertainty in the distance.

### 4.2 Monte-Carlo simulations

The use of photometric redshifts, in general, can introduce systematic errors in the derived galaxy distances and therefore in the luminosity function which we want to derive in what follows. Therefore, we perform Monte-Carlo simulations to test the robustness of the LF estimate under the conditions imposed by the use of photometric redshifts instead of spectroscopic distances. We generate mock MUNICS-like catalogs following a non-evolving LF and compare the output of the LF estimate as described above with the input LF. Note that we can only test the effects of random errors in the distances and not the influence of systematically wrong distances as those shown in Fig. 6.

The simulated LF follows the form of th local 2MASS $`K`$-band LF of Kochanek et al. (2001). We assume that there is no redshift evolution in brightness or number density. The extrapolation of the flux to the rest-frame $`K`$-band was done using the mean $`k`$-correction of the models used for determining photometric redshifts (Fig. 2).

To investigate the influence that errors in the photometric redshifts have on recovering the form of the LF, we simulate four data sets with errors $`\mathrm{\Delta }z`$ drawn from a Gaussian of width $`0.02`$, $`0.055`$ (the value found from the comparison with spectroscopic redshifts in Fig. 4), $`0.1`$, and $`0.2`$. The simulations are repeated 10 times, and the LF is extracted in four redshift bins, $`0.4<z<0.6,0.6<z<0.8,0.8<z<1.0,`$ and $`1.0<z<1.2`$. The results are shown in Fig. 9.

The overall impression from Fig. 9 is that the form is well reproduced, and that, as expected, the influence of the redshift uncertainty is larger at small redshifts where the relative influence on the distance is largest. The form of the LF only becomes significantly biased when the redshift errors exceed $`\mathrm{\Delta }z0.1`$, becoming a significant fraction of the redshift bin size. There, the well-known effect of objects being preferably scattered away from $`L^{}`$, producing too many bright objects and too few lower luminosity objects becomes apparent. For the value $`\mathrm{\Delta }z=0.055`$ (upper right panel in Fig. 9), the form is only slightly affected in the lowest redshift bin, $`0.4<z<0.6`$.

From the last said we conclude that the shape of the luminosity function is not significantly biased as long as the photometric redshifts scatter symmetrically around the true redshifts (see also Subbarao et al., 1996) and their errors are (significantly) smaller than the bin size in $`z`$ over which the luminosity function is averaged.

### 4.3 Influence of the incompleteness corrections

The corrections for incompleteness involve the volume correction, $`1/V_{\mathrm{max}}`$, and the correction for the detection incompleteness at fainter magnitudes, $`1/P(m_{K,i})`$.

Given our magnitude limits, both corrections do not play any role at magnitudes brighter than $`M_K^{}23.5`$. At fainter magnitudes both corrections contribute significantly. The $`V_{\mathrm{max}}`$ term becomes dominant at the faintest levels (as we exclude objects with $`P(m_{K,i})<0.75`$).

In the lowest redshift bin the $`V_{\mathrm{max}}`$ term is ill-determined resulting in a steepening in the LF’s very faint end. The reason is the redshift error associated with photometric redshifts and that the approximation of the visibility function as that of a point source breaks down at low redshifts where the objects are well-resolved. The steepening at of the LF can be explained in terms of low luminosity galaxies having true redshifts just lower than the bin’s lower bound, which are scattered upwards in $`z`$ by the photometric redshift estimate and just make it into the bin. As those galaxies are close to the detection limit, the $`V_{\mathrm{max}}`$ test predicts them to be detectable only in a very small volume of the bin and therefore the $`V_{\mathrm{max}}`$ correction becomes large (factors of 5 and above), producing the observed behavior. Therefore our data do not allow us to place constraints on the faint end slope of the LF, even in our lowest redshift bin.

## 5 Results

In this section we present the results on the evolution of the rest-frame $`K`$-band luminosity function to $`z1.2`$.

In Fig. 10 we plot the rest-frame $`K`$-band luminosity function obtained from the present MUNICS dataset. The plot shows the uncorrected values as well as the $`V_{\mathrm{max}}`$ and incompleteness corrected ones. For comparison, we also show the local $`z=0`$ 2MASS $`K`$-band LF (Kochanek et al., 2001) and a Schechter function evolved with $`z`$ using the best-fit values for the evolution parameters $`\mu =1/\mathrm{\Phi }^{}(0)d\mathrm{\Phi }^{}/dz`$ and $`\nu =dM^{}/dz`$ from Sect. 6. Note that these are not independent Schechter fits in each redshift bin but rather a global estimate of the change in $`\mathrm{\Phi }^{}`$ and $`L^{}`$ with redshift (see below).

Taken at face value, it is apparent that the total rest-frame $`K`$-band LF evolves only mildly to $`z=1.2`$. This is in qualitative agreement with the finding of Cowie et al. (1996) from a much smaller spectroscopic sample, where the authors found that the $`K`$-band LF in the redshift bin $`0.6<z<1`$ was consistent with the one at $`0<z<0.2`$ (using an Einstein-de-Sitter cosmology).

This is further demonstrated by the total $`K`$-band luminosity density shown in Fig. 11. Also here there is only mild evolution with redshift within our sample, with the luminosity density scaling as $`d\mathrm{log}\rho _L/dz=0.24`$.

Note that to compute the previous quantity, we have to assume something about the luminosity in objects at the faint end of the luminosity function which are beyond our magnitude limits. Since we have no handle on the slope of the faint end from our data, we assume that the local slope of $`\alpha =1.09`$ holds. We use the fraction of the total luminosity sampled by the data to the completeness limit $`M_{\mathrm{min}}`$,

$$C_\mathrm{\Phi }=\frac{_{\mathrm{}}^{M_{\mathrm{min}}}\mathrm{\Phi }(M)𝑑M}{_{\mathrm{}}^+\mathrm{}\mathrm{\Phi }(M)𝑑M},$$
(4)
to correct for this effect in Fig. 11. This correction is found to be smaller than a factor of 2 at all redshifts.

Again, we see clearly see the signature of cosmic variance in the luminosity density measurements in the single Mosaic Fields compared with the sample average. The values in the individual fields scatter by $`0.1`$ to $`0.15`$ in $`\mathrm{log}\rho _L`$.

We can compare the luminosity density values obtained by summing over the data with what is expected if the LF is well-represented by a Schechter function, in which case the luminosity density is simply given by $`\rho _L=\mathrm{\Phi }^{}L^{}\mathrm{\Gamma }(2+\alpha )`$. Using the best values for the evolution parameters $`\mu `$ and $`\nu `$ from Sect. 6, these expected values are shown in Fig. 11 as a solid line.

Our result on the evolution of the $`K`$-band LF is in marked contrast to the results obtained for the $`B`$-band LF in the same redshift range from the CFRS and CNOC surveys (Lilly et al., 1995b; Lin et al., 1999) which consistently find an increase in the $`B`$-band luminosity density driven mainly by late spectral types.

At least some of this increase is due to fainter late-type galaxies at (observer’s frame) $`B22`$ mag which first appeared in faint blue number counts (see review by Ellis, 1997). For example, Cowie et al. (1996) found that at $`B24`$, the population of galaxies is a mixture of normal galaxies at modest redshifts and a population of galaxies with a wide range of masses undergoing rapid star formation which are spread out in redshift from $`z=0.2`$ to at least $`z=1.7`$. The remaining part of the increase in blue luminosity is due to luminosity evolution of normal field spirals and early-type systems.

On the other hand, luminosity evolution is inevitable if the $`K`$-band light traces the underlying old stellar population of galaxies, since these stars obviously become younger with look-back time. For the $`K`$ band, pure luminosity evolution predicts brightening by $`0.6`$ mag as the average mass-to-light ratio evolves away from its local value (see MUNICS 3; Pozzetti et al., 1996).

In MUNICS III we have shown a luminosity function based on the VRIJK photometry available at the time. Although we did not present a quantitative analysis of the luminosity function evolution in that paper, it appeared that the LF is consistent with the local shape over the entire range in $`z`$, while we do see departure from the local shape in the present analysis. The explanation is found to be a weak excess of objects at high $`z`$ in the present BVRIJK data (see Fig. 6) discussed above, which is enough to explain the brightening seen at the high-luminosity end. A comparison of the luminosity density in the limits adopted in MUNICS III (Fig. 3) indeed yields only small differences, of the order of 5 to 10 per cent. We wish to emphasize that a non-evolving LF does not imply a non- evolving galaxy population (at least as long as we can only measure the bright ($`L>L^{}`$) part of the LF). Indeed, MUNICS III argued that density evolution is counter balancing luminosity evolution, leading to a constant overall LF. In the current analysis based on broader wavelength coverage and better spectroscopic calibration, the balance between luminosity and density evolution is altered to slightly more luminosity evolution compared to MUNICS III, but the principal results remain the same as will be shown in a detailed quantitative analysis below.

## 6 Likelihood analysis of luminosity function evolution

We define the amount of number density and luminosity evolution in our sample as the change of $`\mathrm{\Phi }^{}`$ and $`M^{}`$ with redshift, respectively, in the Schechter parametrization of the LF. To quantify this, we apply two different estimators for the evolution of the LF with redshift, one operating on the data binned in magnitude and redshift, and one operating solely on individual objects, i.e. unbinned individual magnitude and redshift $`\{M_i,z_i\}`$ pairs.

In both cases we use the Schechter (1976) parametrization of the luminosity function,

$$\mathrm{\Phi }(L)dL=\frac{\mathrm{\Phi }^{}}{L^{}}\left(\frac{L}{L^{}}\right)^\alpha \mathrm{exp}\left(\frac{L}{L^{}}\right)dL$$
(5)
where $`L^{}`$ is the characteristic luminosity, $`\alpha `$ the faint-end slope, and $`\mathrm{\Phi }^{}`$ the normalization of the luminosity function. The corresponding equation in absolute magnitudes reads

$`\mathrm{\Phi }(M)dM`$ $`=`$ $`{\displaystyle \frac{2}{5}}\mathrm{\Phi }^{}\mathrm{ln}10\times \mathrm{\hspace{0.25em}10}^{0.4(M^{}M)(1+\alpha )}`$
$`\times `$ $`\mathrm{exp}\left(10^{0.4(M^{}M)}\right)dM.`$
To estimate the rate of evolution in this parametrization with redshift, we define evolution parameters $`\mu `$ and $`\nu `$ as follows:

$`\mathrm{\Phi }^{}(z)`$ $`=`$ $`\mathrm{\Phi }^{}(0)\left(1+\mu z\right),`$
$`M^{}(z)`$ $`=`$ $`M^{}(0)+\nu z,\mathrm{and}`$ (6)
$`\alpha (z)`$ $`=`$ $`\alpha (0)\alpha .`$
And therefore,

$`\mu `$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Phi }^{}(0)}}{\displaystyle \frac{d\mathrm{\Phi }^{}(z)}{dz}}`$
$`\nu `$ $`=`$ $`{\displaystyle \frac{dM^{}(z)}{dz}}.`$ (7)
Note that the faint end slope of the luminosity function cannot be determined from our data, thus we leave the faint-end slope $`\alpha `$ of the Schechter luminosity function fixed.

For the values of $`M^{}(0)`$ and $`\mathrm{\Phi }^{}(0)`$ we use the values given by Kochanek et al. (2001) for 2MASS,

$`\mathrm{\Phi }^{}(0)`$ $`=`$ $`1.16\times 10^2h^3\mathrm{Mpc}^3`$
$`M^{}(0)`$ $`=`$ $`23.39+5\mathrm{log}h`$
$`\alpha `$ $`=`$ $`1.09.`$
### 6.1 A $`\chi ^2`$-based approach using binned data

First, we apply a $`\chi ^2`$-based estimator of the redshift evolution of the near-infrared luminosity function. The same estimator is applied in MUNICS V to a smaller but purely spectroscopic sample.

To quantify the redshift evolution of $`\mathrm{\Phi }^{}`$ and $`M^{}`$ we compare our luminosity function data in all redshift bins with the local Schechter function evolved according to equation (6) to the appropriate redshift. We do this for a grid of values of $`\mu `$ and $`\nu `$, and calculate the value of $`\chi ^2`$ for each grid point according to

$$\chi ^2(\mu ,\nu )=\frac{1}{n}\underset{i=1}{\overset{N}{}}\frac{\left[\varphi (M_i,z_i)\mathrm{\Phi }_{\mu \nu }(M_i,z_i)\right]^2}{\sigma _i^2},$$
(8)
where the sum is taken over all bins $`i`$, $`\varphi (M_i,z_i)`$ is the measured value of the luminosity function at median redshift $`z_i`$ in the magnitude bin centered on $`M_i`$, and $`\mathrm{\Phi }_{\mu \nu }(M,z)`$ is the value of the local Schechter function at magnitude $`M`$ evolved according to the evolution model defined in equation (6) to the redshift $`z`$. Furthermore, $`\sigma _i`$ is the RMS error of the measured luminosity function value, and $`n`$ is the number of free parameters, i.e. the number of data points used minus the number of parameters derived from the fitting, 23 in our case.

We use the data binned as presented in Fig. 10, spanning $`0.4<z<1.2`$. However, we exclude all data points with a total correction factor (incompleteness and $`V/V_{\mathrm{max}}`$ correction) larger than two. The errors of the measured LF, $`\sigma _i`$, are calculated from the number of objects in the bin and assumed to follow Poisson statistics.

The results of this approach are shown in the left panel of Fig. 12 and the most likely values of $`\mu `$ and $`\nu `$ are listed in Table 2, along with the results from the analysis of the spectroscopic sample in MUNICS V.

### 6.2 A likelihood-based approach using unbinned data

The luminosity function is very steep at its bright end. When the number of very bright (and therefore rare) galaxies in the sample is small, binning of these data in magnitude must be taken with skepticism. As the distribution of objects within any bin will not be uniform and we are dealing with small numbers of objects, the choice of bin centers and widths might influence the outcome of the analysis. Hence we develop a maximum likelihood based method avoiding having to bin the data to determine the values of the evolution parameters $`\mu `$ and $`\nu `$.

Maximum likelihood methods (using binned and unbinned data) have been used in the past by many authors to construct luminosity functions and derive their parameters in various parametric and non-parametric ways (see, e.g., Sandage et al., 1979; Efstathiou et al., 1988; Lin et al., 1996, and references therein). These methods usually suffer from the inconvenience that the normalization has to be determined independently from the shape of the LF. The method presented here allows the shape and normalization (evolution) parameters to be determined simultaneously.

We start with a set of $`N_{\mathrm{obs}}`$ observed galaxies of absolute magnitudes $`M_i`$ and redshifts $`z_i`$. Assuming a luminosity function, $`\mathrm{\Phi }_{\mu \nu }(M,z)`$, and parameterizing its redshift evolution as described before (Eqs. 5 and 6), one can write the probability of observing galaxy $`i`$ having the properties $`M_i`$ and $`z_i`$ as

$$P(M_i,z_i|\mathrm{\Phi }_{\mu \nu })dMdz=\frac{\mathrm{\Phi }_{\mu \nu }(M_i,z_i)(dV(z_i)/dz)}{N_{\mathrm{model}}}dMdz.$$
(9)
The normalization, $`N_{\mathrm{model}}`$, is the double integral over magnitude and redshift of the LF, and hence equal the total number of galaxies expected in the sample if the LF evolves according to our parametrization:

$$N_{\mathrm{model}}=_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}_{M_{\mathrm{min}}(z)}^{M_{\mathrm{max}}}\mathrm{\Phi }_{\mu \nu }(M,z)\frac{dV}{dz}𝑑M𝑑z$$
(10)
Here we need to take the redshift and magnitude limits of the survey into account. The galaxies are taken to be in the redshift range $`z_{\mathrm{min}}<z<z_{\mathrm{max}}`$, and to have absolute magnitudes $`M_{\mathrm{max}}<M<M_{\mathrm{min}}(z)`$.

The lower magnitude limit is clearly a function of $`z`$, and can be written as

$$M_{\mathrm{min}}(z)=m_{\mathrm{min}}2.5\mathrm{log}(1+z)k_c(z)D_m(z),$$
(11)
where $`m_{\mathrm{min}}`$ is the completeness limit of the survey in apparent magnitude, $`k_c(z)`$ is the (type-dependent) $`k`$-correction at redshift $`z`$, and $`D_m(z)`$ is the distance modulus. The $`k`$-corrections for each object are obtained from its best-fitting SED (see Sect. 3).

Since we do not want to deal with incompleteness corrections here, we set

$`m_{\mathrm{min}}`$ $`=`$ $`18.7\mathrm{in}\mathrm{the}K\mathrm{band}`$
$`z_{\mathrm{min}}`$ $`=`$ $`0.4`$
$`z_{\mathrm{max}}`$ $`=`$ $`1.0`$
where MUNICS is to be considered complete (see MUNICS IV for the detailed completeness analysis).

Note that Eq. 9 cannot be used to determine the normalization, $`\mathrm{\Phi }^{}`$, or its evolution with redshift, $`\mu `$, since the dependence on $`\mathrm{\Phi }^{}`$ cancels out (see also the method presented in Lin et al., 1996). A handle on the normalization can be obtained by observing that the total number of objects in the sample, $`N_{\mathrm{obs}}`$, must be reproduced by any successful model, $`\mathrm{\Phi }_{\mu \nu }`$, leading to a boundary condition for Eq. 9. The total number of objects in the sample is by itself subject to uncertainty, obeying at least Poisson statistics (assuming that the total survey volume is large enough such that fluctuations due to large scale structure are averaged out). Instead of a boundary condition, we thus add a second term to Eq. 9 taking the Poisson noise in $`N_{\mathrm{obs}}`$ into account to obtain the final likelihood function

$$\mathrm{ln}=\underset{i=1}{\overset{N_{\mathrm{obs}}}{}}\mathrm{ln}P(M_i,z_i|\mathrm{\Phi }_{\mu \nu })\frac{1}{2}\left(\frac{N_{\mathrm{model}}N_{\mathrm{obs}}}{\sqrt{N_{\mathrm{obs}}}}\right)^2,$$
(12)
The results of this maximum likelihood approach are shown in the right panel of Fig. 12 and the most likely values of $`\mu `$ and $`\nu `$ are listed in Table 2. The dashed line in Fig. 10 shows a Schechter function evolved with $`z`$ using these numbers.

## 7 Discussion

The $`\chi ^2`$ estimates for $`\mu `$ and $`\nu `$ yield a mild decrease in number density of 30% to $`z=1`$ ($`\mu =0.30\pm 0.21`$) accompanied by a brightening of the galaxy population by $`0.64\pm 0.21`$ mag. The uncertainties of the parameters are large, though, because of the relatively large uncertainties of the values of the LF in the individual bins, and because the Schechter LF parameters are almost degenerate if the faint end below $`L^{}`$ cannot be measured (at the steep bright end, negative density evolution can be counterbalanced by brightening).

The results from the likelihood estimator are fully consistent with the $`\chi ^2`$ estimates, with their uncertainties being much smaller, though (see comparison in Table 2). The likelihood estimator gives 25% number density evolution ($`\mu =0.25\pm 0.05`$) and a brightening of $`0.53\pm 0.07`$ mag. $`\mu =0`$ is excluded at high confidence ($`4\sigma `$). The reason for the much steeper contours of the likelihood surface compared to the $`\chi ^2`$ estimator come from the tightness of the constraint in total number of galaxies (Eq. 10). This number is very sensitive to small changes in the Schechter parameters, and therefore provides a very strong constraint on the allowed region in parameter space.

It well worth noting that the parameter estimates obtained by both methods presented here agree with the best estimate obtained in MUNICS V from the spectroscopic sample alone. This makes us confident that the results presented here based on photometric redshifts are, indeed, robust. The formal uncertainties we derive in this work are, however, much smaller than those of the analysis of the purely spectroscopic sample. This is entirely due to the sample size which is a factor of 10 larger than the spectroscopic one in MUNICS V. The results presented here are also consistent with the recent estimate based on a spectroscopic sample by Pozzetti et al. (2003), finding brightening by $`0.54\pm 0.12`$ magnitudes to $`z=1`$ and at most a mild decrease ($`<30\%`$) in number density of red luminous objects.

Obviously, the results of both estimators crucially depend on the parameters assumed for the local ($`z=0`$) values of the luminosity function. We estimate the impact on the parameters due to this by using a direct fit of a Schechter function to the lowest redshift bin ($`0.4<z<0.6`$) of the MUNICS sample itself (upper left panel in Fig. 10) as a reference instead of the $`z=0`$ measurement. The results are effectively unchanged, being consistent within the $`1\sigma `$ uncertainties in both cases, in fact being almost identical for the $`\chi ^2`$ estimator.

Another possible source of systematic error is possible incompleteness in the MUNICS photometric catalogs at high redshift. This would effect the validity of the spectroscopic results in MUNICS V alike since the objects for spectroscopy are selected from the same photometric catalog. The completeness analysis presented in MUNICS IV concluded that detection probabilities are high independent of the sources’ radial profile out to $`z1`$, but that photometry suffers from severe biases for bulge-dominated systems where the lost light fraction can be high for the highest luminosity objects. For lower luminosity objects around $`L^{}`$ the photometry recovers the total magnitudes well. This is a consequence of the relation between surface brightness and total luminosity for elliptical galaxies. This effect scatters very bright objects towards lower luminosities. Since bright objects are very rare, the effect on the shape of the LF around $`L^{}`$ is small, but can be severe at $`LL^{}`$. The total magnitude of disk dominated systems are well recovered regardless of total luminosity. Unfortunately, one cannot correct for this effect without assuming a distribution of morphological types or suitable high resolution imaging which is still unavailable for a sample of this size.

To investigate whether our results are sensitive to these effects, we repeat the analysis described above using a number of restricted redshift ranges. Firstly, lowering the upper limit in redshift down to $`z=0.8`$, we find that the estimate for the density evolution changes by not more than $`\mathrm{\Delta }\mu =0.08`$ with the evolution in luminosity being practically unaffected. The error in $`\mu `$ increases by $`50\%`$, though, because of the smaller number of objects in the restricted-redshift samples. Secondly, we exclude objects at $`z<0.6`$, since the form of the LF in this bin might be altered at the bright end by the use of photometric redshifts (see Sect. 4.2 and Fig. 9). The result of our analysis did not change significantly, since we use all objects to estimate the value of the evolution parameters rather than relying on determinations in individual redshift bins.

Finally, the faint end slope, $`\alpha `$, plays a role in our analysis. The current sample does not go deep enough to constrain $`\alpha `$, which is therefore kept fixed at its local value. The $`\chi ^2`$ estimator is not very sensitive to the actual value of $`\alpha `$, since the faint part of the LF is not sampled much, and the bins are weighted by their $`S/N`$, such that the total $`\chi ^2`$ is dominated by the $`L^{}`$ region (where objects are already numerous and $`V_{\mathrm{max}}`$ corrections are not yet important).

The maximum likelihood estimator is, however, sensitive to the adapted value of $`\alpha `$. Again, this is because of the tightness of the total-number constraint and the fact that faint objects dominate in number (for all likely values of $`\alpha `$). We can estimate the size of the systematic error associated with this by using a value of $`\alpha =1`$ and repeating the analysis. The model predicts fewer galaxies in this case, and therefore the amount of density evolution required is lower, $`\mu =0.15`$. On the other hand, the LF in our lowest redshift (still $`z=0.5`$) bin is inconsistent with $`\alpha =1`$ (fixing $`\mathrm{\Phi }^{}`$ to its local value) at the $`2\sigma `$ level.

## 8 Summary

We present a measurement of the evolution of the rest-frame $`K`$-band luminosity function to $`z1`$ using a sample of $`K`$-selected galaxies drawn from the MUNICS dataset.

Distances and absolute $`K`$-band magnitudes are derived using photometric redshifts from template SED fits to BVRIJK photometry. These are calibrated using $`>500`$ spectroscopic redshifts from a spectroscopic follow-up program described in MUNICS V. A method for obtaining templated SEDs matching the properties of the sample is presented and the photometric redshift estimation is explained. We obtain redshift estimates having a rms scatter of $`0.055`$ and no mean bias.

We use Monte-Carlo simulations to investigate the influence of the errors in distance associated with photometric redshifts on our ability to reconstruct the shape of the LF (Fig. 9).

We construct the rest-frame $`K`$-band LF in four redshift bins spanning $`0.4<z<1.2`$ (Fig. 10) and use two different estimators to derive likely values for the evolution of $`\mathrm{\Phi }^{}`$ and $`M^{}`$ with redshift. While the first estimator relies on the value of the luminosity function binned in magnitude and redshift, the second estimator uses the individually measured $`\{M,z\}`$ pairs alone. Our results are shown in Fig. 12 and are summarized in Table 2. In both cases we obtain a mild decrease in number density by $`25\%`$ to $`z=1`$ accompanied by brightening of the galaxy population by $`0.5`$ to $`0.7`$ mag. These results are fully consistent with an analogous analysis using only the spectroscopic MUNICS sample (MUNICS V).

Among the sources of systematic errors are the local values of the parameters of the LF, especially the faint end slope, and possible incompleteness and biases in the photometry (discussed in detail in MUNICS IV) present in the sample. We conclude that these have some influence on our analysis, sometimes affecting the two estimators in different ways, but that these effects do not influence the general conclusion of this work, namely that there is a mild decrease in total number density of $`20\%`$ to $`30\%`$ to $`z1`$ and brightening by $`0.5`$ to $`0.7`$ mag, taking the systematic errors into account.

To interpret these results in terms of a picture of galaxy evolution is not straight-forward. While the change in characteristic luminosity is consistent with models of passive evolution of the stellar populations, we also have to account for the evolution in number density. Since the $`K`$-band light traces stellar mass, we may attribute a change in $`K`$ band luminosity function either to a change in the mass-to-light ratio alone (passive evolution) or to a combination of a change in $`M/L`$ and in stellar mass. A change in stellar mass can be due to star formation and/or due to merging and accretion. These models cannot be discriminated on basis of the luminosity function alone. Note, however, that in MUNICS III we have assumed a model of passive evolution of the mass-to-light ratio to assign stellar masses to the galaxies in the sample. This model assumed a maximum mass at the given $`K`$-band light by having all stars form at $`z=\mathrm{}`$. Analysis of the number density of objects as a function of stellar mass showed that the total number of galaxies changes only very little (like in this study), but that the number density of systems more massive than $`5\times 10^{10}h^2_{}`$ does decline, and the number density of systems more massive than $`1\times 10^{11}h^2_{}`$ declines even faster, a pattern which is qualitatively consistent with hierarchical clustering, although these models still tend to predict a more rapid decline in number density of massive systems with redshift.

A PLE model for the evolution of $`M/L`$ for the whole galaxy population is clearly too simplistic and provides only an upper envelope to the stellar mass as a function of redshift. A more proper analysis has to model the mass-to-light ratio of each object individually and take star formation and dust extinction into account. We therefore postpone further interpretation to a later paper (Drory et al. 2003, in preparation; MUNICS VI) where we use the photometry and spectroscopy combined with stellar population models to investigate the evolution of the stellar masses as a function of cosmic epoch.

This work was partly supported by the Deutsche Forschungsgemeinschaft, grant SFB 375 “Astroteilchenphysik”. ND acknowledges support by the Alexander von Humboldt Foundation. We would also like to thank the Calar Alto staff for their long-standing support during many observing runs over the last six years. This work is based on observations obtained at (1) the German-Spanish Astronomical Center, Calar Alto, operated by the Max-Planck-Institut für Astronomie, Heidelberg, jointly with the Spanish National Commission for Astronomy, (2) McDonald Observatory, operated by the University of Texas at Austin, and the Hobby - Eberly Telescope, operated by McDonald Observatory on behalf of The University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen, and (3) the European Southern Observatory, Chile, proposal number 66.A-0129 and 66.A-0123.