text:
# Post-collapse dynamics of self-gravitating Brownian particles and bacterial populations

## 1 Introduction

Self-gravitating systems such as globular clusters and elliptical galaxies constitute a Hamiltonian system of particles in interaction that can be supposed isolated in a first approximation . Since energy is conserved, the proper description of stellar systems is the microcanonical ensemble . The dynamical evolution of elliptical galaxies is governed by the Vlasov-Poisson system which corresponds to a collisionless regime. On the other hand, the kinetic theory of globular clusters is based on the Landau-Poisson system (or the orbit averaged Fokker-Planck equation) which describes a collisional evolution. These equations conserve mass and energy. Furthermore, the Landau equation increases the Boltzmann entropy (H-theorem) due to stellar encounters. These equations have been studied for a long time in the astrophysical literature and a relatively good physical understanding has now been achieved . In particular, globular clusters can experience core collapse  related to the “gravothermal catastrophe” concept .

For systems with long-range interactions, statistical ensembles are not equivalent . Therefore, it is of conceptual interest to compare the microcanonical evolution of stellar systems to a canonical one in order to emphasize the analogies and the differences. This can be achieved by considering a gas of self-gravitating Brownian particles  subject to a friction originated from the presence of an inert gas and to a stochastic force (modeling turbulent fluctuations, collisions,…). This system has a rigorous canonical structure where the temperature $`T`$ measures the strength of the stochastic force. Thus, we can precisely check the thermodynamical predictions of Kiessling  and Chavanis  obtained in the canonical ensemble. In the mean-field approximation, the self-gravitating Brownian gas model is described by the Kramers-Poisson system. In a strong friction limit (or for large times) it reduces to the Smoluchowski-Poisson system. These equations conserve mass and decrease the Boltzmann free energy . They possess a rich physical and mathematical structure and can lead to a situation of “isothermal collapse”, which is the canonical version of the gravothermal catastrophe. These equations have not been considered by astrophysicists because the canonical ensemble is not the correct description of stellar systems and usual astrophysical bodies do not experience a friction with a gas (except dust particles in the solar nebula ). Yet, it is clear that the self-gravitating Brownian gas model is of considerable conceptual interest in statistical mechanics to understand the strange thermodynamics of systems with long-range interactions and the inequivalence of statistical ensembles. In addition, it provides one of the first model of stochastic particles with long-range interactions, thereby extending the classical Einstein-Smoluchowski model  to a more general context .

In addition, it turns out that the same type of equations occur in biology in relation with the chemotactic aggregation of bacterial populations . A general model of chemotactic aggregation has been proposed by Keller & Segel  in the form of two coupled partial differential equations. In some approximation, this model reduces to the Smoluchowski-Poisson system . Therefore, there exists an isomorphism between self-gravitating Brownian particles and bacterial colonies. Non-local drift-diffusion equations analogous to the Smoluchowski-Poisson system have also been introduced in two-dimensional hydrodynamics in relation with the formation of large-scale vortices such as Jupiter’s great red spot . These analogies give further physical interest to our Brownian model.

In a recent series of papers , we have studied the dynamics and thermodynamics of self-gravitating Brownian particles confined within a spherical box of radius $`R`$ in a space of dimension $`D`$. In these works, we focused on the pre-collapse regime. In the canonical situation (fixed temperature $`T`$) we showed that a critical temperature $`T_c`$ exists below which the system undergoes a gravitational collapse leading to a finite time singularity at $`t=t_{coll}`$. For $`tt_{coll}`$, the evolution is self-similar in the sense that the density profile evolves as $`\rho (r,t)=\rho _0(t)f(r/r_0(t))`$ where $`f(x)`$ is independent on time. For $`x+\mathrm{}`$, $`f(x)x^\alpha `$ with $`\alpha =2`$. The central density increases as $`\rho _0(t_{coll}t)^1`$ and the core radius decreases as $`r_0(t_{coll}t)^{1/\alpha }`$. For $`T=0`$, the exponent is $`\alpha =2D/(D+2)`$. The scaling profiles can be calculated analytically. These results can be compared with those of Penston  who considered an isothermal collapse modelled by the Euler-Jeans equations. We also introduced a microcanonical description of Brownian particles by letting the temperature $`T(t)`$ evolve in time so as to conserve energy. This can provide a simplified model for the violent relaxation of collisionless stellar systems  or, simply, a numerical algorithm  to determine what is the maximum entropy state at fixed $`E`$ and $`M`$. In the microcanonical situation, there exists a critical energy $`E_c`$ (Antonov energy) below which the system collapses. For $`tt_{coll}`$, there exists a pseudo-scaling regime where $`\alpha `$ passes very slowly from $`\alpha _{max}=2.21\mathrm{}`$ to $`\alpha =2`$. Numerical simulations suggest that $`T(t)`$ remains finite at $`t=t_{coll}`$ so that the true scaling regime corresponds to $`\alpha =2`$, as in the canonical situation .

What happens after $`t_{coll}`$ ? By investigating the case $`T=0`$, we found in  that the evolution continues in the post-collapse regime with the formation of a Dirac peak accreting more and more mass as $`M(t)(tt_{coll})^{D/2}`$ while the density outside the peak evolves self-similarly with decreasing central density $`\rho _0(tt_{coll})^1`$ and increasing core radius $`r_0(tt_{coll})^{\frac{D+2}{2D}}`$. Our aim in this paper is to investigate the post-collapse regime for $`T0`$ in both canonical and microcanonical ensembles. In Sec. 2, we emphasize the analogy between self-gravitating Brownian particles and bacterial populations. In the following, we shall use the astrophysical terminology but we stress that our results apply equally well to biology where the chemotactic model has more concrete physical applications. In Sec. 3, we set the notations and recall the main results concerning the pre-collapse dynamics. In Sec. 4, we study the post-collapse dynamics at $`T=0`$ by a method different from , which can be generalized at finite temperature. The post-collapse dynamics at $`T>0`$ is precisely considered in Sec. 5. In the canonical ensemble, we show that the system forms a Dirac peak whose mass increases as $`M(t)(tt_{coll})^{D/21}`$ while the density profile for $`r>0`$ expands self-similarly with $`\rho _0(tt_{coll})^1`$ and $`r_0(tt_{coll})^{1/2}`$. For large times, the system is made of a Dirac peak of mass $`M`$ surrounded by a light gas of Brownian particles (with negligible self-interaction). Due to thermal motion, complete collapse takes an infinite time (contrary to the case $`T=0`$). For $`t+\mathrm{}`$, the mass contained in the Dirac peak increases as $`1M(t)/M\mathrm{exp}(\lambda t)`$ where $`\lambda `$ is the fundamental eigenvalue of a quantum problem. For $`T0`$, we find that $`\lambda =1/4T+c_D/T^{1/3}+\mathrm{}`$. In the microcanonical ensemble, the post-collapse regime is very pathological. The system tends to create a “Dirac peak of $`0^+`$ mass” surrounded by a uniform halo with infinite temperature. The central structure is reminiscent of a “binary star” containing a weak mass $`2mM=Nm`$ but a huge binding energy comparable to the potential energy of the whole cluster. Since our model is essentially mean-field, the physical formation of binaries is replaced by the type of structures mentioned above. The formation of a “Dirac peak” containing the whole mass in canonical ensemble and the formation of “binaries” in microcanonical ensemble are expected from thermodynamical considerations . These are the structures which provoke a divergence of free energy (at fixed mass and temperature) and entropy (at fixed mass and energy) respectively . All these analytical results are confirmed by numerical simulations of the Smoluchowski-Poisson system in Secs. 6 and 7. We were able in particular to “cross the singularity” at $`t=t_{coll}`$ and describe the post-collapse dynamics.

## 2 Analogy between self-gravitating Brownian particles and bacterial populations

### 2.1 Self-gravitating Brownian particles

We consider a system of self-gravitating Brownian particles described by the $`N`$ coupled stochastic equations

(1) $`{\displaystyle \frac{d𝐫_i}{dt}}=𝐯_i,{\displaystyle \frac{d𝐯_i}{dt}}=\xi 𝐯_i_iU(𝐫_1,\mathrm{},𝐫_N)+\sqrt{2D^{}}𝐑_i(t),`$
where $`\xi `$ is the friction coefficient, $`D^{}`$ is the diffusion coefficient and $`𝐑_i(t)`$ is a white noise satisfying $`𝐑_i(t)=\mathrm{𝟎}`$ and $`R_{a,i}(t)R_{b,j}(t^{})=\delta _{ij}\delta _{ab}\delta (tt^{})`$, where $`a,b=1,\mathrm{},D`$ refer to the coordinates of space and $`i,j=1,\mathrm{},N`$ to the particles. The particles interact via the potential $`U(𝐫_1,\mathrm{},𝐫_N)=_{i<j}u(𝐫_i𝐫_j)`$. In this paper, $`u(𝐫_i𝐫_j)`$ is the Newtonian binary potential in $`D`$ dimensions. The stochastic process (1) defines a toy model of gravitational dynamics which extends the classical Brownian model  to the case of stochastic particles in interaction. In this context, the friction is due to the presence of an inert gas and the stochastic force is due to classical Brownian motion, turbulence, or any other stochastic effect.

Starting from the $`N`$-body Fokker-Planck equation and using a mean-field approximation , we can derive the nonlocal Kramers equation

(2) $`{\displaystyle \frac{f}{t}}+𝐯{\displaystyle \frac{f}{𝐫}}+𝐅{\displaystyle \frac{f}{𝐯}}={\displaystyle \frac{}{𝐯}}\left(D^{}{\displaystyle \frac{f}{𝐯}}+\xi f𝐯\right),`$
where $`𝐅=\mathrm{\Phi }`$ is the smooth gravitational force felt by the particles. The gravitational potential $`\mathrm{\Phi }`$ is related to the density $`\rho =fd^D𝐯`$ by the Poisson equation

(3) $`\mathrm{\Delta }\mathrm{\Phi }=S_DG\rho ,`$
where $`S_D`$ is the surface of the unit $`D`$-dimensional sphere. Equation (2) can be considered as a generalized version of the Kramers-Chandrasekhar equation introduced in a homogeneous medium . In this work, the diffusion and the friction model stellar encounters in a simple stochastic framework. The condition that the Maxwell-Boltzmann distribution is a stationary solution of Eq. (2) leads to the Einstein relation $`\xi =D^{}\beta `$ where $`\beta =1/T`$ is the inverse temperature (we have included the mass of the particles and the Boltzmann constant in the definition of $`T`$). In our case, we do not assume that the medium is homogeneous, so that we have to solve the Kramers-Poisson system. This makes the study much more complicated than usual. Up to date, we do not know any astrophysical application of this model although there could be connexions with the process of planetesimal formation in the solar nebula . Whatever, this model is interesting to develop on a conceptual point of view because it possesses a rigorous thermodynamical structure and presents the same features as more realistic models (isothermal distributions, collapse, phase transitions,…). For this Brownian model, the relevant ensemble is the canonical one since the temperature $`T`$ is fixed. Therefore, the Kramers-Poisson system can be viewed as the canonical counterpart of the Landau-Poisson system. It is interesting to study these two models in parallel to illustrate dynamically the inequivalence of statistical ensembles for systems with long-range interactions.

