We discuss how LHC di-muon data collected to study $B_q \to \mu\mu$ can be used to constrain light particles with flavour-violating couplings to $b$-quarks. Focussing on the case of a flavoured QCD axion, $a$, we compute the decay rates for $B_q \to \mu \mu a$ and the SM background process $B_q \to \mu \mu \gamma$ near the kinematic endpoint. These rates depend on non-perturbative $B_q \to \gamma^{(*)}$ form factors with on- or off-shell photons. The off-shell form factors -- relevant for generic searches for beyond-the-SM particles -- are discussed in full generality and computed with QCD sum rules for the first time. With these results, we analyse available LHCb data to obtain the sensitivity on $B_q \to \mu \mu a$ at present and future runs. We find that the full LHCb dataset alone will allow to probe axion-coupling scales of the order of $10^6$ GeV for both $b\to d$ and $b \to s$ transitions.
Comment: 37 pages, 6 figures, version to appear in JHEP, added multiple subtracted dispersion relation to extend in resonance region including ancillary Mathematica notebook with Form Factors