text:
# Reaction-diffusion dynamics: confrontation between theory and experiment in a microfluidic reactor

## Abstract

We confront, quantitatively, the theoretical description of the reaction-diffusion of a second order reaction to experiment. The reaction at work is $`\text{Ca}^{2+}`$/CaGreen, and the reactor is a T-shaped microchannel, 10 $`\mu `$m deep, 200 $`\mu `$m wide, and 2 cm long. The experimental measurements are compared with the two-dimensional numerical simulation of the reaction-diffusion equations. We find good agreement between theory and experiment. From this study, one may propose a method of measurement of various quantities, such as the kinetic rate of the reaction, in conditions yet inaccessible to conventional methods.

Chemical reaction processes taking place in spatially extended reactors with no active mixing are described by reaction-diffusion equations. In general, these equations are non linear and they possess a rich variety of solutions, giving rise to many different spatio-temporal patterns. Reaction-diffusion equations have been known for more than a century but the first elementary solutions were calculated much more recently by Gálfi and Rácz Gálfi and Rácz (1988). Their analytical solutions were obtained for a reaction of second order in an infinite domain with stepwise initial conditions, in the asymptotic limit of long times. Under these conditions, exact scaling laws were derived showing, for instance, that the position of the maximum of the production rate grows as $`t^{1/2}`$ while the width of the reaction production zone grows as $`t^{1/6}`$. Furthermore, the approach to the asymptotic state was studied by Taitelbaum et al. Taitelbaum et al. (1992) and they observed a slow approach to the asymptotic state, which was modeled by using perturbation analysis and could not be described by the long time asymptotics.

At the experimental level, one is often confronted with the additional presence of advection in the fluid. For this reason, reactions are often studied in small capillary tubes where no large scale flow is present. Such systems have been used Koo and Kopelman (1991); Park et al. (2001) to study a reaction process with initially separated reactants, validating the power law predictions of Gálfi and Rácz at asymptotically long times. Furthermore, recent experiments in porous media Léger et al. (1999) have shown very good agreement between theory and experiment in the asymptotic regime, in the case of one diffusing reactant into a solid substrate.

Here, we address the more general problem of two diffusing, initially separated, reactants in a regime that spans both the initial transients and the asymptotic state. In contrast with previously published work Koo and Kopelman (1991); Léger et al. (1999), we turn our attention to the reaction product and its concentration, instead of the production rate. This analysis takes into account the diffusion of the product as well as the reagents. We compare numerical and experimental measures of the reaction region, beyond the scaling of the different quantities, to study the shape of the concentration profiles.

Microfluidic circuits are used to achieve low Reynolds numbers ($`Re1`$), thus removing any effects of advection without resorting to gels or other media. We believe that the detailed comparison of the shape of the concentration profiles between simulations and experiments allows one to assess to what extent, from the dynamical viewpoint, modeling a reaction by a kinetics of a prescribed order is accurate. This assessment is useful in view of modeling more complex behaviors in laminar, chaotic, or turbulent flows. In this regard, microsystems offer the opportunity to accurately compare theory and experiment since simple geometries can be achieved and, owing to the low Reynolds number at hand, the hydrodynamics is well controlled. Diffusion is the only transport mechanism leading the reagents to react together.

The channel geometry we consider here precisely allows the study of a reaction-diffusion process with stepwise initial conditions. The micro-reactor is T-shaped, similar to that used in Ref. Kamholz et al. (1999). In contrast with these authors, we use shallow channels which favor two-dimensionality as discussed below. We will further infer from the present work a practical method for measuring the kinetics of chemical reactions in the sub-millisecond range, a regime that is out of the reach of traditional stopped-flow techniques Cussler (1997).

Theoretical background: Let us consider a second order irreversible reaction of the type $`A+BC`$, taking place in a two-dimensional space, within a flow driven at uniform velocity U, as shown in Fig. 1.

In the stationary state, after neglecting the derivative along $`X`$ in the Laplacian operator, the reaction-diffusion equations read :

$`U{\displaystyle \frac{A}{X}}`$ $`=`$ $`D_a{\displaystyle \frac{^2A}{Y^2}}kAB`$
$`U{\displaystyle \frac{B}{X}}`$ $`=`$ $`D_b{\displaystyle \frac{^2B}{Y^2}}kAB`$ (1)
where $`A`$ and $`B`$ are the concentrations of the two reactants A and B, $`D_a`$ and $`D_b`$ are the corresponding diffusion coefficients, and $`k`$ is the chemical rate constant. These equations are identical to the typical reaction-diffusion equations, if we replace variable $`X`$ by $`Ut`$, where $`t`$ is the time. This system must be supplemented by an equation governing the evolution of product C, given by :

$$U\frac{C}{X}=D_c\frac{^2C}{Y^2}+kAB$$
(2)
where $`C`$ is the concentration of product C, assumed inert, and $`D_c`$ is the corresponding diffusion coefficient. As initial conditions, we take the reagents to be separated with constant densities for both $`Y<0`$ ($`A=A_0`$ and $`B=0`$) and $`Y>0`$ ($`A=0`$ and $`B=B_0`$). We thus have here a six dimensional parameter space which consists of the three diffusion coefficients of the two reagents and the product, the initial concentration levels of the reagents, and the reaction rate coefficient.

The system of Eqs. 1 can be non-dimensionalized by introducing a characteristic length scale $`\mathrm{}`$, similar to the one introduced in Ref. Gálfi and Rácz (1988)

$$\mathrm{}^2=\frac{\sqrt{D_aD_b}}{k\sqrt{A_0B_0}}.$$
(3)
This length is used to non-dimensionalize $`X`$ and $`Y`$ as:

$$x=\frac{XD}{\mathrm{}^2U},y=\frac{Y}{\mathrm{}}$$

where $`D=\sqrt{D_aD_b}`$. Using these new dimensionless variables and introducing the dimensionless concentrations

$$a=\frac{A}{A_0}\text{and}b=\frac{B}{B_0},$$

the governing equations for the problem become:

$`{\displaystyle \frac{a}{x}}`$ $`=`$ $`\chi {\displaystyle \frac{^2a}{y^2}}{\displaystyle \frac{ab}{\beta }}`$
$`{\displaystyle \frac{b}{x}}`$ $`=`$ $`{\displaystyle \frac{1}{\chi }}{\displaystyle \frac{^2b}{y^2}}\beta ab`$ (4)
in which we define

$$\chi =\sqrt{D_a/D_b},\beta =\sqrt{A_0/B_0}.$$
(5)
This system must also be supplemented by an equation prescribing the evolution of the product concentration. Although the problem can be solved in the general case, we restrict ourselves to the particular case where the third diffusion coefficient, $`D_c`$, is equal to $`D_a`$ (assumed to be the smaller than $`D_b`$). The equation for the evolution of the product concentration thus reads, in dimensionless form:

$$\frac{c}{x}=\frac{1}{\chi }\frac{^2c}{y^2}+ab$$
(6)
with $`c=C/\sqrt{A_0B_0}`$. Note that the reaction rate has disappeared in Eqs. 4 and 6. This means that chemical kinetics does not play any dynamical role: changing $`k`$ rescales the variables of the problem but does not change the structure of the solution. A consequence is that the number of parameters is reduced to two, meaning that the full set of solutions can be represented in a two-dimensional parameter space. A knowledge of both parameters $`\chi `$ and $`\beta `$ is necessary to compare simulations and experiments, if one is interested in both the transient and asymptotic regimes.

Description of the experiment: In our experiments, we used a T-shaped microreactor, whose geometry is represented in Fig 2. This geometry is similar to the one introduced in Ref. Kamholz et al. (1999).

The channels were chemically etched in a glass wafer, which was then anodically bonded to a silicon wafer to create a rectangular cross-section. Throughout the experiments discussed here, the channel depth was 10 $`\mu `$m and the width was 200 $`\mu `$m. The lengths of the T arms were 1 cm while the main channel was 2 cm long. Two reactant species are introduced through the two top branches of the T and they react as they flow side by side in the main channel. We study the reaction of Calcium ions as reactant A and Calcium Green (CaGreen) from Molecular Probes pro  as reactant B (as used, e.g. in Monson and Kopelman (2000)). $`\text{Ca}^{2+}`$ ions bind with the tracer Calcium Green, significantly increasing its fluorescence rea . Epifluorescence microscopy is used to measure the fluorescence emission, yielding spatial images such as Fig. 3. The flows are driven by syringe pumps which also provide measurements of the flow-rates in each separate branch; those ranged between 30 and 100 nl/mn. Owing to the aspect ratio of the micro-channel, the velocity profile is parabolic across the channel depth and is essentially uniform along $`Y`$ direction. Furthermore, the low Reynolds number ensures that the fluid flow is invariant in the $`x`$-direction after some entrance length on the order of the channel thickness.

In the conditions in which we operated (i.e. low flow-rates and thin channels) the tracer concentration can be considered as uniform across the channel depth. As a consequence, the concentration fields can be treated as two-dimensional and Eqs. 4 and 6 can be used to describe the system. We tested the two-dimensionality property by performing a series of diffusion experiments (with no reaction) of fluorescein in water. A fluorescein solution and pure water are driven from the two top arms respectively. In the main channel, they form a diffusing interface which progressively broadens downstream. In thick channels, the broadening was found to display the three-dimensional effects discussed in Refs. Kamholz et al. (1999) and Ismagilov et al. (2000), where the diffusion region scales as $`x^{1/3}`$. In the 10 $`\mu `$m channels, 2D diffusion theory to agrees well with experiment, with a scaling of the diffusion width as $`x^{1/2}`$.

The comparison between the 2D theoretical model and the experiment is further justified by the following reasoning: one may consider that the reacting zone is homogeneous across the channel depth in the major part of the region we analyze, for sufficiently thin channels. This argument relies on estimations of the diffusive time across the channel thickness and on the measurements with fluorescein discussed above. We are thus working mostly in the so-called Taylor-Aris regime where the parabolic velocity profile enhances the apparent diffusion coefficient in the flow direction Pagitsas et al. (1986). One may show however that the streamwise diffusion term is still negligible compared to the transverse term, since the streamwise gradients are very weak. These two features justify the use of the model of Eqs. 4 to compare with our experiment.

The diffusion coefficient of Calcium ions is equal to $`D_{\text{Ca}^{2+}}=7.9\times 10^6`$cm2/s Cussler (1997). The diffusion coefficient of the CaGreen was unknown so we measured it using the method described above, replacing the fluorescein solution with a Calcium Green solution: We checked that the width of the diffused zone increases as the downstream distance raised to the power 1/2, as expected from two-dimensional theory. By fitting the measurements of the width, we obtained an estimate for the Calcium Green diffusion coefficient equal to $`3.2\pm 0.3\times 10^6`$ cm2/s.

Throughout the experiments, we kept the concentration of $`\text{Ca}^{2+}`$ constant at \[$`\text{Ca}^{2+}`$\]=1 mM, while the concentration of CaGreen was varied in the range $`0.1<\text{[CaGreen]}<0.5`$ mM. Therefore, the value of $`\beta `$ varied between $`1.4`$ and $`3.2`$, while $`\chi `$ was fixed at $`\chi =1.57`$.