To simplify the problem further, we shall consider the strong friction limit $`\xi +\mathrm{}`$, or equivalently the limit of large times $`t\xi ^1`$. In that approximation, we can neglect the inertia of the particles. Then, the coupled stochastic equations (1) simplify to

(4) $`\xi {\displaystyle \frac{d𝐫_i}{dt}}=_iU(𝐫_1,\mathrm{},𝐫_N)+\sqrt{2D^{}}𝐑_i(t).`$
Furthermore, to leading order, the velocity distribution is Maxwellian

(5) $`f(𝐫,𝐯,t)=\left({\displaystyle \frac{\beta }{2\pi }}\right)^{D/2}\rho (𝐫,t)e^{\beta \frac{v^2}{2}},`$
and the Kramers equation reduces to the Smoluchowski equation

(6) $`{\displaystyle \frac{\rho }{t}}=\left[{\displaystyle \frac{1}{\xi }}(T\rho +\rho \mathrm{\Phi })\right].`$
It can be shown that the Kramers equation decreases the Boltzmann free energy

(7) $`F[f]=ETS={\displaystyle f\frac{v^2}{2}d^D𝐫d^D𝐯}+{\displaystyle \frac{1}{2}}{\displaystyle \rho \mathrm{\Phi }d^D𝐫}+T{\displaystyle f\mathrm{ln}fd^D𝐫d^D𝐯},`$
i.e. $`\dot{F}0`$ and $`\dot{F}=0`$ at statistical equilibrium (canonical $`H`$-theorem). Similarly, the Smoluchowski equation decreases the free energy $`F[\rho ]`$ which is obtained from $`F[f]`$ by using the fact that the velocity distribution is Maxwellian in the strong friction limit. This leads to the classical expression

(8) $`F[\rho ]=T{\displaystyle \rho \mathrm{ln}\rho d^D𝐫}+{\displaystyle \frac{1}{2}}{\displaystyle \rho \mathrm{\Phi }d^D𝐫}.`$
The passage from the Kramers equation to the Smoluchowski equation in the strong friction limit is classical . It can also be obtained formally from a Chapman-Enskog expansion .

In order to prevent finite time singularities and infinite densities, we can consider a model of self-gravitating Brownian fermions, enforcing the constraint $`f\eta _0`$ (Pauli exclusion principle). The corresponding Kramers equation takes the form

(9) $`{\displaystyle \frac{f}{t}}+𝐯{\displaystyle \frac{f}{𝐫}}+𝐅{\displaystyle \frac{f}{𝐯}}={\displaystyle \frac{}{𝐯}}\left[D^{}{\displaystyle \frac{f}{𝐯}}+\xi f(\eta _0f)𝐯\right].`$
In the strong friction limit, we obtain a Smoluchowski equation of the form

(10) $`{\displaystyle \frac{\rho }{t}}=\left[{\displaystyle \frac{1}{\xi }}(p+\rho \mathrm{\Phi })\right],`$
where $`p(\rho )`$ is the equation of state of the Fermi gas. The fermionic Smoluchowski-Poisson system has been studied in . Generalized Kramers and Smoluchowski equations are introduced in .

### 2.2 The Keller-Segel model

The name chemotaxis refers to the motion of organisms (amoeba) induced by chemical signals (acrasin). In some cases, the biological organisms secrete a substance that has an attractive effect on the organisms themselves. Therefore, in addition to their diffusive motion, they move systematically along the gradient of concentration of the chemical they secrete (chemotactic flux). When attraction prevails over diffusion, the chemotaxis can trigger a self-accelerating process until a point at which aggregation takes place. This is the case for the slime mold Dictyostelium Discoideum and for the bacteria Escherichia coli .

A model of slime mold aggregation has been introduced by Keller & Segel  in the form of two coupled differential equations

(11)
$$\frac{\rho }{t}=(D_2\rho )(D_1c),$$

(12)
$$\frac{c}{t}=k(c)c+f(c)\rho +D_c\mathrm{\Delta }c.$$

In these equations $`\rho (𝐫,t)`$ is the concentration of amoebae and $`c(𝐫,t)`$ is the concentration of acrasin. Acrasin is produced by the amoebae at a rate $`f(c)`$. It can also be degraded at a rate $`k(c)`$. Acrasin diffuse according to Fick’s law with a diffusion coefficient $`D_c`$. Amoebae concentration changes as a result of an oriented chemotactic motion in the direction of a positive gradient of acrasin and a random motion analogous to diffusion. In Eq. (11), $`D_2(\rho ,c)`$ is the diffusion coefficient of the amoebae and $`D_1(\rho ,c)`$ is a measure of the strength of the influence of the acrasin gradient on the flow of amoebae. This chemotactic drift is the fundamental process in the problem.

A first simplification of the Keller-Segel model is provided by the system of equations

(13)
$$\frac{\rho }{t}=D\mathrm{\Delta }\rho \chi (\rho c),$$

(14)
$$\frac{c}{t}=D^{}\mathrm{\Delta }c+a\rho bc,$$

where the parameters are positive constants. An additional simplification, introduced by Jäger & Lauckhaus  consists in ignoring the time derivative in Eq. (14). This is valid in the case where the diffusion coefficient $`D^{}`$ is large. Taking also $`b=0`$, we obtain

(15)
$$\frac{\rho }{t}=D\mathrm{\Delta }\rho \chi (\rho c),$$

(16)
$$\mathrm{\Delta }c=\lambda \rho ,$$

where $`\lambda =a/D^{}`$. Clearly, these equations are isomorphic to the Smoluchowski-Poisson system (6)-(3) describing self-gravitating Brownian particles in a strong friction limit. In particular, the chemotactic flux plays the same role as the gravitational drift in the overdamped limit of the Brownian model. When chemotactic attraction prevails over diffusion, the system is unstable and the bacteria start to aggregate. This blow-up is similar to the collapse of self-gravitating systems in a canonical situation. We note that in the Keller-Segel model, the diffusion coefficient can depend on the density, leading to anomalous diffusion. Such a situation is considered in  where the nonlinear Smoluchowski-Poisson system is studied.

The Keller-Segel model ignore clumping and sticking effects. However, at the late stages of the blow-up, when the density of amoebae has reached high values, finite size effects and stickiness must clearly be taken into account. As a first step, we can propose  to replace the classical equation (15) by an equation of the form

(17)
$$\frac{\rho }{t}=D\mathrm{\Delta }\rho \chi (\rho (\sigma _0\rho )c),$$

which enforces a limitation $`\rho (𝐫,t)<\sigma _0`$ on the maximum density of amoebae. This is the counterpart of the model of self-gravitating Brownian fermions . These type of non-local Fokker-Planck equations also occur in 2D hydrodynamics and astrophysics in relation with the formation of large-scale vortices and galaxies . Their systematic study is clearly of broad interest .

## 3 Collapse dynamics of self-gravitating Brownian particles

### 3.1 The Smoluchowski-Poisson system

At a given temperature $`T`$ controlling the diffusion coefficient, the density $`\rho (𝐫,t)`$ of self-gravitating Brownian particles satisfies the following coupled equations:

(18) $`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`\left[{\displaystyle \frac{1}{\xi }}(T\rho +\rho \mathrm{\Phi })\right],`$
(19) $`\mathrm{\Delta }\mathrm{\Phi }`$ $`=`$ $`S_DG\rho ,`$
where $`\mathrm{\Phi }`$ is the gravitational potential, and $`S_D`$ is the surface of the unit $`D`$-dimensional sphere.

From now on, we set $`M=R=G=\xi =1`$ and we restrict ourselves to spherically symmetric solutions. The equations of the problem become

(20) $`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`(T\rho +\rho \mathrm{\Phi }),`$
(21) $`\mathrm{\Delta }\mathrm{\Phi }`$ $`=`$ $`S_D\rho ,`$
with proper boundary conditions in order to impose a vanishing particle flux on the surface of the confining sphere. These read

(22)
$$\frac{\mathrm{\Phi }}{r}(0,t)=0,\mathrm{\Phi }(1)=\frac{1}{2D},T\frac{\rho }{r}(1)+\rho (1)=0,$$

for $`D>2`$. For $`D=2`$, we take $`\mathrm{\Phi }(1)=0`$ on the boundary. Integrating Eq. (21) once, we can rewrite the Smoluchowski-Poisson system in the form of a single integrodifferential equation

(23)
$$\frac{\rho }{t}=\frac{1}{r^{D1}}\frac{}{r}\left\{r^{D1}\left(T\frac{\rho }{r}+\frac{\rho }{r^{D1}}_0^r\rho (r^{})S_Dr^{{}_{}{}^{}D1}𝑑r^{}\right)\right\}.$$

The total energy is given as the sum of the kinetic and potential contributions

(24)
$$E=\frac{D}{2}T+\frac{1}{2}\rho \mathrm{\Phi }d^D𝐫.$$

The Smoluchowski-Poisson system is also equivalent to a single differential equation

(25)
$$\frac{M}{t}=T\left(\frac{^2M}{r^2}\frac{D1}{r}\frac{M}{r}\right)+\frac{1}{r^{D1}}M\frac{M}{r},$$

for the quantity

(26)
$$M(r,t)=_0^r\rho (r^{})S_Dr^{{}_{}{}^{}D1}𝑑r^{},$$

which represents the mass contained within the sphere of radius $`r`$. The appropriate boundary conditions are

(27)
$$M(0,t)=N_0(t),M(1,t)=1,$$

where $`N_0(t)=0`$, except if the density develops a condensed Dirac peak contribution at $`r=0`$, of total mass $`N_0(t)`$. It is also convenient to introduce the function $`s(r,t)=M(r,t)/r^D`$ satisfying

(28)
$$\frac{s}{t}=T\left(\frac{^2s}{r^2}+\frac{D+1}{r}\frac{s}{r}\right)+\left(r\frac{s}{r}+Ds\right)s.$$

### 3.2 Self-similar solutions of the Smoluchowski-Poisson system

In , we have shown that in the canonical ensemble (fixed $`T`$), the system undergoes gravitational collapse below a critical temperature $`T_c`$ depending on the dimension of space. The density develops a scaling profile, and the central density grows and diverges at a finite time $`t_{coll}`$. The case $`D=2`$ was extensively studied in  and turns out to be very peculiar. Throughout this paper, we restrain ourselves to the more generic case $`D>2`$, although other dimensions play a special role as far as static properties are concerned (see ).

We look for self-similar solutions of the form

(29)
$$\rho (r,t)=\rho _0(t)f\left(\frac{r}{r_0(t)}\right),r_0=\left(\frac{T}{\rho _0}\right)^{1/2},$$

where the King’s radius $`r_0`$ defines the size of the dense core . In terms of the mass profile, we have

(30)
$$M(r,t)=M_0(t)g\left(\frac{r}{r_0(t)}\right),\mathrm{with}M_0(t)=\rho _0r_0^D,$$

and

(31)
$$g(x)=S_D_0^xf(x^{})x^{{}_{}{}^{}D1}𝑑x^{}.$$

In terms of the function $`s`$, we have

(32)
$$s(r,t)=\rho _0(t)S\left(\frac{r}{r_0(t)}\right),\mathrm{with}S(x)=\frac{g(x)}{x^D}.$$

