Exact solutions for differentially rotating galaxies in general relativity

Abstract

A class of stationary axisymmetric solutions of Einstein’s equations for isolated differentially rotating dust sources is presented. The low-energy asymptotic regime is extracted, requiring a self-consistent coupling of quasilocal energy and angular momentum. The Raychaudhuri equation reduces to a balance equation, with two important limits. These limits can be interpreted empirically for rotationally supported configurations such as galaxies. The net energy including quasilocal kinetic contributions vanishes on the inner vortex surface, and the outer rotosurface. These new geometrical objects potentially shed light on virialization. Whether or not abundant collisionless dark matter exists, the new solutions suggest that the phenomenology of galactic rotation curves be fundamentally reconsidered, for consistency with general relativity.

An exact solution of Einstein’s equations for an isolated galaxy sourced by a realistic distribution of stars, treated as a pressureless fluid—dust—is a decades-old open problem that has confounded mathematical relativists for decades. In this Letter we present a new solution to this problem, which unlike previous unsuccessful attempts incorporates the essential physical feature that radial distributions of stars in galaxies rotate by varying amounts. This statement will be precisely defined – but, roughly, each ring of stars in an axially symmetric distribution rotates at slightly different speeds to its neighbours. A consequence of our result is that the amount of collisionless dark matter particles that is conventionally inferred needs to be fundamentally revisited, with profound implications for astrophysics and cosmology.

This problem has been the subject of controversy on account of researchers artificially treating galaxies as rigidly rotating objects [1, 2, 3], which is not only physically unrealistic but is also compounded by mathematical problems [4, 5]. Remarkably, we will find that these issues can be sidestepped by a refinement in understanding the consistent treatment of taking Newtonian limits. In particular, conventionally researchers have naïvely applied nonrelativistic limiting procedures to the gravitational metric before attempting to include nonlinear terms. By contrast, we will include all essential nonlinearities for stationary axisymmetric dust spacetimes before considering the low–energy limit.

The nonlinearity of general relativity arises from the self–interaction of matter and geometry via Einstein’s equations. In place of additive gravitational potentials on a fixed background, the time–averaged motion of matter sources defines regional backgrounds with their own quasilocal energy and angular momentum content [6, 7]. Understanding the hierarchy of regional scales, or the fitting of one regional geometry into another, is an important foundational question in cosmology [8]. Conventionally one assumes the existence of a global asymptotic Minkowski background before applying nonrelativistic limits. Mass and angular momentum are then defined by ideal asymptotic Killing vectors, which obscures their essential quasilocal origin in general relativity. The manner in which quasilocal energy and angular momentum approach asymptotic values in a nonempty universe is little appreciated. It may nonetheless be key to understanding the biggest open problems in physical cosmology.

In this Letter we derive new solutions with self-consistent coupling between quasilocal energy and angular momentum, which can be applied in the low-energy limit to systems at different scales in the fitting problem. Here our primary interest is rotating galaxies. However, our approach will likely further inform other attempts to model the dynamical properties of coarse-grained, cosmological structures in general relativity [9, 10, 11, 12, 13].

Let us consider stationary axisymmetric solutions of the Einstein equations with a dust source. Our solutions will thus apply to numerous astrophysical systems including stellar systems, galaxies and putative dark matter halos on larger scales. As applied to galaxies assuming stationarity—i.e., the presence of a timelike Killing vector, $\partial/\partial t$—means the solutions apply to short time scales relative to those involving galaxy formation and evolution, denoted $t_{\text{evo}}$. Likewise, imposing axisymmetry—i.e., an axial Killing vector, $\partial/\partial\phi$—means that our solutions then apply to an effective dust fluid for time-averaged oscillatory motion, $t_{\text{osc}}$, of individual stars and gas above and below the galactic plane. At the present epoch $t_{\text{osc}}\mathop{\sim}\limits 10\,$Myr and $t_{\text{evo}}\mathop{\sim}\limits 1\,$Gyr. Thus the effective dust fluid applies on time scales, $10^{7}\mathop{\hbox{${\lower 3.8pt\hbox{$<$}}\atop{\raise 0.2pt\hbox{$\sim$}}$}}t\mathop{\hbox{${\lower 3.8pt\hbox{$<$}}\atop{\raise 0.2pt\hbox{$\sim$}}$}}10^{9}\,$yr for nearby galaxies.

The spacetime metrics we consider can thus be written in the Lewis–Papapetrou–Weyl form, viz.

