We have one recurrence for $\widetilde Z(N, \beta, \lambda)$ but it involves a complicated combinatorial sum and a product involving $p$ different terms $\widetilde Z(n_j, \beta, \lambda)$. Here we demonstrate a simpler quadratic recurrence.
Theorem: $\widetilde Z(N, \beta, \lambda)$ satisfies the quadratic recurrence $$\sum_{n=1}^N \left(\frac{N}{p+1} – n \right) p^{-\beta {n \choose 2} – n} \widetilde Z(n, \beta, \lambda) \widetilde Z(N-n, \beta, \lambda)=0, $$ with initial conditions $\widetilde Z(0, \beta, \lambda) = \widetilde Z(1, \beta, \lambda) = 1.$
Proof. We start by doing something surprising. We are going to “insert” another coset of $p \mathbb Z_p$ into $\mathbb Z_p$, and call this new object $\mathcal T$ and we call the extra coset we inserted $\mathcal P$. We label the other cosets as usual. $\mathcal T$ is not a ring, or any type of algebraic object, however it is a well defined measure space if we extend Haar measure in the sensible way. The measure is no longer a probability measure, so we renormalize to produce probability measure $\lambda_{\mathcal T}$. Nothing really changes by adding an extra coset, except we have to keep track of how many particles fall into it and modify the product in the recurrence $Z(N, \beta, \lambda)$ by the appropriate term.
Let us compute the expectation of $N_{\mathcal P}$ with respect to $\lambda_{\mathcal T}$. We do this by summing over the values $n = 0, 1, \ldots, N$ and note that the event $\{N_{\mathcal P} = n\}$ consists of $N$-tuples $\mathbf x$ with $n$ coordinates in $\mathcal P$ and $N – n$ coordinates in $\mathbb Z_p$. It follows, after accounting for the combinatorics of which particles go where, and the renormalization term for $\lambda_{\mathcal T}$, $$\mathbb E[N_{\mathcal P}] = \frac{1}{\widetilde Z(N, \beta, \lambda_{\mathcal T})} \left(\frac{p+1}{p}\right)^N \sum_{n=1}^N n \widetilde Z(n, \beta, \lambda_1) \widetilde Z(N-n, \beta, \lambda).$$ Notice we can expand $\widetilde Z(N, \beta, \lambda_{\mathcal T})$ in exactly the same way,$$\widetilde Z(N, \beta, \lambda_{\mathcal T}) = \left(\frac{p+1}{p}\right)^N \sum_{n=1}^N \widetilde Z(n, \beta, \lambda_1) \widetilde Z(N-n, \beta, \lambda).$$ Consequently, $$\mathbb E[N_{\mathcal P}] \sum_{n=1}^N \widetilde Z(n, \beta, \lambda_1) \widetilde Z(N-n, \beta, \lambda) = \sum_{n=1}^N n \widetilde Z(n, \beta, \lambda_1) \widetilde Z(N-n, \beta, \lambda).$$ The only thing left is to notice that, if what we really wanted was to compute $E[N_{\mathcal P}]$, then we could have just used symmetry, and the fact that there are $p+1$ equal cosets to conclude $E[N_{\mathcal P}] = \frac{N}{p+1}.$ The theorem follows by rearranging terms and using the relationship between $\widetilde Z(n, \beta, \lambda_1)$ and $\widetilde Z(n, \beta, \lambda_1)$ to find $$\sum_{n=1}^N \left(\frac{N}{p+1} – n \right) p^{-\beta {n \choose 2} – n} \widetilde Z(n, \beta, \lambda) \widetilde Z(N-n, \beta, \lambda)=0, $$
Corollary: $Z(N, \beta, \lambda)$ is a rational function in $p^{-\beta}$ with rational coefficients.
Examples
The theorem provides an efficient algorithm to compute examples of $\widetilde Z(N, \beta, \lambda)$. $$\begin{eqnarray} \widetilde Z(0, \beta, \lambda ) &=& 1 \\ \widetilde Z(1, \beta, \lambda) &=& 1 \\ \widetilde Z(2, \beta, \lambda) &=& \frac{(p-1) p^{\beta }}{p^{\beta +1}-1} \\ \widetilde Z(3,\beta, \lambda) &=& \frac{(p-1) p^{3 \beta } \left(-2 p^{\beta +1}+p^{\beta +2}+2 p-1\right)}{\left(p^{\beta +1}-1\right) \left(p^{3 \beta +2}-1\right)} \\ \widetilde Z(4, \beta, \lambda) &=& \frac{-\frac{(p-1)^2 (4-2 (p+1)) p^{\beta -2}}{4 \left(p^{\beta +1}-1\right)^2}-\frac{(3-p) (p-1) \left(-2 p^{\beta +1}+p^{\beta +2}+2 p-1\right) p^{3 \beta -1}}{6 \left(p^{\beta +1}-1\right) \left(p^{3 \beta +2}-1\right)}-\frac{(p-1) (4-3 (p+1)) \left(-2 p^{\beta +1}+p^{\beta +2}+2 p-1\right)}{6 p^3 \left(p^{\beta +1}-1\right) \left(p^{3 \beta +2}-1\right)}}{\frac{1}{24} (4-4 (p+1)) p^{-6 \beta -4}+\frac{1}{6}}. \end{eqnarray}$$
