Giovanni Manfredi
email: giovanni.manfredi@ipcms.unistra.fr

Jean-Louis Rouet

Bruce Miller

Gabriel Chardin

Cosmological structure formation with negative mass

We construct a family of models with negative gravitational mass in the context of Newtonian gravity. We focus in particular on a model that reproduces the features of the so-called Dirac-Milne universe, a matter-antimatter symmetric universe that was recently proposed as an alternative cosmological scenario [A. Benoit-Lévy and G. Chardin, A&A 537, A78 (2012)]. We perform one-dimensional N-body simulations of these negative-mass models for an expanding universe and study the associated formation of gravitational structures. The similarities and differences with the standard cosmological model are highlighted and discussed.

I. INTRODUCTION

It may seem unlikely to study cosmological structure formation using hypothetical negative mass objects. On the other hand, we cannot but note the strangeness of the Standard Cosmological Model which, although impressively concordant on primordial nucleosynthesis, cosmic microwave background (CMB), baryon acoustic oscillations (BAO) and type-1a supernovae (SN1a) luminosity distance, features a very strange composition, with dark matter and dark energy, two unidentified components, supposedly representing approximately 96% of the universe content.

In the standard (ΛCDM) model, it is also difficult to understand the different periods of dominance and insignificance of the various components of the Universe: on one hand, dark energy is remarkably insignificant in the very early stages of the universe, representing less than 10 -120 of the energy density at that epoch, and remaining almost completely negligible until recently. On the other, dark energy, if it exists, is supposed to represent at present about 70% of the universe energy density, and will become strongly predominant in the future. Similarly, matter appears as an insignificant component in the radiation dominated universe from the quark-gluon plasma transition down to temperatures of a few electron-Volts and, after a period of matter domination, will soon become again insignificant. We should therefore be ready to examine the proposition that our description of the universe following the Standard Cosmological Model, although it may represent a fair adjustment to the experimental data, has no connection with the actual composition of the universe.

Looking for descriptions alternative to the Standard Cosmological Model, several authors have noted that our universe is very close to a "coasting" universe, i.e., a universe that neither decelerates nor accelerates, and therefore behaves, at large scales ( 100Mpc), as a gravitationally empty, or Milne [1], universe. For example, Nielsen, Guffanti and Sarkar (NGS) [2] argue that the present SN1a data, twenty years after the discovery of repulsive gravity, is still unable to demonstrate convincingly the acceleration of the universe expansion rate. Similarly, by expanding on the NGS analysis, Blanchard and collaborators [3] have shown that cosmic acceleration by all "local" cosmological probes (at z < 3, where z is the cosmological redshift) is not statistically compelling, while recently the SN1a data of the Hubble Space Telescope [4] show that, at the very least, the Milne universe appears as an excellent fit to these SN1a data. Also, it is well known that the age of our ΛCDM universe is nearly equal to the age of a Milne universe (see [5] and references therein).

These are motivations to look for the possibility that negative mass particles exist in the universe, which would incorporate equal amounts of positive and negative mass particles, thus appearing gravitationally empty. While negative mass solutions seem to violate blatantly the sacrosanct energy conditions, it should be noted that counterexamples have been found to basically all the expressions of energy conditions (for a review, see [6]). Indeed, dark energy itself, with a pressure equal to the negative of its energy density, violates the strong energy condition. On the other hand, recently, Paranjape and collaborators [7] have explicitly constructed perfectly physically acceptable negative-mass "bubbles" sitting in an expanding Einstein-de Sitter universe.

A further motivation to study a universe featuring equal amounts of positive and negative mass particles is provided by the so-called "Dirac-Milne" universe [5], i.e., a symmetric matter-antimatter universe where antimatter particles have a negative gravitational mass, which is remarkably con-cordant without any adjustable parameter. Matter-antimatter universes have been studied by Omnès and his group in the late sixties [8], and later by Cohen, de Rujula and Glashow [9], and it would seem that they are excluded by observation: annihilation at the matter-antimatter domain interfaces produces diffuse gamma-ray flux in the tens of MeV range that contradict existing observational constraints, unless the size of matter and antimatter domains exceeds significantly the gigaparsec scale. However, these two studies suppose, quite understandably, that antimatter respects the Equivalence Principle (EP), which lies at the heart of general relativity.

Is there a simple and natural modification of the EP for antimatter that avoids the annihilation constraint? To answer this question, we note that there exists a physical system, the electronhole gas in a semiconductor, where pseudo-particles and their antiparticles avoid annihilation by creating a charge-free region -i.e., a depletion zone -under either an electric or gravitational field. In this system, electrons gravitate and holes antigravitate [10]. The electron-hole system in a semiconductor provides a physical implementation of Dirac's original idea of the positron as an electron hole in a Dirac sea. A natural redefinition of the EP using two types of negative and positive mass carriers (by analogy with negative and positive charge carriers) was first proposed by Piran [11] following the simulations by Dubinski et al. [12], where they considered underdense regions as negative mass regions, with counterintuitive repulsive behavior. The fact that underdense regions lead to nonlinear structures of progressively larger size, providing discriminating cosmological tests, has been studied in particular by Sheth and collaborators [13,14].

To date, there is no clear-cut direct experimental evidence on the gravitational behavior of antimatter. But this situation may soon change, with several precision experiments currently undertaken at CERN to test the gravitational response of neutral antihydrogen atoms. In particular, the first results of the Gbar [15], ALPHA-g [16], and AEgIS [17] collaborations should be available before the end of this decade, and any deviation from perfect matter-antimatter symmetry will have profound cosmological implications.

In the present work, we first study the various sign combinations of inertial, active and passive gravitational mass that could be used to define negative mass particles in the context of Newtonian gravity. Remarkably, we find that the system corresponding to the analog of an electron-hole system, motivating the Dirac-Milne universe, cannot be represented by any combination of such Newtonian masses, whatever their signs. Instead, the Dirac-Milne scenario must be represented by a matrix Poisson equation that couples unconventionally the two particle species. Subsequently, we investigate gravitational structure formation in such negative-mass universes by means of numerical simulations of an expanding one-dimensional universe, focusing in particular on the Dirac-Milne scenario. Significantly, in this last case, simulations show that structures are produced at relatively early epochs and remain "frozen" in the subsequent history of the universe, without requiring either inflation or dark matter ingredients.

II. NEWTONIAN COSMOLOGY WITH NEGATIVE MASS

The possibility that particles with negative mass exist has long been considered (and sometimes dismissed) in the past [18][19][20], starting from the seminal work of the late Hermann Bondi [21]. Negative mass is not necessarily incompatible with fundamental physical laws, such as the equiv-alence principle, the conservation of energy and momentum, and Newton's third (action/reaction) law. Nevertheless, it does give rise to downright unusual features [20], such as the runaway acceleration of two opposite masses placed near each other, or the fact that, although a positive and a negative mass will individually fall in a gravitational field, the bound system of these two objects will levitate and polarize, with the negative mass levitating above the positive one [18]. In this work, our attention will be focussed on the cosmological consequences of matter that displays a negative mass. This question was raised in the context of Milgrom's Modified Newtonian Dynamics (MOND) [22]. Cosmological simulations using negative Bondi masses were also reported in a recent study [23].