Substituting the ansatz (32) into Eq. (28), we find that

(33)
$$\frac{d\rho _0}{dt}S(x)\frac{\rho _0}{r_0}\frac{dr_0}{dt}xS^{}(x)=\rho _0^2\left(S^{\prime \prime }(x)+\frac{D+1}{x}S^{}(x)+xS(x)S^{}(x)+DS(x)^2\right),$$

where we have set $`x=r/r_0`$. The variables of position and time separate provided that $`\rho _0^2d\rho _0/dt`$ is a constant that we arbitrarily set equal to 2. After time integration, this leads to

(34)
$$\rho _0(t)=\frac{1}{2}(t_{coll}t)^1,$$

so that the central density becomes infinite in a finite time $`t_{coll}`$. The scaling equation now reads

(35)
$$2S+xS^{}=S^{\prime \prime }+\frac{D+1}{x}S^{}+S(xS^{}+DS).$$

The scaling solution of Eq. (35) was obtained analytically in  and reads

(36)
$$S(x)=\frac{4}{D2+x^2},$$

which decays with an exponent $`\alpha =2`$. This leads to

(37)
$$f(x)=\frac{4(D2)}{S_D}\frac{x^2+D}{(D2+x^2)^2},g(x)=\frac{4x^D}{D2+x^2}.$$

Note finally that within the core radius $`r_0`$, the total mass in fact vanishes as $`tt_{coll}`$. Indeed, from Eq. (30), we obtain

(38)
$$M(r_0(t),t)\rho _0(t)r_0^D(t)T^{D/2}(t_{coll}t)^{D/21}.$$

Therefore, the collapse does not create a Dirac peak (“black hole”).

In , we have also studied the collapse dynamics at $`T=0`$ for which we obtained

(39)
$$\rho _0(t)S_D^1(t_{coll}t)^1,$$

as previously, but the core radius is not given anymore by the King’s radius which vanishes for $`T=0`$. Instead, we find

(40)
$$r_0\rho _0^{1/\alpha },$$

with

(41)
$$\alpha =\frac{2D}{D+2}.$$

The scaling function $`S(x)`$ is only known implicitly

(42)
$$\left[\frac{2}{D+2}S(x)\right]^{\frac{D}{D+2}}=Kx^{\frac{2D}{D+2}}S(x),$$

where $`K`$ is a known constant (see  for details), $`S(0)=\frac{2}{D+2}`$, and the large $`x`$ asymptotics $`S(x)f(x)x^\alpha `$. The mass within the core radius is now

(43)
$$M(r_0(t),t)\rho _0(t)r_0^D(t)(t_{coll}t)^{D/2},$$

and it again tends to zero as $`tt_{coll}`$. Comparing Eq. (38) and Eq. (43) suggests that if the temperature is very small, an apparent scaling regime corresponding to the $`T=0`$ case will hold up to a cross-over time $`t_{}`$, with

(44)
$$t_{coll}t_{}T^{D/2}.$$

Above $`t_{}`$, the $`T0`$ scaling ultimately prevails.

## 4 Post-collapse dynamics at $`T=0`$

So far, all studies concerning the collapse dynamics of self-gravitating Brownian particles have concentrated on the time period $`tt_{coll}`$. A natural question arises: what is happening for $`t>t_{coll}`$ ? The first possible scenario is that the state reached at $`t=t_{coll}`$ is in fact a stationary state. However, its is easy to check (see ) that this is absolutely not the case. In addition, the preceding study leads to a sort of paradox . Indeed, we know that the statistical equilibrium state in the canonical ensemble is a Dirac peak . This is not the structure that forms at $`t=t_{coll}`$. This structure is singular at the origin ($`\rho r^2`$) but different from a Dirac peak (in particular the central mass is zero). This means that the evolution must continue after $`t_{coll}`$. In particular, we will show that the Dirac peak predicted by statistical mechanics forms in the post-collapse regime.

The scenario that we are now exploring is the following. A central Dirac peak containing a mass $`N_0(t)`$ emerges at $`t>t_{coll}`$, whereas the density for $`r>0`$ satisfies a scaling relation of the form

(45)
$$\rho (r,t)=\rho _0(t)f\left(\frac{r}{r_0(t)}\right),$$

where $`\rho _0(t)`$ is now decreasing with time (starting from $`\rho _0(t=t_{coll})+\mathrm{}`$) and $`r_0(t)`$ grows with time (starting from $`r_0(t=t_{coll})=0`$). As time increases, the residual mass for $`r>0`$ is progressively swallowed by the dense core made of particles which have fallen on each other. It is the purpose of the rest of this paper to show that this scenario actually holds, as well as to obtain analytical and numerical results illustrating this final collapse stage.

In this section, we present an alternative treatment to that of , where this scenario was analytically shown to hold at $`T=0`$. This new approach is a good introduction to the general $`T0`$ case which is studied in the next section. We refer the reader to  for an explicit solution of the $`T=0`$ post-collapse regime, which we found, leads to a central peak containing all the mass in a finite time $`t_{end}`$.

For $`T=0`$, the dynamical equation for the integrated mass $`M(r,t)`$ reads

(46)
$$\frac{M}{t}=\frac{1}{r^{(D1)}}M\frac{M}{r},$$

with boundary conditions

(47)
$$M(0,t)=N_0(t),M(1,t)=1.$$

We define $`\rho _0`$ such that for small $`r`$

(48)
$$M(r,t)N_0(t)=\rho _0(t)\frac{r^D}{D}+\mathrm{}$$

Up to the geometrical factor $`S_D^1`$, $`\rho _0(t)`$ is the central residual density (the residual density is defined as the density after the central peak has been subtracted). For $`r=0`$, Eq. (46) leads to the evolution equation for $`N_0`$

(49)
$$\frac{dN_0}{dt}=\rho _0N_0.$$

As $`N_0(t)=0`$ for $`tt_{coll}`$, and since this equation is a first order differential equation, it looks like $`N_0(t)`$ should remain zero for $`t>t_{coll}`$ as well. However, since $`\rho _0(t_{coll})=+\mathrm{}`$, there is mathematically speaking no global solution for this equation and non zero values for $`N_0(t)`$ can emerge from Eq. (49), as will soon become clear.

We then define

(50)
$$s(r,t)=\frac{M(r,t)N_0(t)}{r^D},$$

which satisfies

(51)
$$\frac{s}{t}=\left(r\frac{s}{r}+Ds\right)s+\frac{N_0}{r^D}\left(r\frac{s}{r}+Ds\rho _0\right).$$

By definition, we have also $`s(0,t)=\rho _0(t)/D`$.

We now look for a scaling solution of the form

(52)
$$s(r,t)=\rho _0(t)S\left(\frac{r}{r_0(t)}\right),$$

with $`S(0)=D^1`$ and

(53)
$$\rho _0(t)=r_0(t)^\alpha ,$$

where $`r_0`$ is thus defined without ambiguity. Inserting this scaling ansatz in Eq. (51), and defining the scaling variable $`x=r/r_0`$, we find

(54)
$$\frac{1}{\alpha \rho _0^2}\frac{d\rho _0}{dt}\left(\alpha S+xS^{}\right)=S(DS+xS^{})+\frac{N_0}{\rho _0r_0^D}\frac{1}{x^D}(DS+xS^{}1).$$

Imposing scaling, we find that both time dependent coefficients appearing Eq. (54) should be in fact constant. We thus define a constant $`\mu `$ such that

(55)
$$N_0=\mu \rho _0r_0^D,$$

and set

(56)
$$\frac{1}{\alpha \rho _0^2}\frac{d\rho _0}{dt}=\kappa ,$$

with $`\kappa >0`$, as the central residual density is expected to decrease. Equation (56) implies that $`\rho _0(tt_{coll})^1`$, which along with Eq. (53) implies that $`N_0(tt_{coll})^{D/\alpha 1}`$. We thus find a power law behavior for $`N_0`$, which in order to be compatible with Eq. (49), leads to

(57)
$$\rho _0(t)=\left(\frac{D}{\alpha }1\right)(tt_{coll})^1,$$

and then to

(58)
$$\kappa =\frac{1}{D\alpha }.$$

We end up with the scaling equation

(59)
$$\frac{1}{D\alpha }\left(\alpha S+xS^{}\right)+S(DS+xS^{})+\mu x^D(DS+xS^{}1)=0.$$

From Eq. (59), we find that the large $`x`$ asymptotics of $`S`$ is $`S(x)x^\alpha `$. In a short finite time after $`t_{coll}`$, it is clear that the large distance behavior of the density profile ($`rr_0`$) cannot dramatically change. We deduce that the decay of $`S`$ should match the behavior for time slightly less than $`t_{coll}`$ for which $`S(x)x^{\frac{2D}{D+2}}`$. Hence the value of $`\alpha `$ should remain unchanged before and after $`t_{coll}`$. Finally, we obtain the following exact behaviors for short time after $`t_{coll}`$:

(60) $`\rho _0(t)`$ $`=`$ $`{\displaystyle \frac{D}{2}}(tt_{coll})^1,`$
(61) $`r_0(t)`$ $`=`$ $`\left({\displaystyle \frac{2}{D}}\right)^{\frac{D+2}{2D}}(tt_{coll})^{\frac{D+2}{2D}},`$
(62) $`N_0(t)`$ $`=`$ $`\mu \left({\displaystyle \frac{2}{D}}\right)^{\frac{D}{2}}(tt_{coll})^{\frac{D}{2}}.`$
We note the remarkable result that the central residual density $`\rho (0,t)=S_D^1\rho _0(t)`$ displays a universal behavior just after $`t_{coll}`$, a result already obtained in . Moreover, we find that $`N_0(t)`$ has the same form as the mass found within a sphere of radius $`r_0(t)`$ below $`t_{coll}`$, given in Eq. (43).

Moreover, the scaling function $`S`$ satisfies

(63)
$$\frac{D+2}{D^2}\left(\frac{2D}{D+2}S+xS^{}\right)+S(DS+xS^{})+\mu x^D(DS+xS^{}1)=0.$$

The constant $`\mu `$ is determined by imposing that the large $`r`$ behavior of $`s(r,t)`$ (or $`\rho (r,t)`$) exactly matches (not simply proportional) that obtained below $`t_{coll}`$, which depends on the shape of the initial condition as shown in . Equation (63) can be solved implicitly by looking for solutions of the form $`x^D=z[S(x)]`$. After cumbersome but straightforward calculations, we obtain the implicit form

(64)
$$1+\frac{x^D}{\mu }S(x)=\left[1+\frac{x^D}{\mu }\left(S(x)+\frac{2}{D^2}\right)\right]^{\frac{D}{D+2}},$$

which coincides with the implicit solution given in . Note that $`S(x)`$ is a function of $`x^D`$. We check that the above result indeed leads to $`S(0)=D^1`$, and to the large $`x`$ asymptotics

(65)
$$S(x)\mu ^{\frac{2}{D+2}}\left(\frac{2}{D^2}\right)^{\frac{D}{D+2}}x^{\frac{2D}{D+2}}.$$

Note finally that for $`T=0`$, $`N_0`$ saturates to 1 in a finite time, corresponding to the deterministic collapse of the outer mass shell initially at $`r=1`$. Indeed, using Gauss’ theorem, the position of a particle initially at $`r(t=0)=1`$ satisfies

