Stationary solutions of the axially symmetric Einstein–Vlasov system: present status and open problems

The Einstein–Vlasov system consists of the coupled equations for the metric tensor g and density function f on phase space, which in arbitrary coordinates and geometric units ($G = c = 1$) reads

$$\begin{align*} Ric\left(g\right)_{ij} - \frac 12 R\left(g\right) g_{ij} = 8 \pi T\left(g,f\right)_{ij}, \quad p^i \partial_{x^i} f - \Gamma^{k}_{ij}\left(g\right) p^i p^j \partial_{p^k} f = 0, \end{align*}$$

where the second equation is called the Vlasov equation. Here $Ric(g), R(g)$ are the Ricci tensor and scalar of the metric g, $\Gamma^{k}_{ij}(g)$ are the Christoffel symbols of the metric g, and $T_{ij}(g,f)$ is the energy momentum tensor associated with the Vlasov matter. In a coordinate frame $(p^0, p^1, p^2, p^3)$ on the tangent space at $x \in M$ the energy momentum tensor takes the form

$$\begin{equation} T_{ij}\left(x\right) : = \int_{P_x} p_i p_j f\left(x,p\right)\frac{\sqrt{-\text{det} g\left(x\right)}}{- p_0} \mathrm{d} p^1 \mathrm{d} p^2 \mathrm{d} p^3 \end{equation} \tag{ 1 }$$

where P_x is the mass shell at point x and $p_0 = g_{0i} p^i$. For the stationary axisymmetric spacetimes considered in the papers [2, 3] the metric is written in axial coordinates $(t, \rho, z, \phi)$ as

$$\begin{equation} g = - e^{2 \nu} \textrm{d}t^2 + e^{2 \mu} \textrm{d}\rho^2 + e^{2 \mu}\textrm{d}z^2 + \rho^2 B^2 e^{-2 \nu} \left(\textrm{d}\varphi- \omega \textrm{d}t\right)^2, \end{equation} \tag{ 2 }$$

where the metric fields $\nu, \mu, B, \omega$ depend only on the coordinates $\rho, z$. Note that ρ = 0 is the axis of symmetry, and that $(\rho, z)$ are cylindrical coordinates at infinity in the sense that in the appropriate limit ρ is the radius of the symmetry group orbits. The metric field ω identically vanishes for solutions with no net rotation.

In order to construct stationary solutions it is useful to make an ansatz for the density function f on phase space, namely we assume that f depends only on the particle energy E and angular momentum L_z about the axis

$$\begin{equation} f\left(x,p\right) = K\Phi\left(E, L_z\right). \end{equation} \tag{ 3 }$$

Here K is a positive constant, which we call the amplitude, and Φ is a given function. In principle the amplitude could be incorporated into Φ but since it plays an important role in the numerical fixed point methods used in [2, 3, 46, 47] we separate it out. The quantities E and L_z are given by

$$\begin{eqnarray*} && E = -g\left(\partial_t, p^i\right) = e^{2 \nu} p^0 + \omega \left(\rho B\right)^2 e^{-2\nu} \left(p^3 - \omega p^0\right) \\ && L_z = g\left(\partial_\phi, p^i\right) = \left(\rho B\right)^2 e^{-2 \nu} \left( p^3 - \omega p^0\right), \end{eqnarray*}$$

which are constant along the geodesic flow. As a consequence, when f is given by the ansatz (3), the Vlasov equation is satisfied.

In order to work with the integral expressions in the energy momentum tensor it is useful to introduce new momentum variables as follows

$$\begin{equation} v^0 = e^\nu p^0, \quad v^1 = e^\mu p^1, \quad v^2 = e^\mu p^2, \quad v^3 = \rho B e^{-\nu} \left( p^3 - \omega p^0\right). \end{equation} \tag{ 4 }$$

In terms of these coordinates the energy momentum tensor can be written as

$$\begin{equation} T_{ij}\left(x\right) = \int_{P_x} p_i\left(v\right) p_j\left(v\right) f\left(x,v\right) \frac{\mathrm{d}^3 v}{\sqrt{1 + |v|^2}}. \end{equation} \tag{ 5 }$$