$$\mathop{\text{d}\!}s^{2}=-c^{2}e^{2\Phi(r,z)/c^{2}}(\mathop{\text{d}\!}t+A(r,z)\mathop{\text{d}\!}\phi)^{2}+e^{-2\Phi(r,z)/c^{2}}\left[r^{2}\mathop{\text{d}\!}\phi^{2}+e^{2k(r,z)/c^{2}}(\mathop{\text{d}\!}r^{2}+\mathop{\text{d}\!}z^{2})\right]\,,$$ (1)

where $r$ is a radial coordinate, $z$ is an axial coordinate, $\Phi(r,z)$ is related to the conventional Newtonian potential, $A(r,z)$ is related to frame-dragging, and $k(r,z)$ is a conformal factor on the 2–dimensional space of orbits of the isometry group generated by the Killing vectors.

The energy-momentum tensor takes the form

$$T^{\mu\nu}=c^{2}\rho(r,z)\,U^{\mu}U^{\nu}\,,$$ (2)

where $\rho(r,z)$ is the density of particles, each with a 4–velocity $U^{\mu}$, given by

$$U^{\mu}\partial_{\mu}=(-H)^{-1/2}\left(\partial_{t}+\Omega\,\partial_{\phi}\right)\,,$$ (3)

where $d\phi/dt=\Omega(r,z)$ uniquely defines the angular speed of rotation at any point, and $H(r,z)$ is a normalization factor. Since $U^{\mu}U_{\mu}=-c^{2}$, it follows that

$$H=-e^{2\Phi/c^{2}}(1+A\,\Omega)^{2}+e^{-2\Phi/c^{2}}r^{2}\,\Omega^{2}/c^{2}\,.$$ (4)

Einstein’s equations are solved by a hierarchy of quadratures [14, 15, 16]. One finds that $H=H(\eta)$ and $\Omega=\Omega(\eta)$, where $\eta(r,z)$ is the angular momentum of the dust particles measured by Zero Angular Momentum Observers (ZAMOs), and

$$\mathop{\text{d}\!}\Omega=c^{2}\mathop{\text{d}\!}H/({2\eta})\,.$$ (5)

In view of (5) we call $H$ the differential rotation state function, since a choice of $H$ fixes $\Omega$. The norms of the two Killing vectors, $\xi^{\mu}\partial_{\mu}=\partial_{t}$ and $\zeta^{\mu}\partial_{\mu}=\partial_{\phi}$, and their inner product are then respectively given by

$$\xi^{\mu}\xi_{\mu}=-c^{2}e^{2\Phi/c^{2}}=c^{2}\frac{\left(H-\eta\,\Omega/c^{2}\right)^{2}-\left(r\,\Omega/c\right)^{2}}{H}\,,$$ (6)

$$\zeta^{\mu}\zeta_{\mu}=r^{2}e^{-2\Phi/c^{2}}-c^{2}A^{2}e^{2\Phi/c^{2}}=\frac{r^{2}-(\eta/c)^{2}}{(-H)}\,,$$ (7)

$$\xi^{\mu}\zeta_{\mu}=-c^{2}Ae^{2\Phi/c^{2}}=\eta-\frac{(\eta/c)^{2}-r^{2}}{H}\,\Omega\,.$$ (8)

The remaining Einstein equations yield

$$\Xi_{,r}=\frac{1}{2r}\left[g_{tt,z}g_{\phi\phi,z}-g_{tt,r}g_{\phi\phi,r}-(g_{t\phi,z})^{2}+(g_{t\phi,r})^{2}\right]\,$$

$$\Xi_{,z}=\frac{1}{2r}\left[2g_{t\phi,r}g_{t\phi,z}+g_{tt,r}g_{\phi\phi,z}-g_{tt,z}g_{\phi\phi,r}\right],$$ (9)

where $\Xi(r,z)=[\Phi(r,z)-k(r,z)]/c^{2}$. Finally, further requiring that $\Xi_{,rz}=\Xi_{,zr}$ we find that the entire class of solutions is fully determined by a choice of $H$ and a solution of the homogeneous Grad–Shafranov equation [17, 18]

$$\Psi_{,rr}-\frac{1}{r}\Psi_{,r}+\Psi_{,zz}=0,$$ (10)

where

