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

The purpose of this work is to review the status about stationary solutions of the axially symmetric Einstein-Vlasov system with a focus on open problems of both analytical and numerical nature. For the latter we emphasize that the code used to construct stationary solutions in \cite{Ames2016,Ames2019} is open source, see \cite{Ames2023joss}. In the analytical setting the open problems include establishing methods for proving existence of axisymmetric stationary solutions which are far from spherically symmetric, both in the general case and for certain special classes of solutions pointed out in the text. In the numerical setting there are intriguing properties of highly relativistic solutions that demand further attention, such as whether a sequence of such stationary solutions can approach a Kerr black hole, or if they necessarily approach the thin ring limit reminiscent of cosmic strings. The question of whether stationary solutions include states with thin-disk like morphologies as seen in many galaxies is also open. Finally, there are opportunities to extend this research to new settings such as the case of massless particles and coupled black hole-matter systems. We believe that some of the open problems highlighted here are of central importance for the understanding of nature.


Introduction
There are many reasons why studying the axially symmetric Einstein equations coupled to matter is of great interest.In 1955 Wheeler wrote in his seminal work on geons [48]: The simple toroidal geon forms the most elementary object of geon theory much as a simple circular orbit constitutes the first concept of planetary theory.But the simplest physics does not go in the geon case with the simplest mathematics.
Wheeler here expresses the wish to study his ideas on geons in the axially symmetric case but due to the complicated nature of the equations in that case he instead invented a spherically symmetric toy model which he called "Idealized spherically symmetric geons".Although the research field has made tremendous progress since this work was published, analytic results about Einstein's equations coupled to matter for isolated bodies are still very limited in the axially symmetric case.The development of numerical solutions has perhaps been even more dramatic since Wheeler's work from 1955 (where Wheeler in fact solved the Einstein equations numerically in the spherically symmetric case) but numerical studies about Einstein's equations coupled to matter often still concern toy models for the matter such as dust or a scalar field.There are of course exceptions, but at least for the Einstein-Vlasov system only a few investigations have been carried out.Shapiro and Teukolsky were pioneers and initiated an important study in the early 90's on the axisymmetric Einstein-Vlasov system.They considered both the dynamical case [41,43,42] and the stationary case [45,44].More recently this line of research has been continued in [2,3,4,5].Many questions have been answered in these works but several central questions remain open.
The purpose of the present contribution is to review the status about stationary solutions of the axially symmetric Einstein-Vlasov system, including both analytical and numerical work, and in particular to discuss open problems.We hope this discussion will stimulate activity in the field so that progress can be made on central problems.The open source code GECo [6], or further development of this code, may be useful in addressing some of the open problems discussed below.
In order to put the discussion below into context, in the remainder of this section we briefly discuss the spherically symmetric setting from a numerical perspective in Section 1.1, and in Section 1.2 include some equations and concepts from the axially symmetric literature.For a general background on the Einstein-Vlasov system we refer to [10,37,1].In Section 2 we briefly summarize the known results in the field.Further discussion of results in the context of open problems is presented in Section 3. Finally in Section 4 we make a few remarks on the evolutionary setting.

Comparison of numerical methods in the spherically symmetric setting
Spherically symmetric static solutions to the Einstein-Vlasov system are well-studied both analytically and numerically, cf.[10].There is a crucial difference in how solutions in the spherically symmetric case and solutions in the axially symmetric case are obtained.In the former case the equations can be solved as an initial value problem in the sense that given data at r = 0, a non-linear system of integro-differential equations can be solved straightforwardly, at least numerically.In this way a large variety of solutions can be constructed.An example is the class of highly relativistic multi-peak solutions that was discovered in [17].These solutions are believed to be unstable, but still they are easy to obtain numerically.
In the axially symmetric case the situation is different.The numerical methods that have been used in previous works [45,44,2,3] rely on an iteration procedure which must converge in order to produce a solution.Such iterative methods have also been employed in the spherically symmetric setting, where it is observed that convergence is only achieved for a certain range of the parameters.In fact there are indications that this convergent set corresponds with solutions which are dynamically stable.In particular, the highly relativistic multi-peak solutions mentioned above are not accessible with the iterative method.An analysis of the fixed point map discussed above would provide useful insight on this observation -see discussion in Section 3.1.For now, these numerical studies provide some evidence that the solutions obtained by the iteration procedure used in [2,3] are dynamically stable.