where the mass shell condition $g_{ij} p^i p^j = -1$, expressed in the new variables as $(v^0)^2 = 1 + |v|^2$, has been used to eliminate p₀. We take the positive root representing that all particles move forward in time. The particle angular momentum and energy can be expressed respectively as

$$\begin{align*} L_z = \rho B e^{-\nu} v^3 = : \rho s \end{align*}$$

and

$$\begin{align*} E = e^\nu \sqrt{1 + |v|^2} + \omega L_z = : h + \omega \rho s. \end{align*}$$

By plugging in the ansatz for the density function in the expression for the energy momentum tensor the components become integral expressions in the metric fields. The EV system then takes the form a non-linear system of integro-partial differential equations which we state for completeness:

$$\begin{align} \Delta \nu & = 4 \pi \left( \Phi_{00} + \Phi_{11} + \left( 1 + \left(\rho B\right)^2 e^{-4 \nu} \omega^2 \right)\Phi_{33} + 2 e^{-4 \nu} \omega \Phi_{03} \right) \end{align} \tag{ 6 }$$

$$\begin{align} \nonumber & \quad - \frac 1B \nabla B \cdot \nabla \nu + \frac 12 e^{-4\nu} \left(\rho B\right)^2 \nabla \omega \cdot \nabla \omega, \end{align}$$

$$\begin{align} \Delta B & = 8 \pi B \Phi_{11} - \frac 1\rho \nabla \rho \cdot \nabla B, \end{align} \tag{ 7 }$$

$$\begin{align} \Delta \mu & = - 4 \pi \left( \Phi_{00} + \Phi_{11} + \left(\left(\rho B\right)^2 e^{-4 \nu} \omega^2 - 1 \right)\Phi_{33} + 2 e^{-4 \nu} \omega \Phi_{03} \right) \end{align} \tag{ 8 }$$

$$\begin{align} \nonumber & \quad + \frac 1B \nabla B \cdot \nabla \nu - \nabla \nu \cdot \nabla \nu + \frac 1\rho \nabla \rho \cdot \nabla \mu + \frac 1\rho \nabla \rho \cdot \nabla \nu \end{align}$$

$$\begin{align} & \quad + \frac 14 e^{-4\nu} \left(\rho B\right)^2 \nabla \omega \cdot \nabla \omega, \end{align} \tag{ 9 }$$

$$\begin{align} \Delta \omega & = \frac{16 \pi}{\left(\rho B\right)^2} \left( \Phi_{03} + \left(\rho B\right)^2 \omega \Phi_{33} \right) - \frac 3B \nabla B \cdot \nabla \omega + 4 \nabla \nu \cdot \nabla \omega \end{align} \tag{ 10 }$$

$$\begin{align} & \quad - \frac 2\rho \nabla \rho \cdot \nabla \omega, \end{align} \tag{ 11 }$$

where $\Delta u: = \rho^{-1} \partial_\rho(\rho \partial_\rho u) + \partial_z\partial_z u$ and $\nabla u = (\partial_\rho u, \partial_z u)$. The matter components are given by

$$\begin{align} \nonumber \Phi_{00} & = e^{2\mu - 2 \nu} T_{00} \end{align} \tag{ 12 }$$

$$\begin{align} & \quad \,\,\,\, = \frac{2\pi}{B} e^{2\mu - 2\nu} \int_{e^\nu}^\infty \int_{-s_l}^{s_l} E\left(h,s\right)^2 \Phi\left(E\left(h,s\right), \rho s\right) \mathrm{d}s \mathrm{d}h, \end{align}$$

$$\begin{align} \nonumber \Phi_{11} & = T_{\rho \rho} + T_{zz} \end{align}$$

$$\begin{align} & = \frac{2\pi}{B^3} e^{2\mu + 2 \nu} \int_{e^\nu}^\infty \int_{-s_l}^{s_l} \left( s_l^2 - s^2 \right) \Phi\left(E\left(h,s\right), \rho s\right) \mathrm{d}s \mathrm{d}h, \end{align} \tag{ 13 }$$

