Binary neutron star mergers provide insights into strong-field gravity and the properties of ultra-dense nuclear matter. These events offer the potential to search for signatures of physics beyond the standard model, including dark matter. We present the first numerical-relativity simulations of binary neutron star mergers admixed with dark matter, based on constraint-solved initial data. Modeling dark matter as a non-interacting fermionic gas, we investigate the impact of varying dark matter fractions and particle masses on the merger dynamics, ejecta mass, post-merger remnant properties, and the emitted gravitational waves. Our simulations suggest that the dark matter morphology - a dense core or a diluted halo - may alter the merger outcome. Scenarios with a dark matter core tend to exhibit a higher probability of prompt collapse, while those with a dark matter halo develop a common envelope, embedding the whole binary. Furthermore, gravitational wave signals from mergers with dark matter halo configurations exhibit significant deviations from standard models when the tidal deformability is calculated in a two-fluid framework neglecting the dilute and extended nature of the halo. This highlights the need for refined models in calculating the tidal deformability when considering mergers with extended dark matter structures. These initial results provide a basis for further exploration of dark matter's role in binary neutron star mergers and their associated gravitational wave emission and can serve as a benchmark for future observations from advanced detectors and multi-messenger astrophysics.
