A bstract A search for Higgs boson pair production in events with two b -jets and two τ -leptons is presented, using a proton–proton collision dataset with an integrated luminosity of 139 fb − 1 collected at $$ \sqrt{s} $$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:msqrt> <mml:mi>s</mml:mi> </mml:msqrt> </mml:math> = 13 TeV by the ATLAS experiment at the LHC. Higgs boson pairs produced non-resonantly or in the decay of a narrow scalar resonance in the mass range from 251 to 1600 GeV are targeted. Events in which at least one τ -lepton decays hadronically are considered, and multivariate discriminants are used to reject the backgrounds. No significant excess of events above the expected background is observed in the non-resonant search. The largest excess in the resonant search is observed at a resonance mass of 1 TeV, with a local (global) significance of 3 . 1 σ (2 . 0 σ ). Observed (expected) 95% confidence-level upper limits are set on the non-resonant Higgs boson pair-production cross-section at 4.7 (3.9) times the Standard Model prediction, assuming Standard Model kinematics, and on the resonant Higgs boson pair-production cross-section at between 21 and 900 fb (12 and 840 fb), depending on the mass of the narrow scalar resonance.