Exclusive rare decays mediated by b → s 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 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 −10 GeV 2 . We also provide explicitly the contributions proportional to different charge factors, which obey separate dispersion relations. , 0, 0]For diagrams e there are 5 MIs: J e,1 = j[e, 0, 1, 0, 0, 1, 0, 0] J e,2 = j[e, 0, 1, 0, 1, 1, 0, 0] J e,3 = j[e, 0, 1, 1, 0, 1, 0, 0] (B.5) J e,4 = j[e, 0, 2, 0, 1, 1, 0, 0] J e,5 = j[e, 1, 1, 1, 0, 1, 0, 0]