Bayesian constraints on bosonic dark matter in neutron stars within a covariant density-functional equation-of-state ensemble

We investigate neutron stars admixed with repulsively self-interacting bosonic dark matter (DM) using current multi-messenger observations. The dark component is modeled as a complex scalar field coupled to baryonic matter only through gravity, allowing dark-core and dark-halo configurations. To quantify baryonic equation-of-state (EOS) uncertainties, we combine a three-dimensional scan over $(m_χ,λ,F_χ)$ with 27 DDME2-based covariant density-functional EOSs spanning $L_{\rm sym}=30,50,70$ MeV and $-600\le Q_{\rm sat}\le 1000$ MeV. The models are confronted in a Bayesian framework with NICER mass--radius measurements from the Amsterdam and Maryland/Illinois analyses and the GW170817 tidal-deformability constraint. After Bayesian evidence weighting over the EOS ensemble, the EOS-relevant dark-sector scale $μ_χ\equiv m_χ/λ^{1/4}$ is localized around 250 MeV for both NICER data sets. In the $(m_χ,λ)$ parameterization, the marginalized boson-mass posterior peaks at $m_χ\simeq300$--$330$ MeV, while $λ$ remains broadly distributed. Since the bosonic EOS depends on $m_χ$ and $λ$ only through $μ_χ$, their separate marginals are conditional on the adopted independent priors. The inferred DM fraction is more data-set dependent, with posterior modes near 4% for Maryland/Illinois and 12% for Amsterdam. Bosonic DM also shifts the inferred $Q_{\rm sat}$ from negative toward intermediate and positive values, whereas the response of $L_{\rm sym}$ is weaker. This indicates a model-dependent trade-off between DM-induced compactification and the high-density stiffness of baryonic matter. These constraints remain conditional on the selected EOS family and the assumption of a common effective $F_χ$ for all sources.