# An Elo-type rating model for players and teams of variable strength

Düring Bertram(1) University of Warwick, Mathematics Institute, Gibbet Hill
Road, CV47AL Coventry, UK  ,  Fischer Michael(2) University of Vienna,
Fakultät für Mathematik, Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria  and
Wolfram Marie-Therese(1)

###### Abstract.

The Elo rating system, which was originally proposed by Arpad Elo for chess,
has become one of the most important rating systems in sports, economics and
gaming nowadays. Its original formulation is based on two-player zero-sum
games, but it has been adapted for team sports and other settings.
In 2015, Junca and Jabin proposed a kinetic version of the Elo model, and
showed that under certain assumptions the ratings do converge towards the
players’ strength. In this paper we generalise their model to account for
variable performance of individual players or teams. We discuss the underlying
modelling assumptions, derive the respective formal mean-field model and
illustrate the dynamics with computational results.

##  1\. Introduction

Rating systems have become an indispensable tool to rank unobservable
quantities, such as a players’ strength based on observations, for example
outcomes of games. Rating models were originally developed for sports; but are
nowadays also used in gaming and financial markets. The Elo-rating system [1]
is one of the most prominent rating systems – it is used in chess and other
two-player zero sum games. Versions of the Elo-rating have been adopted for
many other sports, for example basketball and football, see [2, 3]. Other
prominent rating systems include the Glicko rating system or Trueskill, see
[4, 5]. Elo and Glicko are based on two-player zero sum games (here a player
can be a single individual or an entire team), while Trueskill is used in
multi-player situations, as for example in online gaming, see [6, 5].

Elo himself tried to confirm the validity of the proposed rating system using
statistical experiments [1]. It was not until 2014 that Jabin and and Junca
[7] showed the convergence of ratings towards the players’ strength for a
continuous kinetic version of the model. Junca [8] later analysed the
convergence of discrete ratings in Robin-round tournament, in which players
compete against all others in a round and discrete ratings are updated after
each such round. However, in this model the players’ strength did not change
in time. Düring et al. proposed a generalisation in [9], in which players
improve and loose skills based on the outcome of games as well as daily
performance fluctuations. A simpler but related learning mechanism was
proposed by Krupp in [10].

Kinetic models have been used very successfully to describe the behaviour of
large interacting agent systems in economics and social sciences. In all these
applications interactions between agents – such as encounters in games, the
trading of goods or the exchange of opinions – are modelled via binary
‘collisions’. Toscani [11] was the first to introduce kinetic models in the
context of opinion formation. His ideas were later generalised for more
complex opinion dynamics [12, 13, 14, 15, 16, 17, 18, 19], or in the context
of wealth distribution [20, 21] or knowledge growth in societies [22, 23]. For
a general overview on interacting multi-agent systems and kinetic equations we
refer to the book of Pareschi and Toscani [24].

The kinetic formulation of the Elo-model by Jabin and Junca assumes that each
player is characterised by a constant strength $\rho$ (being an unobservable
quantity) and a rating $R$, which changes based on the outcome of games. After
each match between player $i$ and $j$ their respective ratings $R_{i}$ and
$R_{j}$ are updated as follows

(1) $\displaystyle\begin{aligned}
{R_{i}^{*}}&={R_{i}}+\gamma(S_{ij}-b(R_{i}-R_{j})),\\\
{R_{j}^{*}}&={R_{j}}+\gamma(-S_{ij}-b(R_{j}-R_{i})).\end{aligned}$

Here $b$ is an odd, monotone, increasing function, usually chosen as
$b(z)=\tanh(\nu z)$ with a scaling constant $\nu\in\mathbb{R}^{+}$. The
parameter $\gamma$ controls the speed of adjustment. The outcome of the game
is given by the random variable ${S_{ij}}$, which takes the values
$\\{-1,1\\}$, corresponding to a win or loss (other, more fine grained
outcomes like a tie can be added in a natural way). It is assumed to equal the
expectation of $b(\rho_{i}-\rho_{j})$, that is

$\displaystyle\langle S_{ij}\rangle=\langle b({\rho_{i}}-{\rho_{j}})\rangle,$

where ${\rho_{i}},{\rho_{j}}$ are the underlying unobservable players’
strength. Note that the interactions (1) are invariant with respect to
translations and both the rating update (1) and the expected game outcome
${S_{ij}}$ depend only on the difference in $\rho$ and $R$, respectively, so
these variables are defined on $\mathbb{R}.$

Jabin and Junca then derived the corresponding macroscopic model for the
distribution of players $f(t,r,\rho)$ with respect to their rating $r$ and
their strength $\rho$:

(2) $\displaystyle\frac{\partial}{\partial
t}f(t,\rho,r)+\frac{\partial}{\partial r}(a[f]f(t,\rho,r))$ $\displaystyle=0$

with

$a[f]=\int\limits_{\mathbb{R}^{3}}w(r-r^{\prime})\big{(}b(\rho-\rho^{\prime})-b(r-r^{\prime})\big{)}f(t,\rho^{\prime},r^{\prime})d\rho^{\prime}dr^{\prime}$

and initial condition $f(0,\rho,r)=f^{0}(\rho,r)$. Here, the even probability
distribution $w$ was introduced, to account for ranking dependent pairings in
tournaments. If $w\equiv 1$ we consider a so-called all-play-all game. If $w$
has compact support only teams with close ratings compete. Possible choices
for $w$ are

(3) $\displaystyle w(r-r^{\prime})$ $\displaystyle=e^{\frac{\log
2}{1+(r-r^{\prime})^{2}}}-1\;\text{ or }\;w(r-r^{\prime})=\chi_{\\{\lvert
r-r^{\prime}\rvert\leq c\\}}.$

where $\chi$ denotes the indicator function (or smoothed variants thereof) and
$c>0$ is the maximal rating difference between paired competitors. If $w>0$
Jabin and Junca [7] showed that solutions to (2) concentrate on the diagonal,
providing the proof that the ratings indeed converge to the underlying
strength.

In this work we propose a generalisation of the Elo-model for teams of players
with fluctuating strengths. Our main contributions are the following

  * •

We propose and analyse an Elo-rating for teams, which includes stochastic
variations in the team strength due to changes in the player setup.

  * •

We formally derive the respective Fokker-Planck equations and analyse their
behaviour for long times.

  * •

We investigate the behaviour of solutions in the special case of competing
teams whose players’ strengths are distributed with a similar variance.

  * •

We illustrate the behaviour of the micro- and macroscopic models with
computational experiments, consolidating and extending the analytical results.

This work is organised as follows: we propose a microscopic generalisation of
the well-known Elo-rating to teams of players and illustrate the behaviour
with microscopic simulations in Section 2. Section 3 focuses on the
corresponding formally derived macroscopic model and its analysis. Next we
investigate the model in the case of homogeneous teams in Section 4 and report
results of computational experiments. Section 5 concludes.

##  2\. A microscopic Elo-type rating for teams

We start by proposing a microscopic version of the Elo-rating for players with
variable strength, which can also be used in the context of teams.

###  2.1. Performance variations in teams and individual players

In the following we consider a microscopic model, which accounts for
performance fluctuations in teams as well as individuals. These fluctuations
may be caused by varying individual or team performance (due to different
line-ups) or for example by card luck. We recall that small performance
fluctuations in the individual strength $\rho$ were modelled in [9] by
stochastic fluctuations in the strength. We will follow a different approach
and replace the constant strength $\rho$ by a random variable $\lambda_{\rho}$
defined on a set of possible outcomes $\Omega_{\rho}$, which can be a finite
set as well as an interval. This then allows us to define the stochastic
process $\\{\lambda_{\rho}(t)\\}_{t\in\mathbb{R}^{+}}$, whose expected value
and variance will be denoted by $\theta:=\langle\lambda_{\rho}\rangle,$ and
$\sigma^{2}:={\operatorname{Var}}[\lambda_{\rho}],$ respectively.

We consider competing teams $T_{i}$, $T_{j}$ instead of individual players.
The corresponding expected value $\theta$ can be interpreted as the mean
strength of the team with a chosen line-up, i.e. a subset of the team’s
players who will be playing in a particular game. We assume similar as in [7],
that the expected outcome of the game between $T_{i},T_{j}$ depends on the
difference of teams’ strengths through $b$:

(4) $\displaystyle\langle S_{ij}\rangle=\langle
b(\lambda_{\rho_{i}}-\lambda_{\rho_{j}})\rangle.$

If two teams $T_{i}$ and $T_{j}$ with ratings ${R_{i}}$ and ${R_{j}}$ meet,
their ratings and strength after the game can be updated using again (1) where
$\gamma$ is a scaling constant controlling the speed of adjustment. It is
usually chosen much smaller than the rating scores, in the hope that a
player’s rating slowly converges to its underlying strength. As discussed in
the introduction we make the following assumption on $b$:

  1. ($B$)

The function $b$ is $C^{3}(\mathbb{R})$, monotonically increasing, bounded,
odd and Lipschitz.

Since $b$ is non-linear the expected value and $b$ cannot be interchanged. To
calculate the expected value of $S_{ij}$ in (4) we use the Law of the
unconscious statistician [25], e.g. in the discrete case we obtain

(5) $\displaystyle\langle
S_{ij}\rangle=\sum_{x_{i}\in\Omega_{\rho_{i}}}\sum_{x_{j}\in\Omega_{\rho_{j}}}b(x_{i}-x_{j})\tilde{p}(x_{i},x_{j})=\sum_{x_{i}\in\Omega_{\rho_{i}}}\sum_{x_{j}\in\Omega_{\rho_{j}}}b(x_{i}-x_{j})p(x_{i})p(x_{j}),$

where $\tilde{p}$ denotes the probability of a possible line up $x_{i}$
playing against a line up $x_{j}$ and the second equality holds if this
happens independently of each other. In the following we always assume this
independence of the stochastic processes for team $T_{i}$ and $T_{j}$.

We can Taylor-expand $\langle S_{ij}\rangle$ in (5) as described, for example,
in [26]. Since we have
$\langle\lambda_{\rho_{i}}-\lambda_{\rho_{j}}\rangle={\theta_{i}}-{\theta_{j}}$
and
${\operatorname{Var}}[\lambda_{\rho_{i}}-\lambda_{\rho_{j}}]={\sigma_{i}}^{2}+{\sigma_{j}}^{2}$,
it follows:

(6) $\displaystyle\begin{aligned} \langle S_{ij}\rangle=\langle
b\big{(}\lambda_{\rho_{i}}-\lambda_{\rho_{j}}\big{)}\rangle&\approx
b({\theta_{i}}-{\theta_{j}})+\frac{1}{2}b^{\prime\prime}({\theta_{i}}-{\theta_{j}})({\sigma_{i}}^{2}+{\sigma_{j}}^{2})\\\
&=:b({\theta_{i}}-{\theta_{j}})+K({\theta_{i}}-{\theta_{j}},{\sigma_{i}},{\sigma_{j}}).\end{aligned}$

Note that the function $K$ is odd in the first argument and even in the other
two. Similarly, the following holds for the variance

(7)
$\displaystyle{\operatorname{Var}}[S_{ij}]\approx\big{(}b^{\prime}({\theta_{i}}-{\theta_{j}})\big{)}^{2}({\sigma_{i}}^{2}+{\sigma_{j}}^{2}).$

###  2.2. Microscopic simulations

In the following we will illustrate the behaviour of the microscopic model
with various simulations. We consider $N$ teams $T_{1}\ldots T_{N}$; each team
has $M$ players with strengths ${\rho_{i_{k}}}$ from which $m<M$ distinct
players are selected as line-up for each match. Let
$\vec{\rho}_{i}=(\rho_{i_{1}},\dots,\rho_{i_{M}})$ denote the vector of all
players in a team $T_{i}$. We assume without loss of generality that the
vector $\vec{\rho}_{i}$ is ordered.

Let us consider first the case that the $m$ players for the line-up are chosen
from the set of $M$ players uniformly. This can be done by generating
$\binom{M}{m}$ normalised vectors $\vec{\lambda}_{i_{k}}\in\\{0,1/m\\}^{M}$
with $|\vec{\lambda}_{i_{k}}|=1$. Then the stochastic process $\lambda^{t}$
selects each vector $\vec{\lambda}_{i_{k}}$ with equal probability and we have

$\displaystyle
P({S_{ij}}=1)=b(\vec{\lambda}_{i_{k}}\cdot\vec{{\rho_{i}}}-\vec{\lambda}_{j_{k}}\cdot\vec{{\rho_{j}}}),$

for a match between $T_{i}$ and $T_{j}$. More realistic line-up selection
would choose players directly proportional to their strength. We recall that
the Elo-rating is translation invariant, hence we shift the expected values
$\theta$ to the interval $[0,10]$ in the following.

As a first example consider the following football-inspired situation of
$N=200$ teams with $M=23$ players each from which $m=11$ players are selected
per match. For any team $T_{i}$ we then have
$\theta_{i}=\frac{11}{23}\sum\rho_{i_{k}}$. We investigate two different
initial setups for the teams:

  1. ($R$1)

For every team $T_{i}$, $i=1,\ldots N$, the players’ strengths $\rho_{i_{k}}$,
$k=1,\dots 23$, is chosen randomly from the interval
$\frac{1}{11}[5-\frac{5}{200}(i-1),5+\frac{5}{200}(i-1)]$. That is
$\rho_{1_{1}}=\dots\rho_{1_{23}}=\frac{1}{11}5$ and
$\rho_{200_{k}}\in\frac{1}{11}[0,10]$ for every $k$. In other words all teams
have an approximate team strength of $\theta_{i}\approx 5$ with increasing
variance $\sigma_{i}^{2}$ as can be seen in Figure 1(a).

  2. ($R$2)

For every team $T_{i}$, $i=1,\ldots N,$ the players’ strengths $\rho_{i_{k}}$
are given by

$\displaystyle\frac{1}{11}(4+\frac{6(i-1)}{197}+\eta_{i_{k}}),\quad
k\in\\{1,\dots,23\\},\quad\eta_{i_{k}}\in\mathcal{N}(0,1).$