It should be added that negative effective mass is a rather old and uncontroversial topic in such disparate fields as materials science [24], semiconductor electronics [25], and cold atomic gases [26]. Usually, it means that a certain wave dispersion relation contains a negative term that can be interpreted, by analogy, as a "negative mass".

In the remainder of this section, we review different models of negative mass in the context of Newtonian gravity. Some of these models have been considered in the past, such as the Bondi masses mentioned above or the standard antigravity scenario whereby unlike masses repel and like masses attract each other. But other cases are possible if one relaxes the hypothesis that Newtonian gravity is fully described by a single gravitational potential, and considers instead two potentials. This approach, described in more detail in Sec. II B, can be viewed as the nonrelativistic limit of a bimetric theory of gravity [27,28]. In this perspective, it is possible to derive a theory that possesses the features required by the Dirac-Milne universe proposed recently by Benoit-Lévy and Chardin [5]. Within the Dirac-Milne cosmology, it is ordinary antimatter that carries a negative mass. In this spirit, we will often refer to particles with positive and negative mass respectively as "matter" and "antimatter" although, generally speaking, negative mass particles need not be identified with antimatter: they are just another kind of matter that responds differently to gravitational fields.

A. Negative active or passive gravitational masses

Let us first consider the cases that can be obtained easily by posing different signs for the various types of masses. As is customary, we distinguish between the active m a and passive m p gravitational masses, in addition to the inertial mass m i . The active gravitational mass is the one that appears on the right-hand side of Poisson's equation:

∆φ = 4πGρ = 4πGm a n, ( 1
)
where n is the number density, ρ the mass density, and φ the gravitational potential. The passive gravitational mass relates the force to the gravitational field through: F = -m p ∇φ. Finally, the inertial mass gives the relation between velocity and momentum: p = m i ṙ (the dot denotes differentiation with respect to time). Therefore Newton's second law of motion ṗ = F becomes r = -(m p /m i )∇φ.

(

) 2
The equivalence principle (EP) demands that the inertial mass be identical to the passive gravitational mass, so that a gravitational field becomes equivalent to an acceleration. But the EP says nothing about the active gravitational mass, which may in principle be different. On the other hand, Newton's third law (the action/reaction principle) requires that the passive and active gravitational masses be the same; otherwise the total force acting on an isolated system of interacting particles would be nonzero, and momentum conservation would be violated.

Since only the ratio m p /m i appears in this Newtonian context, we will assume that the inertial mass is positive and consider the different signs for the active and passive gravitational masses. All possible cases are summarized in Table I Case A is the standard one, with all three masses being positive. For the other cases, we consider a two-component system composed of both matter (type-A particles in Table I) and antimatter (types B, C, or D). Case B is the "antiplasma" scenario, so named because it is the analog of a two-component plasma (made of electrons and ions interacting via the Coulomb force) but with opposite sign for the interaction: like masses attract, whereas unlike masses repel each other. Note that case B does not respect the EP (m p = m i ). Case C is that of the so-called Bondi masses [21]: positive masses attract everything; negative masses repel everything. This case respects the EP, because m p = m i , but gives rise to some surprising features, like the runaway acceleration of a pair of positive and negative masses. Indeed, if a positive and a negative mass are placed near each other, the positive mass will attract the negative mass and the negative mass will repel the positive one: the net result is that both masses accelerate without bound. As only the ratio m p /m i is relevant, the case of Bondi masses may also be viewed as a case where all three masses are negative. This is actually necessary to ensure kinetic energy and momentum conservation during the above-mentioned runaway acceleration.

The case D was never considered before, at least to our knowledge. In this case, both matter and antimatter contribute in the same way to the gravitational field (because their active gravitational masses are both positive). However, while matter is attracted by the field, antimatter is repelled, because it has negative passive gravitational mass. We termed this case "antiinertia" because it can also be viewed (by changing the signs of both m p and m i ) as a case where both gravitational masses are positive while the inertial mass is negative.

Finally, let us note that of the four possible cases of Table I, only A (all masses positive) and C (all masses negative) satisfy at the same time the action/reaction principle and the EP.

In the rest of this work, we will deal with two-component systems made of one positive-mass species (case A, "matter") and a negative-mass species ("antimatter"). For the latter, we will only consider cases B (antiplasma) and C (Bondi). In the Newtonian limit and in the absence of two-body correlations (collisionless approximation), such a system is governed by the following Vlasov-Poisson equations:

∂f ± ∂t + v • ∇f ± - m ± p m ± i ∇φ • ∇ v f ± = 0, (3) ∆φ
= 4πG(m + a n + + m - a n -), (4)
where f ± (r, v, t) are the probability densities (distribution functions) for matter (+) and antimatter (-) in the phase space, and n ± = f ± dv are the number densities. All masses have the same absolute value (|m a | = |m p | = |m i | = m) and their signs are given in Table I.

B. Dirac-Milne scenario

The idea of a matter-antimatter symmetric universe has already been considered in the past [9,29], but was excluded because annihilations at the frontiers of matter and antimatter domains would generate a diffuse gamma ray emission that would be in contradiction with observational data (unless the matter domain we inhabit is virtually the entire visible universe).

Recently, Benoit-Levy and Chardin [5] considered an alternative matter-antimatter symmetric cosmological scenario, where the gravitational interaction between matter and antimatter is repulsive. The latter property is a necessary condition to avoid contact and annihilation between matter and antimatter domains after cosmological recombination. Such matter-antimatter repulsion could in principle be achieved by cases B (antiplasma) and C (Bondi) in Table I. However, the antiplasma scenario does not respect the EP and implies the existence of antimatter gravitational structures resulting from antimatter-antimatter attraction, again leading to unwanted annihilation at the boundaries between matter and antimatter regions. The Bondi case would require a negative inertial mass to ensure energy and momentum conservation, and displays peculiar and counterintuitive features such as the runaway acceleration.

In order to reproduce the Dirac-Milne universe postulated by Benoit-Levy and Chardin [5], we would need a scenario in which matter forms self-gravitating structures such as galaxies and clusters of galaxies, whereas antimatter constitutes a diffuse low-density background. This can be achieved if the following conditions are satisfied: (i) matter attracts matter, (ii) antimatter repels antimatter, and (iii) antimatter and matter repel each other -these conditions are summarized in Table II. The difference with respect to the antiplasma case is that for the Dirac-Milne scenario negative masses repel each other, whereas in the antiplasma case they attract. Also, contrarily to the Bondi case, matter repels antimatter in the Dirac-Milne scenario.