(66)
$$\frac{dr}{dt}=r^{(D1)}.$$

The position of the outer shell is then

(67)
$$r(t)=(1Dt)^{1/D},$$

which vanishes for $`t_{end}=D^1`$.

## 5 Post-collapse dynamics at $`T>0`$

### 5.1 Scaling regime

In the more general case $`T0`$, we will proceed in a very similar way as in the previous section. We define again,

(68)
$$s(r,t)=\frac{M(r,t)N_0(t)}{r^D},$$

where $`N_0`$ still satisfies

(69)
$$\frac{dN_0}{dt}=\rho _0N_0.$$

We now obtain

(70)
$$\frac{s}{t}=T\left(\frac{^2s}{r^2}+\frac{D+1}{r}\frac{s}{r}\right)+\left(r\frac{s}{r}+Ds\right)s+\frac{N_0}{r^D}\left(r\frac{s}{r}+Ds\rho _0\right).$$

By definition, we have again $`s(0,t)=\rho _0(t)/D`$.

We look for a scaling solution of the form

(71)
$$s(r,t)=\rho _0(t)S\left(\frac{r}{r_0(t)}\right),$$

with $`S(0)=D^1`$. As before, we define the King’s radius by

(72)
$$r_0=\left(\frac{T}{\rho _0}\right)^{1/2}.$$

For $`t<t_{coll}`$, we had $`s(r,t)4Tr^2`$ (or $`S(x)4x^2`$). In a very short time after $`t_{coll}`$, this property should be preserved, which implies that the post-collapse scaling function should also behave as

(73)
$$S(x)4x^2,$$

for large $`x`$. Inserting the scaling $`ansatz`$ in Eq. (70), we obtain

(74)
$$\frac{1}{2\rho _0^2}\frac{d\rho _0}{dt}\left(2S+xS^{}\right)=S^{\prime \prime }+\frac{D+1}{x}S^{}+S(DS+xS^{})+\frac{N_0}{\rho _0r_0^D}\frac{1}{x^D}(DS+xS^{}1).$$

Again, this equation should be time independent for scaling to hold, which implies that there exist two constants $`\mu `$ and $`\kappa `$ such that

(75)
$$N_0=\mu \rho _0r_0^D,$$

and

(76)
$$\frac{1}{2\rho _0^2}\frac{d\rho _0}{dt}=\kappa ,$$

with $`\kappa >0`$, as the central residual density is again expected to decrease. Equation (76) implies that $`\rho _0(tt_{coll})^1`$, and then that $`N_0(tt_{coll})^{D/21}`$. We thus find a power law behavior for $`N_0`$, which in order to be compatible with Eq. (69), leads to the universal behavior

(77)
$$\rho _0(t)=\left(\frac{D}{2}1\right)(tt_{coll})^1,$$

and then to

(78)
$$\kappa =\frac{1}{D2}.$$

We end up with the scaling equation

(79)
$$\frac{1}{D2}\left(2S+xS^{}\right)+S^{\prime \prime }+\frac{D+1}{x}S^{}+S(DS+xS^{})+\mu x^D(DS+xS^{}1)=0,$$

where $`\mu `$ has to be chosen so that $`S(x)`$ satisfies the condition of Eq. (73). Its value will be determined numerically in Sec. 6 for $`D=3`$. Note that for small $`x`$, the pre-collapse scaling function satisfies $`S(x)S(0)x^2`$, whereas the post-collapse scaling function behaves as

(80)
$$S(x)S(0)x^D.$$

However, contrary to the $`T=0`$ case, $`S(x)`$ is not purely a function of $`x^D`$.

Finally, we find that the weight of the central peak has a universal behavior for short time after $`t_{coll}`$

(81)
$$N_0(t)=\mu \left(\frac{2}{D2}\right)^{D/21}T^{D/2}(tt_{coll})^{D/21}.$$

Note that $`N_0(t)`$ behaves in a very similar manner to the mass within a sphere of radius $`r_0`$ below $`t_{coll}`$, shown in Eq. (38). In addition, comparing Eq. (81) and Eq. (62), we can define again a post-collapse cross-over time between the $`T0`$ and $`T=0`$ regimes

(82)
$$t_{}t_{coll}T^{D/2},$$

which is similar to the definition of Eq. (44).

### 5.2 Large time limit

Contrary to the $`T=0`$ case, the complete collapse does not take place in a finite time as thermal fluctuations always allow for some particle to escape the central strongly attractive potential. In order to illustrate this point, and obtain more analytic insight on this matter, we will place ourselves in the extreme situation where almost all the mass has collapsed ($`N_01`$), and only an infinitesimal amount remains in the residual profile.

In this limit, the residual density $`\rho (r,t)`$ satisfies the Fokker-Planck equation

(83)
$$\frac{\rho }{t}=T\left(\frac{^2\rho }{r^2}+\frac{D1}{r}\frac{\rho }{r}\right)+\frac{1}{r^{D1}}\frac{\rho }{r},$$

with boundary condition

(84)
$$T\frac{\rho }{r}(1,t)+\rho (1,t)=0.$$

The problem indeed reduces to the study of a very light gas (i.e. with negligible self-interaction) of Brownian particles submitted to the gravitational force $`𝐅=(GM/r^{D1})𝐞_r`$ of a central unit mass. Alternatively, this can also be seen as the probability distribution evolution equation of a system of two Brownian particles moving in their mutual gravitation field.

Equation (83) can be re-expressed as a Schrödinger equation (in imaginary time), thus involving a self-adjoint operator (see Appendix A). The large time behavior is dominated by the first eigenstate. Coming back to the notation of Eq. (83), we find that

(85)
$$\rho (r,t)\mathrm{e}^{\lambda t}\psi (r),$$

where $`\psi `$ satisfies the eigenequation

(86)
$$\lambda \psi (r)=T\left(\psi ^{\prime \prime }+\frac{D1}{r}\psi ^{}\right)+\frac{1}{r^{D1}}\psi ^{},$$

and the same boundary condition as $`\rho `$, i.e.

(87)
$$T\psi ^{}(1)+\psi (1)=0.$$

The eigenvalue $`\lambda `$ will also control the large time behavior of $`\rho _0`$ and $`N_0`$ as Eq. (69) and Eq. (85) both imply that

(88)
$$1N_0(t)\frac{\rho _0(t)}{\lambda }\mathrm{e}^{\lambda t}.$$

We did not succeed in solving analytically the above eigenequation, and for a given temperature, this has to be solved numerically. However, in the limit of very small temperature, we can apply techniques reminiscent from semiclassical analysis in quantum mechanics ($`T\mathrm{}`$). We now assume $`T`$ very small and define $`\varphi `$ such that

(89)
$$\psi (r)=\mathrm{e}^{\frac{\varphi (r)}{T}}.$$

The function $`h=\varphi ^{}`$ satisfies the following non-linear first order differential equation

(90)
$$T\left(h^{}+\frac{D1}{r}h\right)+\frac{h}{r^{D1}}h^2=\lambda T,$$

with the simple boundary condition

(91)
$$h(1)=1.$$

In the limit $`T0`$, the term proportional to $`T`$ in the left-hand side of Eq. (90) can a priori be discarded leading to

(92)
$$h(r)=\frac{2\lambda Tr^{D1}}{1+\sqrt{14\lambda Tr^{2(D1)}}}.$$

If $`4\lambda T<1`$, the above expression is a valid perturbative solution also at $`r=1`$, but cannot satisfies the constraint $`h(1)=1`$. Hence, we conclude that in the limit of small temperature $`4\lambda T1`$, so that the above expression is only valid for $`r`$ not to close to $`r=1`$. The above argument also suggests that $`\lambda T`$ is of order unity and we write

(93)
$$\lambda T=\frac{1}{4}+\mu ^2.$$

To understand how the boundary condition Eq. (91) can be in fact satisfied, one has to come back to Eq. (90), which for $`r=1`$, shows that $`h^{}(1)\lambda T^11`$. This implies that the term $`Th^{}`$ cannot be neglected near $`r=1`$, and that $`h`$ varies in a noticeable way on a length scale from 1 of order $`T`$.

This suggests to define

(94)
$$z(x)=h(1Tx)$$

which satisfies (at order 0 in $`T`$)

(95)
$$z^{}+zz^2=\frac{1}{4}+\mu ^2,$$

and $`z(0)=1`$. This equation has the unique solution

(96)
$$z(x)=\frac{1}{2}+\mu \frac{12\mu \mathrm{tan}(\mu x)}{2\mu +\mathrm{tan}(\mu x)}.$$

For large $`x`$, this function only has a sensible behavior for $`\mu =0`$, which shows that

(97)
$$\underset{T0}{lim}\lambda T=\frac{1}{4},$$

and that Eq. (92) is in fact valid for $`1r(4\lambda T1)0`$. To leading order, we find

(98)
$$h(1Tx)z_0(x)=\frac{1}{2}+\frac{1}{2+x}.$$

Equations (92) and (98), show that $`h(x)`$ goes rapidly from 1 to $`1/2`$ in a small region close to $`r=1`$ where $`h`$ varies on the scale $`T`$. One can even compute the next correction to $`\lambda T`$ by including the next term of order $`T`$ in the equation for $`z`$. Writing

(99)
$$z(x)=z_0(x)+T^{1/3}z_1(xT^{1/3}),$$

we find

(100)
$$z_1^{}+\frac{2}{u}z_1+z_1^2\frac{D1}{2}u=\frac{\mu ^2}{T^{2/3}}=c_D,u=T^{1/3}x.$$

This is again an eigenvalue problem which selects a unique constant $`c_D`$, that we could only solve numerically. Still, this leads to the non trivial result

(101)
$$\lambda =\frac{1}{4T}+\frac{c_D}{T^{1/3}}+\mathrm{}(T0).$$

We now solve the eigenvalue problem (86) (87) in the limit of large temperatures $`T+\mathrm{}`$ (see also Appendix A). We again perform the change of variables (89) and rewrite Eq. (90) in the form

(102)
$$h^{}+\frac{D1}{r}h=\lambda \frac{1}{T}\left(\frac{h}{r^{D1}}h^2\right).$$

Then, we expand the solutions of this equation in terms of the small parameter $`1/T1`$. We write $`h=h_0+\frac{1}{T}h_1+\frac{1}{T^2}h_2+\mathrm{}`$ and $`\lambda =\lambda _0+\frac{1}{T}\lambda _1+\frac{1}{T^2}\lambda _2+\mathrm{}`$. To zeroth order, we have

(103)
$$h_0^{}+\frac{D1}{r}h_0=\lambda _0.$$

The solution of this equation is $`h_0=\lambda _0r/D`$. Using the boundary condition $`h_0(1)=1`$, we obtain

(104)
$$\lambda _0=D,h_0=r,$$

and $`h_n(1)=0`$ for $`n>0`$. To first order, we get

(105)
$$h_1^{}+\frac{D1}{r}h_1=\lambda _1\frac{h_0}{r^{D1}}+h_0^2.$$

Integrating this first order differential equation and using the boundary condition $`h_1(1)=0`$, we obtain

(106)
$$h_1=\frac{D}{2(D+2)}r\frac{1}{2}r^{3D}+\frac{r^3}{D+2},$$

