WIMP Freeze-out dynamics under Tsallis statistics

We generalize thermal WIMP (Weakly Interacting Massive Particle) freeze-out within Tsallis nonextensive statistics. Using Curado-Tsallis $q$-distributions $f_q(E;\mu,T)$ we compute $q$-deformed number and energy densities, pressure, entropy density and Hubble rate, $\{n_q,\rho_q,P_q,s_q,H_q\}$. The Boltzmann equation is generalized accordingly to obtain the comoving abundance $Y_{\chi,q}(x)$ and relic density $\Omega_{\chi,q}h^2$ for a dark-matter candidate $\chi$ in a model-independent setup. The thermally averaged cross section is expanded as $\langle\sigma v\rangle_q \approx a + b\,\langle v_{\rm rel}^2\rangle_q$ up to $p$-wave. The freeze-out parameter $x_f(q)$ is determined from $\Gamma_{{\rm ann},q}(T_f)\simeq H_q(T_f)$ using a $q$-logarithmic inversion, with the expansion rate modified through ultra-relativistic rescalings $R_\rho(q)$ of the effective relativistic degrees of freedom $g_*$ and $g_{*s}$. We show that $x_f$ increases with $q$ and that QCD-threshold features propagate into $Y_{\chi,q}(x)$ and $\Omega_{\chi,q}h^2$. We then perform two $q$-grid scans: fixing $\langle\sigma v\rangle_q$ while varying the dark-matter mass $m_\chi$, and fixing $m_\chi$ while varying the $s$-wave coefficient $a$. For an $s$-wave dominated scenario we construct $\chi^2$ profiles in these planes by comparing $\Omega_{\chi,q}h^2$ with the Planck benchmark $\Omega_c h^2 = 0.120\pm 0.001$. In both cases we find a clear degeneracy in the preferred nonextensive parameter $q_{\rm best}$ along valleys in parameter space. However, fixed-mass scans (varying $\langle\sigma v\rangle_q$) are significantly more constraining than fixed-cross-section scans, reflecting that $\Omega_{\chi,q}h^2$ is mainly controlled by $\langle\sigma v\rangle_q$, so that for realistic cross sections the best-fit $q_{\rm best}$ remains close to the extensive limit $q\to 1$.