$$\begin{align} \nonumber \Phi_{33} & = \left(\rho B\right)^{-2} e^{2\mu + 2 \nu} T_{\varphi \varphi} \end{align} \tag{ 14 }$$

$$\begin{align} & = \frac{2\pi}{B^3} e^{2\mu + 2 \nu} \int_{e^\nu}^\infty \int_{-s_l}^{s_l} s^2 \Phi\left(E\left(h,s\right), \rho s\right) \mathrm{d}s \mathrm{d}h, \end{align}$$

$$\begin{align} \nonumber \Phi_{03} & = e^{2\mu + 2 \nu} T_{0 \varphi} \end{align} \tag{ 15 }$$

$$\begin{align} & \quad \,\,\,\,= - 2\pi \rho B^{-1} e^{2\mu + 2 \nu}\int_{e^\nu}^\infty \int_{-s_l}^{s_l} s E\left(h,s\right) \Phi\left(E\left(h,s\right), \rho s\right) \mathrm{d}s \mathrm{d}h. \end{align}$$

Here

$$\begin{equation} s_l : = B e^{-\nu} \sqrt{e^{-2\nu} h^2 -1} . \end{equation} \tag{ 16 }$$

There are also auxiliary equations which we choose not to write out but refer to [2].

The system above must be complemented with boundary conditions. The following conditions guarantee that the solutions are asymptotically flat

$$\begin{align*} \nu, \mu, \omega \to 0, \quad \textrm{and} \quad B \to 1, \quad \textrm{as} \quad r = |\left(\rho,z\right)| \to \infty. \end{align*}$$

In addition we require that the metric is regular at the axis which implies that

$$\begin{equation} \nu\left(0, z\right) + \mu\left(0,z\right) = \text{ln} B\left(0,z\right) \end{equation} \tag{ 17 }$$

for all z in the solution domain.

Our numerical solutions may be characterized by several quantities. In this review we only use a subset of these quantities and we refer to [2, 3] for a more complete picture.

Two important properties of a solution are the total ADM mass $\mathcal{M}$ and the total angular momentum $\mathcal{J}$. We use the Komar expression for the mass and obtain

$$\begin{equation} \mathcal{M} : = 2\pi \int_{-\infty}^\infty \int_{0}^{\infty} g\left(\rho, z\right)\rho \, \mathrm{d} \rho \,\mathrm{d}z, \end{equation} \tag{ 18 }$$

where

$$\begin{align*} g\left(\rho ,z\right) = e^{2\mu - 2\nu}\rho B T_{00} + \rho B \left(T_{\rho \rho} + T_{zz}\right) + \frac{e^{2\mu + 2\nu}}{\rho B} T_{\varphi \varphi} - e^{2\mu - 2\nu} \rho B \omega^2 T_{\varphi \varphi}. \end{align*}$$

The mass plays an essential role in our iteration scheme. At each step of the iteration the ansatz function is renormalized such that the total mass is unity. We also have a Komar integral expression for the total angular momentum

$$\begin{equation} \mathcal J = -2\pi \int_{-\infty}^\infty \int_{0}^{\infty} e^{2\mu - 2\nu} B \left( T_{0\phi} + \omega T_{\phi \phi} \right)\rho \, \mathrm{d} \rho \,\mathrm{d}z. \end{equation} \tag{ 19 }$$

The total angular momentum is particularly important when we construct highly compact solutions.

An important measure of our solutions is the radius of support of the matter distribution. In spherical symmetry the ratio $2\mathcal M/R_0$, where R₀ is the radius of support in areal coordinates, is a measure of how relativistic a solution is. The ratio $2\mathcal M/R_0$ is often referred to as the compactness ratio.

If we express the metric (2) in spherical coordinates, the radial coordinate r is the isotropic radius. This can be related to the areal radial coordinate R through

$$\begin{align*} R = r \left( 1+ \mathcal M/\left(2 r\right)\right)^2. \end{align*}$$

