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.