"How do theorists actually study rotating black holes?"
When it comes to rotating black holes, the Kerr metric is what becomes central to the analysis. In Boyer-Lindquist coordinates (which are just a generalization of the coordinates used for the metric of a Schwarzschild black hole), it takes the following form:
$$ ds^2 = -\left(1-\frac{2GMr}{\rho^2}\right) dt^2 - \frac{4GMar\sin^2\theta}{\rho^2} dt d\phi + \frac{\rho^2}{\Delta} dr^2 + \rho^2 d\theta^2 + \left(r^2 + a^2 + \frac{2GMa^2 r \sin^2\theta}{\rho^2}\right) \sin^2\theta d\phi^2 $$
with $\rho^2 = r^2 + a^2\cos^2\theta, \Delta = r^2 - 2GMr + a^2$, a being the angular momentum per unit mass of the black hole, G gravitational constant, M mass.
Simply put, physicists analyse the dynamical processes of a spinning black hole collapse by looking at the perturbations of this Kerr metric through the Teukolsky equation, which governs the behaviour of these perturbations. Formally, it is represented as:
$$ \left[\frac{(r^2+a^2)^2}{\Delta} - a^2\sin^2\theta\right] \frac{\partial^2 \psi}{\partial t^2} + \frac{4Mar}{\Delta} \frac{\partial^2 \psi}{\partial t \partial \phi} + \left[\frac{a^2}{\Delta} - \frac{1}{\sin^2\theta}\right] \frac{\partial^2 \psi}{\partial \phi^2} - \Delta^{-s} \frac{\partial}{\partial r} \left(\Delta^{s+1} \frac{\partial \psi}{\partial r}\right) - \frac{1}{\sin\theta} \frac{\partial}{\partial \theta} \left(\sin\theta \frac{\partial \psi}{\partial \theta}\right) - 2s \left[\frac{a(r-M)}{\Delta} + \frac{i\cos\theta}{\sin^2\theta}\right] \frac{\partial \psi}{\partial \phi} - 2s \left[\frac{M(r^2-a^2)}{\Delta} - r - ia\cos\theta\right] \frac{\partial \psi}{\partial t} + (s^2\cot^2\theta - s) \psi = 4\pi \Sigma T,$$
where $\psi$ is the perturbation function, s the "spin weight” of the perturbation (s=0 for scalar fields, $s=\pm 1$ for EM ones, $s=\pm 2$ for gravitational fields), and a (angular momentum), $\Delta$ and $\Sigma$ are functions of r and $\theta$ related to the Kerr metric. You can read more about it here (https://link.springer.com/article/10.1007/s40818-018-0058-8).
Note that solving the latter resorts typically to numerical relativity, which entails discretisation of spacetime into a computational grid due to the complexity and nonlinearity of equation mostly because of the coupling between the black hole's rotation and the perturbations (you can think of this as a manifestation of the frame-dragging effect in rotating spacetimes).
What is known as the Penrose process is also fundamental in the study of rotating black holes. It offers a mechanism to extract energy through the unique geometry of the ergosphere, which is the region exterior to the event horizon where the rotational energy of the black hole becomes accessible. What’s interesting and crucial about this region is that we get the peculiar property that because of frame-dragging (encoded by the Kerr metric’s off-diagonal terms), no physical observer can remain stationary with respect to a distant observer! What the maths then tells us is that a particle entering the ergosphere splits into two, one fragment falling in, while the other escapes to infinity with more energy that the original particle, effectively extracting the black hole’s rotational energy through the Penrose mechanism.
Moreover, by considering a particle falling into a black hole, an important result can be extracted for which I will not go into the detailed mathematics, but simply the key results.
Using the area theorem, one obtains
$$ A = \int \sqrt{|g|} d\theta d\phi = \int \sqrt{(r_+^2 + a^2)^2 \sin^2\theta} d\theta d\phi = 4\pi (r_+^2 + a^2)$$
where the area of the horizon is $S^2$ with $dt=dr=0, r=r_+$
Writing this in terms of mass and angular momentum, then considering their changes for the infalling particle, one obtains the following equation:
$$ \delta M = \left(\frac{\omega_+ \sqrt{M^2 - a^2}}{8\pi a}\right) \delta A + \omega_+ \delta K \equiv \frac{\kappa}{8\pi} \delta A + \omega_+ \delta J,$$
where $\kappa$ is the surface gravity.
Crucially, this shows a correspondence with thermodynamics: through the first law ($ dE=TdS-pdV $), which tells us by analogy that $E\to M$, $S\to A/4$, $T\to \kappa/2\pi$, $pdV\to \omega_+ \delta J$, and directly through the second law that the total area of a black hole event horizon can never decrease, analogously to the total entropy never decreasing in a closed system.
Hopefully, this answer gives you a brief idea of what one may encounter in the study of rotating black holes, which includes the main ideas and tools that physicists use for their study of the latter.
When we look at simple models of gravitational collapse for spherically symmetric matter distribution (such as Oppenheimer-Snyder) we always exploit the fact that the solution to the Einstein field equations just outside any spherically symmetric collapsing matter is Schwarzschild and that this is the only possible solution. This is thanks to Birkhoff's theorem. An analogous theorem exists also for charged matter. Unfortunately, no such theorem exists for rotating matter, as the spherical symmetry is broken. This means that we can't approach the gravitational collapse of rotating matter the same way we would for spherically symmetric charged or uncharged matter collapse.
You are completely correct in this case that unlike the simpler cases of spherically symmetric gravitational collapse, where Birkhoff's theorem provides a clear pathway to understanding the external solutions as Schwarzschild or Reissner-Nordström metrics, rotating systems demand a much more intricate approach.
Briefly put, theorists utilise a blend of methodologies depending on different rotating systems of study. In cases of slow rotating black holes for example, an approach is to treat the rotation as a perturbation to the Schwarzschild solution. However, when you get to highly dynamic and strongly gravitating systems involving rapid rotations, numerical relativity methods are what most likely come into play to study spacetimes whose exact form is unknown - one active area of research currently being the simulation of relativistic binaries and their associated gravitational waves (more about this here: https://arxiv.org/abs/1010.5260), which you may have heard of. Further, when gravitational fields are weak and motion is slow, the post-Newtonian approximation may also sometimes be utilised (especially in the study of binary systems) (read more about it here: https://arxiv.org/abs/1901.08516, https://www.pnas.org/doi/10.1073/pnas.1103127108).