In this paper we denote the areal radius of support by R₀ and the isotropic radius of support by r₀. For a spherically symmetric solution the radius of support can be determined from the cutoff energy E₀ (which is specified in the ansatz function) by the matching condition with a Schwarzschild exterior. The expression in terms of both the areal and isotropic coordinates is

$$\begin{align} E_0 = \sqrt{1 - 2\mathcal M /R_0} = \left(1 - \mathcal M /\left(2 r_0\right)\right)/ \left(1 + \mathcal M /\left(2 r_0\right)\right). \end{align} \tag{ 20 }$$

In spherical symmetry a black hole forms if the mass $\mathcal M$ becomes confined within a Schwarzschild radius of $\mathcal R = 2 \mathcal M$. The compactness $2\mathcal M /R_0$ is thus a useful characterization of the solution. There is no such well-defined criteria in axisymmetry. In our setting, a natural measure of the radius is the length of the axisymmetric Killing vector field which we denote $R_{\mathrm{circ}} : = \rho B e^{-\nu}$. This quantity provides a natural length scale for the solution, in particular when restricted to the reflection plane (z = 0) and evaluated near the boundary of the matter. For Vlasov matter, which typically has an extended atmosphere, it is useful to take the radius at which the compactness inside a cylinder of radius ρ is maximum. We define the compactness parameter $\Gamma : = \max_{\rho \in (0, \infty)} 2 m(\rho)/\bar{R}_{\mathrm{circ}}(\rho)$, where $\bar{R}_{\mathrm{circ}} : = (R_{\mathrm{circ}} )|_{z = 0}$ and where

$$\begin{align} m\left(\rho\right) : = 2\pi \int_{-\infty}^\infty \int_{0}^{\rho} g\left(\tilde{\rho}, z\right)\tilde\rho \, \mathrm{d} \tilde\rho \,\mathrm{d}z. \end{align} \tag{ 21 }$$

Note that $m(\rho) = \mathcal M$ when ρ exceeds the matter support. For the regular solutions we construct $\Gamma \in (0, 1)$.

One reason why this problem is important is obvious, it would give strong support that the numerical solutions obtained by this method are true solutions. Moreover, progress on the fixed point problem could give a new method to generate static solutions in cases where previous strategies have failed. We have not only in mind the axially symmetric Einstein–Vlasov system but also in the case of the Vlasov–Poisson system there is hope that this method can be used to construct new solutions, namely flat steady states. Presently there are only limited results in the literature about such solutions, see [37], and they are of interest as models of disk galaxies, see [18].

A further interesting aspect of this method is related to stability. Namely, as mentioned above, there are indications that solutions obtained by this algorithm are dynamically stable. This observation is particularly exciting in view of the open problem of non-linear stability for the spherically symmetric Einstein–Vlasov system. Hence, if a link between non-linear stability and the static solutions obtained by this algorithm can be established it would certainly be of great interest.

An exciting result in [2, 3] was the discovery that there exist regular stationary toroidal solutions of the Einstein–Vlasov system which admit ergoregions (see for example figures 3 and 4 in [3]). An ergoregion is typically associated with a Kerr black hole but not with a regular stationary solution. The definition of an ergoregion is that the Killing field $\partial_t$ which corresponds to the time translation symmetry becomes spacelike. In our parametrization of the metric it follows that an ergoregion is the set for which

$$\begin{equation} e^{2\nu}-\rho^2 B^2\omega^2 e^{-2\nu}\lt 0. \end{equation} \tag{ 24 }$$

In both [2, 3] in which ergoregions were observed, the following ansatz function was used

$$\begin{align*} \Phi\left(E,L_z\right) = \left(E_0-E\right)_+^0\left(L_z - L_0\right)_+^0. \end{align*}$$

This takes the form of equation (22) with k = 0, and $\psi(L_z)$ taking a similar polytropic form. Note however that the parameter L₀ represents a lower bound, resulting in a solution with net angular momentum. In order to obtain sufficiently relativistic solutions we construct a sequence of solutions and 'gently' decrease the parameter E₀, seeding the solver for each new solution with the previous solution in the sequence.