The mean team strength of the first 198 teams is increasing from values around
4 to values around 10 and the variances $\sigma_{i}^{2}$ are of the same
order. In addition we consider two teams, Germany ($i=199$) and Brazil
($i=200$), whose mean value and variance are motivated by the 2014 FIFA World
Cup results, see [27]. We scale those values to $\theta_{\sf Ger}=10$ as well
as $\theta_{\sf Bra}=9$. However, we do not scale the variance - it is
signficantly higher in the generated data set, as can be seen in Figure 1(b).

(a) The standard deviation $\sigma$ for every team $T_{i}$ based on rule
($R$1).

(b) The standard deviation $\sigma$ regarding $\theta$ based on rule ($R$2).

Figure 1. The standard deviation of the two setups visualised, both calculated
over $10^{4}$ Monte Carlo experiments per team.

We carry out direct Monte Carlo simulation using Bird’s scheme, see [24], for
these two initial setups. We choose time steps of $\Delta t=0.1$ and perform
$25$ matches per time step. The results were then averaged over 50
realisations after $2\cdot 10^{6}$ time-steps. Figure 2 shows the final
distribution of teams for the two different initial setups. In Figure 2(a) we
see clustering around the point $(\theta,R)=(5,5)$ as expected for setup
($R$1). However, an interesting phenomenon is that teams with $\theta<5$
consistently under-perform and conversely $\theta>5$ over-perform, as they lie
above and below the line $\theta=R$, respectively. In Figure 2(b), we see
convergence towards a steady state for the 198 teams created using the rule
($R$2). However, this straight line has a steeper slope than $\theta=R$.
Furthermore, the German and Brazilian team are clear outliers, both are under-
performing relative to their strengths.

(a) Setup ($R$1) with increasing variance in the team performance.

(b) Setup ($R$2) with constant variance $\sigma=1$ in the team performance.

Figure 2. Stationary team distribution for setup ($R$1) and ($R$2).

##  3\. A macroscopic Elo-model for teams

In general, the expected value $\theta$ and variance $\sigma^{2}$ of the
microscopic model (1) are finite since they result from discrete, finite
random processes. Compared to the formal derivation of the macroscopic model
in previous works [7, 10, 1], we need additional assumptions on the moments of
$\sigma$ because of the unboundedness of $K$ in (6) in the 2${}^{\text{nd}}$
and $3^{\text{rd}}$ argument when passing from micro to macro.

Let $f(t,\theta,\sigma,r)$ be the distribution of teams at time $t$ with
expected team performance $\theta$, variance $\sigma^{2}$ and rating $r$. The
derivation of the macroscopic model (9) below is based on the following
assumptions:

  1. ($A$1)

Let $f^{0}\in H^{1}(\mathbb{R}^{3})$ with $f^{0}\geq 0$ and having compact
support. Furthermore, we assume:

$\displaystyle\int_{\mathbb{R}^{3}}f^{0}(\theta,\sigma,r)d\theta d\sigma
dr=1,\;\;\int_{\mathbb{R}^{3}}Rf^{0}(\theta,\sigma,r)d\theta d\sigma
dr=0,\;\;\int_{\mathbb{R}^{3}}\theta f^{0}(\theta,\sigma,r)d\theta d\sigma
dr=0,$ (8) $\displaystyle\int_{\mathbb{R}^{3}}\sigma
f^{0}(\theta,\sigma,r)d\theta d\sigma
dr=1,\;\;\int_{\mathbb{R}^{3}}\sigma^{2}f^{0}(\theta,\sigma,r)d\theta d\sigma
dr=C_{\sigma^{2}}.$

  2. ($A$2)

Let the interaction rate function $w\geq 0$ be an even function with $w\in
C^{2}(\mathbb{R}^{3})\cap L^{\infty}(\mathbb{R}^{3})$.

In Appendix A we derive the following macroscopic Fokker-Planck equation for
the distribution of teams $f=f(t,r,\theta,\sigma)$:

(9) $\displaystyle\begin{split}\frac{\partial}{\partial
t}f(t,\theta,\sigma,r)+\frac{\partial}{\partial
r}(a[f]f(t,\theta,\sigma,r))&=0,\qquad\text{ in }[0,T)\times\mathbb{R}^{3},\\\
f(t=0,\theta,\sigma,r)&=f^{0}(\theta,\sigma,r),\qquad\text{ in
}\mathbb{R}^{3},\end{split}$

with

$\displaystyle a[f]$
$\displaystyle=\int\limits_{\mathbb{R}^{3}}w(r-r^{\prime})\big{(}b(\theta-\theta^{\prime})+\frac{1}{2}b^{\prime\prime}(\theta-\theta^{\prime})(\sigma^{2}+\sigma^{\prime
2})-b(r-r^{\prime})\big{)}f(t,\theta^{\prime},\sigma^{\prime},r^{\prime})\,d\theta^{\prime}d\sigma^{\prime}dr^{\prime}.$

Here $b^{\prime\prime}$ is the second derivative introduced in (6). We note
that the operator $a[\cdot]$ includes an additional correction term resulting
from the variance of the distribution of strengths. This term’s sign depends
on the sign of $\theta-\theta^{\prime}$, either decreasing or increasing the
adjustment of ratings. This is consistent with the under- and over-performance
of teams observed in the microscopic simulations in Section 2.2.2.

###  3.1. Analysis of the Fokker-Planck equation

Since $b$ and $b^{\prime\prime}$ are odd, the total mass is conserved as

$\displaystyle\frac{\partial}{\partial
t}\int_{\mathbb{R}^{3}}f(t,\theta,\sigma,r)\,d\theta d\sigma dr=0.$

and therefore $\int_{\mathbb{R}^{3}}f(t,\theta,\sigma,r)\,d\theta d\sigma
dr=\int_{\mathbb{R}^{3}}f^{0}(\theta,\sigma,r)\,d\theta d\sigma dr=1$ for all
times $t>0$.

Next, we show the existence of a classical solution to (9) following arguments
from [28, 29].

######  Theorem 1.

Assume that the initial datum $f^{0}$

  1. (1)

is compactly supported in the phase space, i.e.
$\operatorname{supp}_{(r,\theta,\sigma)}f^{0}$ is bounded,

  2. (2)

is $C^{1}$-regular and bounded:

$\displaystyle\sum_{0\leq|\alpha|\leq
1}\left\|\nabla_{r}^{\alpha}f^{0}\right\|_{L^{\infty}}<\infty.$

Then, for any $t\in(0,\infty)$, there exists a unique classical solution $f\in
C^{1}\big{(}[0,t)\times\mathbb{R}^{3}\big{)}$ to (9).

###### Proof.

In the following we consider the ’all play all’ setting, that is $w\equiv 1$;
our arguments can, however, be generalised for interaction functions $w$
satisfying (A2). We start by showing that the solution cannot blow up in
finite time. Next we prove local in time existence based on a priori estimates
and a fixed point argument. Global existence follows from a continuation
argument using energy estimates.

First we show that the local solution $f$ remains uniformly bounded. We
rewrite (9) in a non-conservative form,

(10) $\displaystyle\begin{split}\frac{\partial}{\partial
t}f(t,\theta,\sigma,r)+a[f]\frac{\partial}{\partial
r}f(t,\theta,\sigma,r)=-f(t,\theta,\sigma,r)\frac{\partial}{\partial
r}a[f],\end{split}$

where we have

$\displaystyle-f(t,\theta,\sigma,r)\frac{\partial}{\partial
r}a[f]=f(t,\theta,\sigma,r)\int\limits_{\mathbb{R}^{3}}b^{\prime}(r-r^{\prime})f(t,\theta^{\prime},\sigma^{\prime},r^{\prime})\,d\theta^{\prime}d\sigma^{\prime}dr^{\prime},$

which yields

$\displaystyle\left\|\frac{\partial}{\partial r}a[f]\right\|_{L^{\infty}}\leq
L,$

where $L$ is the Lipschitz constant of $b$ and we used that the total mass
equals one. Next we consider a trajectory starting at time
$\tau_{0}\in\mathbb{R}^{+}$ in $(r_{0},\theta_{0},\sigma_{0})$, then the
characteristics are given by

(11) $\displaystyle\frac{\partial}{\partial
t}r(t)=a[f],\qquad\frac{\partial}{\partial
t}\theta(t)=\frac{\partial}{\partial t}\sigma(t)=0.$

Therefore,

$\displaystyle\frac{\partial}{\partial t}f(t,\theta,\sigma,r)\leq
Lf(t,\theta,\sigma,r),$

and Gronwall’s lemma gives

(12) $\displaystyle\left\|f(t)\right\|_{L^{\infty}}\leq
e^{Lt}\left\|f^{0}\right\|_{L^{\infty}}.$

Hence, the solution cannot blow up in finite time.

We continue with the existence of a local solution following Theorem 3.1 in
[29]. To this end we investigate the non-linear transport operator
$\mathcal{H}$:

(13)
$\displaystyle\mathcal{H}(f)=-f(t,\theta,\sigma,r)\frac{\partial}{\partial
r}a[f],\qquad\text{with }\mathcal{H}:=\frac{\partial}{\partial
t}+a[f]\frac{\partial}{\partial r},$

in the following. There exist positive constants $C_{1},C_{2}$ such that

(14) $\displaystyle\begin{aligned}
&|{\mathcal{H}}(f)|\leq\left\|b^{\prime}\right\|_{L^{\infty}}|f|=C_{1}|f|,\\\
&|{\mathcal{H}}(\frac{\partial}{\partial
r}f)|\leq|f|\left\|b^{\prime\prime}\right\|_{L^{\infty}}+\left|\frac{\partial}{\partial
r}f\right|\left\|b^{\prime}\right\|_{L^{\infty}}=C_{2}\big{(}|f|+\left|\frac{\partial}{\partial
r}f\right|\big{)},\end{aligned}$

because of $b$ being Lipschitz and therefore we have $L_{\infty}$ bounds for
its derivatives. Moreover, the map $H$ defined by

$\displaystyle H:C^{1}\big{(}[0,t)\times\mathbb{R}^{3}\big{)}\to
C^{1}\big{(}[0,t)\times\mathbb{R}^{3}\big{)},\quad
f\mapsto(\mathcal{H})^{-1}(-f\frac{\partial}{\partial r}a[f])$

is bounded since we can use the estimates (14) together with the bounded
inverse theorem. Next we consider the solution along the trajectories
$(\theta(t),\sigma(t),r(t))$

$\displaystyle
f(t,\theta(t),\sigma(t),r(t))=f(0,\theta_{0},\sigma_{0},r_{0})+\int_{0}^{t}\frac{\partial}{\partial
r}\bigg{(}a[f(t,\theta(t),\sigma(t),r(t))]f(t,\theta(t),\sigma(t),r(t))\bigg{)}dt.$

We can then use the previous estimates to choose a $t>0$ such that $H$ is a
contraction. Using Banach’s fixed point theorem we obtain a unique local
solution $f\in C^{1}\big{(}[0,t]\times\mathbb{R}^{3}\big{)}$. Let
$\mathcal{F}$ be the $W^{1,\infty}$-norm of $f(t)$,

$\displaystyle\mathcal{F}(t)=\sum_{0\leq|\alpha|\leq
1}\left\|\nabla_{r}^{\alpha}f(t)\right\|_{L^{\infty}}.$

Using (14) and again Gronwall’s lemma we get

$\displaystyle\frac{\partial}{\partial t}\mathcal{F}(t)\leq
C_{3}\mathcal{F}(t)$

and therefore we have the upper bound

$\displaystyle\mathcal{F}(t)\leq\mathcal{F}(0)e^{C_{3}t},\qquad\forall
t\in[0,T).$

The energy bound in $W^{1,\infty}$ allows us to use the standard continuation
principle, giving the global extension of the local solution. ∎

We continue by analysing the behaviour of the moments of $f$. We define the
$s$-th moment for $s\in\mathbb{N}$ with respect to $r$ (and similar the
moments with respect to $\theta,\sigma$),

(15) $\displaystyle m_{s,r}(t)=\int_{\mathbb{R}^{3}}r^{s}f\,d\theta
drd\sigma.$

The evolution with respect to $\sigma$ and $\theta$ is trivial, as the
function does not change with respect to these variables. The evolution of the
second moment w.r.t. to $r$ satisfies:

$\displaystyle\frac{d}{dt}\int\limits_{\mathbb{R}^{3}}r^{2}f\,d\theta
drd\sigma=-\int\limits_{\mathbb{R}^{3}}r^{2}\frac{d}{dr}(a[f]f)\,d\theta
drd\sigma$
$\displaystyle=2\int\limits_{\mathbb{R}^{6}}rw(r-r^{\prime})\Big{(}b(\theta-\theta^{\prime})+\frac{\sigma^{2}+\sigma^{\prime
2}}{2}b^{\prime\prime}(\theta-\theta^{\prime})-b(r-r^{\prime})\Big{)}ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\sigma^{\prime}d\theta
drd\sigma$
$\displaystyle=-\int\limits_{\mathbb{R}^{6}}rb(r-r^{\prime})w(r-r^{\prime})f^{\prime}f\,d\theta^{\prime}dr^{\prime}d\sigma^{\prime}d\theta
drd\sigma-\int\limits_{\mathbb{R}^{6}}r^{\prime}b(r^{\prime}-r)w(r^{\prime}-r)f^{\prime}f\,d\theta^{\prime}dr^{\prime}d\sigma^{\prime}d\theta
drd\sigma$
$\displaystyle\quad+2\int\limits_{\mathbb{R}^{6}}rw(r-r^{\prime})\big{(}b(\theta-\theta^{\prime})+\frac{\sigma^{2}+\sigma^{\prime
2}}{2}b^{\prime\prime}(\theta-\theta^{\prime})\big{)}ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\sigma^{\prime}d\theta
drd\sigma$
$\displaystyle=-\int\limits_{\mathbb{R}^{6}}(r-r^{\prime})b(r-r^{\prime})w(r-r^{\prime})ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\sigma^{\prime}d\theta
drd\sigma$ (16) $\displaystyle<0.$