Comparison between theory and experiment: A typical concentration field is shown in Fig 3. The black region corresponds to the Calcium solution, the grey one to the CaGreen solution, and the brighter region in the center to the bound Ca-CaGreen complex. The background fluorescence, visible in the lower part of the channel, is due to the fluorescence of unbound CaGreen. The asymmetry in the reaction zone is due to a different rate of diffusion of the reagents, because of their different initial concentrations and different diffusion coefficients. This can be physically explained by noting that the species with the higher diffusion coefficient will diffuse further; furthermore the species with the larger initial concentration will diffuse faster due to the large gradient it experiences. In other words, $`\beta `$ and $`\chi `$ both determine the shape of the reaction zone.

We were not able to obtain reliable information on the kinetic coefficient $`k`$ of the $`\text{Ca}^{2+}`$-Calcium Green reaction. We thus took $`k`$ as a free fitting parameter, and adjusted it so as to optimize the agreement between theory and experiment. The best fits were obtained for $`k=1.0\pm 0.47\times 10^6`$ dm3/(mol s), or a length scale $`\mathrm{}=0.84`$ $`\mu `$m for $`\beta =\sqrt{2}`$ and $`\mathrm{}=1.0`$ $`\mu `$m for $`\beta =2`$. These values of $`k`$ and $`\mathrm{}`$ are used throughout this letter in comparing the experiments and numerics.

Figure 4 represents the concentration profiles of the Ca-CaGreen complex (dots) measured at various locations downstream, for $`\beta =\sqrt{2}`$ and 2. The two locations we consider here are $`X=278`$ and 1128 $`\mu `$m (for $`\beta =\sqrt{2}`$ these correspond to $`x=32`$ and $`216`$; for $`\beta =2`$, $`x=23`$ and $`153`$). The plots we show are representative of the evolution of the concentration profiles along the downstream axis. The profiles are asymmetric with respect to $`Y=0`$; in particular the maximum of the concentration is not located at the channel center but is displaced towards the side of the CaGreen solution. For increasing $`X`$, the reaction zone broadens. The solid lines on Fig. 4 represent the value of $`c`$ obtained by numerical solutions of Eqs. 4 and 6, by a standard one dimensional finite-difference scheme with time stepping on a personal computer, using the values of the parameters previously established. The offset on the right hand side of the experimental curves is due to the background fluorescence of the unbound CaGreen (reactant $`a`$). We match this offset in the numerical curves by plotting $`ha+c`$ rather than just $`c`$, where $`h`$ a normalization constant that gives the correct offset. This does not affect the fitting parameters but matches the experimentally measured quantity, i.e. the total fluorescence. We observe good agreement between theory and experiment.

A quantitative measure of the agreement between the simulations and experiments is obtained by looking at some representative quantities for the profiles. In particular, we compare the downstream evolution of three quantities: The position of the maximum, the width at half-height, and the value of the maximum (Fig. 5). These evolutions do not follow simple power laws in the full range of $`x`$ we consider. However, for $`x30`$ dimensionless units, the power laws are better defined. The useful scaling range in the experiments is further limited (at $`X1200`$$`\mu `$m) by lateral wall effects, thus limiting the range over which a power law can be observed; although the scaling range is too small to obtain a clear exponent, our measurements are consistent with the long time limit analytical theory of Gálfi and Rácz Gálfi and Rácz (1988). The maximum concentration shows a saturation effect at large distances, i.e above $`x300`$. These measurements can be compared with the numerics, plotted as full lines. Fig 5 shows there is good agreement between the theory and the experiment. A similar agreement is obtained different values of $`\beta `$ and flowrates $`U`$ we considered.

In conclusion, we demonstrate that two-dimensional reaction-diffusion equations with stepwise initial conditions define not only a theoretical problem, but also a physical situation which can be achieved experimentally. Microfluidic reactors are particularly convenient in order to produce such physical situations. Conversely, 2D reaction diffusion equations can be used to determine, from experimental measurements, the relative concentration of the reactants (micro-titration), the diffusion coefficients or the chemical kinetics. In the present case, we demonstrate a determination of the kinetic coefficient $`k`$ in order to match measured physical quantities (position, width and concentration of the reactant) with numerical predictions. The value of $`k`$ thus obtained for our system corresponds to a reaction half-time of 1 ms, which is faster than the time scales that can be resolved in a stopped-flow apparatus. This indicates that the present approach can be used for analyzing fast chemical kinetics.

The authors acknowledge technical help from Bertrand Lambollez, Jean Rossier and Thibault Colin, and useful discussions with Martin Bazant.


format: text

text:
# DYNAMICS OF A MAGNETIC MONOPOLE IN MATTER, MAXWELL EQUATIONS IN DYONIC MATTER AND DETECTION OF ELECTRIC DIPOLE MOMENTS

(  For a monopole, the analogue of the Lorentz equation in matter is shown to be $`𝐟=g(𝐇𝐯\times 𝐃)`$. Dual-symmetric Maxwell equations, for matter containing hidden magnetic charges in addition to electric ones, are given. They apply as well to ordinary matter if the particles possess T-violating electric dipole moments. Two schemes of experiments for the detection of such moments in macroscopic pieces of matter are proposed. PACS: 14.80.Hv 03.50.De 11.30.Er  )

1. INTRODUCTION

The question of which classical macroscopic fields exert a force on a magnetic monopole of charge $`g`$ in matter is still controversial . For the static force, the formula

$$𝐟=g𝐇,$$
$`(1)`$
instead of $`𝐟=g𝐁`$, is generally accepted. However, for the velocity dependent force, there is no consensus between $`𝐟=g𝐯\times 𝐄`$ and $`𝐟=g𝐯\times 𝐃`$ (we use rationalized equations with $`c=\epsilon _0=\mu _0=1`$). A more general problem is to generalize the macroscopic Maxwell equations to the dual-symmetric matter. The atoms or molecules of such a matter would be made not only of electrically, but also of magnetically charged particles. Thus they can possess

* electric dipole moments coming from decentered electric charges as well as spinning magnetic charges
* magnetic moments comming from spinning electric charges as well as decentered magnetic charges.

After a rederivation and a discussion of Eq.(1), we will present below a consistent solution for the velocity-dependent force and the dual-symmetric Maxwell equations in matter, using simple physical arguments. We will consider only isotropic matter and assume that its electric and magnetic polarizations $`𝐏`$ and $`𝐌`$ are linear in $`𝐃`$ and $`𝐁`$ (or $`𝐄`$ and $`𝐇`$). It will appear that our equations can also take into account the electric dipole moments (e.d.m.) of the ordinary fermions generated by T-violating interactions, and we will propose two kinds of possible measurements of the e.d.m. in macroscopic matter.

2. STATIC FORCE ON A MONOPOLE IN MATTER.

If the force acting on a monopole in matter were $`𝐟=g𝐁`$, a monopole following a closed magnetic line of a permanent magnet could gain energy at each turn, providing a perpetual motion of the first species. This is an argument for chosing $`𝐟=g𝐇`$, the curl of which is zero for a static system.

One might object that the monopole can gain energy at each turn at the expense of the magnetic energy stored by the magnet and will eventually erase the magnetization of the metal. This is indeed what happens when a monopole is circulating through a superconducting loop : the varying flux of the monopole field through the loop produces a counter-electromotive force which damps the supercurrent. However, in the case of a ferroelectric annulus, the magnetized state has the lowest energy and the annulus cannot yield any energy to the monopole.

Another argument for (1) comes from the (gedanken) following experiment : Let us measure the force on a magnetic charge immersed in a ferrofluid. The latter is a practical realization of a liquid magnetizable matter. No static frictional force can perturb the measurement. We protect the monopole from the fluid by a waterproof box. This should not change the result ; anyway the physical monopole is probably dressed by a swarm of ordinary particles. In the absence of the monopole, we denote by $`𝐁\mu 𝐇`$ the field outside the box and by $`𝐁_{box}`$ = $`𝐇_{box}`$ the field inside the box. The fields coming from the monopole will be denoted by a prime. Let us consider two shapes of box (Fig.1) :

a) the box is elongated parallel to $`𝐁`$. Then $`𝐁_{box}`$ = $`𝐇`$ and the measured force is $`𝐟=g𝐇`$.

b) the box is flattened perpendicular to $`𝐁`$. Then $`𝐁_{box}`$ = $`𝐁`$. The force acting on the pole is $`𝐟_1=g𝐁`$, which is different from case a). On the other hand, in front of the box the total field $`𝐁^{tot}=𝐁+𝐁^{}`$ is larger than behind. The magnetic grains of the ferrofluid are therefore attracted toward this region and build a hydrostatic pressure which pushes the box backwards. Quantitatively, the force acting on one grain of magnetic moment $`\stackrel{}{m}`$ in the nonuniform field $`𝐁^{tot}`$ is

$$f_i=m_j_jB_i^{tot}=m_j_iB_j^{tot},$$
$`(2)`$
since $`\times 𝐁^{tot}=0`$ for a static system. Repeated indices are summed over. The resulting macroscopic force by unit volume is

$$\frac{dF_i}{d^3𝐫}=M_j^{tot}_iB_j^{tot}=\frac{\chi }{2}_i(𝐁^{tot}𝐁^{tot})$$

where $`𝐌^{tot}=𝐌+𝐌^{}=\chi (𝐁+𝐁^{})`$ is the magnetization density and $`\chi =\frac{\mu 1}{\mu }`$. This field of force builds the pressure

$$p=\frac{\chi }{2}𝐁^{tot}𝐁^{tot}=\frac{\chi }{2}(𝐁^2+𝐁^2)+𝐌𝐁^{}$$

The first two terms are symmetrical about the box and exert no net force on it. The last term gives

$$𝐟_2=(𝐌𝐁^{})𝑑𝐒=(𝐁^{}d𝐒)𝐌=g𝐌.$$

Here $`d𝐒`$ is the vector representation of a surface element of the box and is directed outward. Permuting $`d𝐒`$ and $`𝐌`$ was allowed because they are parallel in the region where $`𝐁^{}`$ is important. The last equality comes from Gauss theorem for magnetic charges. Adding $`𝐟_1`$ and $`𝐟_2`$ one recovers the result (1) :

$$𝐟=𝐟_1+𝐟_2=g(𝐁𝐌)=g𝐇.$$
$`(3)`$
Most probably, (1) can be generalized to any shape of box. Thus, the relevant field which drives a monopole in matter is $`𝐇`$. It is the field found in a parallel elongated cavity, as for the force $`𝐟=e𝐄`$ driving an electric charge inside a dielectric. Mnemonic : this kind of cavity allows the test charge to follow the force without touching the matter.

Eq.(1) allows trapping a (not too heavy) monopole in the pole of a permanent magnet, where the lines of $`g𝐇`$ converge from all directions. This would not be true for $`g𝐁`$.

Microscopic interpretation

Atomic magnetic dipoles are often pictured as microscopic loops of electrical current. Then $`𝐁`$ appears as the average of the microscopic field $`𝐛`$ over semi-macroscopic volumes sufficiently large compared to the atomic scale. The work of $`𝐛`$ along a straight line $`L`$ is therefore

$$g_L𝐛𝑑𝐥=g_L𝐁𝑑𝐥.$$

In contrast, the work along a line $`L^{}`$ which avoids passing through the loops is

$$g_L^{}𝐛𝑑𝐥=g_L^{}𝐇𝑑𝐥.$$