The solutions admitting ergoregions that we obtain are all highly relativistic and highly rotating, each satisfying the inequality

$$\begin{equation} |\mathcal{J}|>\mathcal{M}^2. \end{equation} \tag{ 25 }$$

We note that the Kerr metric for which (25) holds possesses a naked singularity. Hence, a regular stationary solution satisfying (25) is likely stable (with respect to axially symmetric perturbations) in view of the weak cosmic censorship conjecture. A highly relativistic solution for which (25) does not hold is most likely dynamically unstable; it will collapse to a Kerr black hole if it is perturbed, even in axisymmetry. This is consistent with the observation above in section 1.1 that there are indications that our numerical algorithm only converges to dynamically stable solutions, and hence we are unable to obtain solutions with ergoregions for which $|\mathcal{J}|\lt\mathcal{M}^2$. The above discussion is illustrated in figure 1 where convergence is lost when $E_0\lt 0.72$ for the sequence which does not satisfy the inequality (25). Perhaps a different numerical scheme will be more successful in finding solutions with ergoregions and with $|\mathcal{J}|\lt\mathcal{M}^2$.

Zoom In Zoom Out Reset image size

As alluded to in section 2.1, Shapiro and Teukolsky also looked for solutions containing ergoregions in [46]. They used a delta function based ansatz for both the particle energy and angular momentum ³ , and obtained highly rotating and relativistic solutions. At the limits of the resolution available at the time, Shapiro and Teukolsky were able to compute solutions corresponding to an E₀-value of about 0.745. In the step-function based ansatz studied in [2, 3], solutions with an ergoregion appear in the sequences of solutions studied at around $E_0 \approx 0.66$. Given this, one might speculate that the solutions obtained in [46] were not relativistic enough to contain ergoregions. It would be interesting to look for ergoregions in families of solutions obtained from other (non-step function) ansatzes.

The presence of ergoregions suggests that the sequence of stationary regular solutions may be approaching a family of rotating black hole solutions. Such a quasistationary transition to black hole solutions does not occur for the spherically symmetric Einstein–Vlasov system due to a Buchdahl type inequality, which for a body of mass $\mathcal M$ and radius $\mathcal R$ reads $2 \mathcal M/ \mathcal R \lt 8/9$, see [8]. Hence there is a gap and $2 \mathcal M/ \mathcal R$ cannot be arbitrary close to one. However, if one allows for charge a similar bound relating the mass, radius, and total charge is known [9], and in this case there is no gap; a quasistationary transition to a Reissner–Nordström black hole may thus be possible in spherical symmetry. Indeed, such a transition to the extremal Reissner–Nordström black hole has been shown by Meinel and Hütten [34] in the case of charged dust. Since charge is often considered as the poor man's angular momentum it is natural to ask if a similar transition is possible for rotating solutions. We note that in the case of disk solutions for dust, Meinel [31–33] has answered this question affirmatively by analytic arguments. More general cases have been investigated numerically by Ansorg et al [19, 26, 33].

As discussed in section 2.1 this question is investigated in [3]. The angular momentum parameter L₀ is tuned between highly rotating solutions, and more slowly rotating solutions. For each value a sequence of solutions with gradually decreasing E₀ values is constructed, as shown in figure 1. The high L₀ solution sequences approach the thin-ring limit, while the low L₀ solution sequences eventually terminate, presumably becoming unstable to gravitational collapse. Each such sequence, consisting of tens of solutions, is computationally expensive to complete. The most relativistic solutions we compute are with $L_0 = 0.806\,\,25$. As shown in figure 1, the ratio $\mathcal{J}/\mathcal{M}^2$ becomes very close to 1 as $E_0 \searrow 0.58$. Indeed, the compactness Γ also obtains its maximum value of roughly 0.8 at $E_0 = 0.58$, before decreasing again as the solution sequence bends towards the thin-ring limit (see [3] figure 2(b)). An extremal Kerr black hole has a compactness of $\Gamma = 1$.