$$\Psi=\eta+c^{2}\frac{r^{2}}{2}\int\frac{H^{\prime}}{H}\frac{d\eta}{\eta}-\frac{1}{2}\int\frac{H^{\prime}}{H}\eta\,d\eta\,.$$ (11)

Eq. (11) is an integrability condition for the dust geodesic equations in the absence of pressure in the effective fluid. Furthermore, to model galaxies, we are interested in solutions with reflection symmetry with respect to the $z=0$ plane. The solution to (10) is then given by

$$\Psi(r,z)=\,{\mathpzc C}_{0}\;+\; r\sum_{n=0}^{\infty}{\mathpzc A}_{n}I_{1}({\mathpzc a}_{n}r)\cos({\mathpzc a}_{n}z)$$

$$+r\sum_{m=0}^{\infty}{\mathpzc B}_{m}K_{1}({\mathpzc b}_{m}r)\cos({\mathpzc b}_{m}z),$$ (12)

where $I_{1}$ and $K_{1}$ are the modified Bessel functions of the first and second kind, the constants ${\mathpzc a}_{n},{\mathpzc b}_{m}$ are positive, whilst the constants ${\mathpzc C}_{0}$, ${\mathpzc A}_{n}$ and ${\mathpzc B}_{m}$ are arbitrary.

The trace of the Einstein equations gives

$$8\pi G\rho=\frac{\eta_{,r}^{2}+\eta_{,z}^{2}}{4\eta^{2}r^{2}}\left[\eta^{2}\Delta^{2}(\eta)-c^{4}r^{4}\ell(\eta)^{2}\right]e^{2\Xi}\!,$$ (13)

where $\ell(\eta)=H^{\prime}/H$ and $\Delta(\eta)=2-\eta\,\ell(\eta)$, with $H^{\prime}=dH/d\eta$, while the Raychaudhuri equation reduces to

$$\dot{\Theta}\equiv U^{\nu}\nabla_{\nu}\Theta=-\frac{1}{3}\Theta^{2}+\omega^{2}-\sigma^{2}-R_{\mu\nu}U^{\mu}U^{\nu},$$ (14)

where $\Theta=\nabla_{\mu}U^{\mu}$, $\omega^{2}=\omega^{\mu\nu}\omega_{\mu\nu}$ and $\sigma^{2}=\sigma^{\mu\nu}\sigma_{\mu\nu}$ are the expansion, vorticity and shear scalars respectively. Here $\omega_{\mu\nu}=U_{[\mu;\nu]}$, $\sigma_{\mu\nu}=U_{(\mu;\nu)}-\frac{1}{3}h_{\mu\nu}\Theta$, $h_{\mu\nu}=g_{\mu\nu}+c^{-2}U_{\mu}U_{\nu}$. The dust is non-expanding, $\Theta=0$, $\dot{\Theta}=0$, and $R_{\mu\nu}U^{\mu}U^{\nu}=4\pi\,G\rho$, so that (14) reduces to a balance condition

$$\omega^{2}-\sigma^{2}=4\pi G\rho\,.$$ (15)

Significantly, since the effective dust energy density contains quasilocal kinetic energy, $\rho$ may be negative when shear dominates over vorticity. Indeed,

$$\sigma^{2}=\frac{e^{2\Xi}c^{4}r^{2}(\eta_{,r}^{2}+\eta_{,z}^{2})\,\ell^{2}(\eta)}{8\eta^{2}},$$ (16)

which combined with (13) yields

$$\omega^{2}=\frac{e^{2\Xi}(\eta_{,r}^{2}+\eta_{,z}^{2})\,\Delta^{2}(\eta)}{8r^{2}}\,.$$ (17)

Physical measurements are defined by relevant classes of observers. We first introduce the ZAMO coframe [19]

$${\mathbf{\omega}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}=\Bigl{(}\frac{r\mathop{\text{d}\!}t}{\sqrt{g_{\phi\phi}}}\,,\sqrt{g_{\phi\phi}} \left(\mathop{\text{d}\!}\phi\,-\,\chi\mathop{\text{d}\!}t\right)\,,$$

$$e^{(k-\Phi)/c^{2}}\mathop{\text{d}\!}r\,,e^{(k-\Phi)/c^{2}}\mathop{\text{d}\!}z\Bigr{)}\,,$$ (18)

$g_{\phi\phi}=e^{-2\Phi/c^{2}}r^{2}-c^{2}A^{2}e^{2\Phi/c^{2}}$, and $\chi=-g_{t\phi}/g_{\phi\phi}$ with $g_{t\phi}=-c^{2}\,A\,e^{2\Phi/c^{2}}$. The dual ZAMO tetrad frame is then