Here, we used the short-hand notation
$f^{\prime}=f(t,\theta^{\prime},\sigma^{\prime},r^{\prime})$. Furthermore, we
used that for $r-r^{\prime}<0$ the function $b(r-r^{\prime})<0$ is negative,
since $b$ is odd and monotonically increasing. The latter does not hold in
general for $b+b^{\prime\prime}$, however, the second integral vanishes since
the integrand is still odd in $\theta$ and $r$. Therefore, the second moment
in $r$ decreases over time and we expect convergence towards a stationary
state. Our computational experiments confirm this expected convergence.
However, we are not able to compute these stationary states explicitly as it
was done in [7].

###  3.2. Numerical results for the macroscopic model

We perform several computational experiments illustrating the dynamics of (9)
using a finite difference scheme. It is based on the generalisation of a
finite difference scheme for conservation laws with discontinuous flux
presented by Towers in [30]. This generalisation is straight-forward, as (9)
has only transport in $r$ direction. Let $f^{n}_{j,l,m}$ denote the solution
at the discrete points $(j\Delta r,l\Delta\theta,m\Delta\sigma)$,
$j,l,m\in\mathbb{N},$ and time $t^{n}=n\Delta t$, $n\in\mathbb{N},$ with
discrete positive increments $\Delta r,\Delta\theta$, $\Delta\sigma$ and
$\Delta t$. Then the explicit scheme reads as follows:

$\displaystyle f^{n+1}_{j,l,m}=f^{n}_{j,l,m}-\frac{\Delta t}{\Delta
r}(a^{n}_{j+1/2,l,m}h_{j+1/2,l,m}^{n}-a^{n}_{j-1/2,l,m}h_{j-1/2,l,m}^{n}),$

with cell averages $a^{n}_{j+1/2,l,m}=\frac{1}{\Delta
r}\int_{j}^{j+1}a[f^{n}(r,l\Delta\theta,m\Delta\sigma)]\,dr$. The function $h$
is chosen depending on the sign of the averaged flux (as in the usual Godunov
scheme) that is

$\displaystyle
h^{n}_{j+1/2,l,m}=\begin{cases}f^{n}_{j,l,m},&a^{n}_{j+1/2,l,m}\geq 0,\\\
f^{n}_{j+1,l,m},&a^{n}_{j+1/2,l,m}<0.\end{cases}$

In Figure 3 we visualise the first marginals
$\int_{\mathbb{R}}f(t,r,\theta,\sigma)\,d\sigma$ and
$\int_{\mathbb{R}}f(t,r,\theta,\sigma)\,d\theta$ of the team distribution at
time $t=5$ using timesteps of size $\Delta t=10^{-5}$ and $\nu=1$. The
computational domain is $\Omega=[0,10]\times[0,10]\times[0,1]$ and the spatial
discretisation was set to $\Delta r=\Delta\theta=\Delta\sigma=5\cdot 10^{-2}$,
the initial distribution of teams uniform and normalised. The left plot in
Figure 3 shows that ratings converge towards the mean strength for
$\theta\in[6,8]$, but are blurred for smaller and larger means. We observe a
similar over- and under-performance as in the microscopic simulations in
Figure 2(b). The right plot illustrates the decrease of $m_{2,r}$ in the
direction $\theta$. The larger the uncertainty $\sigma$, the less accurate the
ratings as all teams get a similar rating (around $7$).

Figure 3. Second and third marginal of the team distribution $f$ at time $t=5$
for $f^{0}=1$.

##  4\. Special scaling limits and homogeneous player distributions

Our microscopic computational results suggest that if all teams have the same
player variance, then the ratings converge to the underlying mean team
strength. In this case, however, the integral over $\sigma$ can be seen as a
point evaluation and we can simplify (9) for constant
$\sigma\in\mathbb{R}^{+}_{0}$:

(17) $\displaystyle\begin{split}\frac{\partial}{\partial
t}f(t,\theta,r)+\frac{\partial}{\partial r}(a[f]f(t,\theta,r))&=0\\\
f(t=0,\theta,r)&=f^{0}(\theta,r),\end{split}$

with $a$ changed to

(18) $\displaystyle a[f]$
$\displaystyle=\int\limits_{\mathbb{R}^{2}}w(r-r^{\prime})\big{(}b(\theta-\theta^{\prime})+\sigma^{2}b^{\prime\prime}(\theta-\theta^{\prime})-b(r-r^{\prime})\big{)}f(t,\theta^{\prime},r^{\prime})dr^{\prime}d\theta^{\prime}.$

We discuss the existence of a unique solution and the analysis of the moments.
Furthermore we consider the relative energy to prove convergence of the team
strengths to ratings. The existence of a classical solution itself follows
from the same arguments as in Theorem 1.

######  Theorem 2.

Assume that the initial datum $f^{0}$

  1. (1)

is compactly supported in the phase space, i.e.
${\operatorname{supp}}_{(r,\theta)}f^{0}$ is bounded

  2. (2)

is $C^{1}$-regular and bounded:

$\displaystyle\displaystyle\sum_{0\leq|\alpha|\leq
1}\left\|\nabla_{r}^{\alpha}f^{0}\right\|_{L^{\infty}}<\infty.$

Then, for any $T\in(0,\infty)$, there exists a unique classical solution $f\in
C^{1}\big{(}[0,T)\times\mathbb{R}^{2})$ to (17).

The proof of Theorem 2 can be easily adapted from the proof of Theorem 1 and
is omitted here. Equation (17) is conservative, hence the total mass is
preserved, and the moments with respect to $\theta$ is zero. Again, the second
moment w.r.t. to $r$ is decreasing (using similar arguments as in (3.1)):

$\displaystyle\frac{d}{dt}\int\limits_{\mathbb{R}^{2}}r^{2}f\,d\theta dr$
$\displaystyle=-\int\limits_{\mathbb{R}^{2}}r^{2}\frac{d}{dr}(a[f]f)\,d\theta
drd\theta^{\prime}dr^{\prime}$
$\displaystyle=-\int\limits_{\mathbb{R}^{4}}(r-r^{\prime})\big{(}b(\theta-\theta^{\prime})+\sigma^{2}b^{\prime\prime}(\theta-\theta^{\prime})\big{)}w(r-r^{\prime})ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\theta
dr<0,$

using the short hand notation $f^{\prime}=f(t,\theta^{\prime},r^{\prime})$.

The ratio of $\nu$ and $\sigma$ is important in order to be able to show the
convergence of the team ratings to the average strength. This is the case
under following assumption:

  1. ($B^{\prime}$)

$b+\sigma^{2}b^{\prime\prime}$ is monotonically increasing.

Note that ($B^{\prime}$) holds for example for $b(z)=\tanh(\nu z)$ if
$1+\nu^{2}\sigma^{2}(4-6\operatorname{sech}(z\nu)^{2})>0$. Then we can use
similar arguments as Jabin and Junca [7], who considered the relative energy

$\displaystyle\mathcal{E}(t)=\int\limits_{\mathbb{R}^{2}}(r-\theta)^{2}f(t,\theta,r)\,drd\theta.$

In the following we will show that

(19) $\displaystyle\frac{d\mathcal{E}(t)}{dt}<0.$

We calculate:

$\displaystyle\frac{d}{dt}\int\limits_{\mathbb{R}^{2}}(r-\theta)^{2}f\,d\theta
dr$
$\displaystyle=-\int\limits_{\mathbb{R}^{2}}(r-\theta)^{2}\frac{d}{dr}(a[f]f)\,d\theta
drd\theta^{\prime}dr^{\prime}$
$\displaystyle=-\int\limits_{\mathbb{R}^{4}}(r-r^{\prime})b(r-r^{\prime})w(r-r^{\prime})ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\theta
dr$
$\displaystyle\quad-\int\limits_{\mathbb{R}^{4}}(\theta-\theta^{\prime})\big{(}b(\theta-\theta^{\prime})+\sigma^{2}b^{\prime\prime}(\theta-\theta^{\prime})\big{)}w(r-r^{\prime})ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\theta
dr<0.$

For $r-r^{\prime}<0$ we have $b(r-r^{\prime})<0$, while the opposite holds
true for $(r-r^{\prime})>0$. Because of ($B^{\prime}$) the second term is
positive, yielding the stated energy decay.

Assumption ($B$) together with ($B^{\prime}$) gives us bounds for
$b^{\prime\prime\prime}$, whereas we can deduce $b^{\prime\prime}$ being
Lipschitz, too, with constant $L_{2}$. Following the arguments in Jabin and
Junca, [7], we obtain:

######  Theorem 3.

Let $f^{0}$ be as in Theorem 2 and $w\geq w_{min}>0$ on
${\operatorname{supp}}f^{0}$. Then

$\displaystyle\mathcal{E}(t)\leq\mathcal{E}(0)\exp(-2w_{min}(L+\sigma^{2}L_{2})t),$

where $L,L_{2}$ depend on $b,b^{\prime\prime}$ (and therefore $\nu$) and on
${\operatorname{supp}}f^{0}.$

###### Proof.

The proof is along the lines of [7, 10], adapted for the additional term
related to $b^{\prime\prime}$. We define

$\displaystyle\mathcal{D}(f(t)):=\int\limits_{\mathbb{R}^{4}}(r-r^{\prime}-\theta+\theta^{\prime})w(r-r^{\prime})[b(\theta-\theta^{\prime})+\sigma^{2}b^{\prime\prime}(\theta-\theta^{\prime})-b(r-r^{\prime})]ff^{\prime}\,dr^{\prime}d\theta^{\prime}drd\theta,$

and by our previous calculations we have that
$-\mathcal{D}(f(t))=\frac{d\mathcal{E}(t)}{dt}\leqslant 0$. Using similar
symmetry arguments as before we deduce that

$\displaystyle-\int\limits_{\mathbb{R}^{4}}(r-r^{\prime})w_{min}\sigma^{2}b^{\prime\prime}(r-r^{\prime})ff^{\prime}\,d\theta^{\prime}dr^{\prime}d\theta
dr<0.$

We now split the integrands of $\mathcal{D}$ to obtain

$\displaystyle\mathcal{D}(f(t))\geq$
$\displaystyle\int\limits_{\mathbb{R}^{4}}(r-r^{\prime}-\theta+\theta^{\prime})w_{min}[b(\theta-\theta^{\prime})-b(r-r^{\prime})]ff^{\prime}\,dr^{\prime}d\theta^{\prime}drd\theta$
$\displaystyle+\int\limits_{\mathbb{R}^{4}}(r-r^{\prime}-\theta+\theta^{\prime})w_{min}\sigma^{2}[b^{\prime\prime}(\theta-\theta^{\prime})-b^{\prime\prime}(r-r^{\prime})]ff^{\prime}\,dr^{\prime}d\theta^{\prime}drd\theta.$

Using that $b,b^{\prime\prime}$ are Lipschitz and odd we have

$\displaystyle(r-r^{\prime}-\theta+\theta^{\prime})w_{min}[b(\theta-\theta^{\prime})-b(r-r^{\prime})]$
$\displaystyle\geq Lw_{min}|r-r^{\prime}-\theta+\theta^{\prime}|^{2}$
$\displaystyle(r-r^{\prime}-\theta+\theta^{\prime})w_{min}\sigma^{2}[b^{\prime\prime}(\theta-\theta^{\prime})-b^{\prime\prime}(r-r^{\prime})]$
$\displaystyle\geq\sigma^{2}L_{2}w_{min}|r-r^{\prime}-\theta+\theta^{\prime}|^{2}.$

Therefore, it follows that

$\displaystyle\mathcal{D}(f)\geqslant\int\limits_{\mathbb{R}^{4}}(L+\sigma^{2}L_{2})w_{min}|r-r^{\prime}-\theta+\theta^{\prime}|^{2}f^{\prime}f\,dr^{\prime}d\theta^{\prime}drd\theta.$

We assume w.l.o.g. (due to the translation invariance of the model)

(20) $\displaystyle\int\limits_{\mathbb{R}^{2}}rf(t,r,\theta)\ dr\
d\theta=\int\limits_{\mathbb{R}^{2}}\theta f(t,r,\theta)\,drd\theta=0$

which gives

$\displaystyle\int\limits_{\mathbb{R}^{4}}(r-\theta)(r^{\prime}-\theta^{\prime})\
f^{\prime}f\,dr^{\prime}drd\theta^{\prime}d\theta=0.$

Then we can deduce

$\displaystyle\mathcal{D}(f)\geqslant
2(L+\sigma^{2}L_{2})w_{min}\int\limits_{\mathbb{R}^{2}}|r-\theta|^{2}\ fdr\
d\theta$

and altogether

$\displaystyle-\mathcal{D}(f)=\frac{d\mathcal{E}(t)}{dt}\leqslant-2Lw_{min}\int\limits_{\mathbb{R}^{2}}|r-\theta|^{2}\
f\ dr\ d\theta=-2w_{min}(L+\sigma^{2}L_{2})\mathcal{E}(t).$

Using Gronwall’s lemma we conclude the proof as in [7, 10]. ∎

We conclude by underpinning our analytical results with numerical simulations.

##### Micro- and macroscopic simulations.

For the microscopic simulation we consider $N=500$ players with fixed mean
strengths $\theta_{n}$, chosen uniformly distributed in $[4,10]$. In every
time-step we then choose $\mathcal{N}(\theta_{n},\sigma)$ distributed values
for the evaluation of ${S_{ij}}$. We set $\nu=0.5$ and simulate $10^{6}$ time-
steps of $\Delta t=0.1$ and 25 collisions per time-step over 50 realisations
as in the previous simulations described in Section 22.2. On a macroscopic
level, we use the algorithm presented in Section 33.2 reduced by the dimension
in $\sigma$.

We see a great agreement between the two models in Figure 4 and 5. In
addition, we clearly see the influence of $\sigma$ on the ranking as discussed
in Figure 3. If the variance $\sigma$ is large, all teams are rated equally,
in particular the ratings converge to $7$ for all values of $\theta$. Or
expressed differently: in expectation weaker teams are over-performing and
stronger teams under-performing. We see a similar effect already in our first
microscopic simulations in Figure 2.

