  We present a phenomenological study of the rare double radiative decay $B\to
K \gamma\gamma$ 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 $\sim 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\to K \gamma\gamma$ matrix element is then expressed in terms of the
usual $B\to K$ form factors known from semileptonic rare decays. The integrated
SD branching ratio in the SM in the OPE region turns out to be $\Delta
{\cal{B}}(B \to K \gamma \gamma)_{SM}^{OPE} \simeq 1 \times 10^{-9}$. We work
out the di-photon invariant mass distribution with and without the resonant
background through $B\to K \{\eta_c,\chi_{c0}\}\to K\gamma \gamma$. 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 \to \tau^+ \tau^-$
decays, which constrains the scalar/pseudoscalar Four-Fermi operators with
$\tau^+ \tau^-$, leaves considerable room for new physics in the
one-particle-irreducible contribution to $B\to K \gamma \gamma$ decays. In this
case, we find that the SD $B\to K \gamma \gamma$ 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.
