The Higgs self-coupling is an important parameter in the standard model to clarify the electroweak symmetry breaking mechanism,
which can be measured via Higgs boson pair productions at hadron colliders.
We present a fully differential NNLO QCD calculation of the Higgs boson pair production in association
with a massive gauge boson at hadron colliders.
For $Whh$ production, the NNLO effects reduce the scale uncertainties significantly, and are sizeable in the large transverse momentum region.
For $Zhh$ production, the NNLO effects enhance the NLO total cross sections by $20\%- 50\%$, and change the shape of NLO kinematic distributions. The impact on the extraction of the self-coupling from experiments is significant.
These theoretical results can be utilised in future experimental analysis.