with

(107)
$$\lambda _1=\frac{D^2}{2(D+2)}.$$

Hence, the large temperature behavior of the eigenvalue is

(108)
$$\lambda =D+\frac{D^2}{2(D+2)}\frac{1}{T}+\mathrm{}(T+\mathrm{}).$$

This expansion can be easily carried out to higher orders but the coefficients are more and more complicated. Restricting ourselves to $`D=3`$, we get $`\lambda =3+\frac{9}{10}\frac{1}{T}\frac{477}{700}\frac{1}{T^2}+O(T^3)`$.

## 6 Numerical simulations in the canonical ensemble

In this section, we illustrate the analytical results obtained in the previous section in the case of $`D=3`$. Except when specified otherwise, our simulations have been performed at $`T=1/5<T_c0.397\mathrm{}`$, for which we have obtained $`t_{coll}0.44408\mathrm{}`$.

In order to perform our simulations, we have used a Runge-Kutta algorithm with adaptive step in space and time. We call $`dr`$ the spatial discretization near $`r=0`$ (which we need to take very small as the density profile becomes singular at $`r=0`$). An important numerical problem arises in the numerical integration of Eq. (69), which is crucial in obtaining non zero values for $`N_0(t)`$. As this equation is a first order differential equation with initial condition $`N_0(0)=0`$, any naive integration scheme should lead to a strictly vanishing value for $`N_0(t)`$ for all time, and any $`dr`$. Still, when performing this naive numerical integration, we see that crossing $`t_{coll}`$ generates increasing values for $`M(dr,t)`$, although keeping $`M(0,t)=N_0(t)=0`$ ultimately makes the numerical integration unstable. In order to bypass this problem, we have decided to introduce a numerical scheme where Eq. (69) is replaced by

(109)
$$\frac{dN_0}{dt}=\rho _0^{fit}N_0^{fit}.$$

$`N_0^{fit}`$ and $`\rho _0^{fit}`$ are extracted from a fit of $`M(r,t)`$ to the functional form (we are in $`D=3`$)

(110)
$$M(r,t)N_0^{fit}(t)+\frac{\rho _0^{fit}(t)}{3}r^3+a_5(t)r^5+a_6(t)r^6,$$

in a region of a few $`dr`$, excluding of course $`r=0`$. This functional form is fully compatible with the expected expansion for $`M(r,t)`$, both below ($`a_6=0`$) and above ($`a_5=0`$) $`t_{coll}`$. We find that this numerical scheme allows us to cross smoothly the singularity at $`t_{coll}`$. An effective cut-off is introduced which effectively depends on $`dr`$, and we have checked that the result presented in this section are extremely close to the ones that would be obtained in the ideal limit $`dr0`$. This is illustrated in Fig. 1, where the smoothing effect of our algorithm is shown to act on a very small time region after $`t_{coll}`$. Even more surprisingly, we find that for sufficiently large times (actually very small compared to any physical time scales), our results are essentially independent of $`dr`$, even for unreasonably large values of $`dr`$. We are thus confident that we have successfully crossed the collapse singularity.

In Fig. 1, we plot $`N_0(t)`$ for small time which compares well with the universal form of Eq. (81), where $`\mu =8.38917147\mathrm{}`$ has been determined so as to ensure the proper behavior of $`S(x)`$ for large $`x`$ (see Sec. 5.1). We also illustrate the exponential decay of $`1N_0(t)`$, with a rate in perfect agreement with the value of $`\lambda `$ extracted from solving numerically the eigenvalue problem of Eq. (86). Finally, we show the effect of the numerical spatial discretization $`dr`$ near $`r=0`$. Satisfactorily enough, the value of $`N_0(t)`$ is sensitive to the choice of $`dr`$ only for very small times after the collapse, and we were able to easily reach small enough $`dr`$, in order to faithfully reproduce the post-collapse singularity. In Fig. 2, we convincingly illustrate the post-collapse scaling, and compare the post-collapse scaling function to that obtained analytically below $`t_{coll}`$ (pre-collapse). In Fig. 3, we confirm the validity of our perturbative expansion for $`\lambda `$, in the limit of small temperature. We compare the value of $`c_D`$ extracted from directly solving the full eigenvalue problem to that obtained from Eq. (100), finding a perfect agreement. Finally, in Fig. 4, we compare the numerical value obtained for $`h(r)`$ to the different analytical estimates given in the preceding section, for $`T=0.01`$. The two important regions $`1r\mu ^2T^{2/3}`$ and $`1r\mu ^2T^{2/3}`$ can be clearly identified.

## 7 Post-collapse in the microcanonical ensemble

So far, we have only addressed the post-collapse dynamics in the canonical ensemble. In the microcanonical ensemble, the dynamical equation have to be supplemented with the strict energy conservation condition (see Eq. (24)), which fixes the global temperature $`T(t)`$. For this model, it was shown in  that below a certain energy (Antonov energy), the system collapses with an apparent scaling associated with $`\alpha _{max}2.2`$ for intermediate times (when the temperature still increases in a noticeable way) before entering a scaling regime with $`\alpha =2`$, identical to that obtained in the canonical ensemble. In the limit $`tt_{coll}`$, the temperature and the potential energy both seem to converge to a finite value preserving a constant energy. Closely before the collapse time, the temperature behaves as $`T(t_{coll})T(t)(t_{coll}t)^\gamma `$ with $`\gamma 1/2`$. This section addresses the $`t>t_{coll}`$ time period.

Assuming a spherical mass density, and after integration by parts, the potential energy $`W`$ can be rewritten in the form ($`D>2`$)

(111)
$$W(t)=\frac{1}{2}_0^1\frac{M^2(r,t)}{r^{D1}}𝑑r\frac{1}{2(D2)}.$$

We see immediately that as $`D1>1`$, the occurrence of a finite mass $`N_0(t)0`$ concentrated at $`r=0`$ implies an infinite potential energy, hence an infinite temperature. We thus anticipate that the post-collapse dynamics in the microcanonical ensemble is probably an ill-defined problem. In this extreme regime, let us try to consider the possible flaws of this model in order to describe a consistent dynamics of a reasonable physical self-gravitating system. First, our assumption of uniform temperature is certainly not realistic in a system displaying huge density contrast, and some alternative approaches are needed to incorporate a spatially dependent temperature. This point is certainly crucial and will be addressed in a future work . Furthermore, in this regime, a careful physical analysis predicts that this system of self-gravitating individual particles should lead to the formation of binaries which is probably beyond the description ability of our essentially mean-field approach. In other words, the system may become intrinsically heterogeneous, which probably cannot be captured by our continuous model. Finally, we can think of other physical effects (degeneracy effects of quantum or dynamical origin, finite particle size effects,…) preventing the system from reaching arbitrarily large densities. One way to describe such effects consists in introducing a spatial cut-off $`h`$ or a density cut-off of order $`h^D`$. In such a system, the dynamics first follows the pre-collapse dynamics until the maximum density is approached. Then, the system will ultimately reach a maximum entropy state that we propose to characterize in a simple manner, as in .

We propose to describe the final state as a “core-halo” structure, which for simplicity we modelize as a core of radius $`h1`$ and constant density

(112)
$$\rho _{core}=\frac{DN_0}{S_Dh^D},$$

which mimics a regularized central Dirac peak containing a mass $`N_0`$. In the region $`h<r1`$ stands the halo of constant density

(113)
$$\rho _{halo}=\frac{D(1N_0)}{S_D(1h^D)},$$

containing the rest of the mass. As $`h`$ is small, we can compute the potential (or total) energy and the entropy only including the relevant leading terms. We find

(114)
$$E=\frac{D}{2}T\frac{D}{D^24}\left[\frac{N_0^2}{h^{D2}}+1+\frac{D2}{2}N_0\frac{D}{2}N_0^2\right]+O(h^2),$$

where the first term is the kinetic energy, whereas the entropy (up to irrelevant constants) reads

(115)
$$S=\frac{D}{2}\mathrm{ln}TN_0\mathrm{ln}\left(\frac{N_0}{h^D}\right)(1N_0)\mathrm{ln}(1N_0)+O(h^D).$$

For a given small value of $`h`$, $`S`$ has a local maximum at $`N_0=0`$ provided that $`E>E_c(h)`$, with $`lim_{h0}E_c(h)=\frac{D}{D^24}`$. Below $`E_c`$, the sole entropy maximum resides at $`N_0`$ satisfying the following implicit equation (again in the limit of small $`h`$)

(116)
$$N_0=\frac{D}{\mathrm{ln}\left(\frac{N_0}{h^D}\right)}\frac{1}{\mathrm{ln}h},$$

where the given asymptotics is quantitatively correct only for extremely small values of $`h`$. Hence, we find that the mass included in the core slowly decreases with the core size $`h`$ , resulting in an effective singularity $`\rho (𝐫)=\frac{1}{\mathrm{ln}r}\delta ^D(𝐫)`$. Meanwhile, the temperature diverges as

(117)
$$T\frac{2}{D}W\frac{2}{D^24}\frac{1}{h^{D2}\mathrm{ln}^2h}.$$

and leads rapidly and efficiently to a uniform halo.

In order to relate this result to our actual system, we perform microcanonical post-collapse simulations using the same regularization scheme as in the canonical case in order to describe the evolution of $`N_0(t)`$. In addition to this, we also need to regularize the potential energy which is strictly infinite (as well as $`T`$) when $`N_00`$. Consistently with the previous discussion, we introduce a numerical cut-off $`h`$ by defining

(118)
$$W(t)=\frac{1}{2}_h^1\frac{M^2(r,t)}{r^{D1}}𝑑r\frac{1}{2(D2)}.$$

In Fig. 5, we plot $`N_0(t)`$ and $`T(t)`$ for different values of the cut-off $`h`$. Contrary to the canonical case, the post-collapse dynamics is strongly $`h`$ dependent. We see that in conformity with our result of Eq. (116), the central mass $`N_0`$ clearly decreases as $`h0`$. Therefore, in the microcanonical case, the physical picture is that when the collapse time $`t_{coll}`$ is reached, the temperature increases rapidly, which leads to the rapid homogenization of the system except for a dense and small core, whose mass $`N_0\mathrm{ln}^1T`$ is a decreasing function of the maximum temperature reached. This central structure with weak mass and huge binding energy is similar to a “binary star” structure in stellar dynamics. Binary formation is the physical process that arrests core collapse in globular clusters . This is also the end point of our simple microcanonical Brownian model.

## 8 Conclusion

In this paper, we have investigated the post-collapse dynamics of a gas of self-gravitating Brownian particles in canonical and microcanonical ensembles. Our results also apply to the chemotactic aggregation of bacterial populations in biology. At the collapse time $`t_{coll}`$, the system develops a singular density profile scaling as $`\rho r^2`$. However, the “central singularity contains no mass”, the temperature does not diverge, and the entropy and free energy are finite . Since this profile is not a maximum entropy (resp. minimum free energy) state nor a stationary solution of the Smoluchowski-Poisson system, the collapse continues after $`t_{coll}`$. This solves the apparent paradox reported in .