The axisymmetric Einstein-Vlasov System
In order to review the topic and to describe open problems in a meaningful way we find it useful to formulate the system of equations and to introduce some of the quantities that we use to characterize the solutions.

Equations
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 where the second equation is called the Vlasov equation.Here Ric(g), R(g) are the Ricci tensor and scalar of the metric g, Γ 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 ∈ M the energy momentum tensor takes the form 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, ρ, z, ϕ) as where the metric fields ν, µ, B, ω depend only on the coordinates ρ, z.Note that ρ = 0 is the axis of symmetry, and that (ρ, 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 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 [45,44,2,3] we separate it out.The quantities E and L z are given by 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 In terms of these coordinates the energy momentum tensor can be written as 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 0 .We take the positive root representing that all particles move forward in time.The particle angular momentum and energy can be expressed respectively as 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: where ∆u := ρ −1 ∂ ρ (ρ∂ ρ u) + ∂ z ∂ z u and ∇u = (∂ ρ u, ∂ z u).The matter components are given by Φ 00 = e 2µ−2ν T 00 Here 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 ν, µ, ω → 0, and B → 1, as r = |(ρ, z)| → ∞.
In addition we require that the metric is regular at the axis which implies that ν(0, z) + µ(0, z) = ln B(0, z) for all z in the solution domain.

Solution characteristics
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 M and the total angular momentum J .We use the Komar expression for the mass and obtain where g(ρ, z) = e 2µ−2ν ρBT 00 + ρB(T ρρ + T zz ) + e 2µ+2ν ρB T φφ − e 2µ−2ν ρBω 2 T φφ .
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 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 2M/R 0 , where R 0 is the radius of support in areal coordinates, is a measure of how relativistic a solution is.The ratio 2M/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 R = r(1 + M/(2r)) 2 .
In this paper we denote the areal radius of support by R 0 and the isotropic radius of support by r 0 .For a spherically symmetric solution the radius of support can be determined from the cutoff energy E 0 (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 In spherical symmetry a black hole forms if the mass M becomes confined within a Schwarzschild radius of R = 2M.The compactness 2M/R 0 is thus a useful characterization of the solution.There is no such welldefined 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 circ := ρBe −ν .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 Γ := max ρ∈(0,∞) 2m(ρ)/ Rcirc (ρ), where Rcirc := (R circ )| z=0 and where Note that m(ρ) = M when ρ exceeds the matter support.For the regular solutions we construct Γ ∈ (0, 1).

Brief survey of known results
In this section we briefly summarize the current knowledge regarding stationary axisymmetric solutions of the Einstien-Vlasov system.Relevant details are discussed below in Section 3 in the context of open problems.

Numerical
Stationary solutions to the axially symmetric Einstein-Vlasov system were for the first time constructed numerically in 1993 by Shapiro and Teukolsky in [45,44].In these works several sequences of stationary solutions are investigated with prolate and toroidal spatial density profiles, and include solutions with net total angular momentum.An exciting problem that these works left open was the question whether or not stationary solutions can be constructed which admit ergoregions.In numerical work together with Anders Logg, the authors answered this question affirmatively in [2].
In this work stationary solutions with a variety of morphologies are generated, including disk-like, spindle-like, toroidal, and solutions formed from a composition of ansatz functions (multi-species solutions).See [6] for the numerical code used to generate these solutions.
Further investigations of highly rotating and relativistic sequences of toroidal solutions was carried out in [3].The numerical results suggest two distinct possible limiting spacetimes depending on whether the total angular momentum is larger or less than the mass squared.Sequences for which J becomes less than M 2 eventually terminate, presumably becoming unstable to gravitational collapse.Sequences for which J stays larger than M 2 approach what we dub the "thin-ring limit", for which the solutions appear to have properties similar to those of cosmic strings, in particular a locally conical geometry about the string characterized by a deficit angle.A bisection search is carried out that tunes between these two extremes, and the results suggest a possible quasistationary transition to an extremal Kerr black hole for the critical solution sequence.In all of the numerical studies just cited the reduction scheme presented in Section 1.2.1 and an ansatz of the form Eq. ( 3) is used.The ansatz function Φ(E, L z ) is taken to have a product structure, i.e.

Φ(E, L
In many cases the ansatz function for the energy is taken to be the polytropic one ϕ(E) = (E 0 − E) k + , where (x) + = x if x ≥ 0 and (x) + = 0 if x < 0, and E 0 , k are parameters.The parameter E 0 has a natural interpretation as the cutoff energy for the particles, and its presence is important for obtaining solutions with compact support [36].Different choices of ψ result in different morphologies for the spatial density, among other properties.We note that if ψ is an even function the resulting solutions have zero net angular momentum.If instead the ansatz function ψ(L z ) vanishes for L z < 0 then all particles have angular momentum of the same sign.Solutions generated by such an ansatz have a net angular momentum and we call them rotating.

Analytic
In 2011 the existence of static axially symmetric solutions was shown for the first time by the second author together with Kunze and Rein in [15].In this case the total angular momentum vanishes.This result was generalized in 2014 to the stationary, rotating case [16], and extended to the Einstein-Vlasov-Maxwell case (in which the particles have charge) in 2020 [46].These analytic results rely on a perturbation argument with the consequence that the constructed solutions can only be guaranteed to deviate slightly from spherically symmetric solutions and to have small total angular momentum.Essentially the same system of equations as presented in Section 1.2 is used in the works [15,16], with the technical difference that perturbation parameters corresponding to the angular momentum and relativistic nature of the solutions are introduced.The authors also make use of an equation for ξ := ν + µ.
In addition to the existence results discussed above, there is a recent result by Jabiri [28] which uses a related method to construct stationary solutions.The solutions are obtained as bifurcations from the Kerr spacetime and thus gives a generalization to the case where a black hole is surrounded by Vlasov matter.
We mention also that in the static and cylindrically symmetric setting (in which the solution has unbounded support in one-direction), the existence of solutions and compact support in the radial direction is proved in [26].

Discussion and open problems
In this section we discuss properties of the solutions that were constructed analytically in [15,16] and numerically in [2,3] with the aim of formulating open problems, both analytical and numerical.The first part concerns existence of analytic solutions and ideas on how to extend previous methods to more general ones.The second part concerns highly relativistic compact solutions where ergoregions, the quasistationary transition to black holes and the thin ring limit are discussed.The third part concerns models of galaxies where the aim is to find solutions with the morphology of galaxies as observed in nature.

Existence of far-from spherically symmetric axially symmetric solutions
A common feature of the existence results is that the solutions are obtained as perturbations of known solutions.Hence solutions which are far-from spherically-symmetric in the sense of spatial density and net angular momentum -in particular, solutions containing ergoregions-are not covered by these results.Accordingly, there is extensive room for analytic progress on the existence of stationary solutions.One way to attack this problem is to mimic the numerical approach, i.e. to analytically investigate the iteration scheme of the numerical algorithm and show convergence in some domain of parameter space.In fact there are several reasons why this may be of importance, some of which are discussed below.To this end, let us describe the problem in the simplified case of the spherically symmetric Vlasov-Possion system.This is the Newtonian analogue of the Einstein-Vlasov system.Although existence of static solutions to the Vlasov-Poisson system is well understood, it is nevertheless an interesting question if existence can be shown via the method suggested by the numerical algorithm.
The Vlasov-Poisson system reads This system has the same general structure as the Einstein-Vlasov system but it is much less involved.The aim is to prove the existence of static solutions by the following strategy.Let an ansatz Φ(E, L) for the particle distribution be given, i.e., Here the particle energy E and the modulus of angular momentum L are given by For a given spatial density ρ = ρ(x) of finite mass M > 0 and compact support we define its induced gravitational potential by If U ρ is substituted into the ansatz (23) we obtain a new spatial density We now define an amplitude so that the new spatial density K(ρ) ρ again has mass M .Let us consider the map T : D → D defined by T (ρ) = K(ρ)ρ, where Here D is a suitable domain.If this map has a fixed point ρ * then ) is a static solution of the Vlasov-Poisson system where f * is given by the new ansatz Hence the exact problem that is solved is a priori not known, it is determined once K(ρ * ) is known.
It is an open problem to show that the map T has a fixed point.Similar fixed point problems, based on the same strategy, can be formulated for several related systems; e.g. the spherically symmetric Einstein-Vlasov system and the axially symmetric Einstein-Vlasov system.It is an open problem in each case to show the existence of a fixed point.

Motivations
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, cf.[35], and they are of interest as models of disk galaxies, cf.[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.

Highly relativistic solutions
In this section we discuss properties of highly relativistic solutions, by which we mean that the compactness ratio Γ is large.Such solutions contain ergoregions resembling black hole solutions.A central open problem is whether or not a quasistationary transition to an extremal black hole is possible.This is discussed below.

Ergoregions
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 ∂ 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 In both [2] and [3] in which ergoregions were observed, the following ansatz function was used This takes the form of Eq. ( 22) with k = 0, and ψ(L z ) taking a similar polytropic form.Note however that the parameter L 0 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 0 , 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 We note that the Kerr metric for which (25) holds possesses a naked singularity.Hence, a 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 |J | < M 2 .The above discussion is illustrated in Figure 1 where convergence is lost when E 0 < 0.72 for the sequence which does not satisfy the inequality (25).Perhaps a different nu- As alluded to in Section 2.1, Shapiro and Teukolsky also looked for solutions containing ergoregions in [44].They used a delta function based ansatz for both the particle energy and angular momentum,1 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 0 -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 ≈ 0.66.Given this, one might speculate that the solutions obtained in [44] 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.

Quasistationary transition to black holes
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 M and radius R reads 2M/R < 8/9, cf.[8].Hence there is a gap and 2M/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 [32] 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 [30,29,31] has answered this question affirmatively by analytic arguments.More general cases have been investigated numerically by Ansorg et al. [31,19,25].
As discussed in Section 2.1 this question is investigated in [3].The angular momentum parameter L 0 is tuned between highly rotating solutions, and more slowly rotating solutions.For each value a sequence of solutions with gradually decreasing E 0 values is constructed, as shown in Figure 1.The high L 0 solution sequences approach the thin-ring limit, while the low L 0 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.80625.As shown in Figure 1, the ratio J /M 2 becomes very close to 1 as E 0 ↘ 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 2b).An extremal Kerr black hole has a compactness of Γ = 1.
Given the presumed instability of solutions with J ∼ 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 J is just greater than 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, cf.[31,19,25], 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, cf.[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 0 is successively decreased independently of the value of L 0 .This would be surprising, since as mentioned above, static spherically symmetric charged solutions exist which are arbitrary close to black hole solutions, cf.[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.

Existence of a thin-ring limit to sequences of rotating toroidal solutions
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, 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.

Models of galaxies
One goal in [2] was to construct solutions which resemble galaxies as observed in nature.While the solutions obtained in [2] make no attempt to be full-fledged galaxy models, it is still of interest to compare their properties with real galaxies in order to better understand which properties of galaxies are captured by the fundamental axially symmetric gravitational physics.We point out that models of galaxies are most often obtained within the Newtonian framework, cf.[22] for models based on the Vlasov-Poisson system.One reason why relativistic models could be of interest in this context is that net rotation has no effect on the gravitational field in Newtonian gravity, yet it noticeably alters the geometry in general relativity, even for solutions which are not very relativistic.

Solutions with disk-like morphology
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, cf.e.g.[40,18], 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 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 → ∞ the ansatz becomes independent of L z , thus generating a spherically symmetric spatial density.As L 0 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 0 for the particle energy are fixed, while the parameter L 0 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 < 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] 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 disklike 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.

Spindle-like, toriodal and composite solutions
In [2] we also construct spindle-like and toroidal solutions similar to what was done in [44].The ansatz function we use for spindle solutions is given by and the ansatz function for the toroidal solutions takes the form 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 recieved attention by astrophysicists.In 2017 galaxy surveys revealed that prolate spindle-like galaxies are much more common than previously thought [47].In 2019 it was announced that the VLA telescope had directly imaged a toroidal structure within an active galactic nucleus [23].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 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 [27] and obtained in this way  28) ansatzes, and inspired by Hoag's object [27] 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.

The massless Einstein-Vlasov system
A question related to the opening of this review and to the existence of compact solutions is whether or not massless axially symmetric solutions exist.In the spherically symmetric case massless solutions exist if the solutions are sufficiently compact, cf.[14,11].Indeed, highly compact spherically symmetric massless solutions have the property that the density function vanishes at some finite radius R and the solution can be glued to a vacuum Schwarzschild solution at r = R. Without the gluing procedure the density function f will eventually become strictly positive at a sufficiently large radius, cf.[14].This particular feature is only present in the massless case.Hence, the crucial question is whether or not massless axially symmetric solutions can be constructed such that the density function vanishes outside a compact set.If so, it would be possible to glue the solution to a vacuum solution -although no explicit vacuum solution exists as in the spherically symmetric case.
We find this question exciting in view of the quote by Wheeler given in the introduction.Wheeler wanted to find and investigate the properties of regular solutions of the Einstein-Maxwell system.He named such solutions geons.In spherical symmetry the only solution to the Einstein-Maxwell system is the Reissner-Nordström solution which is not regular.Hence, spherically symmetric geons do not exist strictly speaking, although Wheeler did introduce the concept of idealized spherically symmetric geons [48].His ultimate wish was to study the axially symmetric case.Since an electromagnetic field can be modeled as a photon gas, a possible way to obtain understanding of the properties of solutions of the Einstein-Maxwell system is to study solutions of the massless Einstein-Vlasov system.In the spherically symmetric case highly relativistic solutions of the two systems are indeed very similar, cf.[14].Hence, the question whether or not axially symmetric solutions of the massless Einstein-Vlasov system exist is closely related to the existence of Wheeler's original conception of geons.

Regular solutions about central black holes
An interesting direction in which to extend the study of regular axisymmetric stationary solutions of the Einstein-Vlasov system, is to consider solutions with a central black hole.This imposes inner boundary conditions at the black hole horizon.One motivation is to better understand coupled matterblack hole systems and the effects that matter can have on central black holes.Another reason such a line of investigation is of interest is to provide foundational results for current research on accretion disks.To date and to the author's knowledge, such studies make use of a fixed black hole background solution, as in for example [33], and focus on the plasma physics in this extreme gravitational environment.In the purely gravitational setting, Rioseco and Sarbach have several works studying the dynamics of Vlasov matter on black hole backgrounds, for example [38,39].
The self-gravitating case is much more challenging.In the spherically symmetric setting, Andréasson proves the existence of solutions with a central Schwarzschild black hole in the massless case [11], and in [34] Rein establishes a such existence in the massive case.As mentioned in Section 2.2, Jabiri has proved existence of self-gravitating solutions using a perturbation argument about the Kerr spacetimes [28].This results represents the first proof of existence for axially symmetric self-gravitating solutions about a central black hole.By the nature of the proof however, the matter must be small (close to vacuum), and it remains an open problem to find general axially symmetric self-gravitating solutions of the Einstein-Vlasov system with a central black hole.
A sensible starting place is to probe this problem numerically.One approach would be to modify the code used in [2,3] and documented in [6].Related studies exist in the case of uniformly rotating fluids about black holes, see for example [21,20].

Extension to the Einstein-Vlasov-Maxwell system
Another interesting direction of research, and possible extension of the code [6], is the addition of charge.Allowing the particles to be charged (in addi-tion to having mass) and interact additionally through the electromagnetic field leads to the Einstein-Vlasov-Maxwell system.This system was studied numerically in 2009 by Andréasson and coauthors [13] in the spherically symmetric setting, and in 2020 Thaller used perturbation methods to prove existence of solutions in the axisymmetric setting [46].This result is similar to the work of [16] except that the reference solution is a charged solution of the spherically symmetric Vlasov-Poisson system, and as a result the particle charge is not restricted to be small.The total angular momentum and the strength of the gravitational field are still restricted, and an understanding of general solutions to the axially symmetric Einstein-Vlasov-Maxwell system is lacking.

Remarks on the evolution problem
We have reviewed the status of stationary solutions of the axially symmetric Einstein-Vlasov system.The question of their stability, and more generally, the fate of any axially symmetric initial data requires an evolution code.Shapiro and Teukolsky were pioneers in developing such a code using the particle in cell (PIC) method.More recently, several results on the evolution problem have been carried out, cf.[49,24,4,5].We have in this work left out a discussion about these results.One reason being that, to our knowledge, no evolution code is open source and it is an extensive work to develop such a code.Needless to say there is immense room for improvements and developments of this topic.We will not enter into this here but let us at least bring up one open problem which is related to the discussion above.The highly compact solutions that were constructed in [3] require very high resolution.It is an outstanding problem to determine whether or not these solutions are stable.The high resolution needed to construct the stationary solutions is clearly a severe obstacle.How can the high resolution needed for the stationary solution be carried over to the evolution problem to make a simulation practically feasible?