A phase transition in a Curie-Weiss system with binary interactions

A single-sort continuum Curie-Weiss system of interacting particles is studied. The particles are placed in the space $\mathbb{R}^d$ divided into congruent cubic cells. For a region $V\subset \mathbb{R}^d$ consisting of $N\in \mathbb{N}$ cells, every two particles contained in $V$ attract each other with intensity $J_1/N$. The particles contained in the same cell are subjected to binary repulsion with intensity $J_2>J_1$. For fixed values of the temperature, the interaction intensities, and the chemical potential the thermodynamic phase is defined as a probability measure on the space of occupation numbers of cells, determined by a condition typical of Curie-Weiss theories. It is proved that the half-plane $J_1\,\times\,$\textit{chemical potential} contains phase coexistence points at which there exist two thermodynamic phases of the system. An equation of state for this system is obtained.


Introduction
The mathematical theory of phase transitions in continuum particle systems has much fewer results as compared to its counterpart dealing with discrete underlying sets like lattices, graphs, etc. It is then quite natural that the first steps in such theories are being made by employing various mean field models. In [1], the mean field approach was mathematically realized by using a Kac-like infinite range attraction combined with a two-body repulsion. By means of rigorous upper and lower bounds for the canonical partition function obtained in that paper, the authors derived the equation of state indicating the possibility of a first-order phase transition. Later on, this result was employed in [2,3] to go beyond the mean field, see also [4,5] for recent results. Another way of realizing the mean-field approach is to use Curie-Weiss interactions and then appropriate the methods of calculating the asymptotics of integrals, cf. [6]. Quite recently, this way was formulated as a coherent mathematical theory based on the large deviation techniques, in the framework of which the Gibbs states (thermodynamic phases) of the system are constructed as probability measures on an appropriate phase space, see [7,section II].
In this work, we introduce a simple Curie-Weiss type model of a single-sort continuum particle system in which the space R d is divided into congruent (cubic) cells. For a bounded region V ⊂ R d consisting of N such cells, the attraction between every two particles in V is set to be J 1 /N, regardless of their positions. If such two particles lie in the same cell, they repel each other with intensity J 2 > J 1 . Unlike [1], we deal with the grand canonical ensemble. Therefore, our initial thermodynamic variables are the inverse temperature β = 1/k B T and the physical chemical potential. However, for the sake of convenience we employ the variables p = βJ 1 and µ = β × (physical chemical potential) and define single-phase domains

The model
By N, R and C we denote the sets of natural, real and complex numbers, respectively. We also put N 0 = N ∪ {0}. For d ∈ N, by R d , we denote the Euclidean space of vectors x = (x 1 , . . . , x d ), x i ∈ R. In the sequel, its dimension d will be fixed. By dx we mean the Lebesgue measure on R d . In this section, we use some tools of the analysis on configuration spaces whose main aspects can be found in [9].

The integration
For n ∈ N, let Γ (n) be the set of all n-point subsets of R d . Every subset of this kind, γ ∈ Γ (n) (called configuration), is an n-element set of distinct points x ∈ R d . Let also Γ (0) be the one-element set consisting of the empty configuration. Every Γ (n) is equipped with the topology related to the Euclidean topology of R d . Then, we define that is, Γ 0 is the topological sum of the spaces Γ (n) . We equip Γ 0 with the corresponding Borel σ-field B(Γ 0 ) that makes (Γ 0 , B(Γ 0 )) a standard Borel space. A function G : Γ 0 → R is B(Γ 0 )-measurable if and only if, for each n ∈ N, there exists a symmetric Borel function G (n) : (R d ) n → R such that G(γ) = G (n) (x 1 , . . . , x n ) for γ = {x 1 , . . . , x n }. For such a function G, we also set G (0) = G( ). The Lebesgue-Poisson measure λ on (Γ 0 , B(Γ 0 )) is defined by the relation which should hold for all measurable G : Γ 0 → R + := [0, +∞) for which the right-hand side of (2.1) is finite. For some c > 0, we let ∆ = (−c/2, c/2] d ⊂ R d be a cubic cell of volume υ = c d centered at the origin. Let also V ⊂ R d be the union of N ∈ N disjoint translates ∆ of ∆, i.e., As is usual for Curie-Weiss theories, cf. [6,7], the form of the interaction energy of the system of particles 23502-2 placed in V depends on V. In our model, the energy of a configuration γ ⊂ V is where I ∆ is the indicator of ∆ , that is, I ∆ (x) = 1 if x ∈ ∆ and I ∆ (x) = 0 otherwise. For convenience, in W N above we have included the self-interaction term Φ N (x, x), which does not affect the physics of the model. We also write W N and Φ N instead of W V and Φ V since these quantities depend only on the number of cells in V. The first term in Φ N with J 1 > 0 describes the attraction. By virtue of the Curie-Weiss approach, it is taken equal for all particles. The second term with J 2 > 0 describes the repulsion between two particles contained in one and the same cell. That is, in our model every two particles in V attract one another independently of their location, and repel if they are in the same cell. The intensities J 1 and J 2 are assumed to satisfy the following condition The latter is to secure the stability of the interaction [10], that is, to satisfy Let β = 1/k B T be the inverse temperature. To optimize the choice of the thermodynamic variables we introduce the following ones 4) and the dimensionless chemical potential µ = β × (physical chemical potential). Then, (p, µ) ∈ R + × R is considered to be the basic set of thermodynamic variables, whereas a and υ are model parameters.
The grand canonical partition function in region V is where |γ| stands for the number of points in the configuration γ, and Γ V is the subset of Γ 0 consisting of all γ contained in V. We write Ξ N instead of Ξ V for the reasons mentioned above.