Moreover, the numerical simulations show that $\nu$ can be used to balance the
variance $\sigma^{2}$ and obtain the desired convergence of ratings to the
teams’ average strength as discussed in the previous subsection. If assumption
($B^{\prime}$) holds the long-term behaviour of (17) will be similar to the
original Elo model (2) in [7]. This effect can also be observed on a
microscopic level, see Figure 6, where we compare the long-term behaviour for
$\sigma=2$ and different values of $\nu\in\\{1,0.1,0.01\\}$. For $\nu=1$ all
expected values $\theta\in[4,10]$ converge to $7$ and we get an almost
horizontal line in the long-term run. If we choose $\nu=0.1$ we get the
desired diagonal $\theta=r$ as time increases, as in the case $\nu=0.01$.
However, the smaller value of $\nu$ corresponds to a slower convergence
towards the stationary state, as can be seen in Figure 6, bottom right.

Figure 4. Microscopic team distributions for $\nu=0.5$ and different values of
$\sigma$. Figure 5. Macroscopic results for $\Omega=[4,10]\times[4,10]$ for
different values of $\sigma$ and $\nu=0.5$. Figure 6. A smaller value of $\nu$
leads to slower convergence, however, we can retrieve the desired convergence
of $R$ to the expected value $\theta$ by choosing $\nu$ sufficiently small. In
all simulations we set $\sigma=2$.

##  5\. Conclusion

In this paper we proposed a generalisation of the Elo-rating model for teams
of players with varying strengths, which includes fluctuations in the
performance to account for example for variable line-ups in team sports. Based
on the microscopic interaction rules we then derived the corresponding kinetic
model, proved existence of a solution and analysed different moments of its
solution. These analytical insights indicate the formation of non-trivial
steady states – a hypothesis that is supported by our numerical results.
Furthermore, we considered the special case of similar variance $\sigma^{2}$,
which allowed to formally derive a lower dimensional equation. Under further
smallness assumptions we could then use techniques from Junca and Jabin, see
[7], to show convergence of the rating $r$ to the expected value $\theta$. The
smallness assumption relates to practically relevant parameter values. For
example, in chess the scaling parameter $\nu$ in $b$ is usually quite small,
around $\frac{1}{400}$, as reported in [1, 31, 32]. We were able to show
numerically, both at the microscopic and kinetic level, that a large $\nu$
leads to the loss of convergence of $r\to\theta$. Choosing $\nu$ according to
($B^{\prime}$), we obtain the desired convergence and were able to proof this
analytically. This effect also occurs at the microscopic level.

Nevertheless, the microscopic simulations showed that a large $\sigma$ has a
strong impact on the ratings. The question therefore remains whether
fluctuations in the underlying strength should be included in $\rho$, see [9],
or incorporated in the outcome of the game ${S_{ij}}$ (as proposed in this
paper). Following [9], performance fluctuations could also be included via an
additional random term in the microscopic interactions. This leads to a PDE
with a diffusive term which is of the following form:

$\displaystyle\frac{\partial f(r,\theta,t)}{\partial
t}=-\frac{\partial}{\partial
r}\left(a[f]f(\theta,R,t)\right)+\frac{{\sigma^{2}}}{2}d[f]\frac{\partial^{2}}{\partial\theta^{2}}f(\theta,r,t),$

with

$\displaystyle
a[f]=a[f](r,\theta,t)=\int_{\mathbb{R}^{2}}w(r-r^{\prime})(b(\theta-\theta^{\prime})-b(r-r^{\prime}))f(\theta^{\prime},r^{\prime},t)\,d\theta^{\prime}dr^{\prime},$
$\displaystyle
d[f]=d[f](r,\theta,t)=\int_{\mathbb{R}^{2}}w(r-r^{\prime})f(r^{\prime},\theta^{\prime},t)\,d\theta^{\prime}dr^{\prime}$

where the influence of diffusion is determined by the maximum variance
$\sigma^{2}$ of the team strengths.

Accounting for uncertainty in ratings through an additional functional
dependence and not via diffusion also happens in the Glicko rating [4], an
extension of the Elo rating. However, here the variable $\sigma$ is the
uncertainty of the rating. It is assumed that $\sigma$ increases if players do
not compete and decreases if they participate in tournaments. This microscopic
model could similarly be used to derive a continuous kinetic rating model.
Another interesting direction of future research is the combination of
performance fluctuations with learning effects, as considered in [9, 10].

### Contributions and Acknowledgments

The paper has been conceived and is based on work from all three authors. MF
obtained the analytical results under guidance from BD and MTW. MF carried out
the numerical experiments with support from MTW. All authors worked on and
approved the manuscript.

MF acknowledges partial support via the Austrian Science Fund (FWF) project
F65. MTW acknowledges partial support via the New Frontier’s Project NST-0001
of the Austrian Academy of Sciences ÖAW.

##  Appendix A Appendix

###  A.1. Derivation of the Boltzmann-type equation

We follow the derivation of [9] for the corresponding PDE of Fokker-Planck
type to study the dynamics of the corresponding model. We start with the
evolution equation for the distribution of teams
$f_{\gamma}=f_{\gamma}(\theta,\sigma,R,t)$ with respect to their rating $R$,
intrinsic team-strength $\rho$ and the variance ${\sigma_{i}}$. For a fixed
number of teams, $N$, the interactions (1) induce a discrete-time Markov
process with $N$-particle joint probability distribution
$P_{N}(\theta_{1},\sigma_{1},R_{1},\theta_{2},R_{2},\dots,\theta_{N},\sigma_{N},R_{N},\tau)$.
Then we can state the evolution of the first marginal from

$\displaystyle P_{1}(\theta,\sigma,R,\tau)=\int
P_{N}(\theta,\sigma,R,\theta_{2},\sigma_{2},R_{2},\dots,\theta_{N},\sigma_{N},R_{N},\tau)d\theta_{2}d\sigma_{2}dR_{2}\cdots
d\theta_{N}d\sigma_{N}dR_{N},$

where $\tau$ is the discrete time step using only the one- and two-particle
distribution functions [33, 34] in a single time step,

$\displaystyle\begin{split}&P_{1}(\theta,\sigma,R,\tau+1)-P_{1}(\theta,\sigma,R,\tau)=\\\
&\Bigg{\langle}\frac{1}{N}\Biggl{[}\int_{\mathbb{R}^{6}}P_{2}({\theta_{i}},{\sigma_{i}},R_{i},{\theta_{j}},{\sigma_{j}},R_{j},\tau)w({R_{i}}-{R_{j}})\bigl{(}\delta_{0}(\theta-{\theta_{i}}^{*},R-{R_{i}^{*}})+\delta_{0}(\rho-{\rho_{j}^{*}},R-{R_{j}^{*}})\bigr{)}\cdot\\\
&\cdot
d{\theta_{i}}d{\sigma_{i}}dR_{i}d{\theta_{j}}d{\sigma_{j}}dR_{j}-2P_{1}(\theta,\sigma,R,\tau)\Biggr{]}\Biggr{\rangle}.\end{split}$

Here, $\langle\cdot\rangle$ denotes the mean operator with respect to the
random variables ${S_{ij}}$ and the function $w(\cdot)$ corresponds to the
interaction rate function which depends on the difference of the ratings. This
yields a hierarchy of equations, the so-called BBGKY-hierarchy, see [33, 34],
describing the dynamics of the system of a large number of interacting agents.

A standard approximation is to neglect correlations and assume that

$\displaystyle
P_{2}({\theta_{i}},{\sigma_{i}},R_{i},{\theta_{j}},{\sigma_{j}},R_{j},\tau)=P_{1}({\theta_{i}},{\sigma_{i}},R_{i},\tau)P_{1}({\theta_{j}},{\sigma_{j}},R_{j},\tau),$

By scaling time as $t=2\tau/N$ and performing the thermodynamical limit
$N\to\infty$, we can use standard methods of kinetic theory [33, 34] to show
that the time-evolution of the one-agent distribution function $f_{\gamma}$
(corresponding to $P_{1}$ and $f_{\gamma}f_{\gamma}$ to $P_{2}$) is governed
by the following Boltzmann-type equation:

(21) $\displaystyle\begin{aligned}
&\frac{d}{dt}\int_{{\Omega}}\phi(\theta,\sigma,r)f_{\gamma}(\theta,\sigma,r,t)d\theta
d\sigma dr=\\\
&\frac{1}{2}\Bigg{\langle}\int_{{\Omega}}\int_{{\Omega}}\Big{(}\phi(\theta,\sigma,r^{*})+\phi(\theta^{\prime},\sigma^{\prime},r^{\prime*})-\phi(\theta,\sigma,r)-\phi(\theta^{\prime},\sigma^{\prime},r^{\prime})\Big{)}\cdot\\\
&\cdot
w(r-r^{\prime})f_{\gamma}(\theta,\sigma,r,t)f_{\gamma}(\theta^{\prime},\sigma^{\prime},r^{\prime},t)d\theta^{\prime}d\sigma^{\prime}dr^{\prime}d\theta
d\sigma dr\Bigg{\rangle},\end{aligned}$

where $\phi(\cdot)$ is a (smooth) test function, with support
$\mathrm{supp}({\phi})\subseteq{\Omega}$.

###  A.2. Analysis of the Boltzmann-type equation

#### Conservation of mass

Setting $\phi(\theta,\sigma,r)=1$ in equation (21) we have

$\displaystyle\frac{d}{dt}\int_{\mathbb{R}^{3}}f_{\gamma}(\theta,\sigma,r,t)d\theta
d\sigma dr=0.$

Therefore, the total mass is conserved, that is

$\displaystyle\int_{\mathbb{R}^{3}}f_{\gamma}(\theta,\sigma,r,t)d\theta
d\sigma dr=1,\ \text{ for all times }t\geq 0.$

#### Moments with respect to the rating

We define the $s$-th moment for $s\in\mathbb{N}$ with respect to $r$

(22) $\displaystyle
m_{s,r}(t)=\int_{\mathbb{R}^{3}}r^{s}f_{\gamma}(\theta,\sigma,r,t)d\theta
d\sigma dr.$

Now choose $\phi(\theta,\sigma,r)=r$. Due to ($B$), ($A$1) and the symmetry of
$b(\cdot)$ and $b^{\prime\prime}(\cdot)$ we obtain

(23) $\displaystyle\begin{aligned}
\frac{d}{dt}m_{1,r}(t)=&\frac{1}{2}\gamma\int_{\mathbb{R}^{6}}f_{\gamma}(\theta,\sigma,r,t)f_{\gamma}(\theta^{\prime},\sigma^{\prime},r^{\prime},t)w(r-r^{\prime})\cdot\\\
&\cdot\bigg{(}b(\theta-\theta^{\prime})-\frac{1}{2}b^{\prime\prime}(\theta-\theta^{\prime})({\sigma^{2}}+\sigma^{\prime
2})-b(r-r^{\prime})+\\\
&+b(\theta^{\prime}-\theta)+\frac{1}{2}b^{\prime\prime}(\theta^{\prime}-\theta)({\sigma^{2}}+\sigma^{\prime
2})-b(r^{\prime}-r)\bigg{)}d\theta^{\prime}d\sigma^{\prime}dr^{\prime}d\theta
d\sigma dr=0.\end{aligned}$

Hence the mean value w.r.t. the rating is preserved in time and therefore
$m_{1,r}=0$ for all times $t\geq 0$.

#### Moments with respect to the variance and expected value

We define $m_{s,\theta}$ and $m_{s,{\sigma^{2}}}$ similar to (22) and also
need the boundedness of the second moment $m_{2,{\sigma^{2}}}$.

(24) $\displaystyle\frac{d}{dt}m_{2,{\sigma^{2}}}(t)=0$

which follows directly from (21) when testing with
$\phi(\theta,\sigma,r)={\sigma^{2}}$. Analogue we get

(25) $\displaystyle\frac{d}{dt}m_{1,\theta}(t)=\frac{d}{dt}m_{2,\theta}(t)=0.$

###  A.3. Derivation of the Fokker-Planck equation

We now derive the limiting Fokker-Planck equation in the case
$\gamma\rightarrow 0$. Based on the interaction rules (1), which define the
outcome of a game, we compute the expected values of the following quantities:

$\displaystyle\langle r^{*}-r\rangle=\gamma(\langle S\rangle-b(r-r^{\prime}))$
(26)
$\displaystyle{\operatorname{Var}}[r^{*}-r]=\gamma^{2}{\operatorname{Var}}[{S_{ij}}]$
$\displaystyle\langle(r^{*}-r)^{2}\rangle=\gamma^{2}(\langle
S\rangle-b(r-r^{\prime}))^{2}+{\operatorname{Var}}[r^{*}-r]=\gamma^{2}\big{(}\langle
S\rangle-b(r-r^{\prime}))^{2}+{\operatorname{Var}}[S]\big{)}$

with $S$ analogue to (6) where we used $\langle X^{2}\rangle=\langle
X\rangle^{2}+\operatorname{Cov}[X,X]=\langle
X\rangle^{2}+{\operatorname{Var}}[X]$. Using Taylor expansion of
$\phi(\theta,\sigma,r^{*})$ up to order two around $(\theta,\sigma,r)$, we
obtain

$\displaystyle\langle\phi(\theta,\sigma,r^{*})-\phi(\theta,\sigma,r)\rangle=$
$\displaystyle\langle r^{*}-r\rangle\frac{\partial}{\partial
r}\phi(\theta,\sigma,r)+\frac{1}{2}\langle(r^{*}-r)^{2}\rangle\frac{\partial}{\partial
r^{2}}\phi(\theta,\sigma,r)+$
$\displaystyle+\langle{\mathcal{R}_{\gamma}}(\phi,\theta,\sigma,r,\tau)\rangle,$

where the remainder term ${\mathcal{R}_{\gamma}}$ is given in the Peano-
representation of Taylor’s formula via