Eq.(1) implies that the monopole avoids passing through the microscopic current loops or more likely that the loops move to "dodge" the monopole. This was of course the case with the ferrofluid, but in solid matter the atoms cannot escape from the monopole trajectory. Does it means that (1) is false if the monopole goes through an atom ? Not necessarily. At the approach of the pole, the electron wave functions are deformed and, if the monopole is sufficiently slow, they return adiabatically to the ground states. Thus no energy is exchanged between the monopole and the atom, as if the loop "dodges" the monopole.

3. VELOCITY-DEPENDENT FORCE

In vacuum, the analog of the Lorentz force for a moving monopole is $`𝐟=g𝐯\times 𝐄`$. Accordingly, a piece $`d𝐥`$ of wire carrying a current $`I^{}`$ of magnetic charges is subject to the dual Laplace force $`d𝐟=I^{}d𝐥\times 𝐄`$. Following the ferrofluid example, we consider a wire protected by a waterproof tube in a liquid dielectric (Fig.2) :

a) the tube is flattened perpendicular to $`𝐃=\epsilon 𝐄`$. Then $`𝐄_{tube}`$ = $`𝐃`$ and the measured force is

$$d𝐟=I^{}d𝐥\times 𝐃.$$
$`(4)`$
b) the tube is flattened parallel to $`𝐃=\epsilon 𝐄`$. Then $`𝐄_{tube}`$ = $`𝐄`$. The force acting on the wire is $`d𝐟_1=I^{}d𝐥\times 𝐄`$. On the other hand, on the right of the tube, the total field $`𝐃+𝐃^{}`$ is larger than on the left. The polar molecules are therefore attracted toward this region, building an excess of pressure which pushes the tube toward the left. Calculations like those between Eqs. (2) and (3) give the thrust $`d𝐟_2=I^{}d𝐥\times 𝐏`$, where $`𝐏`$ is the macroscopic electric polarization. In total,

$$d𝐟=d𝐟_1+d𝐟_2=I^{}d𝐥\times (𝐄+𝐏)$$

is equivalent to (4). Thus the field acting on a wire of magnetic current is $`𝐃`$. It is the field found in a perpendicular flattened cavity, as for the Laplace force $`d𝐟=Id𝐥\times 𝐁`$ on an ordinary current in a magnetized matter. Mnemonic : this cavity allows the wire to follow the force without touching the matter.

For a moving monopole, (4) becomes

$$𝐟=g𝐯\times 𝐃.$$
$`(4^{})`$
Microscopic interpretation

In a dielectric, $`𝐄`$ is the average of the microscopic field $`𝐞`$ over volumes sufficiently large compared to the molecular scale. The work of $`𝐞`$ when the wire sweeps a flat surface $`S`$ is

$$I^{}_S𝐞𝑑𝐒=I^{}_S𝐄𝑑𝐒.$$

In contrast, the work of $`𝐞`$ along a surface $`S^{}`$ which avoids cutting the dipole molecules is

$$I^{}_S^{}𝐞𝑑𝐒=I^{}_S^{}𝐃𝑑𝐒.$$

Eq.(4) implies that a moving dual wire avoids cutting the dipole molecules, or that the molecules "dodge" the wire. This has a well-defined topological meaning. Let us recall however that this wire was introduced to make the problem time-independent. For a moving monopole, there is no swept surface and the topological interpretation is lost.

Gathering (1) and (4’), the total force on a magnetic charge is

$$𝐟=g(𝐇𝐯\times 𝐃).$$
$`(5)`$
This result does not take into account dissipation and holds only for sufficiently slow monopoles, such that atoms and molecules evolve adiabatically under the influence of the monopole field.

4. MAXWELL EQUATIONS IN DYONIC MATTER

We consider matter containing magnetic charges $`\pm g`$ bound in magnetically neutral molecules, in addition to ordinary particles. These molecules posses magnetic dipoles of the form $`g𝐫`$, building a macroscopic magnetic polarization $`𝐩^{}`$ ($`𝐫`$ is the north-south charge separation). If they have spin, they also posses electric dipoles of the form $`\gamma ^{}𝐒`$ building a macroscopic electric polarization $`𝐦^{}`$ ($`\gamma ^{}`$ is the "giro-electric" ratio). $`𝐩^{}`$ and $`𝐦^{}`$ are dual respectively to the polarization $`𝐩`$ and the magnetization $`𝐦`$ built by the ordinary particles. The dual-symmetric Maxwell equations for the space average of the microscopic fields are

$$\times \overline{𝐛}_t\overline{𝐞}=𝐣+\delta 𝐣$$

$$\overline{𝐞}=\rho +\delta \rho $$

$$\times \overline{𝐞}_t\overline{𝐛}=𝐣^{}+\delta 𝐣^{}$$

$$\overline{𝐞}=\rho ^{}+\delta \rho ^{}.$$

Here $`(\rho ,𝐣)`$ is the external ordinary charge-current density and $`(\delta \rho ,\delta 𝐣)`$ the induced one, given by

$$\delta \rho =𝐩$$

$$\delta 𝐣=\times 𝐦+_t𝐩$$

Similarly, for the magnetic charge analogues,

$$\delta \rho ^{}=𝐩^{}$$

$$\delta 𝐣^{}=\times 𝐦^{}+_t𝐩^{}$$

From these equations, we can write the dual-symmetric Maxwell equations in matter :

$$\times 𝐇_t𝐃=𝐣$$

$$𝐃=\rho $$

$$\times 𝐄_t𝐁=𝐣^{}$$

$$𝐁=\rho ^{}$$
$`(6)`$
where $`𝐇=\overline{𝐛}𝐦`$, $`𝐃=\overline{𝐞}+𝐩`$ as usual, but $`𝐄=\overline{𝐞}𝐦^{}`$ and $`𝐁=𝐛+𝐩^{}`$. We see that $`𝐄`$ and $`𝐁`$ can no more be interpreted as the spatial averages of the microscopic fields. In that sense they are no more "fundamental" than $`𝐃`$ and $`𝐇`$. In fact the dual of $`𝐄`$ is not $`𝐁`$ but $`𝐇`$, whereas the dual of $`𝐁`$ is $`𝐃`$. The usual relations

$$𝐃=𝐄+𝐏,𝐁=𝐇+𝐌$$
$`(7)`$
are recovered, defining

$$𝐏𝐩+𝐦^{},𝐌𝐦+𝐩^{}.$$
$`(8)`$
It means that the microscopic nature of the dipole is forgotten at the level of the macroscopic Maxwell equations. Only their long range fields in $`1/r^3`$ are relevant. As in ordinary matter, $`𝐄`$ and $`𝐇`$ are the fields found in elongated cavities parallel to the respective fields, whereas $`𝐃`$ and $`𝐁`$ are found in perpendicular flat cavities.

5. THE DYONIC PERMITTIVITY-PERMEABILITY MATRIX

We assume that the polarizations $`𝐏`$ and $`𝐌`$ respond linearly to the macroscopic fields $`𝐃`$ and $`𝐁`$.

$$\left(\begin{array}{c}𝐏\\ 𝐌\end{array}\right)=[\chi ]\left(\begin{array}{c}𝐃\\ 𝐁\end{array}\right)=[\chi ^{}]\left(\begin{array}{c}𝐄\\ 𝐇\end{array}\right),$$
$`(9)`$
with $`[1+\chi ^{}][1\chi ]^1`$. In ordinary matter $`\chi _{11}=\chi _e=\chi _e^{}/\epsilon `$, $`\chi _e^{}\epsilon 1`$ ; $`\chi _{22}=\chi _m=\chi _m^{}/\mu `$, $`\chi _m^{}\mu 1`$, and $`\chi _{12}=\chi _{21}=0`$. In a matter containing only one species of dyon $`(e,g)`$ and antidyon $`(e,g)`$ bound in polar molecules, $`𝐏=𝐩`$ and $`𝐌=𝐩^{}`$ are linked by

$$\frac{𝐩}{e}=\frac{𝐩^{}}{g},\frac{1}{e}\frac{p_i}{D_j}=\frac{1}{g}\frac{p_i}{B_j}$$
$`(10)`$
wherefrom

$$[\chi ]=\mathrm{C}^{te}\left(\begin{array}{cc}e^2& eg\\ eg& g^2\end{array}\right).$$
$`(11)`$
An analogous matrix, with $`eg`$, is obtained with dipoles coming from spinning dyons ($`𝐦/e=𝐦^{}/g`$). We note that $`[\chi ]`$ and $`[\chi ^{}]`$ are symmetrical matrices. This remains true for a mixture of different species of molecules.

Thus, the usual relations $`𝐃=\epsilon 𝐄`$, $`𝐁=\mu 𝐇`$ are replaced by

$$\left(\begin{array}{c}𝐄\\ 𝐇\end{array}\right)=[1\chi ]\left(\begin{array}{c}𝐃\\ 𝐁\end{array}\right),$$
$`(12)`$
The speed of light is

$$c=(det[1\chi ])^{\frac{1}{2}}(c_{vac.}1).$$
$`(13)`$
Whatever they come from, the nondiagonal elements of $`[\chi ]`$ violate P- and T- symmetries, since they connect vector to pseudovectors. However PT is conserved.

6. ENERGY-MOMENTUM TENSOR

The various components of the energy-momentum tensor $`\mathrm{\Theta }^{\mu \nu }`$ can be derived from energy and momentum conservation in simple physical systems. Let us suppose that the whole space is filled with dual-symmetric matter. To get $`\mathrm{\Theta }^{i0}`$ (energy flow) and $`\mathrm{\Theta }^{ij}`$ (momentum flow) one considers a sandwich made of three slab-like regions of the z coordinate, R1 = $`[a,0]`$, R2 = $`[0,b]`$ and R3 = $`[b,b+a]`$. R1 carries uniform electric and magnetic charge-current densities, $`\{\rho ,𝐣;\rho ^{},𝐣^{}\}`$ and R3 carries the opposite densities, such that the fields vanish outside the sandwich. Solving (6) and (12) with appropriate $`\rho ,𝐣,\rho ^{},𝐣^{}`$, any kind of uniform field configuration $`\{𝐄,𝐃;𝐇,𝐁\}`$ can be obtained in R2. These fields are linearly attenuated in R1 and R3. In R3 a power

$$\frac{dW}{dtd^3𝐫}=𝐄𝐣+𝐇𝐣^{}$$
$`(14)`$
is dissipated and a force

$$\frac{d𝐟}{d^3𝐫}=\rho 𝐄+𝐣\times 𝐁+\rho ^{}𝐇𝐣^{}\times 𝐃$$
$`(15)`$
is exerted per unit of volume. The same quantities per unit of area (integrated over $`z`$ in R3) give $`\mathrm{\Theta }^{z0}`$ and $`\mathrm{\Theta }^{zi}`$ in R2.

To get $`\mathrm{\Theta }^{00}`$ (energy density) and $`\mathrm{\Theta }^{0i}`$ (momentum density) one has to "rotate" the sandwich in the 4-dimensional space-time, replacing $`z`$ by $`t`$ and slabs by time-slices or "epoch" T1, T2, T3. During T1 the (3-dimensional) space is filled with uniform current densities $`𝐣`$ and $`𝐣^{}`$, which progressively build uniform fields according (6) and (12). The second epoch is current-free and the uniform fields remain constant. The last epoch destroys the fields with opposite currents. Integrating (14) and (15) over $`t`$ in T3 give $`\mathrm{\Theta }^{00}`$ and $`\mathrm{\Theta }^{0i}`$ in T2. This method is detailled in . One obtains

