  Coexistence of different geometric shapes at low energies presents a
universal structure phenomenon that occurs over the entire chart of nuclides.
Studies of the shape coexistence are important for understanding the
microscopic origin of collectivity and modifications of shell structure in
exotic nuclei far from stability. The aim of this work is to provide a
systematic analysis of characteristic signatures of coexisting nuclear shapes
in different mass regions, using a global self-consistent theoretical method
based on universal energy density functionals and the quadrupole collective
model. The low-energy excitation spectrum and quadrupole shape invariants of
the two lowest $0^{+}$ states of even-even nuclei are obtained as solutions of
a five-dimensional collective Hamiltonian (5DCH) model, with parameters
determined by constrained self-consistent mean-field calculations based on the
relativistic energy density functional PC-PK1, and a finite-range pairing
interaction. The theoretical excitation energies of the states: $2^+_1$,
$4^+_1$, $0^+_2$, $2^+_2$, $2^+_3$, as well as the $B(E2; 0^+_1\to 2^+_1)$
values, are in very good agreement with the corresponding experimental values
for 621 even-even nuclei. Quadrupole shape invariants have been implemented to
investigate shape coexistence, and the distribution of possible
shape-coexisting nuclei is consistent with results obtained in recent
theoretical studies and available data. The present analysis has shown that,
when based on a universal and consistent microscopic framework of nuclear
density functionals, shape invariants provide distinct indicators and reliable
predictions for the occurrence of low-energy coexisting shapes. This method is
particularly useful for studies of shape coexistence in regions far from
stability where few data are available.