$${\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}=\Bigl{(}\frac{1}{r}{\sqrt{g_{\phi\phi}}} \left(\partial_{t}\,+\,\chi\partial_{\phi}\right)\,,\frac{1}{\sqrt{g_{\phi\phi}}}\partial_{\phi}\,,$$

$$\qquad e^{-(k-\Phi)/c^{2}}\partial_{r}\,,e^{-(k-\Phi)/c^{2}}\partial_{z}\Bigr{)}\,.$$ (19)

Starting from the formula for the dust particles’ velocity measured by ZAMOs, $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}=({{{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}}^{1}}\cdot{\mathbf{U}})/({{{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}}^{0}}\cdot{\mathbf{U}})$, we find

$$\eta(r,z)=r\,v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}(r,z)\,.$$ (20)

Furthermore, let us define an effective Lorentz factor $\gamma_{\lower 2.0pt\hbox{$\scriptstyle Z$}}:={{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}}^{0}\cdot{\mathbf{U}}$, so that

$$-H\,\gamma_{\mathrm{Z}}^{2}\,v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}=r\,(\Omega-\chi)\,.$$ (21)

The second congruence of relevance are ideal Stationary Observers (SOs), with coframe [20]

$${\mathbf{\omega}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}=\Bigl{(}e^{\Phi/c^{2}} \left(\mathop{\text{d}\!}t\,+\,A\mathop{\text{d}\!}\phi\right)\,,r\,e^{\Phi/c^{2}}\mathop{\text{d}\!}\phi\,,$$

$$\qquad e^{(k-\Phi)/c^{2}}\mathop{\text{d}\!}r\,,e^{(k-\Phi)/c^{2}}\mathop{\text{d}\!}z\Bigr{)}\,,$$ (22)

The dual SO tetrad frame is then

$${\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}=\Bigl{(}e^{-\Phi/c^{2}} \partial_{t}\,,\frac{e^{\Phi/c^{2}}}{r}\left(\partial_{\phi}-A\,\partial_{t}\right)\,,$$

$$\qquad e^{-(k-\Phi)/c^{2}}\partial_{r}\,,e^{-(k-\Phi)/c^{2}}\partial_{z}\Bigr{)}\,.$$ (23)

By analogy to the ZAMO case for SOs we define,$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}=({{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}^{1}}\cdot{\mathbf{U}})/({{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}^{0}}\cdot{\mathbf{U}})$ and $\gamma_{\lower 2.0pt\hbox{$\scriptstyle S$}}:={{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}}^{0}\cdot{\mathbf{U}}$, so that

$$r\,v_{\lower 2.0pt\hbox{$\scriptstyle S$}}(r,z)=\gamma_{\mathrm{S}}^{-1}e^{-\Phi/c^{2}}r^{2}\Omega\,.$$ (24)

We can define $\eta_{\lower 2.0pt\hbox{$\scriptstyle S$}}=v_{\lower 2.0pt\hbox{$\scriptstyle S$}}r$ by analogy to (20). However, it plays no particular role here. Equivalently, we have

$$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}-v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}=\frac{{{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}^{1}}\cdot{{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}^{0}}}{{{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle S$}}^{0}}\cdot{{\mathbf{e}}_{\lower 2.0pt\hbox{$\scriptstyle Z$}}^{0}}}=\frac{g_{\phi\phi}}{r^{2}}r\chi=\frac{-g_{t\phi}}{r},$$ (25)

so that

$$e^{2\Phi/c^{2}}\gamma_{\lower 2.0pt\hbox{$\scriptstyle S$}}v_{\lower 2.0pt\hbox{$\scriptstyle S$}}=-H\,\gamma_{\mathrm{Z}}^{2}\,(v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}+r\,\chi)\,.$$ (26)

Eq. (25) shows that ZAMOs and SOs measure differentdust velocities at any spacetime event, with a difference directly proportional to the frame-dragging contribution.

In the present framework it is clear how observations have been misapplied to theoretical observables in the past. In particular, several analyses [1, 2, 3] misapply $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$ to the rotation curves of distant galaxies, when $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$ given by (24) is the relevant velocity, since

