Is the fully connected Potts model exactly solvable? Suppose that we have "spins" $\sigma_1,\dots,\sigma_N$, with $\sigma_i\in\{1,\dots,q\}$, for $i=1,\dots,N$, and that our Hamiltonian is
$$
  H = -\frac{J}{N} \sum_\stackrel{i,j=1}{i\ne j}^N \delta(\sigma_i,\sigma_j) - h \sum_{i=1}^N \delta(\sigma_i,1),
$$
in which $\delta$ is a Kronecker delta: $\delta(\sigma_i,\sigma_j)=1$ if $\sigma_i=\sigma_j$, and $\delta(\sigma_i,\sigma_j)=0$ if $\sigma_i\ne \sigma_j$.
This is a fully connected Potts model, defined on the complete graph, where each spin interacts with every other spin in the system.
Define the model partition function as
$$
  Z = \sum_{\sigma_1=1}^q \dots \sum_{\sigma_N=1}^q \exp(-\beta H),
$$
in which $\beta=1/(k_B T)$.
Direct computation of $Z$ is not feasible for any reasonable $N$, since it involves the sum of $q^N$ terms.
Is there in the literature a known way to perform the above summations and find a "manageable" expression $f(J,h,\beta,N$) for the partition function $Z$?
I'm looking for an exact result holding for every $N\geq 2$, and finite $J$, $h$ and $\beta$, which can be computed in polynomial time.
 A: You should not expect close form expressions for the finite-$N$ partition functions.
In fact, this is already the case when $q=2$. The latter is equivalent to the Curie-Weiss model,
in which the spins $\sigma_1,\dots,\sigma_N$ take values in $\{-1,1\}$ and the Hamiltonian takes the form
$$
H = -\frac{J}{N}\sum_{i,j=1}^N \sigma_i\sigma_j - h \sum_{i=1}^N \sigma_i.
$$
Note that I don't impose that $i\neq j$ in the first sum. This only shifts the energy by $J$ and thus plays no role, while slightly simplifying the exposition.
Let me describe the best you can hope for in this model. Introducing the magnetization $M=\sum_{i=1}^N \sigma_i$, the Hamiltonian can be reexpressed as
$$
H = -\frac{J}{N} \Bigl( \sum_{i=1}^N \sigma_i \Bigr)^2 - h M = -\frac{J}{N} M^2 - h M.
$$
From this observation, one can proceed in two different ways, both providing an expression for the partition function.
The first way is combinatorial. Just observe that
$$
Z_N = \sum_{k=0}^N \binom{N}{k} \exp\Bigl( \frac{\beta J}{N} (2k-N)^2 + \beta h (2k-N) \Bigr),
\tag{1}
$$
where the sum is over the number $k$ of spins $\sigma_i$ such that $\sigma_i=1$ (in particular, $M=k-(N-k)=2k-N$).
This is the first "explicit" expression for the partition function. It reduces the partition function from a sum over $2^N$ configurations to a sum over the $N+1$ possible values of the magnetization.
The second approach is via the Hubbard–Stratonovich transformation, which implies that
$$
\exp\bigl( \frac{\beta J}{N} M^2 \bigr) = \sqrt{\frac{N}{\pi\beta J}} \int_{-\infty}^{+\infty} \exp \bigl( - \frac{N}{\beta J} x^2 + 2 M x \bigr) \, \mathrm{d}x.
$$
From this, we can write
\begin{align}
Z_N
&= \sum_{\sigma_1,\dots,\sigma_N} \sqrt{\frac{N}{\pi\beta J}} \int_{-\infty}^{+\infty} \exp \bigl( - \frac{N}{\beta J} x^2 + (2 x + \beta h) M \bigr) \, \mathrm{d}x \\
&= \sqrt{\frac{N}{\pi\beta J}} \int_{-\infty}^{+\infty} \exp \bigl( - \frac{N}{\beta J} x^2 \bigr) \prod_{i=1}^N \underbrace{\sum_{\sigma_i=\pm 1} \exp \bigl( (2 x + \beta h) \sigma_i \bigr)}_{=2\cosh(2x + \beta h)} \, \mathrm{d}x \\
&= \sqrt{\frac{N}{\pi\beta J}} \int_{-\infty}^{+\infty} \exp \bigl( - N \varphi(x) \bigr)  \, \mathrm{d}x ,
\tag{2}
\end{align}
with $\varphi(x) = \frac{1}{\beta J} x^2 - \log\cosh(2x + \beta h) - \log 2$. This is the second "explicit" expression for the partition function. It reduces the partition function from a sum over $2^N$ configurations to an integral.
It does not look likely that one can explicitly evaluate the sum in (1) or the integral in (2). So I believe that this is the best you can hope for. Of course, both (1) and (2) can be used to extract a lot of information on the model. For instance, a saddle-point analysis of (2) would yield sharp approximations of the partition function for finite (but large) values of $N$.
For larger values of $q$, one can obtain expressions similar to those in (1) and (2) (instead of considering the magnetization $M$, one should consider the vector $(N_1,\dots,N_q)$ where $N_k$ is the number of spins taking value $k$).
A: What was proposed in (1) of a previous answer can be generalized to compute $Z$ exactly in polynomial time for $q\ge2$.
The key is that, as in the $q=2$ case, the Hamiltonian only depends on the number of spins taking each of the values.
Define $n_\sigma$ to be the number of spins taking value $\sigma\in\{1,2,\ldots,q\}$. Formally, $n_\sigma=\sum_{i=1}^N\delta(\sigma_i, \sigma)$.
Then observe that
\begin{align}
\sum_{i,j}\delta(\sigma_i,\sigma_j) &= \sum_i \sum_{j}\delta(\sigma_i,\sigma_j)\\
&= \sum_{\sigma=1}^q \sum_{i}\delta(\sigma_i,\sigma)\sum_{j}\delta(\sigma,\sigma_j)\\
&= \sum_{\sigma=1}^q n_\sigma^2
\end{align}
where to get the second equality one groups together spins with the same value $\sigma$.
Therefore the Hamiltonian reads
$$H = -\frac JN\sum_{\sigma=1}^q n_\sigma^2-hn_1$$
Now the partition function be computed by grouping together all configurations sharing the same $\{n_1,n_2,\ldots,n_q\}$ and counting how many there are
\begin{align}
Z_N &= \sum_{\sigma_1,\sigma_2,\ldots,\sigma_N} \exp\left(-\beta H\right)\\
&= \sum_{n_1,n_2\ldots,n_q}  \exp\left(-\beta H\right)\sum_{\sigma_1,\sigma_2,\ldots,\sigma_N} \mathbf{1}\left[n_\sigma=\sum_{i=1}^N\delta(\sigma_i, \sigma)\;\forall\sigma\right]\\
&= \sum_{n_1,n_2\ldots,n_q}\binom{N}{n_1 n_2 \cdots n_q} \exp\left(\frac{\beta J}{N}\sum_{\sigma=1}^q n_\sigma^2+\beta hn_1\right)
\end{align}
which takes $\mathcal O (N^q)$: you need $q$ nested loops, one for each of the $n_\sigma$'s.
It is implicit that the sum only runs over allowed sets $\{n_1,n_2,\ldots,n_q \:|\: n_\sigma \in \{0,1,\ldots,N\} \forall \sigma, \sum_{\sigma=1}^q n_\sigma=N\}$.
