Refined BPS invariants of 6d SCFTs from anomalies and modularity

F-theory compactifications on appropriate local elliptic Calabi-Yau manifolds engineer six dimensional superconformal field theories and their mass deformations. The partition function $Z_{top}$ of the refined topological string on these geometries captures the particle BPS spectrum of this class of theories compactified on a circle. Organizing $Z_{top}$ in terms of contributions $Z_\beta$ at base degree $\beta$ of the elliptic fibration, we find that these, up to a multiplier system, are meromorphic Jacobi forms of weight zero with modular parameter the Kaehler class of the elliptic fiber and elliptic parameters the couplings and mass parameters. The indices with regard to the multiple elliptic parameters are fixed by the refined holomorphic anomaly equations, which we show to be completely determined from knowledge of the chiral anomaly of the corresponding SCFT. We express $Z_\beta$ as a quotient of weak Jacobi forms, with a universal denominator inspired by its pole structure as suggested by the form of $Z_{top}$ in terms of 5d BPS numbers. The numerator is determined by modularity up to a finite number of coefficients, which we prove to be fixed uniquely by imposing vanishing conditions on 5d BPS numbers as boundary conditions. We demonstrate the feasibility of our approach with many examples, in particular solving the E-string and M-string theories including mass deformations, as well as theories constructed as chains of these. We make contact with previous work by showing that spurious singularities are cancelled when the partition function is written in the form advocated here. Finally, we use the BPS invariants of the E-string thus obtained to test a generalization of the Goettsche-Nakajima-Yoshioka $K$-theoretic blowup equation, as inspired by the Grassi-Hatsuda-Marino conjecture, to generic local Calabi-Yau threefolds.