Deciphering the long-distance penguin contribution to $\bar B_{d, s} \to \gamma \gamma$ decays

We compute for the first time the long-distance penguin contribution to the double radiative $B$-meson decays due to the purely hadronic operators acting with the electromagnetic current in the background soft-gluon field from first field-theoretical principles by introducing a novel subleading $B$-meson distribution amplitude. The numerically dominant penguin amplitude arises from the soft-gluon radiation off the light up-quark loop rather than the counterpart charm-loop effect on account of the peculiar analytical behaviour of the short-distance hard-collinear function. Importantly the long-distance up-quark penguin contribution brings about the substantial cancellation of the known factorizable power correction possessing the same multiplication CKM parameters, thus enabling $B_{d, \, s} \to \gamma \gamma$ to become new benchmark probes of physics beyond the Standard Model.