  Exclusive rare decays mediated by $b\to s\ell\ell$ transitions receive
contributions from four-quark operators that cannot be naively expressed in
terms of local form factors. Instead, one needs to calculate a matrix element
of a bilocal operator. In certain kinematic regions, this bilocal operator
obeys some type of Operator Product Expansion, with coefficients that can be
calculated in perturbation theory. We review the formalism and, focusing on the
dominant SM operators ${\cal O}_{1,2}$, we perform an improved calculation of
the NLO matching for the leading dimension-three operators. This calculation is
performed completely analytically in the two relevant mass scales (charm-quark
mass $m_c$ and dilepton squared mass $q^2$), and we pay particular attention to
the analytic continuation in the complex $q^2$ plane. This allows for the first
time to study the analytic structure of the non-local form factors at NLO, and
to calculate the OPE coefficients far below $q^2=0$, say $q^2 \lesssim
-10\,{\rm GeV}^2$. We also provide explicitly the contributions proportional to
different charge factors, which obey separate dispersion relations.