In the canonical ensemble, mass accretes progressively at the center of the system and a Dirac peak forms by swallowing the surrounding particles. Eventually, the Dirac peak contains all the mass. This structure has an infinite free energy $`F=ETS\mathrm{}`$ simply because its binding energy is infinite. This is therefore the most probable structure in canonical ensemble . In the microcanonical ensemble, the maximum entropy state (at fixed mass and energy) consists of a single binary embedded in a hot halo . This is precisely what we see in our numerical simulations. The temperature increases dramatically above $`t_{coll}`$ (resulting in an almost uniform halo) although the mass contained in the core is weak (but finite). We note the “spectacular” fact that almost all the gravitational energy resides in a binary-like core with negligible mass. A similar phenomenon is observed in stellar dynamics for globular clusters having experienced core collapse . This shows that the microcanonical Smoluchowski-Poisson system shares some common properties with kinetic equations usually considered in stellar dynamics (Landau-Fokker-Planck equations), despite its greater simplicity. Clearly, a major drawback of our microcanonical model is to assume that the temperature uniformizes instantaneously, implying an infinite thermal conductivity. We shall relax this simplification in a future work . However, the present study is one of the first dynamical study showing the formation of Dirac peaks and binary-like structures in systems with gravitational interaction.

Our Brownian model is based on the existence of a gaseous medium that generates a friction force. This situation exists in certain astrophysical models such as the transport of dust particles in the solar nebula . Dust particles are submitted to Stokes or Epstein drag. It is clear that when the concentration of particles is important (prior to planetesimal formation), self-gravity has to be taken into account. Thus, our system of self-gravitating Brownian particles could be connected to this astrophysical situation. We just mention this as a possible astrophysical application because it is not our present main motivation to make a precise model of dust-gas-gravity coupling in protoplanetary disks. However, this problem could be considered in future works.

## Appendix A Brownian particle around a “black hole”

We consider the Brownian motion of a particle subject to the gravitational force $`GM𝐫/r^3`$ created by a central mass $`M`$ (“black hole”). We assume that when the particle comes at $`r=0`$, it is captured by the central mass. We denote by $`W(𝐫,t)`$ the density probability of finding the particle in $`𝐫`$ at time $`t`$. It is solution of the Fokker-Planck equation

(119)
$$\frac{W}{t}=\left(TW+W\frac{𝐫}{r^3}\right),$$

where we have set $`G=M=R=\xi =1`$. Let $`W(r,t)`$ denote a spherically symmetric solution of Eq. (119) satisfying the boundary conditions

(120)
$$T\frac{W}{r}(1,t)+W(1,t)=0,$$

(121)
$$W(r,0)=\frac{\delta (rr_0)}{4\pi r_0^2}.$$

We call

(122)
$$𝐉=\left(TW+W\frac{𝐫}{r^3}\right),$$

the current of probability, i.e. $`𝐉dS𝐧`$ gives the probability that the particle crosses an element of surface $`dS`$ between $`t`$ and $`t+dt`$ ($`𝐧`$ is a unit vector normal to the element of surface under consideration).

We introduce the probability $`p(r_0,t)dt`$ that a particle located initially between $`r_0`$ and $`r_0+dr_0`$ arrives for the first time at $`r=0`$ between $`t`$ and $`t+dt`$. We have

(123)
$$p(r_0,t)=_{R_ϵ}𝐉𝑑𝐒=4\pi ϵ^2\left(T\frac{W}{r}+\frac{W}{r^2}\right)_ϵ=4\pi W(0,t),$$

where $`R_ϵ`$ is a ball of radius $`ϵ0`$. The total probability that the particle initially between $`r_0`$ and $`r_0+dr_0`$ has reached the center of the system between $`0`$ and $`t`$ is thus

(124)
$$Q(r_0,t)=_0^tp(r_0,t^{})𝑑t^{}.$$

Finally, we average $`Q(r_0,t)`$ over an appropriate range of initial positions in order to get the expectation $`Q(t)`$ that the particle has been captured at time $`t`$.

With the change of variables

(125)
$$W=\psi e^{\frac{1}{2Tr}},$$

we can transform the Fokker-Planck equation (119) into a Schrödinger equation (in imaginary time) of the form

(126)
$$\frac{\psi }{t}=T\mathrm{\Delta }\psi \frac{1}{4Tr^4}\psi .$$

A separation of the variables can be effected by the substitution

(127)
$$\psi =\varphi (r)e^{\lambda t}.$$

This transformation reduces the Schrödinger equation to a second order ordinary differential equation

(128)
$$\varphi ^{\prime \prime }+\frac{2}{r}\varphi ^{}+\left(\frac{\lambda }{T}\frac{1}{4T^2r^4}\right)\varphi =0,$$

with the boundary condition

(129)
$$\varphi ^{}(1)+\frac{1}{2T}\varphi (1)=0.$$

We note $`\lambda _n`$ the eigenvalues and $`\varphi _n`$ the corresponding eigenfunctions. Since the Schrödinger operator $`H=\mathrm{\Delta }1/4T^2r^4`$ is Hermitian, the eigenfunctions form a complete set of orthogonal functions for the scalar product

(130)
$$fg=_0^1f(r)g(r)4\pi r^2𝑑r.$$

The system can be furthermore normalized, i.e. $`\varphi _n\varphi _m=\delta _{nm}`$. Any function $`f(r)`$ satisfying the boundary condition (129) can be expanded on this basis, as

(131)
$$f(r)=\underset{n}{}f\varphi _n\varphi _n.$$

In particular,

(132)
$$\frac{\delta (rr_0)}{4\pi r_0^2}=\underset{n}{}\varphi _n(r_0)\varphi _n(r).$$

The general solution of the problem (119) (120) can be expressed in the form

(133)
$$W(r,t)=\underset{n}{}A_ne^{\lambda _nt}e^{\frac{1}{2Tr}}\varphi _n(r),$$

where the coefficients $`A_n`$ are determined by the initial conditions (121), using the expansion (132) for the $`\delta `$-function. We get

(134)
$$W(r,t)=e^{\frac{1}{2T}(\frac{1}{r}\frac{1}{r_0})}\underset{n}{}e^{\lambda _nt}\varphi _n(r_0)\varphi _n(r).$$

From this expression, we obtain

(135)
$$p(r_0,t)=4\pi e^{\frac{1}{2Tr_0}}\underset{n}{}e^{\lambda _nt}\varphi _n(r_0)\underset{r0}{lim}\left[\varphi _n(r)e^{\frac{1}{2Tr}}\right].$$

Then, according to Eq. (124), we have

(136)
$$Q(r_0,t)=4\pi e^{\frac{1}{2Tr_0}}\underset{n}{}\frac{1e^{\lambda _nt}}{\lambda _n}\varphi _n(r_0)\underset{r0}{lim}\left[\varphi _n(r)e^{\frac{1}{2Tr}}\right].$$

Finally, averaging over the initial conditions, the probability that the particle has been captured by the central mass at time $`t`$ can be expressed as

(137)
$$Q(t)=\underset{n}{}Q_n(t),$$

where

(138)
$$Q_n(t)=B_n(1e^{\lambda _nt}),$$

and

(139)
$$B_n=4\pi \frac{1}{\lambda _n}\overline{e^{\frac{1}{2Tr_0}}\varphi _n(r_0)}\underset{r0}{lim}\left[\varphi _n(r)e^{\frac{1}{2Tr}}\right].$$

This formally solves the problem. If we consider the large time limit, we just need to determine the first eigenvalue $`\lambda _0(T)`$ of the quantum problem. This has been done analytically in Sec. 5.2 in the limits $`T0`$ and $`T+\mathrm{}`$.

Below, we consider again the high temperature regime where thermal fluctuations prevail over gravity but we do not restrict ourselves to the first eigenvalue. To leading order in the limit $`T+\mathrm{}`$, the Fokker-Planck equation (119) reduces to the pure diffusion equation

(140)
$$\frac{W}{t}=T\frac{1}{r^2}\frac{}{r}\left(r^2\frac{W}{r}\right).$$

However, for consistency (see Sec. 5.2), it is necessary to keep the term of order $`1/T`$ (arising from the gravitational force) in the boundary condition. Hence, we take

(141)
$$\frac{W}{r}(1,t)+\frac{1}{T}W(1,t)=0.$$

The general solution of the diffusion equation (140) with the boundary conditions (141) and (121) can be expressed as

(142)
$$W(r,t)=\underset{n}{}e^{\lambda _nt}\varphi _n(r_0)\varphi _n(r),$$

where $`\varphi `$ is solution of

(143)
$$\varphi ^{\prime \prime }+\frac{2}{r}\varphi ^{}+\frac{\lambda }{T}\varphi =0,$$

(144)
$$\varphi ^{}(1)+\frac{1}{T}\varphi (1)=0.$$

Setting $`\varphi =\chi /r`$, Eqs. (143) and (144) become

(145)
$$\chi ^{\prime \prime }+\frac{\lambda }{T}\chi =0,$$

(146)
$$\chi ^{}(1)=\left(1\frac{1}{T}\right)\chi (1).$$

Equation (145) is readily solved. The eigenvalues can be written $`\lambda _n=Tx_n^2(T)`$, where $`x_n(T)`$ are the solutions of the implicit equation

(147)
$$\mathrm{tan}(x_n)=\frac{x_n}{1\frac{1}{2T}}.$$

The eigenfunctions are

(148)
$$\varphi _n(r)=A_n\frac{\mathrm{sin}(x_nr)}{r}.$$

The general solution of the diffusion equation can thus be written

(149)
$$W(r,t)=\frac{1}{rr_0}\underset{n=0}{\overset{+\mathrm{}}{}}e^{Tx_n^2t}A_n^2\mathrm{sin}(x_nr)\mathrm{sin}(x_nr_0),$$

with

(150)
$$A_n^2=\frac{x_n}{2\pi (x_n\mathrm{sin}x_n\mathrm{cos}x_n)}.$$

If we consider the pure diffusion of a particle in a box, the boundary condition (144) reduces to $`\varphi ^{}(1)=0`$ and the $`x_n`$ are solutions of the implicit equation

(151)
$$\mathrm{tan}(x_n)=x_n.$$

In particular, $`x_0=0`$. This implies that the probability $`W(r,t)`$ converges for large times to a uniform profile $`W(r,+\mathrm{})=3/4\pi `$ which is indeed solution of the diffusion equation in a box. If gravity is taken into account, its first order effect (in the limit $`T+\mathrm{}`$) is to change the boundary condition to Eq. (144). It is as if we had a diffusion across the box  although the true physical process is a capture by the central mass. The eigenvalues are now determined by Eq. (147). The $`x_{n>1}`$ are hardly modified (to first order) with respect to the preceding problem but $`x_0`$ is now different from zero. To first order, we find that $`x_0^2=3/T`$ so that $`\lambda _0=3`$ in agreement with the result of Sec. 5.2. We also note that $`A_0^2=T/4\pi `$ while $`A_{n>0}`$ are independent on $`T`$ (to leading order) and given by Eqs. (150) and (151).

Using Eqs. (123) and (124), the probability that the particle has been captured by the central mass at time $`t`$ is given by

(152)
$$Q(t)=\frac{4\pi }{T}\underset{n=0}{\overset{+\mathrm{}}{}}\frac{1e^{Tx_n^2t}}{x_n}A_n^2\overline{\frac{\mathrm{sin}(x_nr_0)}{r_0}}.$$

