Here rather than using the equation for the Poissonian case (which would imply making sure that the rupture is Poissonian), you can already use the more general equation (as in hazard_curves):
pno = rupture.get_probability_no_exceedance(poes)
curves[imt] *= r_sites.expand(
    pno, total_sites, placeholder=1
)