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 Lsym=30,50,70 MeV and -600≤ Qₛₐₜ≤ 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 μ_χ≡ 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_χ300--$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ₛₐₜ from negative toward intermediate and positive values, whereas the response of Lsym 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.