A search for the production of Higgs boson pairs in proton–proton collisions at a centre-of-mass energy of 13 TeV is presented, using a data sample corresponding to an integrated luminosity of 35.9 fb−1 collected with the CMS detector at the LHC. Events with one Higgs boson decaying into two bottom quarks and the other decaying into two τ leptons are explored to investigate both resonant and nonresonant production mechanisms. The data are found to be consistent, within uncertainties, with the standard model background predictions. For resonant production, upper limits at the 95% confidence level are set on the production cross section for Higgs boson pairs as a function of the hypothesized resonance mass and are interpreted in the context of the minimal supersymmetric standard model. For nonresonant production, upper limits on the production cross section constrain the parameter space for anomalous Higgs boson couplings. The observed (expected) upper limit at 95% confidence level corresponds to about 30 (25) times the prediction of the standard model.