  We investigate the di-Higgs events through $Z_H$-pair production at the CERN
Large Hadron Collider including the pure next-to-leading order (NLO) QCD
correction and the $gg$-fusion contribution in the framework of the littlest
Higgs model with $T$ parity. We employ the diagram subtraction scheme in the
QCD NLO calculations to avoid double counting and keep the convergence of the
perturbative QCD description for the $Z_H$-pair production. We investigate the
dependence of the leading order and QCD corrected integrated cross sections on
the renormalization/factorization scale, and find that the total QCD
corrections slightly reduce the scale uncertainty in the plotted range. By
considering the subsequent decays of the intermediately produced $Z_H$ bosons
and adopting the exclusive four-$b$-jet event selection criterion, the QCD
correction provides considerable enhancement of the kinematic distributions for
final decay products. We find that it is possible to select the signature of
the $Z_H$-pair production from possible standard model background by taking
proper kinematic cuts.