In the Dirac-Milne scenario, matter forms gravitational structures, whereas antimatter, being repelled by everything (including itself), tends to spread across all the available space. Such spread is almost uniform, but not quite: since matter repels antimatter, the latter is expelled from matterdominated overdense regions (galaxies) and forms a low-density almost homogeneous background distributed over the underdense regions in between matter's gravitational structures. This behavior precludes the occurrence of any significant annihilation, in accordance with observation. In the following analyses and computer simulations, we will often make the simplifying hypothesis that antimatter constitutes a low-density homogeneous background uniformly distributed everywhere in space. This approximation is justified for the study of gravitational structure formation, as overdense regions are very much dominated by matter anyway. However, the Dirac-Milne scenario described in the above paragraph does not have a proper Newtonian limit, in the sense that no combination of the signs of the three masses (as in Table I) can reproduce it. The Dirac-Milne scenario can only be accounted for by two gravitational potentials that obey two distinct Poisson equations. In the collisionless approximation, this leads to the following Vlasov-Poisson system:

∂f ± ∂t + v • ∇f ± -∇φ ± • ∇ v f ± = 0, (5)
∆φ

+ = 4πGm(+n + -n -), (6)
∆φ -= 4πGm(-n + -n -) (7)
The above equations cannot be reduced to one of the cases summarized in Table I. This is because there are actually two active gravitational masses for matter, positive in Eq. ( 6) and negative in Eq. ( 7) (for antimatter, the active gravitational mass is always negative). It can be proven that the above equations conserve the total momentum

P tot = mv (f + + f -) drdv
and the total energy

E tot = m 2 (f + + f -)|v| 2 drdv + 1 8πG |∇φ -| 2 -|∇φ + | 2 2 + ∇φ -• ∇φ + dr.
In the remainder of this paper, we will elaborate on some of the consequences of three of the models described so far, namely the antiplasma, Bondi, and Dirac-Milne cases. In particular, we will be concerned with the problem of cosmological structure formation in an expanding universe. But before doing that, we show in the next subsection that all the aforementioned scenarios can be described in the framework of a single general formalism.

C. General formalism

The representation of Newtonian gravity with two Poisson equations, as in Eqs. ( 6)- (7), is not entirely new. For instance, Hohmann and Wohlfart [27] considered a general-relativistic theory in which two classes of particles, respectively with positive and negative mass ratios m p /m i , are described by different metric tensors. They further proved a no-go theorem that rules out all bimetric theories where like particle attract and unlike particles repel each other, i.e., the analog of the antiplasma case that we described in Sec. II. However, the Dirac-Milne scenario considered in the present work should not be affected by this theorem. Hossenfelder [30] also developed a similar bimetric theory which putatively does not suffer from the restriction of the above-mentioned no-go theorem [31].

All the scenarios considered in the preceding sections can be represented in matrix form with the following Vlasov-Poisson equations

∂f ± ∂t + v • ∇f ± -∇φ ± • ∇ v f ± = 0, (8)
∆Φ = 4πGm M n, (9)
where

Φ = φ + φ - , n = n + n - , M = M ++ M +- M -+ M -- , (10)
and M ij = ±1. Thus, M ++ indicates the effect of matter (+) on matter (+), M +-the effect of antimatter (-) on matter (+), M -+ the effect of matter (+) on antimatter (-), and M --the effect of antimatter (-) on antimatter (-). For instance, M ++ = +1 means that matter attracts matter, M +-= -1 means that antimatter repels matter, M --= -1 means that antimatter repels antimatter, etc.

The three cases described in Sec. II (antiplasma, Bondi, and antiinertia) are represented by the following matrices:

M ap = 1 -1 -1 1 , M Bondi = 1 -1 1 -1 , M ai = 1 1 -1 -1 . ( 11
)
Note that in all the above cases the matrix M is non-invertible, i.e. det M = 0. This is expected, as such cases can be described by a single Poisson equation.

Let us now consider the most general case for the matrix M. One must always have M ++ = 1, because this is the standard behavior of ordinary matter (attraction). The other three elements of the matrix can be chosen freely, which leaves us with a total of 2 3 = 8 possibilities. The special case where all elements are positive is trivial, as it corresponds to a situation where matter and antimatter are gravitationally indistinguishable. Then we have the three cases described by Eq. (11). This leaves us with another four cases that cannot simply be represented by a single Poisson equation. One of them is the Dirac-Milne case described in Sec. II B, which yields the matrix:

M DM = 1 -1 -1 -1 . ( 12
)
The remaining three cases can be written as

M 1 = 1 1 1 -1 , M 2 = 1 1 -1 1 , M 3 = 1 -1 1 1 . ( 13
)
Note that in the cases of Eqs. ( 12)-( 13) the matrix M is invertible, i.e. det M = 0, so that they cannot be reduced to one single Poisson equation. All cases where the matrix is invertible should not be affected by the no-go theorem of Hohmann et al. [27].

It is easily shown that the matrix Poisson equation ( 9) is equivalent to the following Lagrangian:

L(φ + , φ -) = ∇Φ T • ∇Φ 8πG + m Φ T M n, (14)
where the superscript T denotes the transpose matrix.

Some final considerations on the above models are in order here. Four of them have the unpleasant property that the interaction between unlike particles is not symmetric -matter attracts antimatter but antimatter repels matter, or viceversa. This group of models include the Bondi, antiinertia, M 2 , and M 3 cases. As a consequence, they all display the peculiar runaway acceleration effect that we described in Sec. II for the Bondi masses. This seems to be a good argument to rule out such models.

Of the remaining models, the antiplasma case leads to the formation of overdense gravitational structures for both matter and antimatter ("antigalaxies"), for which there is no observational evidence. The numerical simulation of this model (see Sec. VI A) also reveals some unobserved features, such as the creation of matter and antimatter streams with opposite velocities. Finally, the case described by the matrix M 1 is similar to the Dirac-Milne scenario, with the significant difference that unlike masses attract each other. Again, this is bound to lead to many annihilation events because matter and antimatter would tend to cluster together.

All in all, the Dirac-Milne case is singled out as the most likely candidate to reproduce the features of the alternative cosmology proposed by Benoit-Lévy and Chardin [5].

III. STEADY STATES AND LINEAR RESPONSE ANALYSIS

It is interesting to note that, out of the eight possible cases mentioned in Sec. II C, spatially homogeneous steady states (i.e., with n + = n -= const. and φ ± = 0) can only exist for the antiplasma and Bondi cases. This is because they are the only cases for which positive and negative mass densities can compensate each other. This is not a problem per se as ordinary gravity itself does not admit steady state solutions.

For a Dirac-Milne universe, one can assume, to an excellent approximation, that antimatter constitutes a homogeneous dilute background with density n 0 , as was discussed in Sec. II B. In that case, one can drop the second Poisson equation ( 7) and a steady state becomes indeed possible for n + = n 0 . This scenario will be justified on more physical grounds when we introduce comoving coordinates in an expanding frame of reference for cosmological applications, see Sec. IV.

Hereafter, we present a linear response analysis of the Vlasov-Poisson equations for the three cases of interest: antiplasma, Bondi, and Dirac-Milne. We expand the distribution functions f ± (r, v, t) = f 0 (v) + f± (r, v, t), the densities n ± = n 0 + ñ± (r, t), and the potentials φ ± (r, t) = φ 0 + φ± (r, t) (where the tilde denotes a small perturbation) and only retain first-order terms. Then we Fourier analyze the perturbations in space and time, by writing for the density fluctuations: ñ± (r, t) = ñ± exp[i(k • r -ωt)], and analogous expressions for the other first-order quantities.