If we average over initial conditions with the weight $`3r_0^2`$ (uniform distribution), we find to leading order in $`T^1`$ that

(153)
$$\overline{\frac{\mathrm{sin}(x_nr_0)}{r_0}}=0,\mathrm{for}n>0,$$

(154)
$$\overline{\frac{\mathrm{sin}(x_0r_0)}{r_0}}=x_0.$$

Hence, the modes $`n>0`$ cancel out. Therefore, in the high temperature regime, the probability of capture is given by

(155)
$$Q(t)=1e^{3t},$$

for all times.

The case $`D=2`$ can be treated by a similar method. Instead of Eqs. (149) (150) and (147), we get

(156)
$$W(r,t)=\underset{n=0}{\overset{+\mathrm{}}{}}e^{Tx_n^2t}A_n^2J_0(x_nr)J_0(x_nr_0),$$

(157)
$$A_n^2=\frac{1}{\pi [J_1^2(x_n)+J_0^2(x_n)]},$$

(158)
$$\frac{x_nJ_1(x_n)}{J_0(x_n)}=\frac{1}{T},$$

where $`J_n`$ is the Bessel function of order $`n`$. For the pure diffusion process, $`x_n=\alpha _{1n}`$ are the zeros of $`J_1`$. If gravity is taken into account, then $`x_{n>1}\alpha _{1n}`$ while $`x_0^2=2/T`$ establishing $`\lambda _0=2`$. The probability that the particle has been captured by the central mass at time $`t`$ is given by

(159)
$$Q(t)=\frac{2\pi }{T}\underset{n=0}{\overset{+\mathrm{}}{}}\frac{1e^{Tx_n^2t}}{x_n^2}A_n^2\overline{J_0(x_nr_0)}.$$

If we average over the initial conditions with a weight $`2r_0`$ (uniform distribution), we get $`\overline{J_0(x_nr_0)}=0`$ if $`n>0`$ and $`\overline{J_0(x_0r_0)}=1`$. Therefore, in the high temperature regime, the probability of capture is given, for all times, by

(160)
$$Q(t)=1e^{2t}.$$

Finally, for $`D=1`$, we obtain

(161)
$$W(r,t)=\underset{n=0}{\overset{+\mathrm{}}{}}e^{Tx_n^2t}A_n^2\mathrm{cos}(x_nr)\mathrm{cos}(x_nr_0),$$

(162)
$$A_n^2=\frac{1}{[1+\frac{\mathrm{sin}(2x_n)}{2x_n}]},$$

(163)
$$x_n\mathrm{tan}(x_n)=\frac{1}{T}.$$

For the pure diffusion process, $`x_n=n\pi `$. If gravity is taken into account, then $`x_{n>1}n\pi `$ while $`x_0^2=1/T`$ establishing $`\lambda _0=1`$. The probability that the particle has been captured by the central mass at time $`t`$ is given by

(164)
$$Q(t)=\frac{2}{T}\underset{n=0}{\overset{+\mathrm{}}{}}\frac{1e^{Tx_n^2t}}{x_n^2}A_n^2\overline{\mathrm{cos}(x_nr_0)}.$$

If we average over the initial conditions with a weight $`1`$ (uniform distribution), we get $`\overline{\mathrm{cos}(x_nr_0)}=0`$ if $`n>0`$ and $`\overline{\mathrm{cos}(x_0r_0)}=1`$. Therefore, in the high temperature regime, the probability of capture is given, for all times, by

(165)
$$Q(t)=1e^t.$$


format: text

text:
# Photon Mass Bound Destroyed by Vortices

(June 13, 2003)

## Abstract

The Particle Data Group gives an upper bound on the photon mass $`m<2\times 10^{16}`$eV from a laboratory experiment and lists, but does not adopt, an astronomical bound $`m<3\times 10^{27}`$eV, both of which are based on the plausible assumption of large galactic vector-potential. We argue that the interpretations of these experiments should be changed, which alters significantly the bounds on $`m`$. If $`m`$ arises from a Higgs effect, both limits are invalid because the Proca vector-potential of the galactic magnetic field may be neutralized by vortices giving a large-scale magnetic field that is effectively Maxwellian. In this regime, experiments sensitive to the Proca potential do not yield a useful bound on $`m`$. As a by-product, the non-zero photon mass from Higgs effect predicts generation of a primordial magnetic field in the early universe. If, on the other hand, the galactic magnetic field is in the Proca regime, the very existence of the observed large-scale magnetic field gives $`m^11`$kpc, or $`m10^{26}`$eV.

INTRODUCTION: The possibility of a non-zero photon mass remains one of the most important issues in physics, as it would shed light on fundamental questions such as charge conservation, charge quantization, the possibility of the charged black holes and magnetic monopoles, etc. The most stringent upper bounds on the photon mass listed by the Particle Data Group PDG , $`m<3\times 10^{27}`$ eV and $`m<2\times 10^{16}`$, are based on the assumption that a massive photon would cause large-scale magnetic fields to be accompanied by an energy density

$$m_A^2\stackrel{~}{A}_\mu \stackrel{~}{A}^\mu $$
(1)
associated with the Proca field $`\stackrel{~}{A}_\mu `$ that describes the massive photongold . This manifests itself in two different ways. The first limit comes from the potential astrophysical effects yam ; chibisov , and the second from an experiment that used a toroidally magnetized pendulum lakes  to measure the magnetic field gradient in a magnetically shielded vacuum. A recent experimentluo  using an improved technique obtained a similar result. Both experiments actually measured the product $`m^2\stackrel{~}{A}`$, where the ambient Proca vector potential is presumably dominated by the field of the galaxy. The value assumed in lakes ; luo , $`\stackrel{~}{A}RB1\mu `$G$`\times `$kpc, is astronomically reasonable as the large-scale, $`R1`$kpc, galactic field has a strength $`B1\mu `$G.

Let us review the standard arguments behind these bounds, which assume that a massive photon at low energies is described by the Proca field gold . (Throughout this paper we denote the Proca field as $`\stackrel{~}{A}_\mu `$, whereas $`A_\mu `$ should be understood as the Maxwellian field.) The Lagrangian density for $`\stackrel{~}{A}_\mu `$ is

$$=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+\frac{1}{2}m_A^2\stackrel{~}{A}_\mu \stackrel{~}{A}^\mu +\stackrel{~}{A}_\mu J^\mu $$
(2)
where $`F_{\mu \nu }`$ is the usual field strength, and $`J^\mu `$ is the conserved current. The $`m_A^2`$-term in Eq. 2, makes $`\stackrel{~}{A}_\mu `$ experimentally observable. Naively, the existence of the Galactic magnetic field $`𝐁_{\mathrm{gal}}=\times \stackrel{~}{𝐀}_{\mathrm{gal}}`$, implies an ambient galactic vector potential

$$\stackrel{~}{A}_{\mathrm{gal}}B_{\mathrm{gal}}R_{\mathrm{gal}},$$
(3)
The associated Proca energy $`m_A^2\stackrel{~}{𝐀}^2`$ can be detected by direct or indirect observations. Indirect observations rely on the effect the Proca energy would have on the galactic plasma chibisov , implying the limit $`m_A<3\times 10^{27}`$eV. The direct detections lakes  are based on measuring the torque on a magnetized ring, which depends on the angle between $`\stackrel{~}{A}_{\mathrm{gal}}`$ and the vector potential of the ring, because the energy density contains a term $`m_A^2\stackrel{~}{𝐀}_{\mathrm{gal}}\stackrel{~}{𝐀}_{\mathrm{ring}}`$ . The null result lakes  implies the limit $`m_A2\times 10^{16}`$eV.

We claim that photon mass bounds cannot be established without specifying the microscopic origin of the mass. In particular, if $`m`$ arises from the commonly accepted Higgs mechanism, the above bounds do not apply over a large portion of the parameter space. It is quite possible for large-scale magnetic fields to be effectively Maxwellian, even if photons are massive. In this case observations of large-scale fields, say from the galaxy or from Jupiter, are not sensitive to $`m`$. This leaves us with the upper bound from laboratory tests of Coulomb’s law, $`m10^{14}`$eV will . It is also possible that the large-scale fields do remain in the Proca regime. But then, the available information about the large-scale magnetic field of the galaxy, and the gas pressure in the galaxy, actually gives a much more stringent bound, $`m^1R1`$kpc, or $`m10^{26}`$eV yam ; chibisov .

THE HIGGS SCENARIO: If photon has a non-zero mass, there are excellent field-theoretic arguments for thinking it should arise from a Higgs-type effect, in which case the above mass limits are invalid. Note that when estimating the Proca field associated with the galactic magnetic field, one must distinguish the actual Proca field $`\stackrel{~}{A}_\mu `$, which can be measured, from its Maxwellian component $`A_\mu `$; $`\stackrel{~}{A}_\mu `$ has 3 physical degrees of freedom (polarizations) as opposed to $`A_\mu `$ , which has only two. So that by giving a mass term to $`A_\mu `$, we are supplementing it with an additional degree of freedom. The Proca field can be written as

$$\stackrel{~}{A}_\mu =A_\mu \frac{1}{g}_\mu \psi ,$$
(4)
where $`\psi `$ is the additional (longitudinal) polarization. Written in this way, the Proca theory is manifestly gauge invariant under:

$$A_\mu A_\mu +\frac{1}{g}_\mu \omega ,\psi \psi +\omega $$
(5)
The new degree of freedom enters the Proca action only in the mass term; it cancels in $`F_{\mu \nu }`$ as well as in the couplings to the conserved current. When computing the galactic Proca field, we must be sure it is not compensated by the additional polarization. This is, in fact, what happens in the Higgs scenario.

Proca theory, Eq. 2, can be extended to the Higgs theory, by promoting $`m_A`$ into a real scalar field $`\varphi =m_A/g`$, with the self-Lagrangian

$$\frac{1}{2}(_\mu \varphi )^2\frac{\lambda ^2}{2}(\varphi ^2\eta ^2)^2.$$
(6)
Then, $`\varphi `$ can be thought of as the modulus of a complex Higgs field, $`H=\varphi \mathrm{e}^{i\psi }`$, and $`\psi `$ is its phase (Goldstone boson), which becomes a longitudinal photon. The static energy in the absence of the electric field is

$$=d^3x\left(\frac{B^2}{2}+\frac{g^2}{2}\varphi ^2\left(𝐀\frac{1}{g}\psi \right)^2+\frac{1}{2}(\varphi )^2+\frac{\lambda ^2}{2}(\varphi ^2\eta ^2)^2\right).$$
(7)
The two important parameters of the theory are the mass of the photon $`m_A=g\eta `$, and the mass of the Higgs particle $`m_\varphi =\lambda \eta `$. With frozen $`\varphi =\eta `$, the theory simply reduces to Proca theory with the energy

$$=d^3x\left(\frac{B^2}{2}+\frac{g^2}{2}\eta ^2\stackrel{~}{𝐀}^2\right),$$
(8)
where now $`\stackrel{~}{A}`$ carries three degrees of freedom, only two of which contribute into $`B`$.

