  We present a detailed theoretical study of nonfactorizable contributions of
the charm-quark loop to the amplitude of the $B_s\to \gamma\,\gamma$ decay.
This contribution involves the $B$-meson three-particle Bethe-Salpeter
amplitude, $\langle 0|\bar s(y)G_{\mu\nu}(x)b(0)|\bar B_s(p)\rangle$, for which
we take into account constraints from analyticity and continuity. The
charming-loop contribution of interest may be described as a correction to the
Wilson coefficient $C_{7\gamma}$, $C_{7\gamma}\to C_{7\gamma}(1+\delta
C_{7\gamma})$. We calculate an explicit dependence of $\delta C_{7\gamma}$ on
the parameter $\lambda_{B_s}$. Taking into account all theoretical
uncertainties, $\delta C_{7\gamma}$ may be predicted with better than 10\%
accuracy for any given value of $\lambda_{B_s}$. For our benchmark point
$\lambda_{B_s}=0.45$ GeV, we obtain $\delta C_{7\gamma}=0.045\pm 0.004$.
Presently, $\lambda_{B_s}$ is not known with high accuracy, but its value is
expected to lie in the range $0.3\le \lambda_{B_s}({\rm GeV})\le 0.6$. The
corresponding range of $\delta C_{7\gamma}$ is found to be $0.02\le \delta
C_{7\gamma}\le 0.1$. One therefore expects the correction given by charming
loops at the level of at least a few percent.