A. Antiplasma

From Eq. ( 3) we get the relationship between the density perturbation and the potential perturbation:

ñ± = φ k • ∇ v f 0 ω -k • v dv. (15)
Using Eq. ( 15) together with Poisson's equation in the antiplasma case, we arrive at the dispersion relation

D(ω, k) ≡ 1 - 2ω 2 J0 k 2 k • ∇ v f 0 ω -k • v dv = 0, ( 16
)
where ω J0 = (4πGmn 0 ) 1/2 is the Jeans frequency. If k • v ω, we can expand the denominator in the above integral

1 ω -k • v ≈ 1 ω + k • v ω 2 + (k • v) 2 ω 3 + (k • v) 3 ω 4 + . . . .
When the equilibrium distribution f 0 is a function of v ≡ |v| only, then the dispersion relation can be expressed in the following way:

D(ω, k) ≈ 1 + 2ω 2 J0 ω 2 + 3k 2 v 2 ω 2 J0 ω 4 , (17)
where v 2 = v 2 f 0 (v)dv. Setting D = 0 then yields:

ω 2 = -2ω 2 J0 + 3k 2 v 2 , ( 18
)
and v 2 = v 2 th = k B T /m for a Maxwell-Boltzmann distribution with temperature T . This expression was derived earlier [32] using a fluid model. Equation (18) is almost identical to the dispersion relation for self-gravitating systems in the presence of a repulsive background, displaying the usual Jeans instability for large wavelengths (small k). The extra factor 2 comes from having two mobile species, so that one should use the reduced mass

m r = m i + m i -/(m i + + m i -) = m/2.

B. Bondi masses

In the Bondi case, the Vlasov-Poisson equations read as follows:

∂f ± ∂t + v • ∇f ± -∇φ • ∇ v f ± = 0, ( 19
) ∆φ = 4πGm (n + -n -) . ( 20
)
We note immediately that the Vlasov equations for f ± are completely identical. It is therefore useful to define g ≡ f + + f -and h ≡ f + -f -, which obey the equations

∂g ∂t + v • ∇g -∇φ • ∇ v g = 0, (21)
∂h ∂t + v • ∇h -∇φ • ∇ v h = 0, ( 22
) ∆φ = 4πGm ∞ -∞ h(r, v, t) dv. ( 23
)
We note that the equation for g is decoupled from the equation for h. The system ( 22)-( 23) evolves self-consistently independently of g, which is simply transported in the phase space like a passive scalar. The system ( 22)-( 23) looks identical to the equations of a single gravitational species evolving in standard Newtonian gravity, but there is an important difference: the distribution function h can take both positive and negative values and in general will have zero mean.

To perform the linear analysis, we consider Eqs. ( 22)-( 23) and expand all quantities to first order, remembering that at zeroth order h 0 = 0 and φ 0 = 0. It follows that the first-order equations reduce to the trivial free transport equation:

∂ h ∂t + v • ∇ h = 0. ( 24
)
There is no field-dependent term at first order and therefore no Jeans instability. Thus, nontrivial effects will only arise in the nonlinear regime.

C. Dirac-Milne

As mentioned earlier on, there is no homogeneous steady-state solution of Eqs. ( 5)-( 7), since antimatter acts repulsively with respect to both matter and antimatter. However, for the very same reason, antimatter tends to spread throughout space thus yielding an approximately uniform background with density n 0 , as was discussed in Sec. II B. In that case, one can drop the second Poisson equation (7) to obtain the following Vlasov-Poisson system for positive-mass matter

∂f + ∂t + v • ∇f + -∇φ + • ∇ v f + = 0, (25)
∆φ

+ = 4πGm(n + -n 0 ). ( 26
)
These are identical to the equations of ordinary gravity in the presence of a neutralizing background. Following the same procedure as done above for the antiplasma case, we obtain the usual Jeans dispersion relation

ω 2 = -ω 2 J0 + 3k 2 v 2 , (27)
which is identical to Eq. ( 18) except for the missing factor 2 in front of ω J0 .

IV. COMOVING COORDINATES AND SCALING

In this section, we introduce the comoving coordinates used to study an expanding universe in Newtonian gravity and adapt them to the case of a matter-antimatter universe with both attractive and repulsive gravity. Let us consider an expanding distribution of matter with spherical symmetry. Its gravitational field has only one component E r (r, t) that depends on time and on a single spatial variable r. This type of system was studied extensively in the past [START_REF] Rouet | Vistas in Astronomy[END_REF][START_REF] Rouet | Applying Fractals in Astronomy[END_REF][START_REF] Miller | [END_REF][36][37][38][39][40][41].

The equation of motion reads as (we ignore, for simplicity of notation, the signs of the masses; these will be reinstated later):

d 2 r dt 2 = E r (r, t), (28)
where E r = -∂ r φ is the gravitational field. We consider an expanding universe with scaling factor a(t) and transform space and time as follows:

r = a(t)r, (29)
dt = b 2 (t)d t, (30)
where comoving coordinates are denoted by an overcaret. Note that we also introduced a scaled time t, which defines a new time-dependent "clock". The velocity transforms as

dr dt = a b 2 dr d t + ȧr, (31)
where the dot stands for time differentiation with respect to t. The scaled equation of motion is then

d 2 r d t2 + 2b 2 ȧ a - ḃ b dr d t + b 4 ä a r = b 4 a 3 Êr , (32)
where Ê(r, t) is the scaled gravitational field. As the density must scale as ρ(r, t) = a 3 (t)ρ(r, t) in order to preserve the total mass, we scale the gravitational field as Ê(r, t) = a 2 (t)E r (r, t), so that Poisson's equation remains invariant in the scaled variables.

Next, we consider a locally planar perturbation embedded in this expanding universe and denote the corresponding comoving coordinate x. In this locally planar system, Poisson's equation can be approximated by its one-dimensional (1D) counterpart: ∂ x Ê = -4πGρ(x, t). The simulations presented in the next sections will all be performed in this 1D planar reference frame.

A. Symmetric universe with positive and negative masses

In a universe filled with equal amounts of positive-mass and negative-mass matter, the scaling factor a(t) should be linear in time, because over large distances attractive and repulsive gravitational fields cancel each other. Therefore, we take

a(t) = t/t 0 , (33)
where t 0 is an "initial" time corresponding to the epoch when the universe becomes transparent to radiation and therefore dominated by matter (in the Dirac-Milne cosmology its value is different from the 3.8 × 10 5 years after the big bang of standard cosmology [5] -more details on this will be given in the next section). In the following, the subscript "0" is used systematically to refer to quantities evaluated at this initial time. For the function b(t) that determines the scaling of time we use: b(t) = (t/t 0 ) 1/2 , which has the advantage of rendering the friction term time-independent in Eq. (32). We stress that, while the choice of a(t) is dictated by cosmological considerations, the choice of b(t) is rather arbitrary and can be made with the purely pragmatic aim of making the scaled equations easier to solve, either numerically or analytically. In this case the relationship between the real time t and the scaled time t is exponential: t = t 0 exp( t/t 0 ), and the equation of motion becomes

