How Many Components? Quantifying the Complexity of the Metallicity Distribution in the Milky Way Bulge with APOGEE

We use data of ~13.000 stars from the SDSS/APOGEE survey to study the shape of the MDF within the region $|l|\leq11^\circ$ and $|b|\leq13^\circ$, and spatially constrained to ${\rm R_{GC}\leq3.5}$~kpc. We apply Gaussian Mixture Modeling (GMM) and Non-negative Matrix Factorization (NMF) decomposition techniques to identify the optimal number and the properties of MDF components in different spatial locations and under different sampling conditions. We find the shape and spatial variations of the MDF (at ${\rm [Fe/H]\geq-1}$~dex) are well represented as a smoothly varying contribution of three overlapping components located at [Fe/H]~=+$0.32$, $-0.17$ and $-0.66$~dex. The bimodal MDF found in previous studies is in agreement with our trimodal assessment once the limitations in sample size and individual measurement errors are taken into account. The shape of the MDF and its correlations with kinematics reveal different spatial distributions and kinematical profiles for the three chemical components co-existing in the bulge region. We confirm the consensus physical interpretation of metal-rich stars as associated with the boxy/peanut X-shape bar, originating from the secular evolution of the early disc. On the other hand, metal-intermediate stars could be the product of in-situ formation at high redshift, in a gas-rich environment characterized by violent and fast star formation. This interpretation would help to link a present-day structure with those observed in formation in the center of high redshift galaxies. We refrain from associating the metal-poor stars with any particular formation mechanism. They seem to be inconsistent with being thick disc or halo stars, but may be the metal-rich tail of the population currently being characterized at lower metallicity from the study of RR Lyrae stars. Conversely, they could be associated with the metal-poor tail of the early thick disc.