A search for Higgs boson pair production in events with two $b$-jets and two $\tau$-leptons is presented, using a proton-proton collision dataset with an integrated luminosity of 139 fb$^{-1}$ collected at $\sqrt{s}=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 $\tau$-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\sigma$ ($2.0\sigma$). 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.
Post-fit distribution of mHH in the hadhad channel.