d 2 x d t2 + 1 t 0 dx d t = exp(-t/t 0 ) Ê. ( 34
)
Note that t = 0 when t = t 0 .

Defining the scaled velocity v, we obtain Hamilton's equations

dx d t = v (35) dv d t = exp(-t/t 0 ) Ê - v t 0 , (36)
from which one can deduce the corresponding 1D Vlasov-Poisson system for the scaled distribution functions F ± (r, v, t). For the antiplasma and Bondi cases, they read as follows:

∂F ± ∂ t + v ∂F ± ∂ x + m ± p m ± i e -t/t 0 ʱ ∂F ± ∂v - 1 t 0 ∂(vF ± ) ∂v = 0, (37)
∂ Ê ∂ x = -4πGm (n + -n-), (38)
where the values of m ± p and m ± i should be read from Table I. For the Dirac-Milne case, we assume that antimatter (negative mass particles) expands constantly under the repulsive action of all other particles (positive and negative) and can thus be modeled as a uniform repulsive background with n -(r, t) = n 0 /a 3 (t), where n 0 is a constant initial density. This is not completely exact, as antimatter is also ejected from overdense matterdominated regions (see Sec. II B) leading to a depletion zone, but still constitutes a very good approximation for the study of structure formation. It is easy to see that, with the scaling defined above, the scaled version of the Vlasov-Poisson equations ( 5)-( 6) for the positive-mass particles becomes:

∂F + ∂ t + v ∂F + ∂ x + e -t/t 0 Ê+ ∂F + ∂v - 1 t 0 ∂(vF + ) ∂v = 0, (39)
∂ Ê ∂ x = -4πGm (n + -n 0 ). ( 40
)
The forthcoming numerical simulations for the antiplasma, Bondi and Dirac-Milne scenarios, will be based on the solution of Eqs. ( 37)-( 38) and ( 39)- (40).

B. Einstein-de Sitter cosmology

For comparison, we briefly review the basic equations for the comoving coordinates in the standard Einstein-de Sitter (EdS) cosmology, without dark energy [38]. For a critical universe, the scale factor behaves as: a(t) = (t/t 0 ) 2/3 . In addition, taking b(t) = (t/t 0 ) 1/2 , the scaled equation of motion (32) becomes autonomous (all coefficients are time-independent):

d 2 x d t2 + 1 3t 0 dx d t - 2 9t 2 0 x = Ê . (41)
Poisson's equation in d dimensions can be solved exactly for a constant density n 0 , yielding the gravitational field: Ê = -ω 2 J0 x/d. At steady state, this field must exactly cancel the third term on the left-hand side of Eq. ( 41), yielding the following relationship between t 0 and ω J0 :

ω 2 J0 t 2 0 = 2 9 d. (42)
The corresponding Vlasov-Poisson equations in comoving coordinates read as [38]:

∂F ∂ t + v ∂F ∂ x + Ê ∂F ∂v - 1 3t 0 ∂(vF ) ∂v = 0, (43)
∂ Ê ∂ x = -4πGm (n -n 0 ), (44)
where the harmonic term in Eq. ( 41) has been incorporated into Poisson's equation [START_REF] Rouet | [END_REF].

The fact that all coefficients in Eqs. ( 41) or ( 43)-( 44) are time-independent signals that a stationary solution of the scaled equations represents a self-similar solution in the real coordinates. Note also that, for an expansion in t 2/3 , the product ω J (t) t is constant in time, i.e.: ω J (t)t = ω J0 t 0 .

Another relevant difference between the Einstein-de Sitter and the Dirac-Milne scenarios is that in the latter the scaled gravitational field decreases exponentially with the scaled time, whereas it is constant in the former. An important consequence of this mathematical fact (confirmed by the forthcoming simulations) is that gravitational structure formation in a Dirac-Milne universe stops after a time that is a small multiple of t 0 , while it continues indefinitely in the Einstein-de Sitter case.

V. TIMESCALES

Before presenting the results of the numerical simulations, it is instructive to identify the principal timescales involved in the problem under study. Two such timescales are the inverse of the initial Jeans frequency ω J0 = √ 4πGmn 0 and the time t 0 , defined as the time when the universe becomes transparent to radiation (end of recombination epoch). The third timescale is related to the importance of nonlinear effects.

To determine the Jeans frequency, we use the current estimate for today's average baryonic matter density of the universe, i.e. 0.25 protons/m 3 (there is no need for dark matter in the Dirac-Milne universe): ρ baryon (now) = 4.18 × 10 -28 kg/m 3 .

Following Benoit-Lévy and Chardin [5], recombination should occur at t 0 = 14 × 10 6 years in a Dirac-Milne universe, compared to 380 000 years for standard cosmology. Thus, t now = 14×10 9 y ≈ 10 3 t 0 , and since the expansion factor is linear in t, today's mass density should be multiplied by 10 9 in order to obtain ρ 0 at t = t 0 . Finally we get:

ω J0 t 0 ≈ 8.3. ( 45
)
At the present epoch, this parameter takes the value: ω J,now t now ≈ 0.26.

A third timescale is given by the initial density fluctuations, which determines the typical timescale over which nonlinear effects become important (known as the bounce time in plasma physics [42]). This nonlinear time t nl can be written as:

ω J0 t nl = n 0 /ñ λ ( 46
)
where ñλ is the amplitude of the initial density fluctuations (which depends on their wavelength λ).

In summary, we have identified three distinctive timescales for this problem:

• The inverse of the Jeans frequency ω J0 , which determines the rate of the Jeans instability, see Eqs. ( 18) and ( 27);

• The recombination time t 0 , which acts as an initial time in our model. Note that if ω J0 t 0 1 the friction term in Eq. ( 34) dominates from the start, so that cosmological structures will never have time to form. In contrast, if ω J0 t 0 1, structures will keep forming for very long times. The value ω J0 t 0 ≈ 8.3 that we find for the Dirac-Milne universe is sufficiently large to allow for structure formation at the beginning of the matter-dominated epoch. Such structure formation should stop at a time (≈ 10 9 y) relatively close to the present epoch, as will be shown in the forthcoming simulations;

• The nonlinear timescale t nl , which plays an important role in the Bondi scenario, for which no Jeans-like instability exists, and gravitational structures can only form via nonlinear couplings. Note that this "parameter" is actually a dynamical quantity, as the level of fluctuations changes in time in the case of an instability. Thus, its initial value may not be the relevant one for the later evolution in the antiplasma and Dirac-Milne scenarios, where the initial stages of the evolution are dominated by the Jeans instability.

In the forthcoming numerical simulations, we will use units in which n 0 = ω J0 = 1, so that the relevant dimensionless parameters are ω J0 t 0 and ω J0 t nl .

