The fiducial cross section for Y(1S) pair production in proton-proton collisions at a center-of-mass energy of 13 TeV in the region where both Y(1S) mesons have an absolute rapidity below 2.0 is measured to be 79 ±11 (stat) ± 6(syst) ± 3(B) pb assuming the mesons are produced unpolarized. The last uncertainty corresponds to the uncertainty in the Y(1S) meson dimuon branching fraction. The measurement is performed in the final state with four muons using proton-proton collision data collected in 2016 by the CMS experiment at the LHC, corresponding to an integrated luminosity of 35.9fb^-1. This process serves as a standard model reference in a search for narrow resonances decaying to Y(1S)μ⁺μ^- in the same final state. Such a resonance could indicate the existence of a tetraquark that is a bound state of two b quarks and two b^- antiquarks. The tetraquark search is performed for masses in the vicinity of four times the bottom quark mass, between 17.5 and 19.5 GeV, while a generic search for other resonances is performed for masses between 16.5 and 27 GeV. No significant excess of events compatible with a narrow resonance is observed in the data. Limits on the production cross section times branching fraction to four muons via an intermediate Y(1S) resonance are set as a function of the resonance mass.