In this work, bubble polarization domains in periodically ordered ferroelectric Pb(Zr0.4Ti0.6)O3 nanodot arrays and their formation mechanisms have been investigated by piezoresponse force microscopy (PFM) and Monte-Carlo simulations. The PFM observations reveal the coexistence of single domain and apparent bubble domain patterns within the same nanodot array, which also exhibit dissimilar polarization reversal processes. The formation of various polarization configurations can be accounted for by the interplay of various factors, such as polarization anisotropy and depolarization field. Using Monte-Carlo simulation, we are able to reproduce bubble and single domains and further predict that these patterns can be tailored by varying the nanodot parameters, including dot height, aspect ratio, etc.