$$\mathrm{\Theta }^{\mu \nu }=\left(\begin{array}{cc}\frac{𝐄𝐃+𝐇𝐁}{2}& 𝐃\times 𝐁\\ 𝐄\times 𝐇& \mathrm{\Theta }^{00}\delta ^{ij}D^iE^jB^iH^j\end{array}\right)$$
$`(16)`$
as in ordinary matter.

The Dirac condition in matter. One way to derive the Dirac condition between an electron and a monopole is to quantize the joint angular momentum of their fields which are

$$𝐃=\frac{e𝐫}{4\pi r^3},𝐁=\frac{g𝐫^{}}{4\pi r^3}$$

where $`𝐫`$ (resp. $`𝐫^{}`$) is the distance from the charge (resp. the pole) to the observation point. According to (15) the momentum density is $`\mathrm{\Theta }^{0i}=(𝐃\times 𝐁)_i`$, from which one gets the angular momentum

$$𝐉=d^3𝐫𝐫\times (𝐃\times 𝐁)=\frac{eg}{4\pi }\widehat{𝐧}$$
$`(17)`$
where $`\widehat{𝐧}`$ is the unit vector from the charge toward the pole. The usual Dirac condition $`eg=2n\pi \mathrm{}`$ is obtained from the quantization rule $`𝐉\widehat{𝐧}=n\mathrm{}/2`$. Note that if the momentum density were $`𝐄\times 𝐇`$, as sometimes advocated (see the discussion in ), the Dirac condition in medium would not be consistent with that in vacuum.

7. APPLICATION TO THE SEARCH FOR AN ELECTRIC DIPOLE MOMENT

The dual-symmetric formalism applies as well to the case were the electron (or the nucleus) possesses an electric dipole moment (e.d.m.) $`\stackrel{}{d}=\gamma ^{}𝐒`$ in addition to the usual magnetic moment $`\stackrel{}{m}=\gamma 𝐒`$. Then we have $`𝐦/\gamma =𝐦^{}/\gamma ^{}`$ and a nondiagonal $`[\chi ]`$ matrix element is generated, like with dyonic molecules (Eqs.10-11) :

$$[\chi ]=\left(\begin{array}{cc}\chi _e+\chi _mr^2& \chi _mr\\ \chi _mr& \chi _m\end{array}\right),$$
$`(18)`$
where $`r\gamma ^{}/\gamma `$. $`\chi _m=\frac{\mu 1}{\mu }`$ comes from the spinning electrons and $`\chi _e=\frac{\epsilon 1}{\epsilon }`$ from polar molecules.

A nonzero $`\chi _{12}`$ may be generated in another way : the e.d.m. tends to align the spin of an electron along the internal electric field of a polar molecule. It couples $`𝐦`$ to $`𝐩`$. Here we consider only the first mechanism.

Eq.(18) suggests two possible measurements of $`r`$ :

a) In Fig.3a, a cylinder of magnetizable, but insulating material is immersed in a large magnetic field $`𝐁_0`$. The inside field $`𝐁`$ induces a small electric polarization $`𝐏=\chi _mr𝐁`$ and an electric field $`𝐄`$. If the cylinder is much broader than high, we have $`𝐁=𝐁_0`$, $`𝐃0`$ and $`𝐄𝐏`$. More generally one has

$$𝐄=x\chi _mr𝐁_0$$
$`(19)`$
where the coefficient $`x<1`$ depends on the container geometry. Let us take a cubic container of size $`L`$. Between the top and the bottom, we can measure a potential difference $`U=EL`$. The ratio between the stored electrostatic energy $`W=\frac{1}{2}\epsilon 𝐄^2L^3`$ and the magnetic one $`W_0=\frac{1}{2\mu }𝐁_0^2L^3`$ is

$$\frac{W}{W_0}=\epsilon \mu (x\chi _mr)^2.$$
$`(20)`$
in terms of common units, we have

$$\frac{U}{\mathrm{volt}}=\mathrm{3\hspace{0.17em}10}^8x\chi _mr\frac{L}{\mathrm{metre}}\frac{B_0}{\mathrm{tesla}}$$

$$\frac{W}{\mathrm{eV}}=\mathrm{2.5\hspace{0.17em}10}^{24}x^2\epsilon \chi _m^2r^2\left(\frac{L}{\mathrm{metre}}\right)^3\left(\frac{B_0}{\mathrm{tesla}}\right)^2.$$

\[useful relations are : 1 tesla = 3 108 volt/metre,

1 (tesla)2 (metre)3 = 0.8 106 joule,

1 eV = 1.6 10-19 joule $`10^4`$ kelvin = $`\mathrm{5\hspace{0.17em}10}^6\mathrm{}`$/metre

and $`\gamma _e=e\mathrm{}/m_e=\mathrm{4\hspace{0.17em}10}^{11}`$ e$`\times `$cm\].

Let us assume $`r=10^{16}`$, which corresponds to an electron e.d.m. of $`\mathrm{2\hspace{0.17em}10}^{27}`$ e$`\times `$cm, $`\epsilon 1`$, $`\chi _m0.5`$ and $`x0.5`$. For a field of 1 tesla, and a cube of 1 meter, a potential difference of about $`\mathrm{0.5\hspace{0.17em}10}^8`$ volt is obtained. $`W10^9`$ eV $`10^5`$ kelvin. The voltmetre has to be cooled at least to this temperature to prevent thermal noise.

b) in Fig.3b the same container is put in a large electric field $`𝐄_0`$. Using the second matrix of (9)

$$[\chi ^{}]\left(\begin{array}{cc}\epsilon 1& \epsilon \chi _m^{}r\\ \epsilon \chi _m^{}r& \chi _m^{}\end{array}\right),$$
$`(18^{})`$
with $`\chi _m^{}\mu 1`$, we predict a small magnetization $`𝐌=\epsilon \chi _m^{}r𝐄`$. If the cylinder is much higher than broad, we have $`𝐄𝐄_0`$, $`𝐇0`$ and $`𝐁𝐌`$. For a cubic container we assume

$$𝐁x\epsilon \chi _m^{}r𝐄_0$$
$`(21)`$
with $`x0.5`$. This field can be measured by a SQUID encircling the container. The phase shift of the wave function in one loop is

$$\phi =eL^2B/\mathrm{}$$

$$=\mathrm{0.5\hspace{0.17em}10}^{12}x\epsilon \chi _m^{}r\left(\frac{L}{\mathrm{metre}}\right)^2\frac{E_0}{10^5\mathrm{volt}/\mathrm{metre}}.$$

This phase can be multiplied by a large number of turns around the cylinder. The ratio between the output (magnetostatic) energy $`W`$ and the input (electric) one $`W_0`$ is still given by (20), but $`W`$ and $`W_0`$ are typically $`10^5`$ times smaller and the temperature must be much lower than in case a).

8. CONCLUSION

We have given arguments that the macroscopic fields acting on magnetic charges and currents are $`𝐇`$ and $`𝐃`$. Comparing with electric charges and currents, one has a unified mnemonic principle : in each case, the acting field is the one found in a parallel-elongated (resp. flat-perpendicular) cavity in which a charge (resp. current wire) can follow the force without touching the medium. In a classical microscopic picture, a monopole avoids passing through the microscopic current loops and a dual current wire avoids cutting the dipole molecules. Quantum mechanically, it means that the perturbation of the atoms and molecules lying on the trajectory of the monopole is adiabatic. This should be the case at low enough velocity in a liquid. The monopole will be presumably accompanied by a swarm of atoms magnetically (or electrically, for a dyon) bound to it. In a solid, such a swarm could forbid the monopole to move without producing cracks.

The dual-symmetric Maxwell equations in matter are formally unchanged, but $`𝐄`$ and $`𝐁`$ can no more be interpreted as the spatial averages of the microscopic fields. The duality correspondance is $`𝐄𝐇`$ and $`𝐃𝐁`$. When dyons are present, or when ordinary particles possess electric dipole moments, $`\epsilon `$ and $`\mu `$ are replaced by a permittivity-permeability matrix $`[1\chi ]`$ whose nondiagonal elements violate the P- and T- symmetries (but not PT). The energy momentum tensor is also unchanged. The usual Dirac condition $`eg=2n\pi \mathrm{}`$ is obtained provided the momentum density is $`𝐃\times 𝐁`$. These results have been obtained under the hypothesis that $`𝐏`$ and $`𝐌`$ are linear in the fields.

As an application of the dual-symmetric formalism, two possible measurements of the electron e.d.m. have been suggested. They are at the limit of the present technological possibilities. However, mechanisms like the $`𝐦𝐩`$ coupling in a polar molecule mentioned in Sect.7 might enhance the signal.

REFERENCES

1. Y. Hara, Electromagnetic force on a magnetic monopole, Phys. Rev., 1985, v. A32, p.1002-1006

2. D. Fayolle, Dynamique d’un monopôle magnétique dans la matière, LYCEN/9961, july 1999, unpublished.

3. F.N.H. Robinson, Electromagnetic stress and momentum in matter, Phys. Rep., 1975, v. C16, p. 313-354.


text:
# The Einstein-Vlasov system

## 1 Introduction

The basic equations of general relativity are the Einstein equations coupled to some other partial differential equations describing the matter content of spacetime. There are many choices of matter model which are of physical interest and Yvonne Choquet-Bruhat has published fundamental results on the Cauchy problem for the Einstein equations coupled to a wide variety of matter models. One of these is collisionless matter described by the Vlasov equation. It is the subject of these lectures.

The Vlasov equation arises in kinetic theory. It gives a statistical description of a collection of particles. It is distinguished from other equations of kinetic theory by the fact that there is no direct interaction between particles. In particular, no collisions are included in the model. Each particle is acted on only by fields which are generated collectively by all particles together. The fields which are taken into account depend on the physical situation being modelled. In plasma physics, where this equation is very important, the interaction is electromagnetic and the fields are described either by the Maxwell equations or, in a quasi-static approximation, by the Poisson equation . In gravitational physics, which is the subject of the following, the fields are described by the Einstein equations or, in the Newtonian approximation, by the Poisson equation. (There is a sign difference in the Poisson equation in comparison with the electromagnetic case due to the replacement of a repulsive by an attractive force.) The best known applications of the Vlasov equation to self-gravitating systems are to stellar dynamics . It can also be applied to cosmology. In the first case the systems considered are galaxies or parts of galaxies where there is not too much dust or gas which would require a hydrodynamical treatment. Possible applications are to globular clusters, elliptical galaxies and the central bulge of spiral galaxies. The ‘particles’ in all these cases are stars. In the cosmological case they might be galaxies or even clusters of galaxies. The fact that they are modelled as particles reflects the fact that their internal structure is believed to be irrelevant for the dynamics of the system as a whole. The Vlasov equation is also used in cosmology to model non-baryonic dark matter (, p. 323). In that case the ‘particles’ are elementary particles.

These lectures are concerned not with the above physical applications but with some basic mathematical aspects of the Einstein-Vlasov system. First the definition and general mathematical properties of this system of partial differential equations are discussed and then the Cauchy problem for the system is formulated. The central theme in what follows is the global Cauchy problem, where ‘global’ means global in time. The known results on this and related problems are surveyed and important methods used are highlighted. Further information on kinetic theory in general relativity may be found in .

