One of the major challenges in cosmology is to understandthe matter content of our universe. According to General Relativity (GR) there are deep connections between what fills the universe and its geometry, its dynamical expansion and eventually its fate. Notably, visible ordinary matter appears to be only a small fraction of the matter in our universe, whereas the bulk comes in the form of non-luminous and non-baryonic particles, dubbed dark matter (DM). Complementary measurements of large scale structures, galaxy formation, gravitational lensing and of the cosmic microwave background (CMB) strongly suggest that more than 80\% of the matter in the universe consists of DM.The most accurate determination for ordinary matter and DM energy density are provided by anisotropies in the CMB and amount to $\Omega_{\hbox{\scriptsize B}} h^2 = 0.02237 \pm 0.00015$ and $\Omega_{\hbox{\scriptsize DM}} h^2 = 0.1200 \pm 0.0012$.On the contrary, there is almost a total lack of information concerning DM fromthe particle physics point of view. Over the last few decades, a Weakly Interacting Massive Particle (WIMP) has become the most studied candidate for DM. One of its key features is the production mechanism in the early universe via the so-called \textit{freeze-out} of a thermal relic. For a DM particle of mass $M$, the freeze-out occurs at a temperature $T \sim M/25$, i.e. the DM particles are \textit{non-relativistic}. However, the freeze-out mechanism is not limited to WIMPs and even applies to cases where interactions are stronger.The key quantity that governs the evolution of the DM abundance is the thermally averaged annihilation cross section $\langle \sigma v \rangle_{{\rm{ann}}}$, encompassed in a Boltzmann equation\begin{equation}(\partial_t + 3H) n = -\langle \sigma v \rangle_{{\rm{ann}}} (n^2-n^2_{{\rm{eq}}}) \, ,\label{Boltzmann_1}\end{equation} where $H$ is the Hubble rate and $n$ is the total number density of DM particles ($n_{{\rm{eq}}}$ is that in equilibrium). The value of $\langle \sigma v \rangle_{{\rm{ann}}}$ depends on the model under consideration, namely its particle content and interactions. It is extremely important to calculate $\langle \sigma v \rangle_{{\rm{ann}}}$ accurately because the predicted present-day DM energy density depends crucially on it through the solution of eq.~(\ref{Boltzmann_1}). The DM mass is in turn fixed as function of the other model parameters to reproduce $\Omega_{\hbox{\scriptsize DM}} h^2$ in the first place (e.g.~couplings, mass splitting with other dark species). However, determining $\langle \sigma v \rangle_{{\rm{ann}}}$ by including the full features of each model and the thermal environment is not an easy task. In a variety of theories, DM interacts with gauge bosons or scalars that induce long-range interactions because of repeated soft exchanges. Remarkably, the inclusion of bound-state effects for DM annihilation has been recently shown to have a large impact on the \textit{overclosure bound} for DM models, namely the largest value of the particle mass compatible with the observed DM energy density. At the same time, it is manifestly subtle and complicated to include bound-state dynamics in a thermal medium due to the intricate interplay between non-relativistic and thermal energy scales. The overall objective of the project is to develop and apply a novel approach for relic density computations in the framework of non-relativistic effective field theories (NREFTs) and open quantum systems (OQSs). Starting from a thermal field theoretic formulation of the problem, the interactions between the DM and the light plasma constituents will systematically be included: \textit{bound-state formation and dissociation}, Sommerfeld effect, DM thermal masses and interaction rates, and the role of \textit{phase transitions} during the freeze-out process. I will provide the DM annihilation cross sections and overclosure bounds for well-motivated classes of models: DM coannihilating with strongly interacting states (QCD-like), DM embedded in spontaneously broken Abelian and non-Abelian gauge gauge theories, composite DM. To this aim cutting-edge effective field theory (EFT) techniques at finite temperature will be used, in particular NREFTs and potential NREFTs (pNREFTs). Moreover, DM particles in medium will be recast in an open-quantum-system framework for the first time, with a focus on DM bound-state formation and their in-medium evolution. By adopting and exploiting such methods, the corresponding overclosure bounds on the DM mass will serve as comprehensive and solid benchmarks to interpret and guide experimental analyses for a wide class of DM models. In the contemporary endeavour of testing and possibly falsifying WIMP-like models with present and upcoming experiments, the outcomes of this project are timely and necessary in order to deliver the most reliable relic density constraints on the models parameters, especially on the DM mass.