Phase transition in Modified Newtonian Dynamics (MONDian) self-gravitating systems

We study the statistical mechanics of binary systems under gravitational interaction of the Modified Newtonian Dynamics (MOND) in three-dimensional space. Considering the binary systems, in the microcanonical and canonical ensembles, we show that in the microcanonical systems, unlike the Newtonian gravity, there is a sharp phase transition, with a high-temperature homogeneous phase and a low temperature clumped binary one. Defining an order parameter in the canonical systems, we find a smoother phase transition and identify the corresponding critical temperature in terms of physical parameters of the binary system.


Introduction
The growth of structures from the initial condition in the early Universe to the galaxies and clusters of galaxies is addressed by the standard model of cosmology. In this scenario, the origin of the structures is quantum fluctuations of a scalar field, the so-called inflaton field. The amplitude of the structures grows after the end of inflation. The standard paradigm for structure formation, ΛCDM [1], has a very good agreement in the early Universe from the CMB observations [2] because the free parameters are tuned to fit the CMB. However, ΛCDM runs into severe difficulties with local and large scale structures [3][4][5]20].
There is another approach to structure formation theory from the statistical mechanics point of view, where the structures form when a phase transition occurs in gravitating systems [6][7][8][9]. One classical example in statistical mechanics is the 2D self-gravitating system with a logarithmic gravitational potential [7]. In this approach, one takes an ensemble of N-body particles as a thermodynamical system, where from the partition function one can derive the thermodynamical quantities. We note that in this system all the particles are under their mutual gravitational interaction; this has to be taken into account in calculating the partition function. It has been shown that for such a system there exist two phases of (a) high-temperature gaseous phase, and (b) clumped low-temperature phase. This approach has an analytical solution only in 2D (logarithmic) gravity. To extend it to 3D, one has to deal with a simpler thermodynamical approach such as studying two-body systems [8].
In the standard theory of structure formation, to have a compatible theoretical result with observation, we need to have a dark matter component in addition to the baryonic component of the cosmic fluid. The effect of dark matter in structure formation is important when the universe was dominated by radiation whose pressure repelled the baryonic matter and prevented the formation of baryonic structures. The perturbation in the dark matter fluid of the Universe (unlike the baryonic and radiation components) would grow through the gravitational instability to form over-dense regions. The consequence of dark matter structures is that after recombination, the gravitational potential of dark matter could accumulate the baryonic matter to form the baryonic part of galaxies within dark matter halos [10,11]. There is another approach to deal with the dynamics of large scale structures by replacing the dark matter with a modification to the gravity law. In some of these models, the extra degrees of freedom such as the scalar or the vector sectors are considered for the gravitational field. A generic formalism of this theory is the Scalar-Vector-Tensor theory which is the so-called Modified gravity model (MOG) [12]. Although this theory can explain the dynamics of galaxies and clusters of galaxies without a need for dark matter [13,14], it predicts unexpected mass-to-light ratio [15]. There are also non-local gravity models, where in one of these theories, the Einstein gravity has been supplemented by the non-local terms, in analogy to the non-local electrodynamics [16]. This theory also educes to the standard Poisson equation in the weak field approximation with an extra term that plays the role of dark matter. The theory provides compatible dynamics for galaxies and clusters of galaxies without the need for a dark matter component [17]. There is also another popular model which is called Modified Newtonian Dynamics (MOND) [18]. In this theory, Newton's second law or usually the Poisson equation is modified in systems with accelerations smaller than a universal acceleration, a 0 . This theory also provides compatible dynamics to the spiral and elliptical galaxies [19]. It is worth mentioning that the main problem of modified gravity models is that they can not explain the observational data on different scales, though see [20] for a possible hybrid solution using both MOND and hot dark matter in the form of sterile neutrinos.
In this work, we study the phase transition for a binary interacting via MONDian gravity, in non-expanding and expanding spaces. For simplicity, we start with an ensemble of two-body objects instead of the N-body system and investigate the phase transition by decreasing the temperature of the system.
The rest of the paper is organized as follows: In section II, we review the statistical mechanics of binary interacting systems under a Newtonian potential in the microcanonical ensemble. In section III, we present the statistical mechanics of a binary system in MOND. Section IV is devoted to studying the statistical mechanics of MONDian systems in the canonical ensemble. We then study and discuss the influence of expanding Universe with scale factor a(τ) on the critical temperature of the phase transition in section V. Section VI summarizes the paper and discusses possible further research.