Let $`(M,g_{\alpha \beta })`$ be a spacetime, i.e. $`M`$ is a four-dimensional manifold and $`g_{\alpha \beta }`$ is a metric of Lorentz signature $`(,+,+,+)`$. Note that $`g_{\alpha \beta }`$ denotes a geometric object here and not the components of the geometric object in a particular coordinate system. In other words the indices are abstract indices. (See , section 2.4 for a discussion of this notation.) It is always assumed that the metric is time-orientable, i.e. that the two halves of the light cone at each point of $`M`$ can be labelled past and future in a way which varies continuously from point to point. With this global direction of time, it is possible to distinguish between future-pointing and past-pointing timelike vectors. The worldline of a particle of non-zero rest mass $`m`$ is a timelike curve in spacetime. The unit future-pointing tangent vector to this curve is the 4-velocity $`v^\alpha `$ of the particle. Its 4-momentum $`p^\alpha `$ is given by $`mv^\alpha `$. There are different variants of the Vlasov equation depending on the assumptions made. Here it is assumed that all particles have the same mass $`m`$ but it would also be possible to allow a continuous range of masses. When all the masses are equal, units can be chosen so that $`m=1`$ and no distinction need be made between 4-velocity and 4-momentum. There is also the possibility of considering massless particles, whose wordlines are null curves. In the case $`m=1`$ the possible values of the four-momentum are precisely all future-pointing unit timelike vectors. These form a hypersurface $`P`$ in the tangent bundle $`TM`$ called the mass shell. The distribution function $`f`$, which represents the density of particles with given spacetime position and four-momentum, is a non-negative real-valued function on $`P`$. A basic postulate in general relativity is that a free particle travels along a geodesic. Consider a future-directed timelike geodesic parametrized by proper time. Then its tangent vector at any time is future-pointing unit timelike. Thus this geodesic has a natural lift to a curve on $`P`$, by taking its position and tangent vector together. This defines a flow on $`P`$. Denote the vector field which generates this flow by $`X`$. (This vector field is what is sometimes called the geodesic spray in the mathematics literature.) The condition that $`f`$ represents the distribution of a collection of particles moving freely in the given spacetime is that it should be constant along the flow, i.e. that $`Xf=0`$. This equation is the Vlasov equation, sometimes also known as the Liouville or collisionless Boltzmann equation.

To get an explicit expression for the Vlasov equation, it is necessary to introduce local coordinates on the mass shell. In the following local coordinates $`x^\alpha `$ on spacetime are always chosen such that the hypersurfaces $`x^0`$=const. are spacelike. (Greek and Latin indices take the values $`0,1,2,3`$ and $`1,2,3`$ respectively.) Intuitively this means that $`x^0`$, which may also be denoted by $`t`$, is a time coordinate and that the $`x^a`$ are spatial coordinates. A timelike vector is future-pointing if and only if its zero component in a coordinate system of this type is positive. It is not assumed that the vector $`/x^0`$ is timelike. One way of defining local coordinates on $`P`$ is to take the spacetime coordinates $`x^\alpha `$ together with the spatial components $`p^a`$ of the four-momentum in these coordinates. Then the explicit form of the Vlasov equation is:

$$f/t+(p^a/p^0)f/x^a(\mathrm{\Gamma }_{\beta \gamma }^ap^\beta p^\gamma /p^0)f/p^a=0$$
(1)
where $`\mathrm{\Gamma }_{\beta \gamma }^\alpha `$ are the Christoffel symbols associated to the metric $`g_{\alpha \beta }`$. Here it is understood that $`p^0`$ is to be expressed in terms of $`p^a`$ and the metric using the relation $`g_{\alpha \beta }p^\alpha p^\beta =1`$. An alternative form of the Vlasov equation which is often useful is obtained by coordinatizing the mass shell using the components of the four-momentum in an orthonormal frame instead of the coordinate components.

The Vlasov equation can be coupled to the Einstein equations as follows, giving rise to the Einstein-Vlasov system. The unknowns are a 4-manifold $`M`$, a (time orientable) Lorentz metric $`g_{\alpha \beta }`$ on $`M`$ and a non-negative real-valued function $`f`$ on the mass shell defined by $`g_{\alpha \beta }`$. The field equations consist of the Vlasov equation defined by the metric $`g_{\alpha \beta }`$ for $`f`$ and the Einstein equation $`G_{\alpha \beta }=8\pi T_{\alpha \beta }`$. (Units are chosen here so that the speed of light and the gravitational constant both have the numerical value unity.) To obtain a complete system of equations it remains to define $`T_{\alpha \beta }`$ in terms of $`f`$ and $`g_{\alpha \beta }`$. It is defined as an integral over the part of the mass shell over a given spacetime point with respect to a measure which will now be defined. The metric at a given point of spacetime defines in a tautological way a metric on the tangent space at that point. The part of the mass shell over that point is a submanifold of the tangent space and as such has an induced metric, which is Riemannian. The associated measure is the one which we are seeking. It is evidently invariant under Lorentz transformations of the tangent space, a fact which may be used to simplify computations in concrete situations. In the coordinates $`(x^\alpha ,p^a)`$ on $`P`$ the explicit form of the energy-momentum tensor is:

$$T_{\alpha \beta }=fp_\alpha p_\beta |g|^{1/2}/p_0𝑑p^1𝑑p^2𝑑p^3$$
(2)
A simple computation in normal coordinates based at a given point shows that $`T_{\alpha \beta }`$ defined by (2) is divergence-free, independently of the Einstein equations being satisfied. This is of course a necessary compatibility condition in order for the Einstein-Vlasov system to be a reasonable set of equations. Another important quantity is the particle current density, defined by:

$$N^\alpha =fp^\alpha |g|^{1/2}/p_0𝑑p^1𝑑p^2𝑑p^3$$
(3)
A computation in normal coordinates shows that $`_\alpha N^\alpha =0`$. This equation is an expression of the conservation of the number of particles. There are some inequalities which follow immediately from the definitions (2) and (3). Firstly $`N_\alpha V^\alpha 0`$ for any future-pointing timelike or null vector $`V^\alpha `$, with equality only if $`f=0`$ at the given point. Hence unless there are no particles at some point, the vector $`N^\alpha `$ is future-pointing timelike. Next, if $`V^\alpha `$ and $`W^\alpha `$ are any two future-pointing timelike vectors then $`T_{\alpha \beta }V^\alpha W^\beta 0`$. This is the dominant energy condition (, p. 91). Finally, if $`X^\alpha `$ is a spacelike vector then $`T_{\alpha \beta }X^\alpha X^\beta 0`$. This is the non-negative pressures condition. This condition, the dominant energy condition and the Einstein equations together imply that the Ricci tensor satisfies the inequality $`R_{\alpha \beta }V^\alpha V^\beta 0`$ for any timelike vector $`V^\alpha `$. The last inequality is called the strong energy condition. These inequalities constitute one of the reasons which mean that the Vlasov equation defines a well-behaved matter model in general relativity. However this is not the only reason. A perfect fluid with a reasonable equation of state or matter described by the Boltzmann equation also have energy-momentum tensors which satisfy these inequalities.

The Vlasov equation in a fixed spacetime is a linear hyperbolic equation for a scalar function and hence solving it is equivalent to solving the equations for its characteristics. In coordinate components these are:

$$dX^a/ds=P^a,dP^a/ds=\mathrm{\Gamma }_{\beta \gamma }^aP^\beta P^\gamma $$
(4)
Let $`X^a(s,x^\alpha ,p^a)`$, $`P^a(s,x^\alpha ,p^a)`$ be the unique solution of (4) with initial conditions $`X^a(t,x^\alpha ,p^a)=x^a`$ and $`P^a(t,x^\alpha ,p^a)=p^a`$. Then the solution of the Vlasov equation can be written as:

$$f(x^\alpha ,p^a)=f_0(X^a(0,x^\alpha ,p^a),P^a(0,x^\alpha ,p^a))$$
(5)
where $`f_0`$ is the restriction of $`f`$ to the hypersurface $`t=0`$. This function $`f_0`$ serves as initial datum for the Vlasov equation. It follows immediately from this that if $`f_0`$ is bounded by some constant $`C`$, the same is true of $`f`$. This obvious but important property of the solutions of the Vlasov equation is used frequently without comment in the study of this equation.

The above calculations involving $`T_{\alpha \beta }`$ and $`N^\alpha `$ were only formal. In order that they have a precise meaning it is necessary to impose some fall-off in the momentum variables on $`f`$ so that the integrals occurring exist. The simplest condition to impose is that $`f`$ has compact support for each fixed $`t`$. This property holds if the initial datum $`f_0`$ has compact support and if each hypersurface $`t=t_0`$ is a Cauchy hypersurface. For by the definition of a Cauchy hypersurface, each timelike curve which starts at $`t=0`$ hits the hypersurface $`t=t_0`$ at a unique point. Hence the geodesic flow defines a continuous mapping from the part of the mass shell over the initial hypersurface $`t=0`$ to the part over the hypersurface $`t=t_0`$. The support of $`f(t_0)`$, the restriction of $`f`$ to the hypersurface $`t=t_0`$ is the image of the support of $`f_0`$ under this continuous mapping and so is compact. Let $`P(t)`$ be the supremum of the values of $`|p^a|`$ attained on the support of $`f(t)`$. It turns out that in many cases controlling the solution of the Vlasov equation coupled to some field equation in the case of compactly supported initial data for the distribution function can be reduced to obtaining a bound for $`P(t)`$. An example of this is given below.

The data in the Cauchy problem for the Einstein equations coupled to any matter source consist of the induced metric $`g_{ab}`$ on the initial hypersurface, the second fundamental form $`k_{ab}`$ of this hypersurface and some matter data. In fact these objects should be thought of as objects on an abstract 3-dimensional manifold $`S`$. Thus the data consist of a Riemannian metric $`g_{ab}`$, a symmetric tensor $`k_{ab}`$ and appropriate matter data, all defined intrinsically on $`S`$. The nature of the initial data for the matter will now be examined in the case of the Einstein-Vlasov system. It is not quite obvious what to do, since the distribution function $`f`$ is defined on the mass shell and so the obvious choice of initial data, namely the restriction of $`f`$ to the initial hypersurface, is not appropriate. For it is defined on the part of the mass shell over the initial hypersurface and this is not intrinsic to $`S`$. This difficulty can be overcome as follows. Let $`\varphi `$ be the mapping which sends a point of the mass shell over the initial hypersurface to its orthogonal projection onto the tangent space to the initial hypersurface. The map $`\varphi `$ is a diffeomorphism. The abstract initial datum $`f_0`$ for $`f`$ is taken to be a function on the tangent bundle of $`S`$. The initial condition imposed is that the restriction of $`f`$ to the part of the mass shell over the initial hypersurface should be equal to $`f_0`$ composed with $`\varphi `$. An initial data set for the Einstein equations must satisfy the constraints and in order that the definition of an abstract initial data set for the Einstein equations be adequate it is necessary that the constraints be expressible purely in terms of the abstract initial data. The constraint equations are:

$`Rk_{ab}k^{ab}+(\mathrm{tr}k)^2`$ $`=16\pi \rho `$ (6)
$`_ak_b^a_b(\mathrm{tr}k)`$ $`=8\pi j_b`$ (7)
Here $`R`$ denotes the scalar curvature of the metric $`g_{ab}`$. If $`n^\alpha `$ denotes the future-pointing unit normal vector to the initial hypersurface and $`h^{\alpha \beta }=g^{\alpha \beta }+n^\alpha n^\beta `$ is the orthogonal projection onto the tangent space to the initial hypersurface then $`\rho =T_{\alpha \beta }n^\alpha n^\beta `$ and $`j^\alpha =h^{\alpha \beta }T_{\beta \gamma }n^\gamma `$. The vector $`j^\alpha `$ satisfies $`j^\alpha n_\alpha =0`$ and so can be naturally identified with a vector intrinsic to the initial hypersurface, denoted here by $`j^a`$. What needs to be done is to express $`\rho `$ and $`j_a`$ in terms of the intrinsic initial data. They are given by the following expressions:

$`\rho `$ $`={\displaystyle f_0(p^a)p^ap_a/(1+p^ap_a)^{1/2}({}_{}{}^{(3)}g)^{1/2}𝑑p^1𝑑p^2𝑑p^3}`$ (8)
$`j_a`$ $`={\displaystyle f_0(p^a)p_a({}_{}{}^{(3)}g)^{1/2}𝑑p^1𝑑p^2𝑑p^3}`$ (9)
If a three-dimensional manifold on which an initial data set for the Einstein-Vlasov system is defined is mapped into a spacetime by an embedding $`\psi `$ then the embedding is said to induce the given initial data on $`S`$ if the induced metric and second fundamental form of $`\psi (S)`$ coincide with the results of transporting $`g_{ab}`$ and $`k_{ab}`$ with $`\psi `$ and the relation $`f=f_0\varphi `$ holds, as above. A form of the local existence and uniqueness theorem can now be stated. This will only be done for the case of smooth (i.e. infinitely differentiable) initial data although versions of the theorem exist for data of finite differentiability.

Theorem 1.1 Let $`S`$ be a 3-dimensional manifold, $`g_{ab}`$ a smooth Riemannian metric on $`S`$, $`k_{ab}`$ a smooth symmetric tensor on $`S`$ and $`f_0`$ a smooth non-negative function of compact support on the tangent bundle $`TS`$ of $`S`$. Suppose further that these objects satisfy the constraint equations. Then there exists a smooth spacetime $`(M,g_{\alpha \beta })`$, a smooth distribution function $`f`$ on the mass shell of this spacetime and a smooth embedding $`\psi `$ of $`S`$ into $`M`$ which induces the given initial data on $`S`$ such that $`g_{\alpha \beta }`$ and $`f`$ satisfy the Einstein-Vlasov system and $`\psi (S)`$ is a Cauchy hypersurface. Moreover, given any other spacetime $`(M^{},g_{\alpha \beta }^{})`$, distribution function $`f^{}`$ and embedding $`\psi ^{}`$ satisfying these conditions, there exists a diffeomorphism $`\chi `$ from an open neighbourhood of $`\psi (S)`$ in $`M`$ to an open neighbourhood of $`\psi ^{}(S)`$ in $`M^{}`$ which satisfies $`\chi \psi =\psi ^{}`$ and carries $`g_{\alpha \beta }`$ and $`f`$ to $`g_{\alpha \beta }^{}`$ and $`f^{}`$ respectively.

The formal statement of this theorem is rather complicated, but its essential meaning is as follows. Given an initial data set (satisfying the constraints) there exists a corresponding solution of the Einstein-Vlasov system and this solution is locally unique up to diffeomorphism. There also exists a global uniqueness statement which uses the notion of the maximal Cauchy development of an initial data set, but this is not required in the following. The first proof of a theorem of the above kind for the Einstein-Vlasov system was given by Yvonne Choquet-Bruhat in .

The problem of extending this local theorem to one which is in some sense global is a very difficult one. In fact with presently available mathematical techniques it is too difficult. One way of making some progress in understanding the general problem is to study the simplified cases obtained by imposing various symmetries on the solutions. Note that if a symmetry is imposed on the initial data for the Cauchy problem this is inherited by the corresponding solutions. (See , section 5.6 for a discussion of this.) This ensures the consistency of restricting the problem to a particular symmetry class.

In the following different symmetry classes will be considered in turn, proceeding from the strongest to the weakest assumptions. First spatially homogeneous solutions are discussed. These are simple enough that it is possible to make statements about the Einstein equations coupled to a general class of matter models. After this has been done, further results which can be obtained in the particular case where the matter is described by the Vlasov equation are presented. The first inhomogeneous solutions to be discussed are those which are static and spherically symmetric. Apart from their intrinsic interest these are of relevance for the study of spherically symmetric collapse which is discussed next. Brief comments are made on dynamical cosmological solutions before concentrating on one question where there are results on solutions of the Einstein-Vlasov system without any symmetry assumptions being required. This concerns the construction of solutions with an isotropic singularity.

## 2 Spatially homogeneous solutions I: general matter models

A solution of the Einstein equations coupled to some matter equations is said to be symmetric under the action of a Lie group $`G`$ if $`G`$ acts by isometries of the metric which also leave the matter fields invariant. A solution of the Einstein-matter equations is called spatially homogeneous if it is symmetric under the action of a Lie group whose orbits are spacelike hypersurfaces. If we think of these as hypersurfaces of constant time then the metric only depends on time and the Einstein equations reduce to ordinary differential equations, an enormous simplification. The equations of motion of matter fields which are defined on spacetime also reduce to ODE’s. Since the Vlasov equation is defined on the mass shell it in general still contains derivatives with respect to the momenta and thus remains a partial differential equation in the spatially homogeneous case.

The spatially homogeneous spacetimes can be classified into various types according to the Lie group involved. The conventional terminology is that there are nine Bianchi types I-IX and one additional type, the Kantowski-Sachs models. The latter will not be discussed further here. We might like to use spatially homogeneous spacetimes as cosmological models. Since there is in a sense no spatial dependence we do need to worry about spatial boundary conditions. Eventually, however, we would like to consider more realistic cosmological models which are inhomogeneous perturbations of the Bianchi spacetimes, and then boundary conditions become important. One simple condition to pose is that the spacetimes involved should contain a compact Cauchy surface. Then there is no danger of extra information coming in from infinity. Supposing that it is desired to impose this condition of spatial compactness the question arises if all Bianchi types are compatible with it. Unfortunately this is not the case. In fact the only ones which are are types I and IX.

A larger class of spatially compact spacetimes where the Einstein equations reduce to ODE’s are the locally spatially homogeneous ones, as defined in . The idea is to require that the spacetime itself be spatially compact while its universal cover is spatially homogeneous. For details see . This allows a much bigger class of Bianchi types to be included.

The idea now is to consider solutions of the Einstein-matter equations which are spatially compact and locally spatially homogeneous only assuming some general conditions on the matter model. These conditions are satisfied in the case of the Vlasov equation but also, for example, in the case of the Euler equation describing a perfect fluid with a physically reasonable equation of state. This generality is a luxury we can only afford due to the assumption of (local) spatial homogeneity. In general solutions of the Euler equations can be expected to form shocks which leads to a breakdown of the solution of the evolution equations. In the homogeneous case this possibility does not arise. It is the absence of shock formation which makes the Vlasov equation particularly convenient to work with when studying inhomogeneous spacetimes.

The matter models to be considered will be defined in terms of some general properties. As usual $`T^{\alpha \beta }`$ denotes the energy-momentum tensor. When a specific matter model has been chosen $`T^{\alpha \beta }`$ will be a functional of some matter variables, denoted collectively by $`F`$, and the spacetime metric $`g_{\alpha \beta }`$. In the following another quantity $`N^\alpha `$ (called the particle current density) will be required. It is also assumed to be a functional of $`F`$ and $`g_{\alpha \beta }`$. Now various properties which will be assumed at appropriate points will be listed.

(1) $`T^{\alpha \beta }V_\alpha W_\beta 0`$ for all future-pointing timelike vectors $`V^\alpha `$ and $`W^\alpha `$ (dominant energy condition)

(2) $`T^{\alpha \beta }(g_{\alpha \beta }+V_\alpha V_\beta )0`$ for all unit timelike vectors $`V^\alpha `$ (non-negative sum pressures condition)

(3) for any $`F`$ and $`g_{\alpha \beta }`$ the conditions $`_\alpha N^\alpha =0`$ and $`_\alpha T^{\alpha \beta }=0`$ are satisfied (conservation conditions)

(4) for any $`F`$ and $`g_{\alpha \beta }`$ the vector $`N^\alpha `$ is future-pointing timelike or zero

(5) for any constant $`C_1>0`$ there exists a positive constant $`C_2`$ such that for any $`F`$ and $`g_{\alpha \beta }`$ with $`N_\alpha N^\alpha C_1`$ and any timelike vector $`V^\alpha `$ the following inequality holds:

$$T^{\alpha \beta }V_\alpha V_\beta C_2(N^\alpha V_\alpha )^2$$

(6) for any constant $`C_1>0`$ there exists a positive constant $`C_2<1`$ such that for any $`F`$ and $`g_{\alpha \beta }`$ with $`N^\alpha N_\alpha C_1`$ and any unit timelike vector $`V^\alpha `$

$$(g_{\alpha \beta }+V_\alpha V_\beta )T^{\alpha \beta }3C_2T^{\alpha \beta }V_\alpha V_\beta $$

(7) if a solution with the given symmetry of the Einstein equations coupled to the given matter model is such that the time coordinate defined above takes all values in the interval $`(t_1,t_2)`$ , if it is not possible to extend the spacetime so as to make this interval longer and if $`t_1`$ or $`t_2`$ is finite then $`\mathrm{tr}k(t)`$ is unbounded in a neighbourhood of $`t_1`$ or $`t_2`$ respectively.

(8) for any constant $`C_1>0`$ there exists a constant $`C_2>0`$ such that $`T_{\alpha \beta }T^{\alpha \beta }C_1`$ implies $`N_\alpha N^\alpha C_2`$

Some comments will now be made concerning the physical motivation of some of these conditions. If a given type of matter can be considered as being made up of particles then a particle current density $`N^\alpha `$ is defined. If the particles have positive rest mass then this vector is future pointing timelike or zero as required by condition (4). If the particles are massless then this condition is still satisfied except for very special types of matter where $`N^\alpha `$ might be null. If particles cannot be created or destroyed then $`N^\alpha `$ is divergenceless as required in condition (3). It is not easy to give an intuitive interpretation of conditions (5) and (6). The meaning of (5) is roughly as follows. If matter is observed from a boosted frame then the particle density is multiplied by a $`\gamma `$-factor arising from the effect of Lorentz contraction on the volume element. The observed energy density is also affected in this way but picks up an additional $`\gamma `$-factor. Hence when a given matter distribution is considered from a boosted frame the multiplicative factor in the observed energy density behaves like the square of that in the observed particle density. As for condition (6), for a perfect fluid it is related to the condition that the speed of sound should be bounded away from the speed of light. The given symmetry referred to in condition (7) will be the Bianchi symmetry being considered.

In a spatially compact spacetime an interesting quantity is the volume of the hypersurfaces of constant time. The Friedman-Robertson-Walker (FRW) models generally used in cosmology have the property that this volume either increases at all times (open models) or increases to maximum after which it decreases again (closed models). The possibility that it is always decreases, a case which is allowed mathematically, is usually ignored since we know that our universe is expanding at the present time. In any case these models can be obtained from those where the volume is always increasing by reversing the direction of time. Thus they present no essentially new mathematical phenomena.