$$\frac{v_{\lower 2.0pt\hbox{$\scriptstyle S$}}}{c}=\frac{\sqrt{1+\left(2e^{2\Phi/c^{2}}r\Omega/c\right)^{2}}-1}{2e^{2\Phi/c^{2}}r\Omega/c}\simeq\frac{r\Omega}{c}e^{2\Phi/c^{2}}+\dots$$ (27)

coinciding, at leading order, with the widely used special relativistic interpretation of the redshift.

To complete our identification of physically relevant velocities in the general case we supplement the velocities $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$ and $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$ by the kinetic and dragging velocities [21, 18]

$$v_{\lower 2.0pt\hbox{$\scriptstyle K$}}:=r\,\Omega\,,$$

$$v_{\lower 2.0pt\hbox{$\scriptstyle D$}}:=r\,\chi\,,$$ (28)

respectively. While $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$ and $v_{\lower 2.0pt\hbox{$\scriptstyle K$}}$ coincide at low velocities, in general their differences can be physically important.

Let us now derive the functional form applicable to systems such as galaxies, with subrelativistic local relative speeds, $v\ll c$, weak pseudo-Newtonian potentials, $\Phi\mathop{\sim}\limits v^{2}/c^{2}$, and nonrelativistic frame-dragging.

From these conditions, we find $A=r\,v_{\lower 2.0pt\hbox{$\scriptstyle D$}}/c^{2}+\mathcal{O}({v^{3}/c^{3}})$. Thus by (4), $H=-1+\mathcal{O}({v^{2}/c^{2}})$. Since $H=H(\eta)$, using (20) and noting that $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}=v_{\lower 2.0pt\hbox{$\scriptstyle K$}}-v_{\lower 2.0pt\hbox{$\scriptstyle D$}}+\mathcal{O}(v^{3}/c^{3})$, we find

$$H(\eta)=-1+\epsilon\,\eta^{2}/(b_{\epsilon}c)^{2}+\mathcal{O}(v^{3}/c^{3})\,,$$ (29)

where $\epsilon\in\{-1,0,1\}$ and $b_{\epsilon}$ are constant lengths which can be identified as impact parameters. We stress that the functional form of $H(\eta)$ in (29) is the only possible choice that allows a dust system to be consistently considered in a low-energy regime. Substituting just the first two terms of (29) into (11) we obtain

$$\eta(r,z)={b_{\epsilon}c}{\mathop{\text{\,tann}\!}}\left({\frac{b_{\epsilon}\,\Psi(r,z)}{b_{\epsilon}^{2}c-\epsilon\,r^{2}c}}\right)\,,$$ (30)

where

$${\mathop{\text{\,tann}\!}}\,(x)=\begin{cases}\tanh(x),&\epsilon=+1\\ x,&\epsilon=0\\ \tan(x),&\epsilon=-1\\ \end{cases}.$$ (31)

Moreover, by (5) it follows that

$$\Omega=\Omega_{0}+\epsilon\frac{\eta}{b_{\epsilon}^{2}}.$$ (32)

Furthermore, to ensure zero rotation on the symmetry axis we choose $\Omega_{0}=0$. We thus find

$$v_{\lower 2.0pt\hbox{$\scriptstyle K$}}=\epsilon\frac{rc}{b_{\epsilon}}\,{\mathop{\text{\,tann}\!}}\left({\frac{b_{\epsilon}\,\Psi(r,z)}{b_{\epsilon}^{2}c-\epsilon\,r^{2}c}}\right)\,,$$ (33)

and

$$v_{\lower 2.0pt\hbox{$\scriptstyle D$}}=c\left(\frac{r}{b_{\epsilon}}\epsilon-\frac{b_{\epsilon}}{r}\right)\,{\mathop{\text{\,tann}\!}}\left({\frac{b_{\epsilon}\,\Psi(r,z)}{b_{\epsilon}^{2}c-\epsilon\,r^{2}c}}\right).$$ (34)

The presence of differential rotation is crucial mathematically since the term $\epsilon r^{2}c$ in (30) regularizes potential divergences at infinity in (Exact solutions for differentially rotating galaxies in general relativity). In particular, we can choose solutions to (10) which are regular about the $z$-axis and diverge as $\Psi\mathop{\sim}\limits r^{\alpha}$, $\alpha<2$ as $r\to\infty$. All these cases are phenomenologically interesting as the sub-quadratic divergences in the modified Bessel function series $I_{1}({\mathpzc a}_{n}r)$ in (Exact solutions for differentially rotating galaxies in general relativity), are thereby regularized. This option was unavailable in the rigidly rotating Balasin-Grumiller (BG) model [3] where only the modified Bessel function series $K_{1}({\mathpzc b}_{n}r)$ in (Exact solutions for differentially rotating galaxies in general relativity) was chosen, being regular as $r\to\infty$ but singular on the $z$-axis. Thus, we find solutions that display self-consistency in the low-velocity regime.

