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 by including the perturbative resummation of enhanced logarithms of m b /Λ QCD at the next-to-leadinglogarithmic accuracy. 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 ξ) by virtue of the operator identities due to the classical equations of motion. The higher-twist contributions to the B d, s → γγ helicity form factors from the two-particle and three-particle bottom-meson distribution amplitudes are evaluated with the perturbative factorization technique, up to the twist-six accuracy. In addition, the subleading power weak-annihilation contributions from both the currentcurrent 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, which cannot be factorized into the light-cone distribution amplitudes of the heavy-meson and the resolved photon as well as the hard-scattering kernel calculable in perturbation theory canonically. Phenomenological explorations of the radiative B d, s → γ γ 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.