Note that in the EdS scenario there are only two relevant timescales, because the product ω J0 t 0 = 2/3 ≈ 0.82 is fixed by Eq. (42). Further, in the EdS case the product ω J (t) t = ω J0 t 0 is an invariant, i.e. it is the same for all epochs, whereas it decreases (as t -1/2 ) in the Dirac-Milne case. This simple fact illustrates a fundamental difference between the EdS and Dirac-Milne universes. In the former, ω J (t) t is of order unit at all epochs, so that structure formation occurs all along the life of this type of universe. In contrast, the Dirac-Milne universe starts with a relatively large value ω J0 t 0 = 8.3, which then decreases and becomes lower than unity near the present epoch, so that structure formation eventually stops.

To put it somewhat figuratively: the EdS universe is a universe that never ages and produces "offsprings" (i.e., structures) during its entire lifetime. The Dirac-Milne universe ages with time and after a certain epoch (corresponding roughly to ω J (t) t ≈ 1) it becomes "sterile" and begets no further gravitational structures.

VI. NUMERICAL RESULTS

In this section we present the results of numerical simulations of the three models with negative mass considered so far, namely the antiplasma, Bondi, and Dirac-Milne scenarios, with particular emphasis on the latter case. For comparison, we will also show some results obtained with a standard Einstein-de Sitter cosmology. The simulations were performed using an N-body code that solves the equations of motion [START_REF] Rouet | Applying Fractals in Astronomy[END_REF] for N interacting particles. For N → ∞ the N-body problem tends to its mean-field limit, described by the Vlasov-Poisson equations discussed in the preceding sections. Typical simulations employed N = 2.5 × 10 5 particles.

As mentioned in the preceding sections, we consider a 3D expanding spherically-symmetric universe and then study planar perturbations in the comoving coordinates. This effectively reduces the problem to one spatial dimension in the local comoving coordinate x, which will be represented in the numerical results. More details on the model can be found in Refs. [START_REF] Miller | [END_REF][36][37][38]. Since we use a 1D approximation, the "particles" are in fact infinite sheets of mass m (in absolute value). Boundary conditions are taken to be spatially periodic with period L. The initial velocity spread is very small (except for the Bondi case which, lacking a linear instability, must be strongly excited in order to observe some interesting dynamics). In all presented results, time is expressed in years elapsed since the Big Bang (t = 0). Density is expressed in terms of n 0 . The units of space and gravitational field are somewhat arbitrary: once we have fixed a certain unit of time (say ω -1 J0 , as was done in the code), space is measured in units of a certain length λ and the gravitational field in units of λω 2 J0 . Thus, λ disappears from the scaled equations of motion ( 32) and ( 41): one chooses an arbitrary value for it, and then all lengths are measured in terms of this unit.

In all cases, the initial density is the sum of a spatially uniform term ρ 0 = mn 0 and a small perturbation ρ with power spectrum |ρ k | 2 ∼ k p , where k is the wavenumber. Initial power spectra of this form, with p ∈ [0, 4], were used in a number of earlier works on structure formation [36,40]. In the present work, we take p = 2, which produces a spectrum that is largest at small wavelengths, and then study the clustering of matter at increasingly larger scales, until the size L of the box is reached.

We start by briefly showing a few results on the antiplasma and Bondi cases, then turn to a more extensive comparison between the Dirac-Milne and Einstein-de Sitter cosmologies.

A. Antiplasma and Bondi cases

As was discussed in Sec. III A, the antiplasma case is subject to a Jeans instability with growth rate of the order of the Jeans frequency ω J0 . Consequently, the system quickly departs from the initial equilibrium and develops nonlinear structures. These structures are shown in snapshots of the matter density and phase-space distributions taken at different instants, as depicted in Fig. 1.

