Massive stars reveal variations of the stellar initial mass function in
  the Milky Way stellar clusters
  We investigate whether the stellar initial mass function (IMF) is universal,
or whether it varies significantly among young stellar clusters in the Milky
Way. We propose a method to uncover the range of variation of the parameters
that describe the shape of the IMF for the population of young Galactic
clusters. These parameters are the slopes in the low and high stellar mass
regimes, $\gamma$ and $\Gamma$, respectively, and the characteristic mass,
$M_{ch}$. The method relies exclusively on the high mass content of the
clusters, but is able to yield information on the distributions of parameters
that describe the IMF over the entire stellar mass range. This is achieved by
comparing the fractions of single and lonely massive O stars in a recent
catalog of the Milky Way clusters with a library of simulated clusters built
with various distribution functions of the IMF parameters. The synthetic
clusters are corrected for the effects of the binary population, stellar
evolution, sample incompleteness, and ejected O stars. Our findings indicate
that broad distributions of the IMF parameters are required in order to
reproduce the fractions of single and lonely O stars in Galactic clusters. They
also do not lend support to the existence of a cluster mass-maximum stellar
mass relation. We propose a probabilistic formulation of the IMF whereby the
parameters of the IMF are described by Gaussian distribution functions centered
around $\gamma=0.91$, $\Gamma=1.37$, and $M_{ch}=0.41$ M$_{\odot}$, and with
dispersions of $\sigma_{\gamma}=0.25$, $\sigma_{\Gamma}=0.60$, and
$\sigma_{M_{ch}}=0.27$ M$_{\odot}$ around these values.