Statistical mechanics of a self-gravitating binary under Newtonian potential: Microcanonical ensemble
The statistical mechanics of self-gravitating systems has been the subject of attention for many years [8,9,24]. Such systems have distinguished physical properties due to the long-range nature of the gravitational force. At thermal equilibrium, these systems are not spatially homogeneous, and intrinsic inhomogeneity character suggests that fractal structures can emerge in a system of gravitationally interacting particles [9]. Here we review the statistical mechanics of a self-gravitating binary system in Newtonian gravity in the microcanonical ensemble [8]. We start with the Hamiltonian of a two-body system, H(P, Q; p, r) = P 2 2M where in 3D the potential is V(r) = −Gm 2 /r, (Q, P) are coordinates and momentum of the center of mass, and (r, p) are the relative coordinates and momentum with the reduced mass. In what follows, for the sake of simplicity we take the mass of these two bodies to be identical, (i.e. M = 2m and µ = m/2). We assume a spherical shape for the two objects with radius of b/2 and that the two-body system is confined in a spherical box of radius R, where r in equation (1) varies within the interval of (b, R). The volume associated with a constant energy of this system (i.e. H = E) in the phase-space [25,26] is given by the density of states First we integrate over Q, which leads to 4πr 3 /3, and inserting the explicit form of the Hamiltonian leads to the following relation: Now we perform integration over p-space. Using the new variable x ≡ p 2 /µ, the square-root term appears from the integration as a result of a property of the Dirac δ-function.
This time we integrate over P-space similar to the previous step in equation (3). The result is where A = 64π 5 m 3 /3. In the case of Newtonian gravity, the potential energy is V(r) = −Gm 2 /r and since the kinetic energy is always positive definite (i.e E − V(r) > 0), the upper limit of the integral (r max ) should be taken in such a way that guarantees the positive sign of the kinetic energy. The upper bound of integration for the following ranges of energy is given by, The second condition has no specific meaning in astrophysics, though it has an analogy with the standard thermodynamics where the size of the box is R. We assume that the container has a fixed volume while we can increase the velocity of the particles, and they are constrained to stay in this volume. Integrating (5) results in From the density of states function g(E) which represents the phase space volume covered by this system with energy E, we can calculate the entropy and the temperature of the system according to the following relations (with Boltzmann constant K B = 1):  Figure 1. Dimensionless temperature versus dimensionless energy (Eq. 9) for a binary system interacting under Newtonian gravity (dashed orange curve) and MONDian gravity (solid blue curve). As expected, there is a flat part for MOND due to entering the DML. Since for constant temperature the energy has a finite change there must be a sharp phase transition.
From equation (7), we obtain the dimensionless temperature for the interval (−Gm where the dimensionless temperature t and energy are defined as Also for the second interval of −Gm 2 /R < E < ∞ in (7), t( ) is given by In Figure (1), we depict t( ) in terms of [25]. It can be seen that the specific heat is positive along AB and CD while it is negative along BC. For a system with energy in the range AB, the two solid spherical objects are in contact, and increasing the energy of the system increases the kinetic energy, or in other words, the temperature of the system. Over the range BC, the two objects detach from each other and start circular motion according to the mutual gravitational force between them. For the CD path in Figure (1) the total energy is larger than zero (i.e. E > 0) and the two objects decouple from each other and behave as free particles.

Statistical mechanics of a self-gravitating binary in MOND: Microcanonical ensemble
In the Modified Newtonian Dynamics (MOND), proposed to solve the dark matter problem on galaxy scales [18], the second law of Newtonian mechanics is modified for small accelerations. The characteristic acceleration for this modification is a 0 = 1.2 × 10 −10 ms −2 [21] where for a ≤ a 0 , the definition of force is modified to F = maµ(a/a 0 ). Here µ(x) is larger than unity for smaller accelerations and is unity for the large accelerations (i.e. a a 0 ). This model can explain the rotation curve of spiral galaxies and predicts the Tully-Fisher relation [27]. However, in this theory, energy and momentum are not well-defined [28].
One of the simple solutions to this problem is that one may interpret MOND as modified gravity rather than modified dynamics. The gravitational acceleration in the deep MOND limit where g a 0 will be [18] The subscript DML stands for the deep-MOND limit, in which the potential from the gravitational acceleration is given by For an isolated binary system in the DML, the gravitational force between the two particles scales as F ∼ 1/r which for objects in circular motion results in a constant velocity, compatible with the flat rotational curves of galaxies. We note that for g a 0 the Newtonian gravity is recovered (i.e. g = g N ) and the potential will behave as φ N = −Gm/r.
Since in MOND the acceleration of a test particle in a gravitating system is stronger than in Newtonian gravity, in addition to studying the dynamics of a system, we can investigate the formation of structures in these two scenarios. For a self-gravitating system, the free-fall time scale represents the strength of clustering of a structure. The ratio of characteristic timescales in Newtonian and MOND gravity is [29,30], For g N < a 0 we would expect a shorter time-scale for the clustering in MOND compare to the Newtonian gravity. We refer interested readers to [23] for a thorough investigation into the free-fall timescales in MOND.
In order to have a continuous transition from the DML to the Newtonian gravity, one needs a transition function µ(a) [31]. There are various transition functions. One of the simple models is µ(x) = x/(1 + x) or a sigmoid function, other transition function can be found in [32,33]. The sigmoid function better fits the observational data [22]. Now let us study the statistical mechanics of two-body objects in MOND in the microcanonical ensemble. The enumeration of the total number of possible states (phase volume) for a binary system with energy E is given by By integrating over P, p and Q, ( similar to the calculations in section (2) for Newtonian gravity, the function g(E) is given by where A = 64π 4 m 3 /3. Taking into account the Newtonian and DML phases for the gravitational potential, the total potential is given in the following two domains, Where the MOND radius r M = √ Gm/a 0 is the scale of transition between the two domains. In order to calculate integral (15), one needs precisely determine r max in terms of energy ranges, i.e., Here the bound of r max results from the positive sign of the kinetic energy. We define the temperature of system from equation (8) and plot t = t( ) for three energy ranges in Fig (1). Here we adopt R/b = 10 10 and e 2 /e 1 = 100 where e 1 = m √ Gma 0 and e 2 = Gm 2 /b. In this figure, the (t, ) diagram for MONDian gravity is almost similar to that for Newtonian, except for just slightly below 0, where for the case of Newtonian gravity t(0) = 0 but for MONDian gravity t(0) > 0. In fact, for MONDian gravity (the solid curve), there is a flat part, in which for constant temperature, the energy has a finite change, indicating a sharp phase transition where the heat capacity diverges.
We can interpret this area as when the binary system enters the deep-MOND where the rotation velocity of the binary objects around their center of mass is v rot = √ Gma 0 /2. In this case, with increasing energy of the system, the orbital size of the binary increases; however, the kinetic energy, which represents the temperature of the system, remains constant. In the next section, we study statistical mechanics and phase transition of a binary system in the MONDian gravity for the canonical ensemble.

Statistical mechanics of a self-gravitating binary in MOND: Canonical ensemble
Let us at first assume an ensemble of binary objects under Newtonian gravity. The system is composed of an ensemble of thermalized binaries with an associated temperature. The partition function associated with a binary system in this ensemble is given by where the parameters and the Hamiltonian are defined in Sec. (2). Integrating over momenta P, p and position Q, equation (18) simplifies to [25,29]: In dimensionless form, using the definition of t as introduced in Sec. (2), the partition can be written as dt and has no known analytical solution, so here we solve it numerically and calculate the mean energy of the system by using: In order to identify a phase transition, we calculate numerically the derivative of energy with respect to the thermodynamical variables. If this quantity diverges or becomes discontinuous, the system undergoes a phase transition [34]. We calculate the derivative of energy with respect to the temperature, which is defined as the specific heat and is shown in Fig. (2). From the numerical calculation of specific heat c v in terms of the temperature, we find that it has a peak at t = t critical . To understand the nature of the detected phase transition, we define an order parameter and study its behavior near the critical temperature. Here we take the mean distance between the companions of a binary system as the order parameter and define it as [7]: For the first-order phase transition, we expect to have an abrupt change in the order parameter. For instance, if we apply it for matter in the liquid and gaseous phases, this parameter sharply changes when a liquid changes to the gaseous state. For the second order phase transition, the order parameter will be a continuous function at the critical temperature. For simplicity, we rewrite Eq. (22) in dimensionless form as: and calculate < x 2 > numerically allowing us to plot the order parameter as a function of temperature in Fig. (3). We notice that there is the phase transition for the order parameter at exactly the same temperature that the specific heat has a peak, whereas the order parameter is a differentiable function at the critical temperature. The order parameter shows that a high-temperature system has a homogeneous phase, while at low temperature < x 2 > vanishes. Now we perform a similar calculation for the temperature dependence of the order parameter in the combination of deep MOND and Newtonian gravity, where for small and large accelerations (i.e a < a 0 and a > a 0 ), DML and Newtonian gravity will have dominant contributions, respectively. Since the potential is a function of distance, we need to make an approximation before calculating the partition function. The potential is defined as φ(r) = − F.dr, we can break this integration into separate three parts with different gravitational potentials. For r r M we have Newtonian gravity, and r r M potential is in DML, and between we are in a regime that is a combination of Newtonian and MONDian, i.e., Three terms in the integration belong to the Newtonian, DML, and mixing of Newtonian and DML, respectively. Here δ is a small constant. The result of integration is, where C is the constant of integration, and we chose it to be ∼ − ln(R). We note that the third term in r.h.s of Eq. (26), that is r M +δ r M −δ F.dr, depends on the interpolating function, (F = mµ(a/a 0 )). A simple approach is by choosing a proper µ with a very fast transition from Newtonian regime to MONDian regime. Then we can ignore this term and the potential simplifies to An alternative approach is to use a simple interpolating function for a single point mass using a hyperbolic substitution [32,33], wherer ≡ 2r r M and r M ≡ Gm a 0 .
We provide the influence of the interpolating function on the critical temperature, see below.
For simplification, we continue with Eq. 26. As a result, the partition function of the ensemble of the binary system in the MONDian gravity (using the MONDian potential in Eq. 18) by integrating over variables P, p and Q is Which in dimensionless representation it simplifies to: Here e 1 = m √ Gma 0 and e 2 = Gm 2 /b. Following the same procedure as the Newtonian case, we calculate the specific heat and order parameter and plot them in Figs. 2 & 3, where the specific heat and the order parameter are given by: < r 2 >= Eq. 31 can also be rewritten in the dimensionless form as: x 2 exp − e 1 ln(bx/R) Comparing Figs. 2 and 3, the phase transition temperature is identical whether we obtain it from the divergence of C v or from a fast change of < x 2 >. Also, we note that C v is always positive for the canonical case, unlike the micro-canonical case.
If we use the specific interpolating function that is introduced in Eq. (27), we observe a similar behavior, but at a different critical temperature (see Fig 4). We see that qualitatively the critical behavior of systems is the same but the interpolating function affects the details of the phase transition.

Statistical mechanics of a self-gravitating binary in MOND: Comoving coordinates
In an expanding universe with characteristic scale factor a(τ), the physical coordinates of r is related to the comoving coordinate q as [35], r ≡ a(τ)q .
The Hamiltonian of a binary system with potential given in Eq. (26) can be written in the comoving coordinate as This Hamiltonian in an expanding universe can also be obtained from the Minkowski-Hamiltonian with the following replacements:  Figure 4. Order parameter as a criterion to detect the phase transition. This plot shows that the variance of the distance between the particles changes fast but remains differentiable near the phase transition.
We can interpret the renormalization process by assigning a dynamics to these parameters. For instance the m → ma(τ) 2 renormalization is related to the kinetic energy of particles that decrease as 1/a 2 with the expansion of the Universe. Here we have the effective gravitational constant (G → Ga(τ) −5 ) as well as the acceleration parameter of MOND (a 0 → a 0 a(τ) −1 ) changing with the scale factor.
To consider a self-gravitating binary at any time in approximate thermal equilibrium, we assume that the characteristic time of the particle motions under their mutual gravitation is shorter than the time variation of the scale factor. This hypothesis is valid for structures that are almost decoupled from the expansion and become virialized. Given this, the Eqs. 28 & 31 in the comoving frame will be as follow: Eqs. 36 & 37 can be rewritten in dimensionless form as: where x ≡ q b . In Fig. (5) , the order parameter < x 2 > for various scale factors a(τ) = 1, 0.1, 0.05 is plotted in the left panel, while the right panel demonstrates critical temperature versus scale factor, showing that the critical temperature decreases as the scale factor increases.

Conclusions
In this work, we have studied thermodynamical phase transition under MONDian gravity. In addition influence of the cosmic expansion on the critical temperature of the detected phase transition has been studied. We have shown that in the microcanonical ensemble of binary systems under MONDian gravity, a sharp phase transition is not present in the Newtonian gravity. Furthermore, we find a smoother phase transition with finite critical temperature by studying the specific heat C v and an order parameter of a binary system in a canonical ensemble. One interesting result is that although both Newtonian and MONDian systems experience a phase transition in the canonical ensemble, they have different critical temperatures. The next steps in our research will be considering all interactions in the N-body system with the MONDian gravity and its connection with cosmological MOND [15]. In other directions, it would also be interesting to study the fractal structure of clusters in the crumpled phase of the system, as well as the equation of state (using the partition function in the grand canonical ensemble) and complexity in the view of [36]. Investigation of the critical temperature of rotating MONDian self-gravitating systems will be another interesting problem.