As particles of like masses attract and particles of unlike masses repel, we observe the segregation of clumps of positive and negative matter, visible in the density plots. Interestingly, there is also a segregation in velocity space, with two counterpropagating streams made of positive and negative particles appearing. Thus, the Jeans instability extracts gravitational potential energy from the system (initially at rest in the comoving frame) and converts it into kinetic drift energy. Nevertheless, by zooming in on one of the streams, one still observes plenty of intricate gravitational structures, as is shown in Fig. 2. Here and in other forthcoming figures, the zoom region changes with time because it follows the trajectories of the particles that are contained initially in a given region of the phase space at a chosen time (here, the small rectangle in the phase space plot in Fig. 1.

This behavior may seem counterintuitive, but performing the linear analysis of an equilibrium consisting of two counterpropagating streams, reveals that such equilibrium is (for large enough drift velocities) stable, in contrast to the one-stream initial conditions which is unstable as was shown in Sec. III A. Thus, the system evolves from an unstable equilibrium (one stream with zero mean velocity) to a stable equilibrium (two streams with finite and opposite velocities). Note also that this is exactly the opposite situation compared to a two-component plasma, for which one-stream equilibria are stable whereas two-stream equilibria are unstable [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF][START_REF] Rouet | [END_REF]. This peculiar behavior in the phase space constitutes one further reason to rule our the antiplasma scenario for cosmological applications.

For the Bondi case, as we saw in Sec. III B, there is no linear instability. Therefore, we need to introduce very large fluctuations in the system's initial condition in order to observe some structure formation. If the initial fluctuations are low, then no gravitational structures are formed in the subsequent evolution. The most salient feature observed in the simulations is that, in contrast to the antiplasma case, positive and negative particles clump together to form "neutral" regions (see Fig. 3, which only shows one eighth of the total computational box for clarity). Also the structures form very early in the evolution and remain basically unchanged thereafter.

B. Dirac-Milne and Einstein-de Sitter universes

We now turn to the more interesting case of the Dirac-Milne universe and its comparison with the Einstein-de Sitter case. As was mentioned in Sec. IV A, one can simplify the treatment by assuming that the negative-mass component remains uniform in the comoving coordinates, so that only positive-mass particles need to be simulated. The initial conditions are identical for both sets of simulations, although the physical initial time (t 0 , representing the recombination epoch) is of course not the same.

The matter density and phase-space distributions for the Dirac-Milne case are provided in Fig. 4 and a zoom is shown in Fig. 5. At first sight, structure formation occurs similarly to the standard (EdS) case, driven by the Jeans instability [START_REF] Miller | [END_REF][36][37][38]. However, some crucial differences are noticeable. In particular, structure formation starts earlier and stops before the present epoch. Indeed, almost no significant changes can be seen in the density and phase-space distributions from the present epoch onwards. Mathematically, this is due to the fact that, in the comoving equations of motion (36), the relative impact of the friction term increases with respect to the gravitational field. Indeed, the scaled velocity in the phase-space plots of Fig. 4 initially increases under the action of the Jeans instability, but later decreases rapidly because of the effect of the friction term in the comoving equations of motion. Structure formation halts when this latter term becomes dominant.

Let us now compare this behavior with that observed for the Einstein-de Sitter universe, which was already simulated in earlier works [START_REF] Miller | [END_REF][36][37][38]. The matter density and phase-space distributions are plotted in Fig. 6, and two consecutive zooms in Figs. 7 and8. In accordance with previous results, the formation of gravitational structures continues for all times, as is particularly clear from the two zooms. The self-similar nature of the matter distribution in the phase space is also apparent: two snapshots taken at different scales (but at the same time) look qualitatively similar, indicating that there is no intrinsic spatial scale for this system. This fact has led us, in earlier works, to compute a fractal dimension for the matter distribution [36,38]. These features are the hallmark of a critical universe.

This different behavior is in accordance with our discussion of Sec. V, where we pointed out that the Dirac-Milne universe "ages" with time, meaning that it ceases to produce gravitational structures. In contrast, the EdS universe does not age and continuously generates new structures at all scales.

Another difference between the Dirac-Milne and EdS cases is that in the former the (scaled) velocity range of the particles decreases after an initial growth (Fig. 4), whereas in the latter it constantly increases (Fig. 6). This signals a change in behavior for the Dirac-Milne universe, which should roughly correspond to the epoch where structure formation ends. In order to understand the relevance of this quantity, let us go back to the definition of the scaled velocity, Eq. ( 31), which can be rewritten as

v = a b 2 v + H(t)r, (47)
where H(t) = ȧ/a is the Hubble parameter. The first term on the right-hand side of Eq. ( 47) represents the so-called peculiar velocity of the masses (residual velocity in the comoving reference frame), while the second is the Hubble redshift term. In Fig. 9, we plot the evolution of the "thermal" peculiar velocity, defined as v th ≡ (v -Hr) 2 1/2 , where the average is taken over all the particles. The difference is clear: for Dirac-Milne, the spread of the peculiar velocities first increases, reaches a peak, and then decreases continuously. Structure formation halts when v th starts decreasing, i.e., around t c ≈ 10 9 years. We also note that this time corresponds roughly to the time when the parameter ω J (t) t reaches unity. Indeed, for the Dirac-Milne case:

ω J (t) t = ω J0 t 0 t t 0 -1/2 ≈ 8.3 t t 0 -1/2 (48)
becomes equal to unity for t c ≈ 69t 0 ≈ 10 9 y. This is in agreement with our earlier discussion in Sec. V.

Peculiar velocities of clusters of galaxies have been estimated by measuring their impact on the cosmic microwave background spectrum [45]. A compilation of several peculiar velocity surveys was published in recent years [46]. These surveys consistently showed that the observed velocities exceed the prediction of the ΛCDM model by a factor of ≈ 5. Although the discrepancy between the Dirac-Milne and EdS predictions is even higher, our result (Fig. 9) goes in the right direction. A proper comparison between the ΛCDM and Dirac-Milne cosmologies will be studied in a forthcoming work. Finally, we show in Fig. 10 the power spectra of the matter density for the EdS and Dirac-Milne cases. Both spectra begin (at t = t 0 ) as a power law |ρ k | 2 ∼ k p with p = 2 up to k max = 2π/d, where d ≈ 1 is the average initial interparticle distance. For larger wavenumbers the spectrum is flat (white noise), simply because there are not enough particles to resolve such short distances.

Initially, the low wavenumber amplitudes (leftmost part on Fig. 10) are amplified but keep their shape intact: this is the linear stage of the instability. The linear growth rate of gravitational structures observed in these simulations is significantly larger for the Dirac-Milne case compared to the EdS one. However, it is also well-known that the linear growth is faster for the concordant ΛCDM model than for the EdS case -see e.g. Ref. [47] -so that our Dirac-Milne result again appears to go, at least qualitatively, in the right direction.

As time increases, nonlinearities start to play a predominant role and the spectra take a characteristic shape that is qualitatively similar both in the Dirac-Milne and the EdS cases, as can be seen from Fig. 10. At low wavenumbers (region I), the spectrum is still of the type |ρ k | 2 ∼ k 2 , a remnant of the initial condition. At larger wavenumbers, an intermediate region (II) appears, characterized by a power-law spectrum with negative exponent, slightly steeper for the Dirac-Milne (p = -0.78) compared to the EdS (p = -0.67) case. This power-law region is rather robust (it extends over roughly four decades) and constitutes the signature of hierarchical clustering in the phase space, as was seen in the previous plots (Figs. 4 and6). For even larger wavenumbers, the spectra are again flat (region III).

The spectra display a peak separating regions I and II, which can be taken as the typical cluster size at a certain epoch. In both cases, the peak initially moves towards smaller and smaller wavenumbers (larger scales), although quite noticeably faster for the Dirac-Milne case (Fig. 11). But for longer times (> 10 9 y) the behaviors diverge: for EdS, the position of the peak k peak (t) keeps moving to larger and larger scales (eventually reaching the size of the computational box, at which point the simulation would not be valid anymore), whereas for Dirac-Milne it saturates at a constant value. A similar saturation effect is also expected for ΛCDM. At the present epoch, the typical cluster size in comoving coordinates ( ∼ k -1 peak ) is almost two orders of magnitude larger for the Dirac-Milne universe compared to the EdS case. We also stress that such cluster size is determined primarily by nonlinear effects occurring during the matter-dominated epochs.

In summary, Figs. 10-11 clearly show that, compared to EdS, structure formation begins earlier and proceeds initially faster in the Dirac-Milne scenario. However, the production of new gravitational structures does not continue forever (as in EdS), but rather stops a few billion years after the Big Bang. This confirms the intuition that structure formation is linked to the dimensionless parameter ω J (t)t, as mentioned in Sec. V. For the Dirac-Milne case, this quantity is large at recombination (ω J0 t 0 ≈ 8.3) so that structures initially form quickly and efficiently; then it decreases (as t -1/2 ) and, when it falls below unity, structure formation halts. In contrast, for EdS, ω J (t)t = ω J0 t 0 = 2/3 is constant and of order unity for all times, so that structure formation occurs at a moderate rate for all epochs.

As a closing remark to this section, we note that here we compared the Dirac-Milne universe to the EdS case mainly for reasons of computational simplicity. This should be accurate enough for matter-dominated epochs (i.e., between recombination and the vacuum epoch), after which the full ΛCDM model should be used instead of the simpler EdS. Both the Dirac-Milne and the ΛCDM models describe a universe that "ages" with time, meaning that structure formation halts after a certain time (see Sec. V). Interestingly, the Dirac-Milne cosmology can reproduce this feature of the observed universe without resorting to dark energy. The timescales of structure formation in the Dirac-Milne and ΛCDM universes are not too different either. Indeed, for ΛCDM, structure formation ends when: a(t)/a(t now ) = (t/t now ) 2/3 = (Ω m /Ω Λ ) 1/3 , which (using the current estimations for the normalized densities: Ω m = 0.3, Ω Λ = 0.7) yields t = t Λ ≈ 0.65 t now ≈ 9×10 9 y. This is larger than the equivalent "freezing" time t c ≈ 10 9 y defined above for the Dirac-Milne cosmology, but in both cases structure formation ended before the present epoch. Also keep in mind that this is a rather crude estimate of t c . A closer inspection of the phase-space portraits (Fig. 4, and other snapshots not shown here) suggests that structure formation stops some time between 10 9 and 10 10 years after the Big Bang, which is approximately compatible with ΛCDM.

VII. DISCUSSION AND CONCLUSIONS

The scope of this work was twofold. Firstly , we developed a general formalism that allows for the existence of negative masses in a Newtonian framework. The family of models that we obtained includes not only those that can be defined by the different signs of the active, passive, and inertial masses, but also a more general class that can only be characterized by a set of two Poisson's equations for the gravitational potential. These models could be viewed as the Newtonian limit of some bimetric extension of General Relativity. There is a total of seven nontrivial models that can be constructed in this way.

We produced circumstantial evidence that only one of these models possesses the right features that make it compatible with the alternative "Dirac-Milne" cosmology recently proposed by Benoit-Levy and Chardin [5]. This scenario assumes a universe composed of equal amounts of matter and antimatter, in which the latter antigravitates. In the model that we propose to simulate the Dirac-Milne scenario, all gravitational interactions between matter and antimatter are repulsive, except for matter-matter interactions, which must of course be attractive. These features imply that antimatter is ejected from overdense matter-dominated regions (galaxies), while spreading out in underdense regions to form a uniform dilute repulsive background, with very few annihilation events occurring, in accordance with observations. In a cosmological context, it is this repulsive background that induces an expansion rate linear in time, a(t) ∼ t, for the Dirac-Milne universe, faster than the expansion rate of the matter-dominated epochs in the standard cosmological model, which scales as a(t) ∼ t 2/3 . Hence, the Dirac-Milne universe does not suffer from the horizon problem and does not need primordial inflation to explain the current homogeneity at large scales. Indeed, the dilute repulsive background could be viewed as a cosmological constant that decreases in time, with the corresponding vacuum energy decreasing volumetrically as a -3 , see Eqs. ( 39)- (40).

The second purpose of this work was to study the implications of the Dirac-Milne cosmology on gravitational structure formation. To do this, we devised a local 1D model embedded in a 3D spherically expanding universe and performed N-body simulations of both a Dirac-Milne and an Einstein-de Sitter (Ω m = 1, Ω Λ = 0) universe. In both cases we observed gravitational structure formation, with clusters and subclusters developing from an almost uniform initial condition. After the nonlinear stage of the evolution is reached, both models display a power-law behavior in the wavenumber spectrum of the matter density -with, however, one crucial difference: whereas for EdS the formation of structures continues indefinitely, for the Dirac-Milne universe it stops a few billion years after the Big Bang, after which time the gravitational structures stay frozen in comoving coordinates. Today, the typical size of the structures should be over an order of magnitude larger in the Dirac-Milne universe.

In this work, for the sake of computational simplicity, we compared the Dirac-Milne universe to the Einstein-de Sitter case, which should be accurate for matter-dominated epochs before the effect of the cosmological constant has become dominant. Nevertheless, for a proper analysis of structure formation, the Dirac-Milne scenario should be compared to the standard ΛCDM model. It should be possible to do this in the framework of the present 1D approach, thus keeping the computational cost relatively low. The results of such comparison will be reported in future works. For the time being, we can conclude that the Dirac-Milne universe generates hierarchical structures compatible with those observed in our universe, and whose formation halts at an epoch earlier than the present one, as is also the case for the standard cosmological model. All in all, these preliminary results are encouraging and deserve to be further verified and quantitatively compared both with observations and with the predictions of the standard ΛCDM scenario.

FIG. 1 :FIG. 2 :

 12 FIG.1: Antiplasma case. Left frames: Matter density ρ(x, t) for positive mass (red) and negative mass (blue) particles at four different times. Right frames: Corresponding distributions in the scaled phase space. The small rectangle (right column, third panel) indicates the region that is zoomed in on Fig.2. Time is measured in years elapsed after the Big Bang.

FIG. 3 :

 3 FIG.3: Bondi case. Left frames: Matter density ρ(x, t) for positive mass (red) and negative mass (blue) particles at three different times. Right frames: Corresponding distributions in the scaled phase space. Time is measured in years elapsed after the Big Bang.

FIG. 4 :

 4 FIG. 4: Dirac-Milne case. Left frames: Matter density ρ(x, t) for the positive-mass particles at five different times. Right frames: Corresponding distributions in the scaled phase-space. The fourth set of frames from the top corresponds roughly to the present epoch.

FIG. 5 :

 5 FIG.5: Dirac-Milne case. Left frames: Matter density ρ(x, t) for the positive-mass particles at three different times. Right frames: Corresponding distributions in the scaled phase-space. The figures are zooms of the region denoted by a small rectangle in Fig.4(bottom right frame).

FIG. 6 :FIG. 7 :

 67 FIG.6: Einstein-de Sitter case. Left frames: Matter density ρ(x, t) for the positive-mass particles at five different times. Right frames: Corresponding distributions in the scaled phase-space. The third set of frames from the top corresponds roughly to the present epoch.

FIG. 8 :FIG. 9 :

 89 FIG.8: Einstein-de Sitter case. Further zooms of the region indicated by a rectangle on Fig.7.

FIG. 11 :

 11 FIG.11: Wavenumber corresponding to the peak in the power spectrum for the Dirac-Milne and Einstein-de Sitter universes as a function of time, in comoving coordinates.

TABLE I :

 I . Signs of the three types of mass for the four cases considered in the text. The absolute values of the masses are supposed to be the same. As only the ratio m p /m i is relevant, we assume for all cases m i > 0, with no loss of generality.

		Case	Name	Active grav. mass Passive grav. mass Inertial mass
	Matter	A	Standard	+	+	+
		B Antiplasma	-	-	+
	Antimatter C	Bondi	-	+	+
		D Antiinertia	+	-	+

TABLE II :

 II Interactions between matter (+) and antimatter (-) particles in the Dirac-Milne universe.

	Type of matter Type of matter Interaction
	+	+	Attraction
	-	-	Repulsion
	-	+	Repulsion
	+	-	Repulsion

  Power spectra of the matter density for the EdS (top panel) and Dirac-Milne (bottom panel) universes. Time is expressed in years after the Big Bang. The thick lines correspond to the present epoch (t = 14 × 10 9 y).

		10 4
		10 2						t=3.8e+05 t=1.4e+07
		10 0						t=1.2e+08 t=1.1e+09
	k 2 >/N <ρ	10 -2						t=1.4e+10 t=8.5e+10 t=7.6e+11
		10 -4						t=4.6e+12 t=2.9e+13
		10 -6						slope : 2.01 slope : -0.65
		10 -8	10 -3	10 -2	10 -1	10 0	10 1	10 2
					k
		10 4						t=1.4e+07
								t=2.0e+07
		10 2						t=2.9e+07 t=4.2e+07
	k >/N 2 <ρ	10 -2 10 0						t=6.0e+07 t=1.2e+08 t=3.7e+08 t=1.4e+10
		10 -4						t=7.6e+11 t=2.9e+13
		10 -6						slope : 2.00 slope : -0.79
		10 -8	10 -3	10 -2	10 -1	10 0	10 1	10 2
					k
	FIG. 10:

Acknowledgments

The numerical simulations were partly performed on the computer cluster at the Centre de Calcul Scientifique en région Centre-Val de Loire (CCSC).