  We present a Bayesian analysis of renormalization-group invariant
nucleon-nucleon interactions at leading order in chiral effective field theory
($\chi$EFT) with momentum cutoffs in the range 400--4000 MeV. We use history
matching to identify relevant regions in the parameter space of low-energy
constants (LECs) and subsequently infer the posterior probability density of
their values using Markov chain Monte Carlo. All posteriors are conditioned on
experimental data for neutron-proton scattering observables and we estimate the
$\chi$EFT truncation error in an uncorrelated limit. We do not detect any
significant cutoff dependence in the posterior predictive distributions for
two-nucleon observables. For all cutoff values we find a multimodal LEC
posterior with an insignificant mode harboring a bound $^1{S}_0$ state. The
$^3P_0$ and $^3P_2$ phase shifts emerging from the Bayesian analysis are less
constrained and typically more repulsive compared to the results of a phase
shift optimization. We expect that our inference will impact predictions for
nuclei. This work demonstrates how to perform inference in the presence of
limit-cycle-like behavior and spurious bound states, and lays the foundation
for a Bayesian analysis of renormalization-group invariant $\chi$EFT
interactions beyond leading order.