Given the presumed instability of solutions with $\mathcal{J} \sim \mathcal{M}^2$, and the conjectured relationship between stability and the iterative solution method, it is not surprising that obtaining solutions very near the extremal Kerr solution, even when $\mathcal{J}$ is just greater than $\mathcal{M}^2$, is challenging. With more computational effort and resolution can one continue this bisection search and reach all the way to a black hole solution? It is not clear that the answer to this question is affirmative. The numerical method used in the fluid case, see [19, 26, 33], is different from the one used in [3] and it is not straightforward to adapt that method to the Einstein–Vlasov system. Hence, from this discussion an essential question arises: is a transition of stationary solutions to black hole solutions possible? Although such a transition has been established in the case of dust there is no a priori requirement that a quasistationary transition to black hole solutions must occur for solutions of the Einstein–Vlasov system.

Since solutions of the Einstein–Vlasov system share many properties of solutions of field theoretical models such as the Einstein–Dirac system, see [12], we find it to be a question of fundamental importance to understand whether or not such a transition is possible. Perhaps there is a maximum value of the compactness parameter Γ which is strictly below one as in the spherically symmetric case. Perhaps solutions necessarily approach the thin ring limit when E₀ is successively decreased independently of the value of L₀. This would be surprising, since as mentioned above, static spherically symmetric charged solutions exist which are arbitrary close to black hole solutions, see [9, 13]. It may be that a different numerical algorithm is required to answer these questions or it may be that higher resolution is sufficient.

The thin-ring limit discussed in the section above has interesting properties of its own. As shown in [3] the limiting members of solution sequences approaching this limit appear to display a local conical geometry around the matter, reminiscent of cosmic strings. Due to the near Dirac nature of the spatial matter distribution in this limit (i.e. it is highly focused on a ring), proving existence of such solutions may be analytically tractable. Indeed, in the spherically symmetric case it was utilized in [7] that the highly compact solutions approached an infinitely thin shell where the matter sources become Dirac distributions. This feature was essential for deriving upper and lower bounds on the compactness of the solutions, i.e. upper and lower bounds on $2m/r$. Perhaps a similar study is possible in the axially symmetric case by utilizing the Dirac nature of the thin ring solutions.

An important class of galaxies are disk galaxies such as the Milky Way. Given the prevalence of disk galaxies in the observable Universe, it is reasonable to ask whether such a morphology is represented in the space of stationary solutions of the Einstein–Vlasov system. There exist models of disk galaxies which are confined to the plane, see e.g. [18, 42], but it would be desirable to obtain three-dimensional solutions which are disk-like in the sense that the core of the density is as close to planar as possible.

An ansatz in which the z-component of the angular momentum is taken to have a Gaussian profile generates oblate spheroidal solutions. The ansatz is given by

$$\begin{equation} \psi\left(E,L_z\right) = \frac 1L_0 \text{exp}\left(L_z^2/L_0^2\right). \end{equation} \tag{ 26 }$$

This ansatz is symmetric in L_z and therefore generates non-rotating solutions of the Einstein–Vlasov system. A rotating version can be constructed by additionally imposing that all particles have L_z of the same sign. In the limit $L_0 \to \infty$ the ansatz becomes independent of L_z , thus generating a spherically symmetric spatial density. As L₀ is decreased, particles with higher angular momentum are more heavily weighted compared to those with low angular momentum. This yields a more flattened shape.

In [2] the flattening of such non-rotating oblate spheroids is studied. Parameters k and E₀ for the particle energy are fixed, while the parameter L₀ is gradually reduced from 10 to 1.5. Figure 2 shows the spatial density for a selection of solutions. Within the parameter range shown the spatial density distribution stretches to its most flattened form, while for parameters $L_0 \lt 1.5$ our numerical algorithm does not converge. We interpret this lack of convergence as a failure of the configuration to remain gravitationally bound. The final solution in the sequence does not have a very flattened form, though higher density contours appear far more disk-like than lower density contours—see also figures 9 and 10 in [2].

Zoom In Zoom Out Reset image size