$\displaystyle\langle{\mathcal{R}_{\gamma}}(\phi,\theta,\sigma,r,\tau)\rangle$
$\displaystyle=\frac{1}{2}\langle(r^{*}-r)^{2}\rangle\frac{\partial}{\partial
r^{2}}(\phi(\theta,\sigma,\bar{r})-\phi(\theta,\sigma,r))$
$\displaystyle=\frac{1}{2}\gamma^{2}\big{(}\langle
S\rangle-b(r-r^{\prime}))^{2}+{\operatorname{Var}}[S]\big{)}\frac{\partial}{\partial
r^{2}}(\phi(\theta,\sigma,\bar{r})-\phi(\theta,\sigma,r))$

for some $0\leq c\leq 1$ with $\bar{r}$ defined as

$\displaystyle\bar{r}=cr+(1-c)r^{*}.$

Next we rescale time as $\tau=\gamma t$ and insert the expansion in (21). This
yields

$\displaystyle\begin{split}\frac{d}{d\tau}\int_{\mathbb{R}^{3}}&\phi(\theta,\sigma,r)f_{\gamma}(\theta,\sigma,r,\tau)d\theta
d\sigma
dr=\frac{1}{2\gamma}\int_{\mathbb{R}^{3}}{\tilde{\mathcal{R}}_{\gamma}}(\phi,\theta,\sigma,r,\tau)f_{\gamma}(\theta,\sigma,r,\tau)d\theta
d\sigma dr+\\\ +\int_{\mathbb{R}^{6}}&\bigg{(}\frac{\partial}{\partial
r}\phi(\theta,\sigma,r)\big{(}b(\theta-\theta^{\prime})+K(\theta-\theta^{\prime},\sigma,\sigma^{\prime})-b(r-r^{\prime})\big{)}\cdot\\\
&\cdot
w(r-r^{\prime})f_{\gamma}(\theta,\sigma,r,\tau)f_{\gamma}(\theta^{\prime},\sigma^{\prime},r^{\prime},\tau)\bigg{)}d\theta^{\prime}d\sigma^{\prime}dr^{\prime}d\theta
d\sigma r\end{split}$

whereas the remainder ${\tilde{\mathcal{R}}_{\gamma}}$ is given by

(27)
$\displaystyle\begin{split}{\tilde{\mathcal{R}}_{\gamma}}(\phi,\theta,\sigma,r,\tau)=\int_{\mathbb{R}^{3}}\langle{\mathcal{R}_{\gamma}}&(\phi,r^{\prime*},\theta^{\prime},\sigma^{\prime},r^{\prime},\tau)\rangle
w(r^{\prime}-r)f_{\gamma}(\theta^{\prime},\sigma^{\prime},r^{\prime},\tau)d\theta^{\prime}d\sigma^{\prime}dr^{\prime}\\\
+\gamma^{2}\int_{\mathbb{R}^{3}}&\frac{\partial^{2}}{\partial r^{\prime
2}}\phi(\theta^{\prime},\sigma^{\prime},r^{\prime})\bigg{(}\big{(}\langle
S\rangle-b(r^{\prime}-r)\big{)}^{2}+{\operatorname{Var}}[S]\bigg{)}\cdot\\\
&\cdot
w(r-r^{\prime})f_{\gamma}(\theta^{\prime},\sigma^{\prime},r^{\prime},\tau)d\theta^{\prime}d\sigma^{\prime}dr^{\prime}.\end{split}$

All summands will vanish for $\gamma\to 0$ with similar arguments as in [9].
Let us assume that $\phi(\theta,\sigma,r)$ belongs to the space
$\mathcal{C}_{2+\delta}(\mathbb{R}^{3})=\\{h:\mathbb{R}^{3}\rightarrow\mathbb{R},\
\|D^{\zeta}h\|_{\delta}<+\infty\\}$, where $0<\delta\leq 1$, $\zeta$ is a
multi-index with $|\zeta|\leq 2$ and the seminorm $\|\cdot\|_{\delta}$ is the
usual Hölder seminorm

$\displaystyle\|f\|_{\delta}=\sup_{x,y\in\mathbb{R}^{3}}\frac{|f(x)-f(y)|}{|x-y|^{\delta}}.$

Equations (6),(7) together with conservation laws (24) and (25) guarantee the
boundedness of both expectation $\langle S\rangle$ and variance
${\operatorname{Var}}[S]$. Then with this choice of $\phi(\theta,\sigma,r)$,
both summands containing $\frac{\partial}{\partial r^{2}}\phi$ vanish using
the same arguments as in [11, 35].

Therefore, the density $f_{\gamma}(\theta,\sigma,r,\tau)$ converges to
$f(\theta,\sigma,r,\tau)$ which solves

(28) $\displaystyle\begin{aligned}
\frac{d}{d\tau}\int_{\mathbb{R}^{3}}&\phi(\theta,\sigma,r)f(\theta,\sigma,r,\tau)d\theta
d\sigma
dr=\int_{\mathbb{R}^{3}}f(\theta,\sigma,r,\tau)\frac{\partial}{\partial
r}\phi(\theta,\sigma,r)\cdot\\\
&\cdot\bigg{[}\int_{\mathbb{R}^{3}}w(r-r^{\prime})(b(\theta-\theta^{\prime})+K(\theta-\theta^{\prime},\sigma,\sigma^{\prime})-b(r-r^{\prime}))f(\theta^{\prime},\sigma^{\prime},r^{\prime},\tau)\\\
&\qquad\qquad d\theta^{\prime}d\sigma^{\prime}dr^{\prime}\bigg{]}d\theta
d\sigma dr\end{aligned}$

It remains to show that for suitable boundary conditions equation (28) gives
the desired weak formulation of the Fokker-Planck equation. We calculate

$\displaystyle\int_{\mathbb{R}}\bigg{(}$ $\displaystyle
f(\theta,\sigma,r,\tau)\phi(\theta,\sigma,r)\big{(}\int_{\mathbb{R}^{3}}w(r-r^{\prime})\cdot$
$\displaystyle\cdot(b(\theta-\theta^{\prime})+K(\theta-\theta^{\prime},\sigma,\sigma^{\prime})-b(r-r^{\prime}))f(\theta^{\prime},\sigma^{\prime},r^{\prime},\tau)d\theta^{\prime}d\sigma^{\prime}dr^{\prime}\big{)}\bigg{)}_{r=-\infty}^{r=+\infty}drd\sigma$

This term is zero, if

(29) $\displaystyle\lim_{|r|\rightarrow+\infty}f(\theta,\sigma,r,\tau)=0$

These boundary condition are guaranteed for the Boltzmann equation
$f_{\gamma}(\theta,\sigma,r,\tau)$ by mass conservation and preservation of
the first moment $m_{1,r}$, see (23). Then (28) is the weak form of the
Fokker-Planck equation

$\displaystyle\frac{d}{d\tau}\int_{\mathbb{R}^{3}}$
$\displaystyle\phi(\theta,\sigma,r)f(\theta,\sigma,r,\tau)d\theta d\sigma dr=$
(30) $\displaystyle-\int_{\mathbb{R}^{3}}$
$\displaystyle\phi(\theta,\sigma,r)\frac{\partial}{\partial
r}\bigg{[}f(\theta,\sigma,r,\tau)\int_{\mathbb{R}^{3}}w(r-r^{\prime})\cdot$
$\displaystyle\cdot\big{(}b(\theta-\theta^{\prime})+K(\theta-\theta^{\prime},\sigma,\sigma^{\prime})-b(r-r^{\prime}))f(\theta^{\prime},\sigma^{\prime},r^{\prime}\big{)}d\theta^{\prime}d\sigma^{\prime}dr^{\prime}\bigg{]}d\theta
d\sigma dr.$

## References

  * [1] A. Elo  “The Rating of Chessplayers, Past and Present”  Ishi Press, 1986
  * [2] S. Price  “How FIFA’s New Ranking System Will Change International Soccer” Accessed: 2021-10-22  In _Forbes_ , 2018  URL: https://www.forbes.com/sites/steveprice/2018/06/11/how-fifas-new-ranking-system-will-change-international-soccer/?sh=7cc5f8536c41
  * [3] Nate Silver and Reuben Fischer-Baum  “How We Calculate NBA Elo Ratings” Accessed: 2021-10-22, 2015  URL: https://fivethirtyeight.com/features/how-we-calculate-nba-elo-ratings
  * [4] Mark Glickman  “Parameter Estimation in Large Dynamic Paired Comparison Experiments”  In _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ 48, 1998  DOI: 10.1111/1467-9876.00159
  * [5] R. Herbrich, T. Minka and Thore Graepel  “TrueSkill: A Bayesian skill rating system”  In _NIPS_ , 2006, pp. 569–576
  * [6] Tom Minka, Ryan Cleven and Yordan Zaykov  “TrueSkill 2: An improved Bayesian skill rating system” Accessed: 2021-10-22  In _Technical Report_ , 2018  URL: https://www.microsoft.com/en-us/research/uploads/prod/2018/03/trueskill2.pdf
  * [7] Pierre-Emmanuel Jabin and Stephane Junca  “A Continuous Model For Ratings”  In _SIAM Journal on Applied Mathematics_ 75.2 Society for Industrial and Applied Mathematics, 2015, pp. 420–442  DOI: 10.1137/140969324
  * [8] Stéphane Junca  “Contractions to update Elo ratings for round-robin tournaments” Preprint, 2021  URL: https://hal.archives-ouvertes.fr/hal-03286591
  * [9] B. Düring, M. Torregrossa and MT Wolfram  “Boltzmann and Fokker-Planck Equations Modelling the Elo Rating System with Learning Effects”  In _J. Nonlinear Sci._ 29, 2019, pp. 1095–1128  DOI: 10.1007/s00332-018-9512-8
  * [10] Katja Krupp  “Kinetische Modelle für die Rangeinstufung von Spielern”, 2016
  * [11] Giuseppe Toscani  “Kinetic models of opinion formation”  In _Commun. Math. Sci._ 4, 2006
  * [12] Laurent Boudin and Francesco Salvarani  “Opinion dynamics: kinetic modelling with mass media, application to the Scottish independence referendum”  In _Physica A_ 444 Elsevier (North-Holland), Amsterdam, 2016, pp. 448–457  DOI: 10.1016/j.physa.2015.10.014
  * [13] Lorenzo Pareschi, Giuseppe Toscani, Andrea Tosin and Mattia Zanella  “Hydrodynamic models of preference formation in multi-agent societies”  In _J. Nonlinear Sci._ 29.6 Springer US, New York, NY, 2019, pp. 2761–2796  DOI: 10.1007/s00332-019-09558-z
  * [14] Bertram Düring, Ansgar Jüngel and Lara Trussardi  “A kinetic equation for economic value estimation with irrationality and herding”  In _Kinet. Relat. Models_ 10.1 American Institute of Mathematical Sciences (AIMS), Springfield, MO, 2017, pp. 239–261
  * [15] Emiliano Cristiani and Andrea Tosin  “Reducing complexity of multiagent systems with symmetry breaking: an application to opinion dynamics with polls”  In _Multiscale Model. Simul._ 16.1 Society for IndustrialApplied Mathematics (SIAM), Philadelphia, PA, 2018, pp. 528–549  DOI: 10.1137/17M113397X
  * [16] Giacomo Albi, Lorenzo Pareschi, Giuseppe Toscani and Mattia Zanella  “Recent advances in opinion modeling: control and social influence”  In _Active Particles, Volume 1. Modeling and Simulation in Science, Engineering and Technology_ Birkhäuser, Cham, 2017
  * [17] Giacomo Albi, Lorenzo Pareschi and Mattia Zanella  “Boltzmann games in heterogeneous consensus dynamics”  In _J. Stat. Phys._ 175.1 Springer, 2019, pp. 97–125
  * [18] Bertram Düring, Peter Markowich, Jan-Frederik Pietschmann and Marie-Therese Wolfram  “Boltzmann and Fokker-Planck equations modelling opinion formation in the presence of strong leaders”  In _Proc. R. Soc. A._ 465, 2009, pp. 3687–3708  DOI: 10.1098/rspa.2009.0239
  * [19] Bertram Düring and Marie-Therese Wolfram  “Opinion dynamics: inhomogeneous Boltzmann-type equations modelling opinion leadership and political segregation”  In _Proc. R. Soc. A._ 471, 2015, pp. 20150345  DOI: 10.1098/rspa.2015.0345
  * [20] Bertram Düring, Daniel Matthes and Giuseppe Toscani  “Kinetic equations modelling wealth redistribution: A comparison of approaches”  In _Phys. Rev. E_ 78 American Physical Society, 2008, pp. 056103  DOI: 10.1103/PhysRevE.78.056103
  * [21] Lorenzo Pareschi and Giuseppe Toscani  “Wealth distribution and collective knowledge: a Boltzmann approach”  In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 372.2028 The Royal Society Publishing, 2014, pp. 20130396
  * [22] Martin Burger, Alexander Lorz and Marie-Therese Wolfram  “On a Boltzmann Mean Field Model for Knowledge Growth”  In _SIAM Journal on Applied Mathematics_ 76.5, 2016, pp. 1799–1818
  * [23] Martin Burger, Alexander Lorz and Marie-Therese Wolfram  “Balanced growth path solutions of a Boltzmann mean field game model for knowledge growth”  In _Kinetic and Related Models_ 10, 2016  DOI: 10.3934/krm.2017005
  * [24] Lorenzo Pareschi and Giuseppe Toscani  “Interacting multiagent systems. Kinetic equations and Monte Carlo methods”  Oxford University Press, 2013
  * [25] J.K. Blitzstein and J. Hwang  “Introduction to Probability”, Chapman & Hall/CRC Texts in Statistical Science  CRC Press/Taylor & Francis Group, 2014
  * [26] Haym Benaroya and Seon Han  “Probability models in engineering and science”  Taylor & Francis, 2013
  * [27] “Brazil vs. Germany: Goalimpact of Lineups” Accessed: 2021-08-30, http://www.goalimpact.com/blog//2014/07/brazil-vs-germany-goalimpact-of-lineups.html
  * [28] Eitan Tadmor Seung-Yeal Ha  “From particle to kinetic and hydrodynamic descriptions of flocking”  In _Kinetic and Related Models_ 1, 2008
  * [29] Dario Benedetto, Emanuele Caglioti and Mario Pulvirenti  “A kinetic equation for granular media”  In _ESAIM: Mathematical Modelling and Numerical Analysis_ 31.5 EDP Sciences, 1997, pp. 615–641
  * [30] John Towers  “Convergence Of A Difference Scheme For Conservation Laws With A Discontinuous Flux”  In _SIAM Journal on Numerical Analysis_ 38, 1999
  * [31] Federation Internationale Des Echecs  “The Official Laws of Chess and Other Fide Regulations”  Collier Books, 1990
  * [32] K. Harkness  “Official Chess Handbook”  D. McKay Co., 1967
  * [33] C. Cercignani  “The Boltzmann equation and its applications”  Springer, New York, NY, 1988
  * [34] Carlo Cercignani, Reinhard Illner and Mario Pulvirenti  “The mathematical theory of dilute gases”  Springer Science & Business Media, 2013
  * [35] Stephane Cordier, Lorenzo Pareschi and Cyrille Piatecki  “Mesoscopic Modelling of Financial Markets”  In _Journal of Statistical Physics_ 134, 2009, pp. 161–184


