This paper reports a search for Higgs boson pair (hh) production in association with a vector boson (W or Z ) using 139 fb−1 of proton–proton collision data at √s = 13 TeV recorded with the ATLAS detector at the Large Hadron Collider. The search is performed in final states in which the vector boson decays leptonically (W → lν, Z → ll, νν with l = e, μ) and the Higgs bosons each decay into a pair of b-quarks. It targets V hh signals from both non-resonant hh production, present in the Standard Model (SM), and resonant hh production, as predicted in some SM extensions. A 95% confidence-level upper limit of 183 (87) times the SM cross-section is observed (expected) for non-resonant V hh production when assuming the kinematics are as expected in the SM. Constraints are also placed on Higgs boson coupling modifiers. For the resonant search, upper limits on the production cross-sections are derived for two specific models: one is the production of a vector boson along with a neutral heavy scalar resonance H , in the mass range 260–1000 GeV, that decays into hh, and the other is the production of a heavier neutral pseudoscalar resonance A that decays into a Z boson and H boson, where the A boson mass is 360–800 GeV and the H boson mass is 260–400 GeV. Constraints are also derived in the parameter space of two-Higgs-doublet models.