In [2] models with net angular momentum are also studied. It is clear that rotation has an effect on the geometry, and hence the morphology, but it does not give rise to the very flattened form that we look for (see for example figure 9 of [2]). It should be stressed that a systematic investigation of this effect was not carried out in [2] but it would be surprising if merely an increase of the angular momentum would result in sufficiently flattened solutions.

In conclusion, our solutions do not satisfyingly resemble the extremely flattened galaxies which are observed in nature. We find it to be a central open problem to answer the question whether or not it is possible to obtain solutions of the Einstein–Vlasov system which resemble realistic disk galaxies. Perhaps a different choice of ansatz function will work out, or a different numerical method, or that it is simply not possible. It would be a great contribution to find out the answer to this question.

A different aspect is that it is often claimed that halos of dark matter may be crucial to the pronounced flattened shape of disk galaxies. It is a question of great importance to investigate if such halos surrounding disk-like solutions do improve the convergence of the numerical algorithm. If this turns out to be the case, it would give significant support to the established claim that disk galaxies are embedded in dark matter halos. To our knowledge this has not been investigated for the Einstein–Vlasov system and we find it to be an essential open problem.

In [2] we also construct spindle-like and toroidal solutions similar to what was done in [46]. The ansatz function we use for spindle solutions is given by

$$\begin{equation} \psi_\textrm{spindle}\left(L_z\right) = \begin{cases} \left(1 - Q |L_z| \right)^l, |L_z| < 1/Q \\ 0, |L_z| \unicode{x2A7E} 1/Q , \end{cases} \end{equation} \tag{ 27 }$$

and the ansatz function for the toroidal solutions takes the form

$$\begin{equation} \psi_\textrm{torus}\left(L_z\right) = \begin{cases} \left(|L_z| - L_0\right)^l, |L_z| > L_0 \\ 0, |L_z| \unicode{x2A7D} L_0. \\ \end{cases} \end{equation} \tag{ 28 }$$

Here $Q,L_0$ and l are parameters. As a remark we mention that the ansatz (27) was used for the Vlasov–Poission system in [18] to investigate the rotation velocities of stars in galaxies. In that setting, this ansatz gives rise to flat rotation curves similar to those found in observations.

Interestingly, astrophysical objects with spindle-like and toroidal structures have recently received attention by astrophysicists. In 2017 galaxy surveys revealed that prolate spindle-like galaxies are much more common than previously thought [49]. In 2019 it was announced that the VLA telescope had directly imaged a toroidal structure within an active galactic nucleus [24]. In view of the observational evidence of these types of objects we find that a more careful study of spindle solutions and toroidal solutions is motivated, where the features of the numerical solutions should be compared with the features of the observed objects.

In the context of galaxy morphology, let us also discuss composite objects which are obtained by combining different ansatz functions. Examples of composite astrophysical objects are numerous, and include disk galaxies with a central bulge, galaxies with dark matter halos, as well as ring galaxies. In the case of the Vlasov–Poisson system there are several results in the literature about composite solutions, but for the Einstein–Vlasov system they were first obtained in [2].

It turns out to be sensitive to combine ansatz functions. Our numerical algorithm does not converge for an arbitrary non-trivial linear combination of ansatz functions even if they individually give rise to solutions. For instance, we are not able to obtain convergence by combining an ansatz for a polytropic central bulge with a toroidal ansatz. The spindle solutions turned out to be useful in constructing composite objects. In [2] we used the following form of ansatz function for a composite solution

$$\begin{align*} \Phi\left(E, L_z\right) = K_s \Phi_\textrm{spindle}\left(E,L_z\right) + K_t \Phi_\textrm{torus}\left(E,L_z\right), \end{align*}$$

where K_s and K_t are amplitudes, or weights, of the two ansatz functions. An example of a solution inspired by Hoag's object [28] and obtained in this way is shown in figure 3. An open issue is to understand which combinations of ansatz functions give rise to composite solutions and to find out if these solutions resemble the morphology of composite objects found in nature.

Zoom In Zoom Out Reset image size