In The canonical ensemble we have relaxed the constraint of having a fixed value of the energy, and thus defined the canonical ensemble. However, we still supposed that the number of particles of our system was fixed: we can remove also this constraint, and so consider systems that exchange not only energy but also particles.
The statistical ensemble of such systems is called *grand canonical ensemble*.

We therefore consider a system like the one represented below:

We again have a heat bath at temperature $T$ and a small subsystem that this time can also exchange particles with the bath; let us call $N_{\text{b}}$ the number of particles of the heat bath, $N_{\text{s}}$ the particles of the subsystem (which can both vary from instant to instant), $N=N_{\text{s}}+N_{\text{b}}$ the total number of particles of the system (which is instead fixed) and $V$ the total volume.
As we have done for the canonical ensemble, we rely on the tools of the microcanonical ensemble, namely we must determine the phase space volume of the system taking also into account that the number of particles is not fixed. This volume is:

$\Omega (E,N,V)=\sum _{N_{\text{s}}=0}^{N}\int d\Gamma _{\text{s}}d\Gamma _{\text{b}}\delta ({\mathcal {H}}_{\text{s}}+{\mathcal {H}}_{\text{b}}-E)$

where

$E$ is the total energy of the system, and the sum is taken exactly because the number of particles in the system is not fixed and in principle can run from zero to the largest possible value, i.e.

$N$; furthermore all the quantities with the subscript "s" are relative to the subsystem while those with "b" to the heat bath, and it is implied that the former are computed for

$N_{\text{s}}$ particles while the latter for

$N-N_{\text{s}}$.
Integrating over

$d\Gamma _{\text{b}}$ we can write

$\Omega$ as:

$\Omega (E,N,V)=\sum _{N_{\text{s}}=0}^{N}\int d\Gamma _{\text{s}}d\Gamma _{\text{b}}\delta [{\mathcal {H}}_{\text{b}}-(E-{\mathcal {H_{\text{s}}}})]=\sum _{N_{\text{s}}=0}^{N}\int d\Gamma _{\text{s}}\Omega _{\text{b}}(E-{\mathcal {H}}_{\text{s}},V_{\text{b}},N-N_{\text{s}})$

If we now divide both sides by

$\Omega (E,V,N)$:

$\sum _{N=0}^{N}\int d\Gamma _{\text{s}}{\frac {\Omega _{\text{b}}(E-{\mathcal {H}}_{\text{s}},V_{\text{b}},N-N_{\text{s}})}{\Omega (E,V,N)}}=1$

we recognise, in analogy to what we have seen in

The microcanonical ensemble, that now the phase space probability density is:

$\rho _{N_{\text{s}}}(\mathbb {Q} ,\mathbb {P} )={\frac {\Omega _{\text{b}}(E-{\mathcal {H}}_{\text{s}},V_{\text{b}},N-N_{\text{s}})}{\Omega (E,V,N)}}$

where the subscript

$N_{\text{s}}$ has been added in order to remember that this probability density depends on the number of particles that are in the subsystem, and the normalization condition requires that we must not only integrate

$\rho$ in

$(\mathbb {Q} ,\mathbb {P} )$, but also sum over all the possible numbers of particles.
Now, just like we have done for the canonical ensemble we can use the fundamental postulate of statistical mechanics to write:

$\Omega _{\text{b}}(E-{\mathcal {H}}_{\text{s}},V_{\text{b}},N-N_{\text{s}})=e^{{\frac {1}{k_{B}}}S_{\text{b}}(E-{\mathcal {H}}_{\text{s}},V_{\text{b}},N-N_{\text{s}})}$

and since the subsystem we are considering is again macroscopic but much smaller than the heat bath we have

${\mathcal {H}}_{\text{s}}\ll E$ and

$N_{\text{s}}\ll N$, so we can expand

$S_{\text{b}}$ in the exponential:

$\Omega _{\text{b}}(E-{\mathcal {H}}_{\text{s}},V_{\text{b}},N-N_{\text{s}})=e^{{\frac {1}{k_{B}}}\left(S_{\text{b}}(E,V_{\text{b}},N)-{\mathcal {H}}{\frac {\partial }{\partial E}}S_{\text{b}}(E,V_{\text{b}},N)-N_{\text{s}}{\frac {\partial }{\partial N}}S_{\text{b}}(E,V_{\text{b}},N)+\cdots \right)}$

From thermodynamics we know that:

${\frac {\partial S_{\text{b}}}{\partial E}}={\frac {1}{T}}\quad \qquad {\frac {\partial S_{\text{b}}}{\partial N}}=-{\frac {\mu }{T}}$

where

$\mu$ is the chemical potential, and in the end we get (removing the subscript "s" to have a more general notation):

$\rho _{N}(\mathbb {Q} ,\mathbb {P} )={\frac {e^{-\beta ({\mathcal {H}}(\mathbb {Q} ,\mathbb {P} )-\mu N)}}{\sum _{N}\int d\Gamma e^{-\beta ({\mathcal {H}}(\mathbb {Q} ,\mathbb {P} )-\mu N)}}}$

where:

${\mathcal {Z}}=\sum _{N}\int d\Gamma e^{-\beta ({\mathcal {H}}(\mathbb {Q} ,\mathbb {P} )-\mu N)}$

is called

*grand partition function*.

Just like we have done for the canonical ensemble, we could wonder what relations does this newly found ensemble have with the ones that we already know, in particular the canonical one; what we are now going to show is that they are equivalent when the system considered is macroscopic, since this time the (relative) fluctuations of the number of particles will be incredibly small.