# The Nascent Red-Sequence at ${\mathbf{z}\sim 2}$11affiliation: Based on
observations with the NASA/ESA Hubble Space Telescope, obtained at the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555

Andrew W. Zirm22affiliation: Johns Hopkins University, 3400 N. Charles Street,
Baltimore, MD, 21218<EMAIL_ADDRESS>, S.A. Stanford33affiliation:
Institute of Geophysics and Planetary Physics, LLNL, L-413, Livermore, CA
94550 44affiliation: Physics Dept., University California at Davis, One
Shields Avenue, Davis, CA 95616 , M. Postman55affiliation: Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218 , R.A.
Overzier22affiliation: Johns Hopkins University, 3400 N. Charles Street,
Baltimore, MD, 21218<EMAIL_ADDRESS>, J.P. Blakeslee66affiliation:
Department of Physics and Astronomy, Washington State University, Pullman, WA
99164 77affiliation: Herzberg Institute of Astrophysics, NRC of Canada, 5071
W. Saanich Road, Victoria, BC, V9E 2E7, Canada , P. Rosati88affiliation: ESO-
European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748, Garching
bei München, Germany , J. Kurk99affiliation: Max-Planck-Institut für
Astronomie, D-69117 Heidelberg, Germany , L. Pentericci1010affiliation:
Osservatorio Astronomico di Roma, I-00040 Monte Porzio Catone, Italy , B.
Venemans1111affiliation: IOA, Cambridge University, Madingley Road, Cambridge,
CB3 0HA, UK , G.K. Miley1212affiliation: Leiden Observatory, Leiden
University, Postbus 9513, NL-2300 RA Leiden, The Netherlands , H.J.A.
Röttgering, M. Franx, A. van der Wel22affiliation: Johns Hopkins University,
3400 N. Charles Street, Baltimore, MD, 21218<EMAIL_ADDRESS>, R.
Demarco22affiliation: Johns Hopkins University, 3400 N. Charles Street,
Baltimore, MD, 21218<EMAIL_ADDRESS>, W. van Breugel33affiliation:
Institute of Geophysics and Planetary Physics, LLNL, L-413, Livermore, CA
94550 1313affiliation: University of California, Merced, PO Box 2039, Merced,
CA 95344

###### Abstract

We present new constraints on the evolution of the early-type galaxy color-
magnitude relation (CMR) based on deep near-infrared imaging of a galaxy
protocluster at $z=2.16$ obtained using NICMOS on-board the Hubble Space
Telescope. This field contains a spectroscopically confirmed space-overdensity
of Lyman-$\alpha$ and H-$\alpha$ emitting galaxies which surrounds the
powerful radio galaxy MRC 1138-262. Using these NICMOS data we identify a
significant surface-overdensity ($=6.2\times$) of red $J_{110}-H_{160}$
galaxies in the color-magnitude diagram (when compared with deep NICMOS
imaging from the HDF-N and UDF). The optical-NIR colors of these prospective
red-sequence galaxies indicate the presence of on-going dust-obscured star-
formation or recently formed ($\mathrel{\hbox to0.0pt{\lower
3.0pt\hbox{$\mathchar 536\relax$}\hss}\raise 2.0pt\hbox{$\mathchar
316\relax$}}1.5$ Gyr) stellar populations in a majority of the red galaxies.
We measure the slope and intrinsic scatter of the CMR for three different red
galaxy samples selected by a wide color cut, and using photometric redshifts
both with and without restrictions on rest-frame optical morphology. In all
three cases both the rest-frame $U-B$ slope and intrinsic color scatter are
considerably higher than corresponding values for lower redshift galaxy
clusters. These results suggest that while some relatively quiescent galaxies
do exist in this protocluster both the majority of the galaxy population and
hence the color-magnitude relation are still in the process of forming, as
expected.

###### Subject headings:

galaxies: evolution — galaxies: formation — galaxies: high-redshift —
galaxies: stellar content — galaxies: clusters: individual: MRC 1138-262

††journal: Accepted for publication in the Astrophysical Journal

##  1\. Introduction

The color-magnitude diagram is a powerful diagnostic of galaxy evolution and
formation. The presence, as early as $z\sim 1.5$, of a prominent and low-
scatter, ‘red-sequence’ (RS) in galaxy clusters places useful constraints on
the possible evolutionary pathways in galaxy color and luminosity (Mullis et
al., 2005; Stanford et al., 2005, 2006; Bell et al., 2004; Faber et al.,
2007). The red colors of the primarily early-type RS galaxies are due to the
observed filters spanning the 4000Å spectral break. The universality and
prominence of the RS in appropriately chosen filters have been used to
discover high-redshift clusters (e.g., Gladders & Yee, 2005). Moreover, the
defining characteristic of galaxy clusters, i.e., the large numbers of
galaxies all at the same redshift, allows the slope and intrinsic scatter of
the RS to be measured with great precision. Based on studies of galaxy
clusters at $z<1.3$, the slope of the RS does not appear to evolve and
therefore is more likely the by-product of the mass-metallicity relation as
observed in local galaxy samples (e.g., Tremonti et al., 2004) rather than the
result of a mass-age trend. The scatter, however, is likely due to the
fractional age differences between the RS galaxies (e.g., Blakeslee et al.,
2003). By constructing a set of model galaxies with different star-formation
histories and timescales it is possible to fit for the mean epoch of last
significant star-formation by matching the intrinsic scatter of the RS. Such
studies at $z\sim 1$ have derived formation redshifts of $z_{\rm form}\sim
2.0-2.5$ (e.g., Blakeslee et al., 2006; van Dokkum & van der Marel, 2007). At
redshifts beyond $z\sim 1.5$, however, the 4000Å-break moves into the near-
infrared and galaxy clusters, and therefore the RS, have not been observed
closer to the inferred epoch of formation for early-type galaxies. Hence, to
uncover the younger or forming red-sequence at higher redshifts requires deep
near-infrared imaging of suspected (or, preferably, confirmed) protocluster
fields.

Figure 1.— $J_{110}-H_{160}$ vs. $H_{160}$ color-magnitude diagram for the MRC
1138-262 NICMOS field (open black circles). The large yellow star is the radio
galaxy itself. The blue background points are from the NICMOS data of the
Ultra Deep Field and Hubble Deep Field North. The deep fields cover 2.5x the
area of our observations. Also shown are the spectroscopically confirmed
$H\alpha$ (orange filled circles) and Lyman-$\alpha$ (green filled circles)
emitters. The three dot-dash lines show where the CMRs of lower redshift
clusters would lie under different assumptions. The top line is the Coma
cluster with no evolution, simply redshifted to $z=2.16$ and observed through
the NICMOS filters. The next line down is the same but for the $z=1.24$
cluster RDCS1252. Finally, if we passively de-evolve RDCS1252 to redshift
$z=2.16$ (almost exactly two Gyr), assuming a median age for the 1252 galaxies
of about 3 Gyr (or $z_{\rm form}\sim 3$), we get the third line.

We have undertaken a NICMOS imaging program to study the red galaxy population
in a protocluster at $z=2.16$. Broad and narrow-band imaging, both in the
optical and near-infrared, of the field surrounding the powerful radio galaxy
MRC 1138-262 ($z=2.16$) have identified more than 100 candidate companion
galaxies. This target served as the proof-of-concept for the successful VLT
Large Program summarized in Venemans et al. (2007). There are surface-
overdensities of both line-emitting candidates (Lyman-$\alpha$ and H$\alpha$),
X-ray point sources, sub-mm selected galaxies and red optical–near-infrared
galaxies (Pentericci et al., 2002; Kurk, 2003; Kurk et al., 2004a, b; Croft et
al., 2005; Stevens et al., 2003). Fifteen of the Ly$\alpha$ and 9 of the
H$\alpha$ emitters have now been spectroscopically confirmed to lie at the
same redshift as the radio galaxy. The $I-K$-selected Extremely Red Objects
(EROs; $I-K>4.3$ Vega) seem concentrated around the RG but have no
spectroscopic redshifts at this time. However, by obtaining deep images
through the NICMOS $J_{110}$ and $H_{160}$ filters, which effectively span the
4000Å-break at $z=2.16$, accurate and precise colors and basic morphological
parameters can be measured for the red galaxy population. In this paper we
present the first results from this project. The article is organized as
follows: in Section §2 we describe the data and their reductions, in Section
§5 we present the comparison between the red galaxy counts in this field and
in deep field data, in Section §6 we present the full color-magnitude diagram
and our fits to the “red sequence.” We use a
$(\Omega_{\Lambda},\Omega_{M})=(0.73,0.27)$, $H_{0}=71$ ${\rm km}$ ${\rm
s^{-1}}$ ${\rm Mpc^{-1}}$ cosmology throughout. At $z=2.16$ one arcsecond is
equivalent to 8.4 kpc. All magnitudes are referenced to the AB system (Oke,
1974) unless otherwise noted.

##  2\. Observations, Data Reductions and Photometry

Figure 2.— $g_{435}-I_{814}$ vs. $J_{110}-H_{160}$ color-color diagram using
the ACS and NICMOS data. Arrows represent limits where the galaxy is only
detected in a single band for that color. Filled circles indicate
spectroscopically-confirmed Lyman-$\alpha$ (green) and H-$\alpha$ (orange)
emitting protocluster members. The yellow star is the radio galaxy. The blue,
green and red grids indicate the regions occupied by galaxies with an
exponentially-decaying star-formation rate $\tau=0.15$ Gyr (red), $\tau=0.4$
Gyr (green) and a $\tau=1000$ Gyr (blue) at $z=2.16$ for three ages (0.1, 1
and 3 Gyr) and two different extinctions ($E(B-V)=0.0,1.0$).

The NICMOS instrument on-board HST is capable of deep near-infrared imaging
more quickly than from the ground but with a relatively small field-of-view
($51\arcsec\times 51\arcsec$). In the case of MRC 1138-262, we know that
galaxies are overdense on the scale of a few arcminutes (Kurk et al., 2004a;
Croft et al., 2005) and are thus well-suited for observations with NICMOS
camera 3 on HST. We used 30 orbits of HST time to image 10 of the 24 confirmed
members and $\sim 70$ of the candidate (narrow-band excess sources and EROs)
protocluster members in both the $J_{110}$ and $H_{160}$ filters. We used
seven pointings of NICMOS camera 3 in both filters and one additional pointing
in $H_{160}$ alone. This single ‘outrigger’ $H_{160}$ pointing was included to
obtain rest-frame optical morphological information for a small concentration
of candidate members. These observations reach an AB limiting magnitude
($m_{10\sigma}$; 10$\sigma$, $0\farcs 5$ diameter circular aperture) of
$m_{10\sigma}=24.9$ mag in $J_{110}$ and $m_{10\sigma}=25.1$ mag in $H_{160}$.
The same field was imaged in the $g_{475}$ ($m_{10\sigma}=27.5$ mag) and
$I_{814}$ ($m_{10\sigma}=26.8$ mag) filters using the Wide-Field Channel of
the Advanced Camera for Surveys on HST as part of a Guaranteed Time program (#
10327; Miley et al. 2006). These optical data are useful for their higher
angular resolution and their coverage of the rest-frame far-UV, thus extending
the observed SEDs of candidate protocluster members to shorter wavelengths
where young stars and on-going star-formation dominate the emitted spectrum.
In particular, the $g_{475}$ and $I_{814}$ data allow us to partially
differentiate obscured star-formation from evolved stellar populations in the
candidate RS galaxies.

The NICMOS images were reduced using the on-the-fly reductions from the HST
archive, the IRAF task PEDSKY and the dither/drizzle package to combine the
images in a mosaic. The dither offsets were calculated using image cross-
correlation and were refined with one further iteration of cross-correlation.
Alignment of the pointings relative to each other was accomplished using a
rebinned version of the ACS $I_{814}$ image as a reference. The final mosaic
has a pixel scale of $0\farcs 1$. Galaxies were selected using the
$H_{160}$-band image for detection within SExtractor (Bertin & Arnouts, 1996).
We used a $2.2\sigma$ detection threshold with a minimum connected area of 10
pixels. We also corrected the NICMOS data for the count-rate dependent non-
linearity (de Jong, 2006). Total galaxy magnitudes were estimated by using the
MAG_AUTO values from SExtractor.

The $J_{110}-H_{160}$ colors were determined by running SExtractor (Bertin &
Arnouts, 1996) in two-image mode using the $H_{160}$ image for object
detection and isophotal apertures. The $J_{110}$ image was PSF-matched to the
$H_{160}$ band. The resulting colors and magnitudes are shown in Figure 1. For
galaxies which are not detected at $2\sigma$ significance in the
$J_{110}$-band (those to the right of the thick dashed line, representing
$J_{110,tot}>26.7$, in Fig. 1) we consider the color to be a lower limit.