In local spatially homogeneous spatially compact spacetimes with reasonable matter the same pattern is found. It can be proved that they also share some other significant physical properties with the FRW models. These spacetimes can be parametrized by a Gaussian time coordinate based on a (locally) homogeneous hypersurface. Suppose that the solution is maximal in the sense that it cannot be extended to a larger interval of Gaussian time. Then we would like to know two things. Firstly, if the volume is always increasing then the time of existence in the future is infinite and the spacetime is future geodesically complete. Secondly, if the volume is increasing at some time then the time of existence in the past is finite and as the time of breakdown is approached some geometrical invariant of the spacetime geometry diverges. This rules out the possibility of extending the spacetime in some way which is not globally hyperbolic. This is desirable form the point of view of the strong cosmic censorship hypothesis.

In  theorems of the desired type were proved. The first says that if conditions (1), (2) and (7) above are satisfied by some matter model then inextendible locally spatially homogeneous spatially compact solutions of the Einstein-matter equations where the volume is always increasing satisfy the first conclusion mentioned above. In particular, they are future geodesically complete. The second says that if (1)-(8) are satisfied and the spacetime is not vacuum then the curvature invariant $`G^{\alpha \beta }G_{\alpha \beta }`$ diverges as the finite past limit of the domain of existence is approached. The fundamental intuitive reason for this is that a finite amount of matter is being squeezed to zero volume as the singularity is approached.

These theorems apply to the Einstein-Vlasov system since in that case, as will now be discussed, conditions (1)-(8) are satisfied. In that case $`T^{\alpha \beta }`$ and $`N^\alpha `$ have been defined above. It has also been stated that (1) and (3) hold. Condition (2) is a consequence of the non-negative pressures condition mentioned above. To check condition (4) it is merely necessary to observe that if $`V^\alpha `$ is a future-pointing timelike vector then $`N_\alpha V^\alpha <0`$ unless $`N^\alpha =0`$ at a given spacetime point. Condition (7) is an existence theorem which says, roughly speaking, that as long as the geometry does not break down too badly, the solution of the matter equations (in this case the Vlasov equation) cannot break down in finite time. Here the mean curvature $`\mathrm{tr}k`$ plays the role of a controlling quantity whose boundedness ensures the continued existence of the solution. It remains to check conditions (5), (6) and (8) and to do this we may choose a frame whose timelike member is $`V^\alpha `$ in order to do the calculation. Then the inequalities of (5) and (6) become $`\widehat{T}^{00}C_2(\widehat{N}^0)^2`$ and $`\delta ^{ij}\widehat{T}_{ij}3C_2\widehat{T}^{00}`$. (The hats here indicate the use of indices associated to an orthonormal frame.)

$`N_\alpha N^\alpha `$ $`=\left({\displaystyle f(p^a)p_\alpha /p^0𝑑p^1𝑑p^2𝑑p^3}\right)\left({\displaystyle f(q^a)q^\alpha /q^0𝑑q^1𝑑q^2𝑑q^3}\right)`$
$`={\displaystyle f(p^a)f(q^a)p^\alpha q_\alpha /(p^0q^0)𝑑p^1𝑑p^2𝑑p^3𝑑q^1𝑑q^2𝑑q^3}`$
$`{\displaystyle f(p^a)f(q^a)/(p^0q^0)𝑑p^1𝑑p^2𝑑p^3𝑑q^1𝑑q^2𝑑q^3}`$
$`=\left({\displaystyle f(p^a)/p^0𝑑p^1𝑑p^2𝑑p^3}\right)^2`$
Hence

$`\widehat{N}^0`$ $`={\displaystyle f(p^a)𝑑p^1𝑑p^2𝑑p^3}`$
$`\left({\displaystyle f(p^a)p^0𝑑p^1𝑑p^2𝑑p^3}\right)^{1/2}\left({\displaystyle f(p^a)/p^0𝑑p^1𝑑p^2𝑑p^3}\right)^{1/2}`$
$`(\widehat{T}^{00})^{1/2}(N_\alpha N^\alpha )^{1/4}`$
This shows that (5) holds. It follows directly from the definitions that the inequality of (6) holds with $`C_2=1/3`$ even without restricting $`N_\alpha N^\alpha `$ to be bounded. Finally

$`T_{\alpha \beta }T^{\alpha \beta }`$ $`=\left({\displaystyle f(p^a)p_\alpha p_\beta /p^0𝑑p^1𝑑p^2𝑑p^3}\right)\left({\displaystyle f(q^a)q^\alpha q^\beta /q^0𝑑q^1𝑑q^2𝑑q^3}\right)`$
$`={\displaystyle f(p^a)f(q^a)(p^\alpha q_\alpha )^2/(p^0q^0)𝑑p^1𝑑p^2𝑑p^3𝑑q^1𝑑q^2𝑑q^3}`$
$`{\displaystyle f(p^a)f(q^a)(p^\alpha q_\alpha )/(p^0q^0)𝑑p^1𝑑p^2𝑑p^3𝑑q^1𝑑q^2𝑑q^3}`$
$`=N_\alpha N^\alpha `$
Thus the general theorems apply to give results on the dynamics of locally spatially homogeneous solutions of the Einstein-Vlasov system. It is known that Bianchi type IX solutions cannot expand forever, while models of the other types force the volume to be monotone. Thus we can say the following about inextendible non-vacuum spatially compact solutions of the Einstein-Vlasov system with Bianchi symmetry. If the Bianchi type is IX they have curvature singularities after finite proper time both in the past and in the future. For all other Bianchi types models which are expanding at some time have a curvature singularity at a finite time in the past and are future geodesically complete. The statement about geodesic completeness also holds in the vacuum case. The statement about curvature singularities, however, does not. In some cases there is a Cauchy horizon. This has been discussed in detail in . See also .

## 3 Spatially homogeneous solutions II: application of dynamical systems

We have now obtained a crude picture of the dynamics of spatially homogeneous cosmological models of the Einstein-matter equations for a variety of matter models. It is reasonable to hope that in the case of the Vlasov equation this can be considerably refined, in order to get a detailed picture of the asymptotics of the models near an initial singularity or in a phase of unlimited expansion. This has not yet been achieved in general but for certain cases results were obtained in . These extended theorems concerning the simpler case of massless particles obtained in . The models considered were of the simplest Bianchi types I, II and III. It was assumed that a further symmetry is present so that the spacetimes have a total of four Killing vectors. These are the so-called LRS models (locally rotationally symmetric). The reason for making this assumption is that then the Vlasov equation can be solved explicitly and only the Einstein equations remain to be handled. The equations can be reduced to a system of ODE’s in contrast to the general case, although the coefficients of the system involve one function which is not known explicitly. Fortunately it suffices to know certain qualitative features of this function in order to determine the asymptotic behaviour of the solutions of the ODE’s at early and late times. In fact to do this analysis it was also necessary to assume invariance under certain reflections, a piece of information which will be suppressed in the following for simplicity.

To describe the results it is useful to introduce the generalized Kasner exponents $`p_i`$. Each of the homogeneous hypersurfaces has an induced metric $`g_{ab}`$ and a second fundamental form $`k_{ab}`$. Let $`\lambda _i`$ denote the eigenvalues of the second fundamental form with respect to the metric. By definition this means that they are the solutions of the eigenvalue equation $`det(k_{ab}\lambda g_{ab})=0`$. Suppose now that $`\mathrm{tr}k`$ is non-zero. Then we can define $`p_i=\lambda _i/(_j\lambda _j)`$. The quantities $`p_i`$ are functions of $`t`$ and their sum is equal to one. The Kasner solution of the vacuum Einstein equations is given by

$$ds^2=dt^2+t^{2p_1}dx^2+t^{2p_2}dy^2+t^{2p_3}dz^2$$
(10)
for the Kasner exponents $`p_i`$, which are constants satisfying $`_jp_j=1`$ and $`_jp_j^2=1`$. These equations are called the first and second Kasner relations. In this case the generalized Kasner exponents are equal to the quantities $`p_i`$ in this metric form and this explains their name.

Consider a solution of the Einstein-Vlasov system of Bianchi type I which is LRS. Choose the time of the initial singularity to be $`t=0`$. The following statements hold. For each $`i`$ the quantity $`p_i(t)`$ converges to $`1/3`$ as $`t\mathrm{}`$. This is the value of the generalized Kasner exponents in a spatially flat FRW model. This means that the spacetime isotropizes at late times. At early times there are three possibilities. The first is that the $`p_i`$ are identically equal to $`1/3`$ at all times. This is the FRW case. The second is that they have limits $`(1/2,1/2,0)`$. The third, which is the generic case, is that they have the limits $`(2/3,2/3,1/3)`$. In this last case the second Kasner relation is satisfied asymptotically and so, in a certain sense, the solution with matter is approximated near the singularity by a vacuum solution.

In type III the initial singularity is similar to that in type I but in type II it is oscillatory. As $`t0`$ the generalized Kasner exponents approach both the values $`(2/3,2/3,1/3)`$ and $`(0,0,1)`$ as closely as desired in any neighbourhood of $`t=0`$. In the expanding direction the type II solution approaches an explicitly known self-similar dust spacetime whose generalized Kasner exponents have the values $`(3/8,3/8,1/4)`$. In type III the generalized Kasner exponents approach $`(0,0,1)`$ in the expanding direction. It is possible to get more detailed information about the asymptotics of these models in the limits $`t0`$ and $`t\mathrm{}`$. The proofs of all these statements make use of the fact that for ODE’s the highly developed theory of dynamical systems is available. Getting out the finer features of the expanding phase of Bianchi III models is more delicate that the other cases just discussed and was carried out using centre manifold theory in .

A key aspect of the work on the asymptotics of types I, II and III was writing the dynamical system in cleverly chosen variables. As we go to a singularity or to infinity in a phase of unlimited expansion the most obvious variables go to zero or infinity. If dimensionless variables can be found which, at least for some of the solutions, converge to finite non-zero limits in these regimes then this is a great help in analysing the asymptotics. Often the original dynamical system can be extended to a smooth dynamical system on a compact region. This avoids the problem, found with many choices of variables, that the solutions expressed in the given variables run off to infinity in a way which is hard to control. The strategy of dimensionless variables and compactification has been carried much further in the case of spatially homogeneous solutions of the Einstein-Euler system. For an account of this see .

## 4 Static solutions

This section is concerned with static spherically symmetric solutions of the Einstein-Vlasov system. These may play a role in describing the long-time behaviour of solutions of the full dynamical equations and so they have a natural place in these lectures. Before coming to the Einstein-Vlasov system it is worth spending a little time thinking about the corresponding non-relativistic problem, where a lot more is known. There are results on static solutions of the Vlasov-Poisson system which are not spherically symmetric and stationary solutions which are not static . Nothing comparable has yet been done in the case of the Einstein-Vlasov system. This is a gap which should be filled. From now on we restrict consideration to the static spherically symmetric case.