By (13) and (29) in each limit the effective density is

$$8\pi G\rho_{\epsilon}=\frac{\eta_{,r}^{2}+\eta_{,z}^{2}}{2r^{2}}\left[\frac{b_{\epsilon}^{4}-\epsilon^{2}r^{4}}{\left(b_{\epsilon}^{2}-\eta^{2}/c^{2}\right)^{2}}\right]e^{2\Xi}\,,$$ (35)

which vanishes on the cylinders $r=b_{\epsilon}$, independent of the choice of solution to (10). The global structure of the spacetime is now determined by: (i) choosing $\rho=\rho_{+}-\rho_{-}$, $b_{+}>b_{-}>0$ with an $H(b_{+},b_{-},\eta)=-1+f(b_{+},b_{-})\eta^{2}$ such that $\lim_{b_{-}\to 0}f=1/(b_{+}c)^{2}$ and $\lim_{b_{-}\to\infty}f=-1/(b_{-}c)^{2}$; (ii) requiring consistent coupling of quasilocal energy and angular momentum in the ultraviolet and infrared limits that bound the interior and exterior respectively. The ultraviolet limit involves the relevant physics of the coarse-grained matter source, whilst the infrared limit involves the approach to asymptotic infinity far from the compact source. While $r=b_{-}$ and $r=b_{+}$ are precisely defined boundaries, a smooth transition from the virial interior to exterior involves some phenomenological freedom. See Table 1.

Spacetime region		$\vphantom{\Big{|}}v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$	$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$	$\rho$
UV limit: Vortex	$\vphantom{\Big{|}}\epsilon=-1,r<b_{-}$	$v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}<0$	$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}>0$	$\rho<0$
Virial interior	$\vphantom{\Big{|}}\epsilon=-1,r>b_{-}$	$v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}<0$	$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}>0$	$\rho>0$
Virial exterior	$\vphantom{\Big{|}}\epsilon=+1,r<b_{+}$	$v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}>0$	$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}>0$	$\rho>0$
IR limit: Virialized	$\vphantom{\Big{|}}\epsilon=+1,r>b_{+}$	$v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}<0$	$v_{\lower 2.0pt\hbox{$\scriptstyle S$}}<0$	$\rho<0$

Table 1: ZAMO and SO velocities, $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$ and $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$, and effective compact quasilocal energy density, $\rho$.

There are only two cases in which $\Theta=0$, $\omega=\sigma$ and $\rho=0$; i.e., the dust congruence has zero expansion, vorticity and shear exactly cancel, and the net quasilocal energy density, $4\pi G\rho=R_{\mu\nu}U^{\mu}U^{\nu}$, vanishes.

In the interior ($\epsilon=-1$), $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$ and $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$ are antiparallel, since to maintain zero angular momentum a ZAMO must counter the dragging of the central rotation. The inner vortex surface arises when maintain this equilibrium requires the net quasilocal energy density to change sign, becoming shear-dominated. On the vortex boundary

$$g_{\phi\phi}=b_{-}^{2}\left(1-\frac{1}{2}\frac{\Psi(r,z)^{2}}{c^{2}b_{-}^{2}}\right)+\mathcal{O}(v^{4}/c^{4})\,,$$ (36)

by (7), (30), and other metric components are also regular, all relative velocities, $v$, being subrelativistic. The inner vortex surface is found at $r=b_{-}\gg L$, where $L$ is a typical scale length of the coarse-grained dust. E.g., for the Milky Way, $L\mathop{\sim}\limits 50\,$kpc. This scale is not unique, and depends crucially on phenomenological interpretation. Crosta et al. [22, 23] used the BG model to fit a rotation curve to the Milky Way, as inferred from the proper motion and redshifts of a sample of $7.2\times 10^{6}$ young stars from the GAIA–DR3 survey, lying within $r<19\,$kpc of the galactic centre and $-1<z<1\,$kpc about the galactic plane. By (30) for large values of the impact parameter, $b$, almost rigid rotation is obtained. This explains why the $K_{1}({\mathpzc b}_{n}r)$ series alone can be fit to observations close to the galactic plane and far from its centre. However, the BG model will inevitably fail globally on account of its axial singularity, which we sidestep via the differential rotation regularization term.