We also measured similarly PSF-matched, isophotal colors using the two ACS
bands and have used them to construct a $g_{475}-I_{814}$ versus
$J_{110}-H_{160}$ color-color diagram (Figure 2). We compared these colors to
those of model SEDs for different ages, star-formation histories and dust
extinctions. Using the 2007 Charlot & Bruzual population synthesis models we
have constructed spectral energy distributions for galaxies with an
exponentially-decaying star-formation rate with time constants of
$\tau=0.15,0.4,1000.0$ Gyr (the red, green and blue grids in Fig. 2
respectively). Each model’s colors are calculated for ages of 0.1, 1 and 3 Gyr
and for $E(B-V)=0.0$ and $1.0$. Aging of the population moves primarily the
$J_{110}-H_{160}$ color to the red while the dust extinction significantly
reddens the $g_{475}-I_{814}$ color. From this analysis it appears that the
$\tau=0.4$ Gyr model represents well the colors of a majority of the red
$J_{110}-H_{160}$ galaxies.

To extend the wavelength coverage for the protocluster galaxies we also
incorporated ground-based $U_{n}$-band data from LRIS-B on the Keck telescope,
$K_{s}$-band imaging from VLT/ISAAC and three band IRAC imaging (the 3.6, 4.5
and 5.8 $\mu$m bands) from the Spitzer Space Telescope. The Keck $U$-band data
(PI W. van Breugel) were obtained in late January and early February of 2003.
The ISAAC data (PI G. Miley) were taken in Period 73 in service mode. The
Spitzer data are from the IRAC Guaranteed Time program (PI G. Fazio, Program
#17). We have smoothed the imaging data for all bands, apart from the IRAC
data, to match the resolution of the $U_{n}$-band image (approximated by a
FWHM $\sim 1\arcsec$ Gaussian). We then used SExtractor to measure galaxy
magnitudes within a $0\farcs 5$ radius circular aperture for each smoothed
image. To incorporate the IRAC data, which has much poorer angular resolution,
we derived aperture magnitudes which were then corrected to match the smoothed
data. These aperture corrections were derived using the photometric curves-of-
growth for 20 stars in the field. The resulting catalog was used to generate
photometric redshift estimates as described below in Section § 3.

Figure 3.— upper panel: Distribution of high-confidence ($>95\%$) photometric
redshifts and their selection function assuming a uniform $N(z)$ for our model
template galaxies (yellow curve) for $1.1\leq(J_{110}-H_{160})\leq 2.1$
galaxies in the MRC 1138-262 NICMOS field. The peak between $z=2.1$ and
$z=2.4$ is statistically highly significant. middle panel: Sum of redshift
probability distributions for all the galaxies in the upper panel. 38.5% of
the total probability is contained in the redshift interval from 2 to 2.3.
lower 2 panels: two examples of the probability distribution function for
individual galaxies

##  3\. Photometric Redshifts

We have used the ACS ($g_{475}$, $I_{814}$), NICMOS ($J_{110}$, $H_{160}$),
ground-based $U_{n}$-band from Keck/LRIS-B, $K_{S}$-band imaging from
VLT/ISAAC and Spitzer/IRAC imaging to estimate photometric redshifts for our
$H_{160}$-band selected sample. We input a catalog of aperture galaxy
magnitudes, based on the matched, smoothed images described above, into the
Bayesian photometric redshift code (BPZ) of Benítez (2000) using a uniform
prior. We felt that the default prior, based on optical galaxy selection and
spectroscopy in the HDF-N, would not necessarily represent the redshift
distribution for our near-infrared selected galaxies. We generated our own
extensive set of template spectral energy distributions using the models of
Charlot & Bruzual (2007). All these SEDs are $\tau$ models with values for
$\tau=[0.15,0.4,1.0,2.0,1000.0]$ Gyr and ages $=[0.05,0.1,0.5,1.0,2.0,3.0]$
Gyr. We also included models with internal dust extinction ranging from
$E(B-V)=[0.0,0.1,0.3,0.5,0.75,1.0]$ mag and metallicity of
$(Z/Z_{\odot})=[0.3,1.0,2.5]$. We focused particular attention on the
$J_{110}-H_{160}$ selected surface-overdensity. In the upper panel of Figure 3
we present the high confidence ($>95\%$) photo-$z$ distribution for the NIR-
color selected ($1.1\leq(J_{110}-H_{160})\leq 2.1$) subsample. We ran
extensive simulations by redshifting our template set, adding appropriate
photometric errors and using BPZ to recover the redshifts. The yellow curve
represents the redshift selection function for this color cut, template set
and filters assuming that these model galaxies follow a uniform $N(z)$ over
this redshift interval. The simulation results were free of significant
systematic errors and the random errors are estimated to be $\delta z/z\sim
0.1$. Based on these SED fits, the approximate luminosity-weighted ages of the
red galaxies lie between 1 and 2.5 Gyrs and their stellar masses are of order
a few $\times 10^{10}M_{\odot}$. These stellar masses are reasonable as are
the absolute magnitudes (see Figure 6). More detailed SED modeling is deferred
to a future paper.

There is a clear excess of galaxies between $z=2.0$ and $z=2.5$. For each
galaxy fit by BPZ we have generated the full redshift probability
distribution. In the lower panel of Fig. 3 we show the $H_{160}$-band
weighted-average of these probability distributions. There is a clear peak
(containing 38.5% of the total probability compared to only 17% of the total
selection function in the same redshift interval) between $z=2.0$ and $z=2.3$,
consistent with the significant peak in the redshift histogram itself.

##  4\. NICMOS Galaxy Morphologies

Figure 4.— Distribution of inferred stellar ages (in terms of $\tau$) for both
the concentrated ($n\geq 2.5$, red line) and diffuse ($n<2.5$, blue line)
galaxies which are well-resolved in the NICMOS data. Constant star-formation
models, for which the $e$-folding time is infinite, are placed at the left-
hand edge of the plot. The blue and red distributions are quite different. Of
particular note is that the most evolved galaxies generally have high $n$
while the low $n$ galaxies dominate the star-forming population.

NICMOS camera 3 provides good angular resolution over its entire field-of-
view. The FWHM of the PSF in our final mosaic is $\approx 0\farcs 27$. To
exploit this resolution we have used the GALFIT code (Peng et al., 2002) to
fit analytic Sérsic surface-brightness profiles (Sersic, 1968) to all the
$H_{160}\leq 24.5$ sources in our $H_{160}$-band mosaic. A model point-spread
function was created for each of these galaxies individually by generating a
TinyTim simulated PSF (Krist, 1993) at the galaxies’ positions in each
exposure and then drizzling these PSFs together in exactly the same fashion as
for the data themselves (see Zirm et al. 2007). We restricted the Sérsic
index, $n$, to be between 1 and 5. We will present a full analysis of the
morphologies of these galaxies in a future paper. For the current work, we use
these derived sizes and profile shapes to assist us in selecting the
morphological “early-type” members of the red galaxy population.

Figure 5.— Histogram of the color distributions for the 1138 and deep fields
(blue). The deep field data has been normalized by total area to the 1138
data. Note the clear excess of red galaxies in the 1138 field. At
$1.1\leq(J_{110}-H_{160})\leq 2.1$ (horizontal dotted lines) for galaxies
brighter than the 2$\sigma$ $J_{110}$-band limit (dashed line) there is an
overdensity of a factor of 6.2 in the 1138 field.

In Figure 4 we show the distribution of galaxy ages derived via these SED fits
as parametrized by the $\tau$ value for the best-fitting model for those
galaxies with high and low Sérsic index ($n\geq 2.5$, red line, and $n<2.5$,
blue line). It is clear that while there is substantial overlap between these
distributions they are not identical and that they differ in the sense that
one might expect, namely, that the concentrated galaxies appear to be
comprised of older stellar populations. This trend gives us some confidence in
trying to select the “early-type” galaxies using these data which is important
for our discussion of the color-magnitude relation in Section 6.

##  5\. Surface Overdensity of Red Galaxies

Figure 6.— Linear fits to the rest-frame $U-B$ (Vega) color-magnitude diagrams
for three different sub-samples of the $H_{160}$-band selected NICMOS sample.
Panel (A) shows the fit (solid line) and intrinsic scatter ($1\sigma$; dotted
lines) for a sample selected to have $1.1\leq(J_{110}-H_{160})\leq 2.1$. The
crossed-out point are those which are rejected as outliers in more than half
of the realizations (see §6). Both the observed and intrinsic scatter are
smaller than the initial color cut. Panel (B) shows the fit and intrinsic
scatter for a photometric redshift selected sample with $2.0<z_{\rm
phot}<2.5$. The stars indicate galaxies whose preferred photometric template
has an age $<4\tau$, while circles represent galaxies with older than this
limit. Panel (C) shows the fit and scatter for those galaxies which meet the
same redshift cut but also are well-resolved with a high Sérsic index
($n>2.5$) and best-fit by an age $\geq 4\tau$ template.

To compare this protocluster field to more generic ‘blank’ field data we have
compiled catalogs for the public NICMOS data in both the Hubble Deep Field
North (HDF-N) and the Ultra Deep Field (UDF). Figure 1 shows the
$J_{110}-H_{160}$ color-magnitude diagram (open black circles) and the color
distributions for both the MRC 1138-262 and the combined HDF-N and UDF galaxy
catalogs (blue circles). The deep field data were also $H_{160}$-band
selected. The area of the two deep fields is roughly 2.5 times the area of our
protocluster observations. We have applied no correction to the deep field
number counts to account for clustering in those fields. The color histogram
in Figure 5 shows the area-normalized galaxy counts from the two deep fields
(blue line) and from the 1138 field to the same ($2\sigma$) limiting magnitude
of $J_{110}=26.7$ mag (AB). The red dashed line shows the difference between
the two color distributions. It is clear that the radio galaxy field is
overdense in red galaxies by a large factor. For sources with colors between
$1.1\leq(J_{110}-H_{160})\leq 2.1$, the horizontal(vertical) dotted lines in
Fig.1(5), and brighter than our $J_{110}$-band $2\sigma$ limit ($26.7$), we
calculate an area-normalized overdensity of $6.2$ when compared to the deep
fields data. We note that the exact value of the measured overdensity is
rather sensitive to systematic color offsets between the protocluster and deep
field data. A redward shift of 0.05 for the deep field galaxies would lower
the measured overdensity to $5.0$. However, we are confident that these
systematic offsets remain small ($<0.05$ mag) since we have used the same
instrument, filters, selection technique and photometric code with very
similar input parameters for both the deep field and 1138 datasets. Looking
back at Figure 2 we can see that many of the spectroscopically-confirmed line
emitters (filled blue circles) and red galaxies in the overdensity are well-
represented by the $\tau=0.4$ Gyr model (green lines) at different ages and
extinctions.

This current work is not the first to observe red galaxies in this field. Kurk
et al. (2004) identified a small ($\sim 1.5\times$) surface overdensity of
extremely red objects (EROs; $I-K>4.3$ Vega magnitudes) using ground-based $I$
and $K$ band data. Many of these EROs are also identified as red in the NICMOS
$J_{110}-H_{160}$ color. More recently Kodama et al. (2007) observed this
field using the wide-field NIR imager, MOIRCS, on the Subaru telescope. These
authors found several bright (presumably massive) red galaxies over a wider
field-of-view but to shallower depths than the NICMOS data presented here. 24
of their color-selected protocluster candidates are within our NICMOS mosaic.
23 of the 24 are identified in our data as being red in $J_{110}-H_{160}$.
Furthermore, 18 of the 94 galaxies which satisfy our color criteria (and have
$J_{110}<26.67$) are also identified by Kodama et al. as protocluster
candidates. The much larger number of red galaxies in the NICMOS data is
primarily due to fainter galaxies detected at high significance in our deeper
data.

##  6\. An Emergent Red-Sequence?

To study the colors and magnitudes of these galaxies in more detail and to
possibly identify a red-sequence in the 1138 field we have split the galaxies
into three sub-samples defined by $J_{110}-H_{160}$ color, photometric
redshift and morphology (Sérsic index). The first sample (Sample A) comprises
all 56 galaxies with $1.1\leq(J_{110}-H_{160})\leq 2.1$ and $H_{160}<24.5$ and
includes the radio galaxy itself. Sample B is made up of all 28 galaxies with
a robust photometric redshift between 2.0 and 2.5 and $J_{110}-H_{160}>0.75$
and $H_{160}<26.0$. This liberal color cut is included to select galaxies
which comprise the large observed surface-overdensity. Finally, sample C
contains seven galaxies with the same photometric redshift cut but which also
have well-resolved $H_{160}$-band surface-brightness profiles with Sérsic
index $n>2.5$. All of these galaxies’ SEDs are also best-fit by models with
relatively little on-going star-formation. We use a limit of (age $\geq
4\times\tau$, cf. Grazian et al 2007 ). Therefore, sample C mimics the color,
morphological and photometric redshift selection of early-type galaxies in
clusters at $z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{$\mathchar
536\relax$}\hss}\raise 2.0pt\hbox{$\mathchar 316\relax$}}1$. The photometry,
photo-$z$s and sizes of the sample C galaxies are listed in Table 1, their
rest-frame color-magnitude distributions are shown in Figure 6 and the two-
dimensional spatial distribution of the Sample A galaxies is plotted in Figure
7. We note that because the measured overdensity is a factor of 6, we
statistically expect one of every seven sample A galaxies to be a field
galaxy. However, this should not effect our results significantly.

For these three sample selections we have fit a line and measured the
intrinsic scatter about that best-fit line (see Fig. 6). For comparison to
lower redshift galaxy clusters we have transformed our observed
$J_{110}-H_{160}$ color and $H_{160}$ magnitudes into rest-frame $U-B$ and $B$
(Vega), respectively, using the following expressions:

$(U-B)_{\rm rest}=0.539\times(J_{110}-H_{160})_{\rm obs}-0.653$ (1)

$M_{B,{\rm rest}}=H_{160,{\rm obs}}-0.170\times(J_{110}-H_{160})_{\rm
obs}-43.625$ (2)

