  Employing the systematic framework of soft-collinear effective theory (SCET)
we perform an improved calculation of the leading-power contributions to the
double radiative $B_{d, \, s}$-meson decay amplitudes in the heavy quark
expansion. We then construct the QCD factorization formulae for the subleading
power contributions arising from the energetic photon radiation off the
constituent light-flavour quark of the bottom meson at tree level. Furthermore,
we explore the factorization properties of the subleading power correction from
the effective SCET current $J^{(A2)}$ at ${\cal O} (\alpha_s^0)$. The
higher-twist contributions to the $B_{d, \, s} \to \gamma \gamma $ helicity
form factors from the two-particle and three-particle bottom-meson distribution
amplitudes are evaluated up to the twist-six accuracy. In addition, the
subleading power weak-annihilation contributions from both the current-current
and QCD penguin operators are taken into account at the one-loop accuracy. We
proceed to apply the operator-production-expansion-controlled dispersion
relation for estimating the power-suppressed soft contributions to the double
radiative $B_{d, \, s}$-meson decay form factors. Phenomenological explorations
of the radiative $B_{d, \, s} \to \gamma \, \gamma$ decay observables in the
presence of the neutral-meson mixing, including the CP-averaged branching
fractions, the polarization fractions and the time-dependent CP asymmetries,
are carried out subsequently with an emphasis on the numerical impacts of the
newly computed ingredients together with the theory uncertainties from the
shape parameters of the HQET bottom-meson distribution amplitudes.