In the exterior ($\epsilon=+1$), $v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$ and $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}$ are parallel. At large spatial distances $v_{\lower 2.0pt\hbox{$\scriptstyle S$}}-v_{\lower 2.0pt\hbox{$\scriptstyle Z$}}$ does not vanish in general, even for asymptotically flat spacetimes. Their difference embodies the quasilocal angular momentum that integrates to a global asymptotic charge in the case of the Kerr geometry. However, the observed universe has zero global angular momentum. To enable direct embedding into a nonrotating cosmological background, differentially rotating distributed sources should consistently couple quasilocal energy and angular momentum.

Our class of exact solutions does indeed admit such geometries, with a bounding surface which we hereby name the rotosurface. It has the following novel geometric properties. In the limit $r\to b_{+}$ at finite $z$, even though the Killing vector norms diverge, the vorticity, shear, Ricci and Kretschmann scalars converge to zero: $\omega^{2}\mathop{\sim}\limits\sigma^{2}\mathop{\sim}\limits R_{\mu\nu\alpha\beta}R^{\mu\nu\alpha\beta}\mathop{\sim}\limits e^{2\Xi}\Psi^{2}/(8|b_{+}-r|^{4}])\to 0$, as $\Psi$ is finite and $e^{2\Xi}\mathop{\sim}\limits e^{-\Psi^{2}/(bc^{2}|b_{+}-r|)}$ as $r\to b_{+}$.

The rotosurface is a new mathematical object, whose topology, geometry, fundamental physics and cosmological implications remain to be fully explored. Analogously to an acceleration horizon for uniform linear acceleration in Minkowski space, it is a limiting rotating surface where a fictitious observer would need to rotate with $v\to c$ in order to have zero angular momentum with respect to the vanishing local dust. However, it is a timelike surface. Indeed, it provides a precise mathematical definition of one type of finite infinity surface [6, 24, 25]. Furthermore, since relevant curvature scalars vanish it provides a new mathematical setting for the concept of virialization, and the nomenclature adopted in Table 1.

Beyond the rotosurface other quasilocal energy sources will dominate: the effective thermal pressure of galaxies in clusters and the kinetic energy of expansion in a void dominated universe [26]. The fact that the quasilocal energy of the isolated system is negative beyond the rotosurface is consistent with: (i) positive spatial curvature energy and binding energy being negative relative to the background; (ii) a small kinematical backreaction, $\Omega_{\cal Q}$, of virialized structures that is negative in the Buchert averaging scheme [27, 28].

Since the new solutions successfully exhibit the essential physics of differential rotation they naturally apply to rotationally supported systems, e.g., disc galaxies. A sample exact solution for a Milky Way-like galaxy with a disc mass of $10^{11}{\rm M}_{\odot}$ and rotosurface at $700\,$kpc (Fig. 1) produces a realistic rotation curve (Fig. 2) consistent with our quasilocal Newtonian limit. The interplay of central bulges and/or dark matter halos is more complex. Nonetheless, having uncovered a consistent incorporation of the low-energy regime in macroscopic general relativity, extensions to pressure-supported systems follow naturally [29]. While the scales of globular clusters, dwarf galaxies and virialized galaxy clusters are different—each involving different particles—the underlying equations are universal.

Finally, we note that since the rotosurface is regular but timelike it may provide a novel mathematical arena for understanding the asymptotic charges of the Bondi–Metzner–Sachs group.

Figure 1: Angular momentum density for a Milky Way-like galaxy. The rotosurface is found at 700 kpc.

Figure 2: Rotation velocity curve for the Milky Way-like galaxy of Fig. 1. The virial interior and exterior of Table 1 are found beyond the scales shown.

Acknowledgments DLW is supported by Marsden Fund grant M1271 administered by the Royal Society of New Zealand, Te Apārangi. We thank Roy Kerr,Federico Re, and Chris Stevens for important physical and mathematical insights that helped us refine our understanding. We thank John Forbes, Christopher Harvey-Hawes, Morag Hills, Emma Johnson, Zachary Lane, Antonia Seifert, Shreyas Tiruvaskar and Michael Williams for useful discussions, and Frederic Hessman for correspondence.