Statistical mechanics of the $N$-queens problem

We investigate the $N$-queens problem as a lattice gas -- a model in which $N$ queens are placed on an $N \times N$ chessboard with pairwise repulsive interactions along shared rows, columns, and diagonals -- from the perspective of statistical mechanics. The ground states are exactly the $Q(N)$ solutions of the classical $N$-queens problem, with entropy per queen $s_0 \approx \ln N - \gamma$ ($\gamma \approx 1.944$). This entropy reflects a characteristic constraint hierarchy: each successive geometric constraint -- columns, then diagonals -- reduces the entropy from the free-placement value $\ln N$ by a definite constant. We derive the exact high-temperature energy $E/N \to 5/3$ as $N \to \infty$. Extensive Monte Carlo simulations with $10^8$ sweeps per temperature point for $N = 8$--$1024$ reveal that the specific heat per queen $C_v/N$ converges to a universal function of $T$ as $N \to \infty$. The converged curve features a non-divergent peak $C_v^{\max}/N \approx 1.63$ at $T^* \approx 0.235\,J$, establishing the absence of a thermodynamic phase transition. Combined with the trivially exact high-temperature entropy $S(\infty)/N = (1/N) \ln \binom{N^2}{N}$, the convergence of $C_v/N$ enables a thermodynamic integration of $C_v/T$ from $T = \infty$ to $T = 0$ that recovers the ground-state entropy -- and hence the Simkin constant $\gamma$ -- purely from Monte Carlo data. This provides an independent thermodynamic route to a fundamental combinatorial constant. Thermodynamic integration yields $\gamma_{\rm MC} = 1.946 \pm 0.003$ at $N = 1024$, within $0.1\%$ of the precise combinatorial value $\gamma = 1.94400(1)$. We further present a transfer-matrix-based tensor network formulation that encodes the non-attacking constraints into a rank-9 site tensor with 17 nonzero elements, providing a complementary exact-enumeration route.