In the end, what worked was those 3 lines of code:

expr = V_r_integral_minus_V_theta_integral.subs(sp.sin(theta)**2, (1-sp.cos(2*theta))/2)
expr = sp.simplify(expr)
separated_expr = sp.collect(expr, [sp.cos(2*theta), sp.log(r), r**2])

So yeah, manually collect all terms, and manually use a trigonometric identity. I don't know if there's a more straightforward way to do it.