  A search for decays of the Higgs and Z boson to pairs of J/$\psi$ or
$\Upsilon$(nS) (n = 1, 2, 3) mesons, with their subsequent decay to
$\mu^+\mu^-$ pairs, is presented. The analysis uses data from proton-proton
collisions at $\sqrt{s} =$ 13 TeV, collected with the CMS detector at the LHC
in 2017 and corresponding to an integrated luminosity of 37.5 fb$^{-1}$. While
an observation of such a decay with this sample would indicate the presence of
physics beyond the standard model, no significant excess is observed. Upper
limits at 95% confidence level are placed on the branching fractions of these
decays. In the J/$\psi$ pair channel, the limits are 1.8 $\times$ 10$^{-3}$ and
2.2 $\times$ 10$^{-6}$ for the Higgs and Z boson, respectively, while in the
combined $\Upsilon$(nS) pair channel, the limits are 1.4 $\times$ 10$^{-3}$ and
1.5 $\times$ 10$^{-6}$, respectively, when the mesons from the Higgs and Z
boson decay are assumed to be unpolarized. When fully longitudinal and
transverse polarizations are considered the limits reduce by about 22-29% and
increase by about 10-13%, respectively.
