Mapping population synthesis event rates on model parameters II: Convergence and accuracy of multidimensional fits

Binary population synthesis calculations and associated predictions, especially event rates, are known to depend on a significant number of input model parameters with different degrees of sensitivity. At the same time, for systems with relatively low formation rates, such simulations are heavily computationally demanding and therefore the needed explorations of the high-dimensional parameter space require major -- often prohibitive -- computational resources. In the present study, to better understand several key event rates involving binary evolution and binaries with two compact objects in Milky Way-like galaxies and to provide ways of reducing the computational costs of complete parameter space explorations: (i) we perform a methodical parameter study of the \emph{StarTrack} population synthesis code ; and (ii) we develop a formalism and methodology for the derivation of {\em multi-dimensional fits} for event rates. We significantly generalize our earlier study, and we focus on ways of thoroughly assessing the accuracy of the fits. We anticipate that the efficient tools developed here can be applied in lieu of large-scale population calculations and will facilitate the exploration of the dependence of rate predictions on a wide range binary evolution parameters. Such explorations can then allow the derivation of constraints on these parameters, given empirical rate constraints and accounting for fitting errors. Here we describe in detail the principles and practice behind constructing these fits, estimating their accuracy, and comparing them with observations in a manner that accounts for their errors.