  We study the rare decay $Z_H \rightarrow H^0 \gamma$ at
$\mathcal{O}(\alpha_{{\rm ew}}^{3}\alpha_s)$ accuracy including two-loop QCD
corrections in the context of the littlest Higgs model (LHM) without
$T$-parity. We revisit analytically and numerically the leading-order (LO)
contributions of the one-loop diagrams induced by massive fermions, scalars and
charged gauge bosons in the LHM, and further study the NLO QCD correction to
this decay process. We perform the numerical calculation by taking the LHM
input parameters $f=3,~4~\text{TeV}$ and $0.1< c <0.6$, and discuss the
numerical results of the decay width up to the QCD NLO within the recent
experimentally constrained LHM parameter space region. Our results show that
the two-loop QCD correction always reduces the LO decay width and the
top-induced QCD correction is the dominant contribution at the QCD NLO. For
$f=4~\text{TeV}$ and $c=0.3$, the NLO QCD corrected decay width reaches
$75.099~\text{keV}$ and the NLO QCD relative correction is about $-11.0\%$.