Naively, one may think that as long as $`\lambda ^2\eta ^4`$ is bigger than the energy density of the galactic magnetic field $`B^2`$ (with $`B\mu `$G $`10^{25}`$ GeV2), one can ignore the fluctuations in the Higgs field, and the theory should effectively reduce to Proca. This is not correct, because the system has a choice of lowering a huge Proca energy, stored over a large volume, by locally exciting the Higgs field. Even if the Higgs field is much heavier than the value of the magnetic field in question for a big portion of parameter space, the local Higgs energy is in fact less costly than the alternative Proca gauge-field energy. Crudely speaking, if $`m`$ is due to a Higgs effect, then the Universe is effectively a type II superconductor where magnetic fields create Abrikosov vortices.

This effect can be readily demonstrated by taking an example with a constant magnetic field $`B_z=2B`$. The corresponding Maxwellian vector potential (up to gauge transformation) is $`A_\theta =Br`$, and in the absence of a third polarization could naively contribute a divergent Proca energy. However, this energy is canceled by non-trivial winding of longitudinal photon. The Proca energy density is

$$g^2|\varphi |^2\left(Br\frac{1}{gr}_\theta \psi \right)^2,$$
(9)
which cancels if on average $`_\theta \psi =gBr^2`$ . This configuration of $`\psi `$ is impossible in the Proca theory, but can occur in the Higgs case if there is an uniform density of zeros of Higgs field $`\varphi `$, around which the phase winds non-trivially, producing vortices. The integral around a closed circle of radius $`r`$

$$\frac{1}{2\pi }_\theta \psi =N(r)$$
(10)
defines a winding number $`N(r)`$ which is equal to the number of vortices located inside the circle. Around each vortex $`\psi `$ changes by $`2\pi `$. The system cancels the Proca energy by creating uniform density of vortices $`n=gB/\pi `$ . The cancellation cannot be exact because $`N(r)`$ is a discrete number, so the residual Proca energy density is $`gB\eta ^2`$ . Equating this to the Higgs energy, $`gB\eta ^2\lambda ^2\eta ^4`$, we get a critical value of the magnetic field $`B_c`$

$$B_c\lambda ^2\eta ^2/g.$$
(11)
For $`B>B_c`$ it is energetically favorable for $`\varphi `$ to vanish everywhere, and the theory becomes Maxwellian. The same value of $`B_c`$ can be obtained by requiring that the Higgs cores overlap. That is, the inter-vortex distance becomes equal to the inverse Higgs mass:

$$\frac{1}{\sqrt{n_c}}\frac{1}{\sqrt{B_cg}}\frac{1}{\lambda \eta }$$
(12)
CRITICAL VALUES OF B: There are several interesting critical values of the system parameters. The first is the magnetic field given by Eq. 11. For $`B>B_c`$ the Higgs VEV vanishes and the photon becomes massless. Even if the galactic magnetic field is above $`B_c`$ this limit can still be of interest, because the extra-galactic magnetic field can be below $`B_c`$. Then the photon will be massive outside the galaxy, but massless inside. In such a regime, the information about $`m`$ can only come from extra-galactic observations.

For $`B<B_c`$, the photon is massive everywhere, and there are two regimes: Proca and non-Proca.

The system can only be in the Proca regime (i.e. vortex formation is unfavorable) when its size satisfies $`R\mathrm{\hspace{0.17em}1}/\sqrt{gB}`$. For the galaxy, assuming $`m_A10^{14}`$eV and $`\lambda =1`$, this limit requires $`\eta 10^{22}`$GeV!

If $`R\mathrm{\hspace{0.17em}1}/\sqrt{gB}`$, vortices are energetically favored, but two sub-regimes are possible. The first occurs when $`R\lambda \eta /gB`$. Then the system classically creates vortices out of vacuum and neutralizes the Proca energy. For the galaxy, assuming $`m_A10^{14}`$eV and $`\lambda =1`$, this requires $`g10^{16}`$ and $`\eta 10^2`$eV. Such a light, weakly-charged Higgs is compatible with all existing experimental data and naturalness bounds. At each point, the typical number of magnetically ovelapping vortices is $`\stackrel{~}{N}m_A^2gB`$. If $`\stackrel{~}{N}\mathrm{\hspace{0.17em}1}`$, the field is effectively Maxwellian. For the galaxy, assuming $`m_A10^{14}`$eV and $`g10^{16}`$, we get $`\stackrel{~}{N}\mathrm{\hspace{0.17em}10}^5`$. (Situations with $`\stackrel{~}{N}\mathrm{\hspace{0.17em}1}`$ could provide experimental signatures that would be smoking guns for the Higgs scenario.)

The opposite case occurs when $`R\lambda \eta /gB`$. Vortices are still energetically favorable, but the system cannot create them classically so that their existence will depend on pre-existing conditions such as phase transitions in the early Universe. Due to its very small charge, the Higgs field decoupled from ordinary matter very early so that a phase transition with vortex formation could have preceded formation of the magnetic field. The evolution of such vortices is not yet understood, but is expected to be different from more conventional cosmic-string networks vilenkin .

Although we focused on Proca-Higgs cases, our analysis can be extended to alternative gauge-invariant, ghost-free theories of the photon mass dgs  in which $`E`$-field of a point charge for $`rm^1`$ is not screened but rather modified to a higher inverse power law $`1/r^3`$. In such cases the constraints may be even milder.

PRIMORDIAL MAGNETIC FIELD: As an interesting by-product, the non-zero photon mass naturally predicts generation of a self-sustained primordial magnetic field in the early universe. Indeed if photon aquires mass by the Higgs mechanism, than the thermal phase transition in $`\varphi `$ would inevitably form vortices by the Kibble mechanism Kibble .

Due to small charge, $`\varphi `$ is never in thermal equilibrium with the Standard Model species, but because of the large self-coupling, it is in thermal equilibrium with itself. Thus, due to the usual high-temperature symmetry restoration, the expectation value of $`\varphi `$ had to vanish at early times. The only situation in which $`\varphi `$ would not vanish, would be if it never was in a thermal equilibrium since inflation. This is unlikely – even if $`\varphi `$ had no direct coupling with inflaton, it would still be produced gravitationally with a Gibbons-Hawking temperature ($`T_{GH}10^{14}GeV`$ for the standard inflation), unless $`m_\varphi >T_{GH}`$.

The vorticies are produced when the temperature of the $`\varphi `$ field drops to $`T_\varphi m_\varphi `$. The standard big bang nucleosynthesis requires that the temperature in $`\varphi `$ quanta be smaller than the temperature in the standard model sector, $`T_\varphi <T_{SM}`$, and the vortex network would form before galaxies if

$$m_\varphi >\frac{T_\varphi }{T_{SM}}10^3\mathrm{eV}.$$
(13)
Thus, formation of a primordial magnetic field is a direct consequence of the photon mass. The evolution of this magnetic field after formation, and its role for the galactic dynamo will be discussed elsewhere progress .

THE PROCA REGIME:  If the galaxy is in the Proca regime the averaged magnetic pressure is (see below)

$$p_{\mathrm{magnetic}}=\frac{𝐁^2}{24\pi }\frac{m^2\stackrel{~}{𝐀}^2}{24\pi }.$$
(14)
(The electric pressure is much smaller because the interstellar medium is a good conductor). In a stable system, this magnetic pressure must be counterbalanced by the plasma pressure and/or the plasma kinetic energy. The interstellar medium of our galaxy is in approximate equipartition assuming conventional electrodynamics shu ; the kinetic energy density of the plasma, the plasma pressure, and the standard magnetic energy density $`B^2/(8\pi )`$ are comparable to each other. Therefore, the “massive” part of the full magnetic pressure cannot exceed the standard part $`m^2\stackrel{~}{𝐀}^2𝐁^2`$, which, together with the estimate $`\stackrel{~}{A}RB`$, gives our bound

$$mR^110^{26}\mathrm{eV}.$$
(15)
In fact, this upper bound, derived from energy equipartition, was already discussed by Yamaguchi yam , but then dismissed because the energy source of the magnetic field is unknown gold . But, no matter what the source of the energy, the interstellar medium of the galaxy must provide the pressure support against the anomalous negative magnetic pressure, and the Yamaguchi estimate is, in fact, correct.

We can repeat the above analysis in terms of the Lorentz force. In Proca theory, one can still calculate magnetic fields from Ampere’s law, if a new current density $`m^2\stackrel{~}{𝐀}`$ is added to the usual electric current $`𝐣`$. Then, approximately, $`BRj+m^2R^2B`$. Assume for the moment that $`m`$ is equal to the PDG upper bound, $`10^{16}`$eV. Then $`m^2R^2>>1`$, and the usual current must precisely balance the “massive current”, $`4\pi 𝐣=m^2\stackrel{~}{𝐀}`$. But this current results in a huge Lorentz force density, $`jB`$, which cannot be possibly counterbalanced by the large-scale pressure gradients or accelerations.

Finally, the same upper bound can be presented as a virial theorem chibisov . The virial theorem relates mean values of different forms of energy for a system of particles and fields executing a bound motion. For the Proca theory, the energy-momentum tensor is

$$T_\nu ^\mu =\frac{1}{4\pi }\left(F^{\mu \alpha }F_{\nu \alpha }+\frac{1}{4}\delta _\nu ^\mu F^2\right)+\frac{m^2}{4\pi }\left(\stackrel{~}{A}^\mu \stackrel{~}{A}_\nu \frac{1}{2}\delta _\nu ^\mu \stackrel{~}{A}^2\right).$$
(16)
The magnetic pressure is defined as the magnetic part of $`T_i^i`$, and gives Eq. 14. The virial theorem, proved along standard lines ll , is

$$=T\frac{m^2}{4\pi }d^3r\stackrel{~}{A}^2,$$
(17)
where $``$ is the energy of the system, $`T`$ is the mean kinetic energy of the particles, $`\stackrel{~}{A}^2=\mathrm{\Phi }^2\stackrel{~}{𝐀}^2`$ is the mean squared 4-potential. Assuming that $`B`$-fields are much larger than the $`E`$-fields, we obtain

$$U_g+\frac{m^2}{8\pi }d^3r\stackrel{~}{𝐀}^2=2T+\frac{1}{8\pi }d^3r𝐁^2,$$
(18)
which shows that, for virialized motion, kinetic energy (which includes plasma pressure and kinetic energy of the macroscopic motion) plus Maxwellian magnetic energy is equal to the gravitational energy plus the Proca part of the magnetic energy. The bound in Eq. 15 assumes that the virial theorem can be applied, approximately, to the random part of the galactic motion, after the mean rotation of the galaxy has been excluded.

CONCLUSIONS:  When trying to measure $`m`$ one must distinguish between measurements performed on large and small scales. If the photon acquires mass by the Higgs mechanism, the large-scale behavior of the photon might be effectively Maxwellian. If, on the other hand, one postulates the Proca regime for all scales, the very existence of the galactic field implies $`m<10^{26}`$ eV, as correctly calculated by Yamaguchi and Chibisov yam ; chibisov .

We thank S. Dimopoulos, G. Farrar, M. Shaposhnikov, A. Vilenkin for useful discussions. E.A. thanks the NYU Cosmology and Particle Physics Group for a stimulating visit. G.D. and A.G. are supported by the David and Lucile Packard Foundation, and G.D. is also supported by the Alfred P. Sloan Foundation and by NSF PHY-0070787.