Grand canonical partition function of hypothetical particles

I have to calculate the grand canonical partition function of a system of hypothetical particles, wherein each single-particle quantum state can be occupied by up to 3 particles.

Obviously, this is a sort of joke, referring to fermions (with a maximum of 2 particles per state) and bosons (unlimited particles per state). It is assumed that these hypothetical particles do not interact with each other.

So I tried viewing each single-particle quantum state as a separate grand canonical ensemble, following the approach on https://en.wikipedia.org/wiki/Fermi%E2%80%93Dirac_statistics

At chemical potential $mu$ and temperature $T$, where the energy of the state is $epsilon$, I get:
begin{equation}
mathcal{Z} = sum_{n=0}^{3}{exp{left(frac{n(mu-epsilon)}{k_B T}right)} = frac{1-exp{left(4frac{mu-epsilon}{k_B T}right)}}{1-exp{left(frac{mu-epsilon}{k_B T}right)}}}
end{equation} where I used the finite geometric progression.

Now I also have to determine the average occupation number $langle n_i rangle$ for a state with energy $epsilon_i$ at temperature $T=0$.

In general, we have
begin{equation}
langle n_i rangle = k_B T frac{partial ln{mathcal{Z}}}{partial mu}
end{equation}

which yields me $langle n_i rangle =2-frac{1}{1+exp(x)}+tanh(x)$ where I defined $x=frac{mu-epsilon_i}{k_B T}$. (I used Wolfram Mathematica for simplifying the algebra.)

Clearly at $T=0$ this expression is ill-defined, but by taking the limit $Trightarrow 0$ we see that $langle n_irangle=0$ if $epsilon_i>mu$, $langle n_irangle=3/2$ if $epsilon_i=mu$ and $langle n_irangle=3$ if $epsilon_i<mu$, correct?