  In the upcoming years, present and next-generation gravitational wave
observatories will detect a larger number of binary neutron star (BNS) mergers
with increasing accuracy. In this context, improving BNS merger numerical
simulations is crucial to correctly interpret the data and constrain the
equation of state (EOS) of neutron stars (NSs). State-of-the-art simulations of
BNS mergers do not include muons. However, muons are known to be relevant in
the microphysics of cold NSs and are expected to have a significant role in
mergers, where the typical thermodynamic conditions favour their production.
Our work is aimed at investigating the impact of muons on the merger remnant.
We post-process the outcome of four numerical relativity simulations of BNS
mergers performed with three different baryonic EOSs and two mass ratios
considering the first $15$ milliseconds after merger. We compute the abundance
of muons in the remnant and analyse how muons affect the trapped neutrino
component and the fluid pressure. We find that depending on the baryonic EOS,
the net fraction of muons is between $30 \%$ and $70 \%$ the net fraction of
electrons. Muons change the flavour hierarchy of trapped (anti-)neutrinos such
that deep inside the remnant, muon anti-neutrinos are the most abundant,
followed by electron anti-neutrinos. Finally, muons and trapped neutrinos
modify the neutron-to-proton ratio, affecting the remnant pressure by up to
$7\%$ when compared with calculations neglecting them. This work demonstrates
that muons have a non-negligible effect on the outcome of BNS merger
simulations, and they should be included to improve the accuracy of a
simulation.