The small color corrections used in these relations were derived using a
family of $\tau$-models with a range of ages (0.1-12 Gyr), $\tau$ (0.1-5 Gyr)
and three metallicities ($0.4,1$ and $2.5Z_{\odot}$).

Figure 7.— Spatial distribution of galaxies (red circles) relative to the
radio galaxy MRC 1138-262 (yellow star) for the color-selected sample defined
in the text (A) and shown in the first panel of Fig. 6. The irregular black
outline encloses the coverage of the 7 NICMOS camera 3 pointing with both
$J_{110}$ and $H_{160}$ imaging.

To fit the “CMR” we used a bootstrap re-sampling technique to estimate the
error on the fitted slope. Then, by assuming that all the red galaxies lie on
this fit line, we ran Monte Carlo realizations of the contribution of the
photometric errors to the observed color scatter about the fit line, i.e., by
fixing a color-magnitude relation we calculate the measurement scatter with
zero intrinsic scatter. We then calculate the intrinsic scatter by subtracting
(in quadrature) the estimated measurement scatter from the observed scatter.
We show these fits (solid line) and the intrinsic scatters (dotted lines) for
the three samples (A, B and C) in Figure 6. The fits to both Sample’s A and B
have nearly identical rest-frame $U-B$ slopes, $0.027$ and $0.026$
respectively, and intrinsic scatters ($0.10$ and $0.12$). While these slopes
are comparable to those found for the well-populated lower redshift cluster
CMRs, the intrinsic scatters are considerably higher. However, the scatter
measured for the eight galaxy Sample C is comparable to that of the lower
redshift samples but with a much steeper slope ($0.130$). When these scatters
are compared to model predictions based on lower redshift clusters
(specifically RDCS 1252.9-2927 at $z=1.24$; Gobat et al. 2008) we find that
the 1138 protocluster has lower than predicted scatter. This may suggest that
the 1138 protocluster is in a more advanced evolutionary state than RDCS 1252
was at $z=2.2$.

We have calculated three representative color-magnitude relations for
comparison to the colors and magnitudes of the red galaxies (three dot-dash
lines in Fig. 1). We have taken two lower redshift clusters, Coma at $z=0.023$
and RDCS 1252.9-2927 at $z=1.24$, and transformed them to the observed filters
and $z=2.16$ under the assumption that the colors do not evolve. In this no
evolution case (the two dot-dash lines in Fig. 1) the CMRs appear at the red
edge of the observed overdensity. There is almost exactly 2 Gyr of cosmic time
between $z=2.16$ and $z=1.24$ in our adopted cosmology. From Blakeslee et al.
(2003) we know that the median redshift of last significant star-formation for
the RDCS 1252 galaxies is between $z=2.7-3.6$. Therefore, if we observe those
galaxies at $z=2.16$ they will be significantly younger and hence bluer. In
fact, this passively de-evolved line (bluest dot-dash line in Fig. 1, labeled
‘$z_{\rm form}\sim 3$’) does fall blueward of the red galaxy overdensity. We
discuss the implications for these comparisons in Section 7.

We have also translated the Kodama et al. ground-based $J-K$ colors to our
NICMOS filters assuming all the red galaxies lie at $z=2.16$. These bright
galaxies also fall along the passively de-evolved line with the radio galaxy.
We have used our suite of SED models to estimate the color transformation from
their ground-based $J-K_{S}$ to our NICMOS $J_{110}-H_{160}$ color. Roughly,
the Kodama et al. bright red galaxies fall where the RDCS 1252 passive line
crosses our color cut at $J_{110}-H_{160}=1.1$. This result hints at a
possible bi-modality in the red galaxy population of this protocluster.
Namely, that there are faint red galaxies that are inconsistent with
passively-evolving cluster members either due to large amounts of dust, or due
to higher redshifts of formation but that the more luminous protocluster
members may have already finished forming and seem consistent with passive
evolution to the present-day.

##  7\. Discussion

We have identified a (6.2$\times$) surface-overdensity and a corresponding
photometric redshift ‘spike’ of red $J_{110}-H_{160}$ galaxies which are
likely associated with a known protocluster at $z=2.16$. The optical-NIR
spectral energy distributions of these sources suggest that they comprise both
evolved galaxies as well as dust-obscured star-forming galaxies. Based on our
SED fits from the photometric redshift determinations, the approximate
luminosity-weighted ages of these sources lie between 1 and 2.5 Gyrs and their
stellar masses are of order a few $\times 10^{10}M_{\odot}$. Detailed modeling
of the SEDs for the protocluster population, along with their morphologies, is
reserved for a future paper.

Comparison with the CMRs of lower redshift clusters shows that the red galaxy
overdensity primarily lies blueward of the no-evolution predictions. That the
red galaxies in 1138 are also redder than the $z_{\rm form}\sim 3$ case
suggests both that there are galaxies with significant dust content, an
assertion supported by the SED fits, and also that they were perhaps formed at
higher redshift than the RDCS1252 galaxies. Of course, without a classical,
low-scatter red-sequence to use as a baseline there remains considerable
uncertainty in the age of the population as a whole. The results of Steidel et
al. (2005) suggest that protocluster galaxies are older than their “field”
counterparts at $z\sim 2.3$ and that these ages and stellar masses were
broadly consistent with evolution to lower redshift cluster galaxies. However,
their protocluster members were all UV-selected and star-forming. With future
spectroscopy of our red galaxy sample it will be possible to see if these
differences persist when looking at a more varied galaxy sample.

For three samples of galaxies drawn from the full $H_{160}$-band selected
dataset we have fit a color-magnitude relation and estimated the intrinsic
scatter. The CMR at $z=2.16$ is not as well-defined as at $z\sim 1$ or 0. For
sample C, made up of 8 galaxies, the color, best-fit spectral template,
morphology and photo-$z$ all point towards them being (proto-)elliptical
galaxies within the protocluster. For this small sample, the estimated
intrinsic scatter is rather low and may suggest that these galaxies represent
the forming red-sequence in this protocluster. The slope of this relation is
extremely steep compared to lower redshift clusters. The slope of the CMR is
generally assumed to be a manifestation of the mass-metallicity relation and
would therefore flatten at higher redshift. The major caveat regarding the
steep slope of Sample C is that none of these galaxies are spectroscopically
confirmed protocluster members. Therefore, this “relation” may just be a
random, although somewhat unlikely, coincidence rather than a nascent CMR.
However, further deep NIR imaging coverage of this field is required to
identify additional members of this proto-elliptical galaxy class.

Support for program # 10404 was provided by NASA through a grant
(GO-10404.01-A) from the Space Telescope Science Institute, which is operated
by the Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-26555. The work of SAS was performed in part under the auspices
of the U.S. Department of Energy, National Nuclear Security Administration by
the University of California, Lawrence Livermore National Laboratory under
contract No. W-7405-Eng-48. JK is financially supported by the DFG, grant SFB
439\. WvB acknowledges support for radio galaxy studies at UC Merced,
including the work reported here, with the Hubble Space Telescope and the
Spitzer Space Telescope via NASA grants HST # 10127, SST # 1264353, SST #
1265551, SST # 1279182.

## References

  * Bell et al. (2004) Bell, E. F., et al. 2004, ApJ, 608, 752
  * Benítez (2000) Benítez, N. 2000, ApJ, 536, 571
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Blakeslee et al. (2003) Blakeslee, J. P., et al. 2003, ApJ, 596, L143
  * Blakeslee et al. (2006) Blakeslee, J. P., et al. 2006, ApJ, 644, 30
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Charlot & Bruzual (2007) Charlot, S. & Bruzual, G. 2007, in prep
  * Croft et al. (2005) Croft, S., Kurk, J., van Breugel, W., Stanford, S. A., de Vries, W., Pentericci, L., & Röttgering, H. 2005, AJ, 130, 867
  * de Jong (2006) de Jong, R. 2006.  Correcting the nicmos count-rate dependent non-linearity.  Technical report, STScI
  * Faber et al. (2007) Faber, S. M., et al. 2007, ApJ, 665, 265
  * Gladders & Yee (2005) Gladders, M. D. & Yee, H. K. C. 2005, ApJS, 157, 1
  * Gobat et al. (2008) Gobat et al., 2008, submitted to A&A
  * Grazian et al. (2007) Grazian, A., et al. 2007, A&A, 465, 393
  * Kodama et al. (2007) Kodama, T., Tanaka, I., Kajisawa, M., Kurk, J., Venemans, B., De Breuck, C., Vernet, J., & Lidman, C. 2007, MNRAS, 377, 1717
  * Krist (1993) Krist, J. 1993, in ASP Conf. Ser. 52: Astronomical Data Analysis Software and Systems II 536
  * Kurk (2003) Kurk, J. D. 2003, May.  The cluster environments and gaseous halos of distant radio galaxies.  Ph. D. thesis, Leiden University
  * Kurk et al. (2004a) Kurk, J. D., Pentericci, L., Röttgering, H. J. A., & Miley, G. K. 2004a, A&A, 428, 793
  * Kurk et al. (2004b) Kurk, J. D., Pentericci, L., Overzier, R. A., Röttgering, H. J. A., & Miley, G. K. 2004b, A&A, 428, 817
  * Miley et al. (2006) Miley, G. K., et al. 2006, ApJ, 650, L29
  * Mullis et al. (2005) Mullis, C. R., Rosati, P., Lamer, G., Böhringer, H., Schwope, A., Schuecker, P., & Fassbender, R. 2005, ApJ, 623, 85
  * Oke (1974) Oke, J. B. 1974, ApJS, 27, 21
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266
  * Pentericci et al. (2002) Pentericci, L., Kurk, J. D., Carilli, C. L., Harris, D. E., Miley, G. K., & Röttgering, H. J. A. 2002, A&A, 396, 109
  * Sersic (1968) Sersic, J. L. 1968, Atlas de galaxias australes Cordoba, Argentina: Observatorio Astronomico, 1968
  * Stanford et al. (2005) Stanford, S. A., et al. 2005, ApJ, 634, 129
  * Stanford et al. (2006) Stanford, S. A., et al. 2006, ApJ, 646, 13
  * Steidel et al. (2005) Steidel, C. C., Adelberger, K. L., Shapley, A. E., Erb, D. K., Reddy, N. A., & Pettini, M. 2005, ApJ, 626, 44
  * Stevens et al. (2003) Stevens, J. A., et al. 2003, Nature, 425, 264
  * Tremonti et al. (2004) Tremonti, C. A., et al. 2004, ApJ, 613, 898
  * van Dokkum & van der Marel (2007) van Dokkum, P. G. & van der Marel, R. P. 2007, ApJ, 655, 30
  * Venemans et al. (2007) Venemans, B. P., et al. 2007, A&A, 461, 823
  * Zirm et al. (2007) Zirm, A. W., et al. 2007, ApJ, 656, 66

Table 1Properties of Red Galaxies in the MRC 1138-262 Field: $H_{160}<24.5\wedge 1.1\leq(J_{110}-H_{160})\leq 2.1$ ID |  $z_{\rm phot}$aaThe errors on these photometric redshifts are approximately 0.1-0.2 | OddsbbThe ODDS parameter approaches unity when the probability distribution function has a single narrow peak | SED TypeccBest-fit broad-band template type: 0: dusty and star-forming 1: quiescent | Sérsic | $r_{e}$ | $r_{e}$ |  $H_{160}$ddTotal magnitudes based on the SExtractor MAG_AUTO parameter | $U_{n}-g_{475}$ | $g_{475}-I_{814}$ | $I_{814}-J_{110}$ | $J_{110}-H_{160}$
---|---|---|---|---|---|---|---|---|---|---|---
|  | (BPZ) |  | Index | ($\arcsec$) | (kpc) | (AB) | (AB) | (AB) | (AB) | (AB)
|  |  |  | ($n$) |  |  |  |  |  |  |
606 | 2.04 | 0.98 | 1.00 | 3.4 $\pm$ 0.4 | 0.21 $\pm$ 0.01 | 1.79 $\pm$ 0.10 | 23.09 $\pm$ 0.03 | 1.10 | 0.79 | 1.18 | 1.06
507 | 2.09 | 1.00 | 1.00 | 5.0 $\pm$ 0.4 | 0.24 $\pm$ 0.01 | 1.99 $\pm$ 0.08 | 22.50 $\pm$ 0.01 | 28.88 | -27.16 | 2.22 | 1.52
314 | 2.10 | 1.00 | 1.00 | 5.0 $\pm$ 0.7 | 0.15 $\pm$ 0.01 | 1.28 $\pm$ 0.05 | 23.03 $\pm$ 0.01 | 1.92 | 1.38 | 1.89 | 1.29
586 | 2.10 | 1.00 | 1.00 | 5.0 $\pm$ 0.3 | 0.40 $\pm$ 0.03 | 3.34 $\pm$ 0.23 | 21.77 $\pm$ 0.01 | 1.80 | 1.89 | 2.22 | 1.61
547 | 2.23 | 1.00 | 1.00 | 5.0 $\pm$ 0.6 | 0.12 $\pm$ 0.00 | 1.02 $\pm$ 0.04 | 22.70 $\pm$ 0.01 | 1.30 | 1.13 | 1.67 | 1.36
493 | 2.25 | 1.00 | 1.00 | 3.5 $\pm$ 0.3 | 0.56 $\pm$ 0.04 | 4.72 $\pm$ 0.36 | 22.75 $\pm$ 0.02 | 0.63 | 1.28 | 2.19 | 1.64
312 | 2.42 | 1.00 | 1.00 | 5.0 $\pm$ 0.3 | 0.38 $\pm$ 0.02 | 3.20 $\pm$ 0.19 | 22.00 $\pm$ 0.01 | -28.02 | 1.45 | 1.94 | 1.55
127 | 2.44 | 1.00 | 1.00 | 2.7 $\pm$ 0.6 | 0.12 $\pm$ 0.02 | 1.01 $\pm$ 0.18 | 22.64 $\pm$ 0.01 | 0.34 | 1.04 | 1.37 | 1.40