We present a phenomenological study of the rare double radiative decay B → Kγγ in the Standard Model (SM) and beyond. Using the operator product expansion (OPE) technique, we estimate the short-distance (SD) contribution to the decay amplitude in a region of the phase space which is around the point where all decay products have energy ∼ m b /3 in the rest frame of the B-meson. At lowest order in 1/Q, where Q is of order m b , the B → Kγγ matrix element is then expressed in terms of the usual B → K form factors known from semileptonic rare decays. The integrated SD branching ratio in the SM in the OPE region turns out to be ∆B(B → Kγγ) OP E SM ≃ 1 × 10 −9 . We work out the di-photon invariant mass distribution with and without the resonant background through B → K{η c , χ c0 } → Kγγ. In the SM, the resonance contribution is dominant in the region of phase space where the OPE is valid. The present experimental upper limit on B s → τ + τ − decays, which constrains the scalar/pseudoscalar Four-Fermi operators with τ + τ − , leaves considerable room for new physics in the one-particle-irreducible contribution to B → Kγγ decays. In this case, we find that the SD B → Kγγ branching ratio can be enhanced by one order of magnitude with respect to its SM value and the SD contribution can lie outside of the resonance peaks. *