  Chiral effective field theory (chiEFT) provides a systematic approach to
describe low-energy nuclear forces. Moreover, chiEFT is able to provide
well-founded estimates of statistical and systematic uncertainties -- although
this unique advantage has not yet been fully exploited. We fill this gap by
performing an optimization and statistical analysis of all the low-energy
constants (LECs) up to next-to-next-to-leading order. Our optimization protocol
corresponds to a simultaneous fit to scattering and bound-state observables in
the pion-nucleon, nucleon-nucleon, and few-nucleon sectors, thereby utilizing
the full model capabilities of chiEFT. We study the effect on other observables
by demonstrating error-propagation methods that can easily be adopted by future
works. We employ mathematical optimization and implement automatic
differentiation to attain efficient and machine-precise first- and second-order
derivatives of the objective function with respect to the LECs. We use
power-counting arguments to estimate the systematic uncertainty that is
inherent to chiEFT and we construct chiral interactions at different orders
with quantified uncertainties. Statistical error propagation is compared with
Monte Carlo sampling showing that statistical errors are in general small
compared to systematic ones. In conclusion, we find that a simultaneous fit to
different sets of data is critical to (i) identify the optimal set of LECs,
(ii) capture all relevant correlations, (iii) reduce the statistical
uncertainty, and (iv) attain order-by-order convergence in chiEFT. Furthermore,
certain systematic uncertainties in the few-nucleon sector are shown to get
substantially magnified in the many-body sector; in particlar when varying the
cutoff in the chiral potentials. The methodology and results presented in this
Paper open a new frontier for uncertainty quantification in ab initio nuclear
theory.