There are two methods which have been used to construct static spherically symmetric solutions of the Vlasov-Poisson system. The first may be called the ODE method. In a spherically symmetric static spacetime there are two constants of motion of the particles, namely the energy $`E`$ and the modulus of the angular momentum $`L`$, which are useful in constructing solutions of the Vlasov equation. In fact $`E`$ and $`L`$, like any quantity conserved along geodesics, are solutions of the Vlasov equation. The same is true of any function $`\mathrm{\Phi }(E,L)`$ of these quantities. Jeans’ theorem says that in a spherically symmetric static solution of the Vlasov-Poisson system the distribution function is a function of $`E`$ and $`L`$. Thus a natural procedure is to make an ansatz for the distribution function by choosing a particular function $`\mathrm{\Phi }`$. The Einstein equations then reduce to a system of integrodifferential equations for the metric coefficients as functions of a radial coordinate. What remains to be done is to analyse the global properties of solutions of this system.

The second method is a variational one. The Vlasov-Poisson system can be expressed as an infinite-dimensional Hamiltonian system. (Cf. .) It is degenerate in the sense that instead of a symplectic structure there is only a Poisson structure. This leads to a large class of conserved quantities known as Casimir invariants. In a Hamiltonian system a minimum of the Hamiltonian is a time-independent solution of the equations of motion. This suggests a variational route to finding static solutions. When Casimir invariants are present there are more general possibilities. If $`C`$ is a Casimir invariant then a minimum of $`H+C`$ is a time-independent solution of the equations of motion. An advantage of this method is that apart from giving results on the existence of static solutions it can also provide information on the stability of the solutions obtained. This method has been applied extensively to the Vlasov-Poisson system by Guo and Rein (see  and references therein).

Returning to the Einstein-Vlasov system, there is a paper  where the energy-Casimir method has been applied but it seems to be much harder than the Vlasov-Poisson case and the results are much more limited. More straightforward is the ODE method. One cautionary note is in order. The direct analogue of Jeans’ theorem is not true in the case of the Einstein-Vlasov system. Counterexamples were constructed by Schaeffer . Nevertheless we can still assume a distribution function of the form $`\mathrm{\Phi }(E,L)`$ and procede from there. A theorem on global existence in the radial coordinate for a rather general choice of $`\mathrm{\Phi }`$ was obtained in . There is a difficulty concerning the physical relevance of these solutions. If we would like to use them to model globular clusters, for instance, then we would like to obtain configurations of finite total mass. The easiest way to prove this is if the spatial density has compact support. The general existence theorem does not give any information on this. In fact whether it is true or not depends on the choice of the function $`\mathrm{\Phi }`$ in quite a delicate way. A criterion for the finiteness of the mass in a large class of functions $`\mathrm{\Phi }`$ was given in .

All known static solutions of the Einstein-Euler system with a physically reasonable equation of state are spherically symmetric and the density is a monotone decreasing function of the radius. In the case of the Einstein-Vlasov system another kind of configuration is possible where the support of the density is a thick shell, i.e. the region between two concentric spheres. In order to achieve this the function $`\varphi `$ must depend on the angular momentum. If it only depends on the energy then the system is equivalent to a solution of the Einstein-Euler system with an equation of state which is in general not explicitly known. The existence of shell solutions was proved in .

## 5 Spherically symmetric collapse

An interesting situation to consider is that of an isolated system consisting of matter undergoing gravitational collapse. The traditional model for this, following Oppenheimer and Snyder, is the collapse of a homogeneous spherical cloud of dust. Unfortunately when inhomogeneities are introduced into the Oppenheimer-Snyder model they often lead to pathologies such as shell-crossing singularities. The advantage of using collisionless matter is that it avoids some of (and perhaps all) the problems associated with dust.

A natural first step towards understanding spherical collapse is to fully understand the case where there is no collapse. If we have only a small amount of matter then it is to be expected that its self-gravitation will not suffice to keep it together and that it will spread out and disperse to infinity. For collisionless matter this has been proved, as described in more detail later. For dust it is not true since even small amounts of matter can develop shell-crossing singularities. Even dust without gravitation can do so and so this effect has nothing to do with gravity at all.

Consider initial data for the Einstein-Vlasov system which are spherically symmetric and asymptotically flat. In a suitable coordinate system only data for the distribution function need be given since the metric can then be determined by solving equations on the hypersurfaces of constant time. This is a reflection of the familiar statement that there is no gravitational radiation in spherical symmetry. Now let us make the initial data small in the following sense. In the presence of fixed bounds on the extent of the support of the initial data for $`f`$ in position and velocity space we require the maximum of $`f`$ to be small. For small data it can be shown that the solution exists globally in a suitable time coordinate and, more importantly, that it is geodesically complete. Moreover, various quantities such as the energy density of the matter decay to zero as $`t\mathrm{}`$ . Thus it can be seen that for small data the solution disperses and the situation is completely under control.

What happens for large data? It is known that for large data a trapped surface (and presumably a black hole) can form . We might nevertheless get global existence in a singularity-avoiding time slicing like maximal or polar slicing. The latter is also sometimes called a Schwarzschild time coordinate. There is a theorem  which says that if a singularity forms in a solution of the spherically symmetric Einstein-Vlasov system then the first singularity (as measured in Schwarzschild time) occurs at the centre of symmetry. Note that the shell-crossing singularities of dust occur away from the centre. A corresponding result for maximal slicing has been proved in .

A general mathematical result on the behaviour of spherically symmetric asymptotically flat solutions of the Einstein-Vlasov system has not yet been obtained. In the absence of further analytical progress, attempts have been made to study the problem numerically. One theme which plays an important role is that of critical collapse. Suppose that we have a family of data depending on a parameter $`\lambda `$ for the spherically symmetric Einstein-Vlasov system which interpolates between weak and strong data, with $`\lambda =0`$ corresponding to data for flat space. For $`\lambda `$ sufficiently small the theorem already mentioned tells us that the matter disperses. For $`\lambda `$ sufficiently large we might expect collapse to a black hole and this is indeed seen numerically. More precisely, it is seen that for $`\lambda `$ smaller than a certain value $`\lambda _{}`$ the matter disperses while for $`\lambda >\lambda _{}`$ it collapses to a black hole. In general some of the matter falls into the black hole while some escapes. Let $`M(\lambda )`$ be the mass of the black hole formed when the initial data corresponds to the parameter value $`\lambda `$. If no black hole is formed $`M(\lambda )`$ is defined to be zero. One of the questions which comes up in the study of critical collapse is whether the function $`M(\lambda )`$ is continuous at $`\lambda _{}`$ or not.

In  numerical evidence was presented that $`M(\lambda )`$ is not continuous. In other words, the limit of $`M(\lambda )`$ as $`\lambda \lambda _{}`$ from above is strictly positive. This is different from what is found for some other matter models, such as the massless scalar field. Olabarrieta and Choptuik  confirmed this finding and were able to present a more detailed picture of what happens. It is convenient for the numerical calculations to take initial data where there are no particles at the centre and no particles on purely radial orbits. In that case as long as the solution remains regular no particle can reach the centre due to conservation of angular momentum. Thus we have a dynamical configuration of matter with a hole in the middle. This allows difficulties with the singularity of polar coordinates at $`r=0`$ to be avoided. It is found in  that the solution evolves towards an unstable static shell solution before turning away again and dispersing or collapsing. The mass of the shell solution sets the mass gap in the graph of $`M(\lambda )`$. The connection between the shell solutions observed numerically in collapse calculations and those whose existence has been shown rigorously is not clear.

## 6 Isotropic singularities and Fuchsian methods

In the last section results for certain asymptotically flat solutions of the Einstein-Vlasov system were described. These are spherically symmetric and hence have three Killing vectors. There seem to be no other symmetry assumptions on asymptotically flat spacetimes which can be usefully studied at the present time. The obvious symmetry class which comes to mind is axial symmetry. In that case, however, there is only one Killing vector, which is very little, and even that has fixed points, which leads to singularities of the equations obtained when the symmetry is factored out. At the moment it seems to offer no advantage over the general case. In the case of spacetimes evolving from data on a compact Cauchy surface (cosmological spacetimes) there is a variety of interesting symmetry types with two or three Killing vectors and a number of papers on solutions of the Einstein-Vlasov system with these symmetries. They will not be reviewed here since some choices had to be made in order to limit the volume of the lectures. A good point of entry into the literature is .

There is one mathematical result on the Einstein-Vlasov system which does not require any symmetry assumptions and it will be the subject of the remainder of this section. It concerns solutions of the Einstein-Vlasov system with massless particles and it would be interesting to know if analogous results hold for massive particles. The idea is to construct large classes of solutions of the equations whose singularities have a particular structure, the isotropic singularities.

Given a spacetime with a foliation by spacelike hypersurfaces we can define the generalized Kasner exponents as in section 3 in terms of the eigenvalues of the second fundamental form. In the inhomogeneous case these are functions $`p_i`$ on spacetime which in general depend on both the time and space coordinates. The condition for an isotropic singularity (at least intuitively) is that all generalized Kasner exponents should tend to $`1/3`$ as the singularity is approached. Thus the solution looks like an isotropic FRW model near the singularity. The actual definition used in the theorem is a different one. A spatially flat FRW model is conformally flat. In the definition of an isotropic singularity it is assumed that the given metric is conformal to a metric which is smooth at the singularity.

The requirement of an isotropic singularity is a restriction on the spacetimes considered but it has been proved by Anguige  that there is a very large class of solutions of the Einstein-Vlasov system with isotropic singularities. In particular, he does not have to make any symmetry assumptions. The solutions can be parametrized by certain data on the singularity which can be given freely. The method of proving this is to use Fuchsian methods (although Anguige does not use this terminology). This technique is of wider importance in studying singularities of solutions of the Einstein-matter equations and will now be discussed in a more general context.

Let a system of partial differential equations with smooth coefficients be given and suppose we would like to investigate the existence of solutions which become singular on a certain hypersurface. For simplicity assume that this is the coordinate hypersurface $`t=0`$. Suppose that in some way it was possible to guess the asymptotic behaviour of the solutions in the approach to the singularity. This might be done by studying explicit solutions or by using trial and error to get a formally consistent asymptotic expansion. Then express the solution $`u`$ to be constructed in terms of an explicit function $`u_0`$ having the expected asymptotics near the singularity and a remainder $`v`$ which is expected to be regular and vanish at $`t=0`$. Now rewrite the original equation for the unknown $`u`$ as an equation for $`v`$ whose coefficients depend on $`u_0`$. Since $`u_0`$ is singular it is to be expected that the equation for $`v`$ is singular at $`t=0`$. Thus the problem of finding a singular solution of a regular equation has been replaced by that of finding a regular solution of a singular equation.

In favourable cases the singular equation obtained by this method is a Fuchsian equation of the form

$$t_tv+N(x)v=tf(t,x,v,v_x)$$
(11)
where the matrix-valued function $`N`$ has some positivity property. There are theorems which guarantee that an equation of this kind has a unique solution $`v`$ which is regular and vanishes at $`t=0`$. One theorem of this kind was proved in , where it was applied to study singularities in Gowdy spacetimes. Since then there have been a number of other applications. (See , section 6.2, for more information on this.) Unfortunately the Einstein-Vlasov system does not fit into the framework of this theorem due to the fact that it is an integrodifferential equation rather than a differential equation. For this reason Anguige had to prove his theorem by doing a direct iteration.

## 7 Outlook

In these lectures a selection of work on the Einstein-Vlasov system has been surveyed. Although the results were only discussed on a very general level, without getting into details, it was still necessary to leave out a lot of interesting topics. Some of them were mentioned briefly, some not at all. It should be clear that this is an area of research where there are many open problems and many promising directions to be explored. The references given here should provide a good starting point for those wanting to follow this road.