jats:titleAbstract</jats:title>jats:pThis paper reports a search for Higgs boson pair (jats:italichh</jats:italic>) production in association with a vector boson (jats:inline-formulajats:alternativesjats:tex-math$$W; {\text {o}r}; Z$$</jats:tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> mml:mrow mml:miW</mml:mi> <mml:mspace /> mml:mrow mml:mtexto</mml:mtext> mml:mir</mml:mi> </mml:mrow> <mml:mspace /> mml:miZ</mml:mi> </mml:mrow> </mml:math></jats:alternatives></jats:inline-formula>) using 139 fbjats:inline-formulajats:alternativesjats:tex-math$$^{-1}$$</jats:tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> mml:msup <mml:mrow /> mml:mrow mml:mo-</mml:mo> mml:mn1</mml:mn> </mml:mrow> </mml:msup> </mml:math></jats:alternatives></jats:inline-formula> of proton–proton collision data at jats:inline-formulajats:alternativesjats:tex-math$$\sqrt{s}=13,\text {TeV}$$</jats:tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> mml:mrow mml:msqrt mml:mis</mml:mi> </mml:msqrt> mml:mo=</mml:mo> mml:mn13</mml:mn> <mml:mspace /> mml:mtextTeV</mml:mtext> </mml:mrow> </mml:math></jats:alternatives></jats:inline-formula> recorded with the ATLAS detector at the Large Hadron Collider. The search is performed in final states in which the vector boson decays leptonically (jats:inline-formulajats:alternativesjats:tex-math$$W\rightarrow \ell \nu ,, Z\rightarrow \ell \ell ,\nu \nu $$</jats:tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> mml:mrow mml:miW</mml:mi> mml:mo→</mml:mo> mml:miℓ</mml:mi> mml:miν</mml:mi> mml:mo,</mml:mo> <mml:mspace /> mml:miZ</mml:mi> mml:mo→</mml:mo> mml:miℓ</mml:mi> mml:miℓ</mml:mi> mml:mo,</mml:mo> mml:miν</mml:mi> mml:miν</mml:mi> </mml:mrow> </mml:math></jats:alternatives></jats:inline-formula> with jats:inline-formulajats:alternativesjats:tex-math$$\ell =e, \mu $$</jats:tex-math><mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> mml:mrow mml:miℓ</mml:mi> mml:mo=</mml:mo> mml:mie</mml:mi> mml:mo,</mml:mo> mml:miμ</mml:mi> </mml:mrow> </mml:math></jats:alternatives></jats:inline-formula>) and the Higgs bosons each decay into a pair of jats:italicb</jats:italic>-quarks. It targets jats:italicVhh</jats:italic> signals from both non-resonant jats:italichh</jats:italic> production, present in the Standard Model (SM), and resonant jats:italichh</jats:italic> 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 jats:italicVhh</jats:italic> 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 jats:italicH</jats:italic>, in the mass range 260–1000 GeV, that decays into jats:italichh</jats:italic>, and the other is the production of a heavier neutral pseudoscalar resonance jats:italicA</jats:italic> that decays into a jats:italicZ</jats:italic> boson and jats:italicH</jats:italic> boson, where the jats:italicA</jats:italic> boson mass is 360–800 GeV and the jats:italicH</jats:italic> boson mass is 260–400 GeV. Constraints are also derived in the parameter space of two-Higgs-doublet models.</jats:p>