Transforming the partition function
Now we use a concrete form of the energy as in (2.2) to bring (2.5) to a more convenient form. For a given = 1, . . . , N and a configuration γ ∈ Γ V , we set γ = γ ∩ ∆ , that is, γ is the part of the configuration contained in ∆ . Then, |γ | will stand for the number of points of γ contained in ∆ . Note that Then, cf. (2.2) and (2.5), x, y ∈γ

23502-3
To rewrite the integrand in (2.5) in a more convenient form we set where ∈ N N 0 is a vector with nonnegative integer components , = 1, 2, . . . , N. Then, (2.5) takes the form 8) where F N is as in (2.7) and ν(γ) ∈ N N 0 is the vector with component |γ |, = 1, . . . , N. For n, m ∈ N 0 , the Kronecker δ-symbol can be written Applying this in (2.8) we get Here, In getting the second line of (2.10), we use (2.6), and then the integral with λ is written according to (2.1). Note that the expression under the integral in the last line of (2.10) factors in j, which allows for writing it in the form Now we apply this in (2.9) and obtain where p is as in (2.4) and Note that, for p = 0, π turns into the (non-normalized) Poisson distribution with parameter υe µ . Hence, alternating the cell size amounts to shifting µ.

Single-phase domains
By a standard identity we transform (2.11) into the following expression and, cf. (2.4) and (2.12), Note that E is an infinitely differentiable function of all its arguments. Set By the following evident inequality we obtain from (2.15) and (2.14) that By virtue of Laplace's method [11], to calculate the large N limit in (2.16) we should find the global maxima of E(y, p, µ) as a function of y ∈ R.
Remark 2.1. From the estimate in (2.17) it follows that: (a) the integral in (2.13) is convergent for all p > 0 and µ ∈ R since a > 1, see (2.3) and (2.4); (b) for fixed p and µ, as the bounded from the above function E(y, p, µ) has global maxima, each of which is also its local maximum.
To get (b) we observe that (2.17) implies lim | y |→+∞ E(y, p, µ) = −∞; hence, each pointȳ of global maximum belongs to a certain interval (ȳ − ε,ȳ + ε), where it is also a maximum point. Since E is everywhere differentiable in y, thenȳ is the point of global maximum only if it solves the following equation By (2.14) and (2.15) this equation can be rewritten in the form Remark 2.2. As we will see from the proof of theorem 2.1 below, the equation in (2.19) has at least one solution for all p > 0 and µ ∈ R. Since both K 1 and K take only strictly positive values, these solutions are also strictly positive.
Definition 2.1. We say that (p, µ) belongs to a single-phase domain if E(y, p, µ) has a unique global maximumȳ ∈ R such that That is, not every point of global maximum corresponds to a point in a single-phase domain.
The condition in (2.19) determines the unique probability measure Q p,µ on N 0 such that which yields the probability law of the occupation number of a single cell. Then, the unique thermodynamic phase of the model corresponding to (p, µ) ∈ R is the product of the copies of the measure defined in (2.21). It is a probability measure on the space of all vectors n = (n ) ∞ =1 , in which n ∈ N 0 is the occupation number of -th cell. The role of the condition in (2.20) is to yield the possibility to apply Laplace's method for asymptotic calculation of the integral in (2.13). By direct calculations, it follows that In dealing with the equation in (2.19) we fix p > 0 and consider E 1 as a function of y ∈ R and µ ∈ R. Then, for a given µ 0 , we solve (2.19) to findȳ 0 and then check whether it is the unique point of global maximum and (2.20) is satisfied, i.e., whether (p, µ 0 ) belongs to a single-phase domain. Then, we slightly vary µ and repeat the same. This will yield a function µ →ȳ(µ) defined in the neighbourhood of µ 0 , which dependends on the choice of p and satisfiesȳ(µ 0 ) =ȳ 0 . In doing so, we use the analytic implicit function theorem based on the fact that, for each fixed p > 0, the function R 2 (y, µ) → E 1 (y, p, µ) can be analytically continued to some complex neighbourhood of R 2 , see (2.15) and (2.19). For the reader's convenience, we present this theorem here in the form adapted from [12, section 7.6, page 34]. For some p 0 > 0, let B ⊂ C 2 be a connected open set containing R 2 such that the function (y, µ) → E 1 (y, p 0 , µ) should be analytic in B.
(2.24) Remark 2.3. In the sequel, we also use the version of the implicit function theorem in which we do not employ the analytic continuation of E 1 to complex values of p. Let the conditions of proposition 2.1 regarding E 1 and E 2 be satisfied. Then, there exist open sets D i ⊂ R, i = 1, 2, 3, and a continuous functionȳ : D 3 × D 2 → D 1 such that p 0 ∈ D 3 , µ 0 ∈ D 2 and the following holds The partial derivative ofȳ(p, µ) over µ ∈ D 2 is given by the right-hand side of (2.24). In the sequel, by writingȳ(µ) we assume both the function as proposition 2.1, defined for a fixed p known from the context, and that as in remark 2.3 with the fixed value of p.
For a fixed p 0 > 0, assume that (p 0 , µ 0 ) belongs to a single-phase domain. By proposition 2.1 there exists ε > 0 such that the function (µ 0 − ε, µ 0 + ε) µ →ȳ(µ) can be defined by the equation E 1 (y, p 0 , µ) = 0. Its continuation from the mentioned interval is related to the fulfilment of the condition E 2 (ȳ(µ), p 0 , µ) < 0, cf. (2.20), which may not be the case. At the same time, by (2.23) we have that holding for all y ∈ R, p > 0 and µ ∈ R. In view of this and (2.24), it might be more convenient to use the inverse function y →μ(y) since the µ-derivative of E 1 is always nonzero. Its properties are described by the following statement obtained from the analytic implicit function theorem mentioned above. Recall that only positive y solves the equation in (2.19). Proof. For a single-phase domain, R, take (p 0 , µ 0 ) ∈ R. By remark 2.3 the function (p, µ) →ȳ(p, µ), defined by the equation E 1 (y, p, µ) = 0 is continuous in some open subset of R containing (p 0 , µ 0 ). By the continuity of E 2 (ȳ(p, µ), p, µ) and the fact that E 2 (ȳ(p 0 , µ 0 ), p 0 , µ 0 ) < 0 (since (p 0 , µ 0 ) ∈ R) we get that E 2 (ȳ(p, µ), p, µ) < 0 for (p, µ) in some open neighbourhood of (p 0 , µ 0 Note that, up to the factor υ −1 ,n(p, µ) is the particle density in phase Q p,µ . For a fixed p, by proposition 2.3 n(p, ·) in an increasing function on I p , which can be inverted to giveμ(p,n). By Laplace's method we then get the following corollary of proposition 2.3.

(2.27)
Let N p be the image of I p under the map µ →n(p, µ). Then, the inverse mapn →μ(p,n) is continuously differential and increasing on N p . By means of this map, for a fixed p, the pressure given in (2.27) can be written as a function ofn P =P(n) = υ −1 E(pn, p,μ(p,n)),n ∈ N p , (2.28) which is the equation of state.

The phase transition
Recall that the notion of the single-phase domain was introduced in definition 2.1, and each domain of this kind is an open subset of the open right half-plane {(p, µ) : p > 0, µ ∈ R}, see proposition 2.3. With this regard we have the following possibilities: (i) the whole half-plane {(p, µ) : p > 0, µ ∈ R} is such a domain; (ii) there exist more than one single-phase domain. In case (i), for all (p, µ) there exists one phase (2.22). In the context of this work, a phase transition is understood as the possibility of having different phases at the same value of the pair (p, µ). If this is the case, (p, µ) is called a phase coexistence point. Clearly, such a point should belong to the common topological boundary of at least two distinct single-phase domains. That is, to prove the existence of phase transitions we have to show that possibility (ii) takes place and that these single-phase domains have a common boundary. We do this in theorems 2.1 and 2.2 below.
Let R be a single-phase domain. Take (p 0 , µ 0 ) ∈ R and consider the line l p 0 = {(p 0 , µ) : µ ∈ R}. If the whole line lies in R, by proposition 2.3,ȳ(µ) is a continuously differentiable and increasing function of µ ∈ R. In our first theorem, we prove that this is the case for small enough p 0 . Proof. In view of remark 2.1, we have to show that, for fixed p p 0 and all µ ∈ R, E(y, p, µ) has exactly one local maximum such that (2.20) holds. For x ∈ R, we set, cf. (2.15), Similarly to (2.17), we get Note also that, for the functions defined in (2.29), we have By means of (2.29) we rewrite the equation in (2.19) in the following form Our aim is to show that there exists p 0 > 0 such that, for each p ∈ (0, p 0 ], the following holds: (a) the second line in (2.32) defines an increasing unbounded functionμ(x), x ∈ R; (b) pφ 2 (x, p) 1 − δ for some δ ∈ (0, 1) and all x ∈ R. Indeed, the function mentioned in (a) can be inverted to give an unbounded  increasing functionx(µ), µ ∈ R, such that the solution of (2.19) isȳ(µ) =x(µ) − µ. Then, by (2.29) and (2.23) we get where the latter inequality follows by (b). Thus, to prove both (a) and (b), it is enough to show that there exists positive p 0 such that for all x ∈ R and p ∈ (0, p 0 ]. (2.33) By (2.29) we have Since Φ(x, p) 1, we get from the latter Fix some x 0 > 0 and then set continuously depends on p > 0 and ψ(p) → 0 as p → 0. For each x > 0, one finds ξ ∈ (0, 1], dependent on x and p, such that From this and (2.30) we get (2.40)

23502-9
For the function defined in (2.39) and x 0 as in (2.37), we pick p 02 such that ψ(p) x 0 for all p p 02 .
For such values of p, this yields We apply this in (2.40) and get for all x x 0 and p p 02 . In (2.42),x(y) is the inverse of the function R x → y = pφ 1 (x, p). Note that dx(y) dy where y p = pφ 1 (x p , p). For the same p, let x p be such that aυp = exp(−x p − 2υe x p ), cf. (2.37). Then, by (2.36) and (2.44), we conclude that the inequality in (2.45) holds also for y ∈ (0, y p ], y p := pφ 1 (x p , p).
As we have seen in the proof of theorem 2.1, the mentioned two intervals may overlap, i.e., it may be that y p > y p , if p is small enough. Let us prove that this is not the case for big p. That is, let us show that there exists p 1 > 0 such that, cf. (2.44) and (2.41), for all p p 1 , there exists x ∈ R such that the following holds p 1 φ 2 (x, p 1 ) 1, and pφ 2 (x, p) > 1, for p > p 1 .
To this end, we estimate the denominator of (2.34) from the above and the numerator from the below. For x = ap/2, by (2.35) it follows that where we used the fact that p > 0. In the sum in the numerator of (2.34), we take just two summands corresponding to n 1 = 1, n 2 = 0 and n 1 = 0, n 2 = 1, and obtain, by (2.47), the following estimate

Numerical results
Here, we present the results of numerical calculations of the functions which appear in the preceding part of the paper.
We begin by considering the extremum points of the functions that appear in section 2.4. According to definition 2.1, the line l p = {(p, µ) : µ ∈ R} lies in a single-phase domain, if the function R + y → E(y, p, µ), see (2.14), has a unique non-degenerate global maximum for all µ ∈ R. The corresponding condition in (2.19) determines an increasing functionȳ(µ), see proposition 2.3, which can be inverted to giveμ(y), see (2.42). In theorem 2.1, we show that this holds for small enough p. Figure 1  For p = p c , the function y → E(y, p c , µ) still has a unique global maximum, which gets degenerate, i.e.,   Let us now turn to the maximum points of E(y, p, µ). For a as in (3.1) and p = 6, figure 3 presents the dependence of E on y for µ 1 = −2.3080 (curve a), and µ 2 = −1.4700 (curve b). This provides an illustration to (2.52). The curve plotted in figure 4 corresponds to the critical value of µ defined by the condition D(µ) = 0. That is, (p, µ c ) with p = 6 and µ c = −1.890291 is the phase coexistence point whose existence was proved in theorem 2.2. Figure 5 presents the dependence of E(ȳ 1 (µ), p, µ) (line 1, red) and E(ȳ 2 (µ), p, µ) (line 2, blue) on µ ∈ M p , p = 6, cf. (2.51). Their intersection occurs at µ = µ c = −1.890291.
The curves plotted in figure 6 present the isotherms -the dependence of the pressure onȳ, which is equivalent to the dependence on the densityn, see (2.26), calculated from (2.28) for a number of fixed values of p.

Concluding remarks
In this work, we proved the existence of multiple thermodynamic phases at the same values of the extensive model parameters -temperature and chemical potential. In contrast to the approach of [1], we deal directly with thermodynamic phases in the grand canonical setting. To the best of our knowledge, this is the first result of this kind.