Also at ]Centre of Excellence on Advanced Mechanics of Materials, Indian
Institute of Science.

Also at ]Centre of Excellence on Advanced Mechanics of Materials, Indian
Institute of Science.

# Geometric Thermodynamics of Collapse of Gels

Asif Raza<EMAIL_ADDRESS>[  Debasish Roy<EMAIL_ADDRESS>[
Computational Mechanics Laboratory, Department of Civil Engineering, Indian
Institute of Science, Bangalore  Sanhita Das<EMAIL_ADDRESS>Department
of Civil and Infrastructure Engineering, Indian Institute of Technology
Jodhpur

###### Abstract

Stimulus-induced volumetric phase transition in gels may be potentially
exploited for various bio-engineering and mechanical engineering applications.
Since the discovery of the phenomenon in the 1970s, extensive experimental
research has helped in understanding the phase transition and related critical
phenomena. Yet, little insight is available on the evolving microstructure. In
this article, we aim at unravelling certain geometric aspects of the
micromechanics underlying discontinuous phase transition in polyacrylamide
gels. Towards this, we use geometric thermodynamics and a Landau-Ginzburg type
free energy functional involving a squared gradient, in conjunction with
Flory-Huggins theory. We specifically exploit Ruppeiner’s approach of
Riemannian geometry-enriched thermodynamic fluctuation theory that has been
previously employed to investigate phase transitions in van der Waals fluids
and black holes. The framework equips us with a scalar curvature that relates
to the microstructural interactions of a gel during phase transition and at
critical points. This curvature also provides an insight into the universality
class of phase transition and the nature of polymer-polymer interactions.

††preprint: APS/123-QED

##  I Introduction

Stimulus-sensitive gels – a class of polymers, find extensive applications in
biomedical and mechanical engineering [1, 2]. Exposed to external stimuli such
as temperature, pH etc., they may undergo volumetric changes, which allows one
to control their volume through stimuli manipulations. Biomedical applications
where such volume control is crucial include drug delivery systems, tissue
engineering, implants etc. [3].

Gels may undergo both continuous and discontinuous phase transitions [4, 5].
The discontinuous phase transition, which is of primary interest, is similar
in many ways to the phase transition in van der Waals fluids or magnetic
systems [3, 6]. Both are volumetric phase transitions in systems with
disordered microstructures. The swollen and shrunken phases in gels correspond
respectively to the gaseous and liquid phases in van der Waals fluids.
Transition in gels may however be triggered not only by temperature, but
several other means such as composition of the solvent, pH and ionic changes,
irradiation by light and electric fields [7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 14, 18, 19, 20, 21]. These features make such materials attractive for
biomedical applications.

Understanding the phase behavior of gels, particularly the processes of phase
separation and spinodal decomposition, is crucial for optimizing their
performance in diverse applications. Phase equilibrium based on Flory-Huggins
theory offers the foundational context necessary for an insightful
understanding of phase transition [22]; but it fails to anticipate the
kinetics and the evolution of morphology with time in the separation process.
To address this, the notion of gradient energy is introduced through the
Landau-Ginzburg functional [23, 24]. The Landau-Ginzburg theory is often
applied to understand the behavior of order parameters near critical points
and to track the transition among different phases in materials [25, 26, 27,
28, 29, 30].

Phase transition and the associated critical phenomena are in need of a more
nuanced treatment affording an insight into the underlying molecular
interactions [31]. Experiments such as dynamic light scattering, friction
measurement, calorimetry [32, 4, 33, 10, 31, 34] have established the
existence of a critical point and led to an understanding of the various
response functions at the critical point. Several observations have come to
light, viz. increase and subsequent divergence of the intensity and
correlation of scattered light, lowering and subsequent divergence of the
swelling and collapsing relaxation times, and decrease of the bulk modulus to
zero [31, 32]. Singularities in the specific heat, osmotic compressibility
have been characterized using critical exponents [35] and the behaviour of the
gel at the critical point is observed to conform to a 3D Ising system.

Similar to other systems belonging to the same universality class, there are
attractive and repelling driving forces that cause the gel to vacillate
between collapsed and swollen phases. Based on the chemical composition of
gels, Tanaka surmised without sufficient evidence that the interactions could
be hydrophobic, electrostatic, ionic or through van der Waals forces [31].

An alternative approach to an accounting of microstructural interactions is
the thermodynamic fluctuation theory. The classical fluctuation theory however
fails where fluctuations are of the order of the system size [36]; thus it
cannot be directly exploited for phase transitions. Ruppeiner proposed a
geometric approach by introducing a Riemannian structure to the thermodynamic
manifold [36], which eliminates the lack of covariance associated with large
fluctuations. Once equipped with a Riemannian structure that naturally
enforces a certain class of frame indifference, the Ricci curvature of the
thermodynamic manifold may be used to gather information on microstructural
interactions.

Ruppeiner’s approach has been implemented to understand the microstucture of
van der Waals fluids, black holes, magnetic substances and elastomers
undergoing strain-induced crystallization [37, 38, 39, 40, 41, 42, 43, 44,
45]. The sign of the scalar curvature indicates the type of interaction among
the mesoscopic entities [44, 38]. This quantity has also been used in [44] to
establish a similarity between phase transitions in van der Waals fluids and
AdS black holes.

In this article, we analyze the volumetric phase transition in polyacrylamide
gels using Ruppeiner’s approach, the aim being to understand the underlying
interactions. We start with a brief description of the Riemannian geometry,
defining the relevant quantities in Sec. II, followed by an analysis in Sec.
III of phase transition using an approach that combines the classical mean
field theory with Landau-Ginzburg. We then determine the curvature for the
thermodynamic manifold of present interest in Sec. IV. Variations of curvature
with temperature and polymer volume fraction are discussed and some
predictions drawn on the microstructure and the nature of interactions during
phase transition. We also analyze the divergence of curvature in the vicinity
of the critical point to ascertain the universality class of phase transition.
Finally, we conclude the article in Sec. V with a summary of observations.

##  II Riemannian Thermodynamic Manifold

The thermodynamic fluctuation theory states that the probability ($P$) of
fluctuations of independent variables $(x_{1},x_{2},...x_{n})$ in a
thermodynamic system is proportional to the number of microstates that is in
turn related to the distance ($\Delta l$) among neighbouring fluctuation
states in the thermodynamic configuration space [36] as shown in Eq. (1).
Thus, the closer the two states, the higher the probability of fluctuation
between them.

$P(x_{1},x_{2},....x_{n})\propto\exp-\big{(}{\frac{1}{2}\Delta l^{2}}\big{)}$
(1)

$\Delta l$ has a positive definite quadratic form given by,

$\begin{split}(\Delta
l)^{2}=&-\dfrac{1}{k_{B}T}\left(\dfrac{\partial^{2}{\Psi}}{\partial{T}^{2}}\right)(\Delta
T)^{2}\\\
&+\dfrac{1}{k_{B}T}\sum_{i,j=2}^{n}\left(\dfrac{\partial^{2}{\Psi}}{\partial{x_{i}}\partial{x_{j}}}\right)\Delta
x_{i}\Delta x_{j}\end{split}$ (2)

where $\Psi$ is the free energy density of the thermodynamic system, $k_{B}$
the Boltzmann constant and $x_{1}=T$ the temperature of the system.

Let us consider a Riemannian manifold with metric $\mathbf{g}$ (written as a
matrix $[g]=g_{ij}$) and consisting of points representing the thermodynamic
states. The distance $\Delta l$ may then be written in terms of $\mathbf{g}$
as follows,

$\begin{split}\Delta l^{2}=\sum^{n}_{i,j=1}g_{ij}\Delta x_{i}\Delta
x_{j}\end{split}$ (3)

Combining Eqs. (2) and (3), $\mathbf{g}$ may be extracted in terms of $\Psi$
as follows.

$\mathbf{g}=\dfrac{1}{k_{B}T}\left[\begin{array}[]{ccccc}-\dfrac{\partial^{2}{\Psi}}{\partial{T}^{2}}&0&0&\cdots&0\\\
0&\dfrac{\partial^{2}{\Psi}}{\partial{x_{2}}^{2}}&\dfrac{\partial^{2}{\Psi}}{\partial{x_{2}}\partial{x_{3}}}&\cdots&\dfrac{\partial^{2}{\Psi}}{\partial{x_{2}}\partial{x_{n}}}\\\
0&\dfrac{\partial^{2}{\Psi}}{\partial{x_{3}}\partial{x_{2}}}&\dfrac{\partial^{2}{\Psi}}{\partial{x_{3}}^{2}}&\cdots&\dfrac{\partial^{2}{\Psi}}{\partial{x_{3}}\partial{x_{n}}}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&\dfrac{\partial^{2}{\Psi}}{\partial{x_{n}}\partial{x_{2}}}&\dfrac{\partial^{2}{\Psi}}{\partial{x_{n}}\partial{x_{3}}}&\cdots&\dfrac{\partial^{2}{\Psi}}{\partial{x_{n}}^{2}}\end{array}\right]$
(4)

We further define an affine connection with the Christoffel symbols
$\mathbf{\Gamma}_{ij}^{k}$ given by,

$\mathbf{\Gamma}_{ij}^{k}=g^{km}\frac{1}{2}\left(\frac{\partial{g_{im}}}{\partial{x^{j}}}+\frac{\partial{g_{jm}}}{\partial{x^{i}}}-\frac{\partial{g_{ij}}}{\partial{x^{m}}}\right)$
(5)

where $g^{ij}$ are components of the inverse of $[g]$.
$\mathbf{\Gamma}_{ij}^{k}$ may then be used to determine the components of the
fourth order Riemannian curvature tensor $\tilde{\mathbf{R}}$,

$\tilde{R}_{ijk}^{l}=\frac{\partial}{\partial
X^{i}}\mathbf{\Gamma}_{jk}^{l}-\frac{\partial}{\partial
X^{j}}\mathbf{\Gamma}_{ik}^{l}+\mathbf{\Gamma}_{im}^{l}\mathbf{\Gamma}_{jk}^{m}-\mathbf{\Gamma}_{jm}^{l}\mathbf{\Gamma}_{ik}^{m}$
(6)

$\tilde{\mathbf{R}}$ may be contracted to obtain the following expressions for
the second order symmetric Ricci tensor $\hat{\mathbf{R}}$ and the Ricci
scalar curvature $R$ respectively:

$\hat{R}_{ij}=\tilde{R}_{ikj}^{k}$ (7)

$R=g^{ij}\tilde{R}^{k}_{ikj}$ (8)

Recall that $g^{ij}$, being components of $[g]^{-1}$, constitute the
contravariant form of the metric tensor $\mathbf{g}$.

Specifically, for a thermodynamic system with two fluctuating variables,
$x_{1}$ and $x_{2}$, the Ricci scalar curvature $R$ may be determined as:

$R=\dfrac{1}{\sqrt{g}}\left[\dfrac{\partial}{\partial
x_{1}}\left(\dfrac{1}{\sqrt{g}}\dfrac{\partial g_{22}}{\partial
x_{1}}\right)+\dfrac{\partial}{\partial
x_{2}}\left(\dfrac{1}{\sqrt{g}}\dfrac{\partial g_{11}}{\partial
x_{2}}\right)\right]$ (9)

where $g\equiv g_{11}g_{22}$ is the determinant of $[g]$.

##  III Thermodynamics and phase diagrams of gels

###  III.1 Free energy

The total free energy density $\Psi$ of a gel has a form as in [4] and is
hypothesized to comprise of three parts as shown below.

$\Psi=\Psi_{m}+\Psi_{el}+\Psi_{surr}$ (10)

The first term $\Psi_{m}$ corresponds to the free energy of mixing while the
second term $\Psi_{el}$ denotes the elastic free energy due to the expansion
of the network structure. The last term $\Psi_{surr}$ is the thermal energy
exchanged with the reservoir.

A modified form of Flory’s formula [22] is used to express $\Psi_{m}$,

$\Psi_{m}=k_{B}T\left(\dfrac{\Theta}{2T}n_{1}v_{2}+n_{1}\ln{(1-v_{2})}\right)$
(11)

where $n_{1}$ is the number of solvent molecules per unit reference volume,
$v_{1}$ is the volume fraction of the solvent, $v_{2}$ is the volume fraction
of the polymer (solute) and the latter is related to $v_{1}$ through
$v_{2}=1-v_{1}$. Also, $\frac{\Theta}{2T}$ is a dimensionless quantity
characterizing the interaction between polymer and solvent molecules. Here,
$\Theta$ is the Flory temperature at which no interactions between polymer
chains and solvent molecules exist.

Further, if $\bar{v}$ denotes the volume of an individual solvent molecule,
then $v_{1}=\frac{n_{1}\bar{v}}{1+n_{1}\bar{v}}$ and
$v_{2}=\frac{1}{1+n_{1}\bar{v}}$. Thus $\Psi_{m}$ can be expressed in terms of
$n_{1}$ and $\bar{v}$ as follows,

$\Psi_{m}=k_{B}T\left[\dfrac{\Theta}{2T}\left(\dfrac{n_{1}}{1+n_{1}\bar{v}}\right)+n_{1}\ln{\left(\dfrac{n_{1}\bar{v}}{1+n_{1}\bar{v}}\right)}\right]$
(12)

We also assume that the gel is not completely devoid of the solvent in its
reference configuration. If the volume fraction of polymer in the reference
configuration is $v_{2_{0}}$ with the number of solvent molecules being
$n_{1_{0}}$, we can express the elastic free energy density $\Psi_{el}$ in
terms of the volumetric deformation $v_{2}/v_{2_{0}}$ as follows,

$\Psi_{el}=\dfrac{k_{B}T\nu_{e}}{2}\left[3\left(\dfrac{v_{2}}{v_{2_{0}}}\right)^{-2/3}-3+\ln{\left(\dfrac{v_{2}}{v_{2_{0}}}\right)}\right]$
(13)

where $\nu_{e}$ is the number of freely-moving chains between successive
entanglements or cross-links. Substituting the expressions for $v_{2}$ and
$v_{2_{0}}$ in Eq. (13), we obtain,

$\begin{split}\Psi_{el}=\dfrac{3}{2}k_{B}T\dfrac{S}{\bar{v}(1+n_{1_{0}}\bar{v})^{2}}\Bigg{[}\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{-2/3}-1\\\
+\ln{\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{1/3}}\Bigg{]}\end{split}$
(14)

where $S=\nu_{e}\bar{v}(1+n_{1_{0}}\bar{v})^{2}$.

To determine $\Psi_{surr}$, we first adopt the experimentally obtained
expression for the specific heat capacity $C_{v}$, reproduced below from [35].

$C_{v}(T)=A|t|^{-\alpha_{\pi}}\left[1+D|t|^{\Delta}\right]+B+C(T-T_{c})$ (15)

where $|t|=\frac{|T-T_{c}|}{T_{c}}$ and $T_{c}$ is the critical temperature.
Also, $\alpha_{\pi}$ is the specific heat critical exponent with a value of
$-0.05$ and $\Delta$ is the correction-to-scaling exponent with a value of
$0.5$. $A,B,C$ and $D$ are constants determined through curve fitting [35].
Their values are $-241.581,182,-1.193$ and $-1.115$ respectively. The
surrounding free energy $\Psi_{surr}$ is obtained as follows,

$\Psi_{surr}=\Psi_{c}+BT(1-\ln{T})-\dfrac{1}{2}CT(T+2T_{c}-2T_{c}\ln{T})$ (16)

so as not to violate $C_{v}=-T\dfrac{\partial^{2}\Psi}{\partial T^{2}}$.
Specifically, we have $C_{v}=-T\dfrac{\partial^{2}\Psi_{surr}}{\partial
T^{2}}$ as $\Psi=\Psi_{surr}+\Psi_{m}+\Psi_{el}$ and
$\dfrac{\partial^{2}\Psi_{m}}{\partial T^{2}}$ and
$\dfrac{\partial^{2}\Psi_{el}}{\partial T^{2}}$ are both identically $0$.

Characteristic to any discontinuous phase transition, $C_{v}(T)$ must diverge
at the critical point. From Eq. (15), it is noted that $C_{v}(T)$ has a point
of discontinuity at $T=T_{c}$. Consequently, the expression of $\Psi_{surr}$
obtained through an integration of $C_{v}$ is decomposed into two cases. The
first case is when $T<T_{c}$ and the second case arises when $T>T_{c}$.

Specifically, for $T<T_{c}$, $\Psi_{c}$ in Eq. (16) is given by,

$\begin{split}\Psi_{c}=&AT\times
10^{-4}\\{(29.98+8.42D)\tilde{T}^{5}+(5.69+18.31D)\tilde{T}^{4}\\\
&+(12.86+49.84D)\tilde{T}^{3}+(39.58+206.25D)\tilde{T}^{2}\\\
&+10^{4}[(0.025+0.275D)\tilde{T}+(1+D)-(1+D)\ln{T}]\\}\end{split}$ (17)

where $\tilde{T}=T/T_{c}$ is the reduced temperature. For $T>T_{c}$,
$\Psi_{c}$ assumes the form,

$\begin{split}\Psi_{c}=&AT\Bigg{\\{}\dfrac{26.83\times
10^{-4}}{\tilde{T}^{2.95}}+\dfrac{1.68\times
10^{-2}D}{\tilde{T}^{2.45}}+\dfrac{1.28\times 10^{-2}}{\tilde{T}^{1.95}}\\\
&+\dfrac{0.19D}{\tilde{T}^{1.45}}-\dfrac{1.05}{\tilde{T}^{0.95}}-\dfrac{2.22D}{\tilde{T}^{0.45}}-19.05\tilde{T}^{0.05}-1.17D\tilde{T}^{0.55}\Bigg{\\}}\end{split}$
(18)

Combining Eqs. (12), (13) and (16), we obtain the complete expression for
$\Psi$ as,

$\begin{split}\Psi(T,n_{1})=k_{B}T\left[\dfrac{\Theta}{2T}\left(\dfrac{n_{1}}{1+n_{1}\bar{v}}\right)+n_{1}\ln{\left(\dfrac{n_{1}\bar{v}}{1+n_{1}\bar{v}}\right)}\right]+\dfrac{3}{2}k_{B}T\dfrac{S}{\bar{v}(1+n_{1_{0}}\bar{v})^{2}}\left[\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{-2/3}-1-\ln{\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{-1/3}}\right]\\\
+\Psi_{c}+BT(1-\ln{T})-\dfrac{1}{2}CT(T+2T_{c}-2T_{c}\ln{T})\end{split}$ (19)

Figure 1: Plot of pressure ($\Pi$) vs. volume fraction ($v_{2}$) at various
temperatures. The black dot represents the critical point.

###  III.2 Pressure and phase diagrams

To characterize the phase transition and the critical phenomenon, we must
obtain the relevant phase diagram. The phase diagram is a graphical
representation of the variations of pressure, temperature, and volume
fraction. The relation among the three quantities is mathematically
represented by the equation of state for the osmotic pressure. The osmotic
pressure is the thermodynamic force conjugate to the volume fraction $v_{1}$
of the solvent related [4]. The equation of state relating osmotic pressure
$\Pi$ with the volume fraction $v_{2}$ and temperature $T$ is given by the
following constitutive relation,

$\Pi=-N_{A}\left(\dfrac{\partial\Psi}{\partial n_{1}}\right)$ (20)

where $\Psi$ is the total free energy. Using (19) in (20), we obtain the
equation of state as,

$\Pi=-RT\left\\{\dfrac{\Theta}{2T}\dfrac{1}{(1+n_{1}\bar{v})^{2}}+\dfrac{1}{1+n_{1}\bar{v}}+\ln{\left(\dfrac{n_{1}\bar{v}}{1+n_{1}\bar{v}}\right)}+\dfrac{S}{(1+n_{1_{0}}\bar{v})^{3}}\left[\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{1/3}-\dfrac{1}{2}\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)\right]\right\\}$
(21)

where $N_{A}$ is Avogadro’s number and $R$ is the gas constant. Eq. (21) may
also be written in terms of the volume fraction $v_{2}$,

$\Pi=-RT\left\\{\dfrac{\Theta}{2T}(v_{2})^{2}+v_{2}+\ln{\left(1-v_{2}\right)}+S(v_{2_{0}})^{3}\left[\left(\dfrac{v_{2}}{v_{2_{0}}}\right)^{1/3}-\dfrac{1}{2}\left(\dfrac{v_{2}}{v_{2_{0}}}\right)\right]\right\\}$
(22)

Fig. 1 shows the $\Pi-v_{2}$ projection of the phase diagram for a
polyacrylamide gel which is obtained using Eq. (22) (the material constants
for which are available in [5]): $k_{B}=1.38\times 10^{-23}$ J K-1,
$\Theta=400$ K, $S=600$, $v_{2_{0}}=0.01$, $R=8.3145$ J K-1 mol-1. The volume
of the individual solvent molecule, $\bar{v}=3.7\times 10^{-28}$ m3, and this
value is adopted from [46].

Fig. 1 shows the isotherms at temperatures $T=292\text{ K},294\text{
K},295\text{ K},296\text{ K},298\text{ K and }303\text{ K}$. In Fig. 1, we
observe all the salient features typical of a P-V diagram in a discontinuous
phase transition. Of all the isotherms, one is the critical isotherm on which
lies the critical point, a point of inflection which can be characterized by
the critical osmotic pressure $\Pi_{c}$, the critical volume fraction
$v_{2_{c}}$ and the critical temperature $T_{c}$. To determine $\Pi_{c}$,
$v_{2_{c}}$ and $T_{c}$, we simultaneously solve $\frac{\partial\Pi}{\partial
v_{2}}=0$ and $\frac{\partial^{2}\Pi}{\partial v_{2}^{2}}=0$. Due to the
equations being strongly nonlinear, they need to be solved numerically. For
polyacrylamide gels, we use the values of constants as mentioned earlier and
obtain $\Pi_{c},v_{2_{c}}$ and $T_{c}$ as 0.82 Pa, 0.1344 and 296 K
respectively, which is shown by the black dot in Fig. 1.

Above the temperature corresponding to the critical isotherm, i.e. for
$T>T_{c}$, distinct swollen or shrunken phases do not exist. Hence, no phase
transition takes place and a single stable homogeneous ”supercritical gel”
phase exists. Below the critical temperature $T<T_{c}$, there is always a
certain zone of instability, marked by negative compressibility during which
phase transition from a swollen to a shrunken phase takes place. Thus, the
critical point serves as the boundary line. Below it, there is a chance of
phase transition taking place; but above it, no phase transition occurs.

Resuming our analysis on the stability of different phases, all isotherms
below $T_{c}$ contain a zone of negative compressibility. The zone begins with
a local maximum and ends with a local minimum. The locus of these extrema is
the spinodal curve. It is important to note that negative compressibility is
thermodynamically infeasible and is a theoretical artifact that arises from
the form of the equation of state. In fact, this is also present in the P-V
diagram for van der Waals fluids [5, 47]. Physically no such state of negative
compressibility exists. The gel abruptly transitions from a swollen to a
shrunken phase at constant pressure. This unphysical outcome might be
attributed to the absence of gradient dependent terms in the expression for
free energy, thereby accounting for the heterogeneity occurring during phase
transition. Experimental evidence [48, 49, 50, 51] indicates that pressure is
constant in this region even though it is in non-equilibrium. Nonetheless, the
point is that, due to a thermodynamically non-equilibrium nature in this
region, the conventional definition of pressure is called into question and
requires further research.

###  III.3 Modified free energy and pressure

The gel in the zone of so-called negative compressibility undergoes
spontaneous phase separation by the process of spinodal decomposition [27, 52,
30]. This is inadequately represented in the free energy expression in Eq.
(19). To address this limitation, one way would be to use a gradient-enhanced
term in the Landau-Ginzburg functional [25, 29, 28, 23, 24]. This functional
leads to a generalized chemical potential that drives the diffusion in the
system. Specifically, we consider a free energy function that includes a
squared gradient term, which should be meaningful under the nonequilibrium
conditions underlying phase transition.

Let $\Psi_{grad}$, the free energy term that incorporates gradient dependence,
be given by,

$\Psi_{grad}=\frac{\eta
k_{B}T}{2}\left(\dfrac{1+n_{1}\bar{v}}{1+n_{1_{0}}\bar{v}}\right)(\nabla
n_{1})^{2}$ (23)

where, $\eta=\eta_{1}\tilde{T}^{2}+\eta_{2}\tilde{T}+\eta_{3}$ is the scaling
parameter of the gradient term. We also have $\eta_{1}=5.8484\times 10^{4}$,
$\eta_{2}=-12.3426\times 10^{4}$ and $\eta_{3}=6.4762\times 10^{4}$.

Now, the final form of the free energy ($\Psi_{mod}$) of the system after
incorporating the gradient-dependent term is given by,

$\Psi_{mod}=(1-\kappa)\Psi+\kappa\Psi_{grad}$ (24)

where $\kappa$ is the gradient energy parameter which plays a crucial role in
the region of phase-separation and weighs the contribution of $\Psi$ against
$\Psi_{grad}$ [30, 25, 29, 28, 23, 24]. It is evaluated as,

$\kappa=\\{\varphi(1-\varphi_{1})+\varphi_{1}\\}\varphi_{2}$ (25)

where $\varphi=(\gamma_{1}+\gamma_{2})$,
$\varphi_{1}=H\left(-\dfrac{\partial^{2}\Psi}{\partial n_{1}^{2}}\right)$,
$\varphi_{2}=H\left(-\dfrac{\partial\Psi}{\partial n_{1}}\right)$ with
$\gamma_{1}=\tanh{(\omega_{1})}$ and
$\gamma_{2}=0.5+0.5\tanh{100(\tanh{\omega_{2}}-0.1)}$. $\gamma_{1}$ and
$\gamma_{2}$ are the parameters corresponding to the local maximum and local
minimum of the spinodal region respectively. The constituent parameters are
given as,

$\omega_{1}=-2\times 10^{5}[0.5+0.5\tanh{(-\zeta)}]\vartheta\zeta_{1}$ (26)
$\omega_{2}=14[0.5+0.5\tanh{(\zeta)}]\vartheta\zeta_{2}$ (27)

where,
$\zeta=\tanh{\left(-\dfrac{1}{(90h)^{2}}\dfrac{\partial^{3}\Psi}{\partial
n_{1}^{3}}\right)}$,
$\vartheta=\tanh{\left(-\dfrac{1}{(90h)^{0.75}}\dfrac{\partial^{3}\Psi}{\partial
n_{1}^{3}}\right)}$,
$\zeta_{1}=\exp{\left\\{-\dfrac{1}{h^{\phi_{1}}}\left(\dfrac{\partial^{2}\Psi}{\partial
n_{1}^{2}}\right)^{2}\right\\}}$ and
$\zeta_{2}=\exp{\left\\{-\dfrac{1}{h^{\phi_{2}}}\left(\dfrac{\partial^{2}\Psi}{\partial
n_{1}^{2}}\right)^{2}\right\\}}$ with $h^{\phi_{1}}$ and $h^{\phi_{2}}$ being
the length scales associated with the local maximum and the local minimum of
the spinodal region respectively. Currently adopted expressions for $\phi_{1}$
and $\phi_{2}$ are given by
$\phi_{1}=\exp{\left\\{-841.6669+838.1726\tilde{T}\right\\}}+0.9955(\tilde{T})^{-0.0008}$,
$\phi_{2}=\exp{\left\\{-427.4063+424.0418\tilde{T}\right\\}}+0.9758(\tilde{T})^{-0.0043}$,
We also use $h=0.1\times 10^{-14}k_{B}^{4}$ and $\nabla n_{1}=10^{12}$.

Consequently, via Eq. (24), the modified pressure $(\Pi_{mod})$ may be
expressed as,

$\Pi_{mod}=-N_{A}\left(\dfrac{\partial\Psi_{mod}}{\partial n_{1}}\right)$ (28)

The above modification is consistent with experimental observations [48, 49,
50, 51] and the Maxwell tie-line construction [53] which gives the pressure
and the region where the transition from swollen to shrunken phase takes
place. In other words, both phases co-exist in this region. In line with [53],
a Maxwell tie-line is drawn subtending equal areas between the local maximum
and minimum and the tie-line. At any point on the tie-line, both swollen and
shrunken phases co-exist in a phase-separated microstructure and the relative
volume fraction of each phase may be determined through the lever-arm rule.
The locus of intersection of the tie-lines with the isotherms is known as the
coexistence curve. Between the coexistence curve and the spinodal curve, the
phases are metastable. At lower volume fractions, the metastable phase is
predominantly swollen (MSW) and at higher volume fractions, it is
predominantly shrunken (MSH). At volume fractions lower than any point on the
swollen limb of the coexistence curve, a homogeneously swollen (SW) phase
exists, while at volume fractions above the shrunken limb of the coexistence
curve, a homogeneously shrunken phase (SH) exists. The coexistence and
spinodal curves meet at the critical point.

It is important to characterize the critical point and determine $\Pi_{c}$,
$v_{2_{c}}$ and $T_{c}$ so that all three projections of the phase diagram may
be represented in terms of reduced dimensionless quantities. One direct
advantage of such a representation is that the critical exponents [4]
associated with the phase transition may be readily calculated. Another
significant advantage is the clear identification of the stable, unstable, and
metastable phase boundaries, i.e. the spinodal and coexistence curves.

The phase diagram projection $\Pi$ \- $v_{2}$ is redrawn using the reduced
pressure $\tilde{\Pi}=\Pi_{mod}/\Pi_{c}$, reduced volume fraction
$\tilde{v}_{2}=v_{2}/v_{2_{c}}$ and reduced temperature $\tilde{T}=T/T_{c}$
and shown in Fig. 2. The critical point is more distinctly illustrated in Fig.
2 and so are the stable and unstable zones. For isotherms $\tilde{T}<1$, the
part of the isotherm where it is parallel to the horizontal axis, is the
region where transition from swollen to shrunken phase takes place. In Fig. 2,
the coexistence curve is represented by the black dashed curve while the
spinodal curve is represented by the black dash-dotted curve.

Figure 2: $\tilde{\Pi}$-$\tilde{v_{2}}$ phase diagram denoting the various
phases and containing the various isotherms. The black dot denotes the
critical point and the black dashed curve represents the coexistence curve.
The spinodal curve is represented by the black dash-dotted curve. Figure 3:
$\tilde{T}-\tilde{v}_{2}$ phase diagram denoting the various phases. The red
and blue solid curves represent the swollen and shrunken limbs of the
coexistence curve respectively. The black dot denotes the critical point and
the black solid curve represents the isobar at $\tilde{\Pi}=0$. The dotted
horizontal line at $\tilde{T}=0.976$ corresponds to the phase transition
temperature and the black dash-dotted curve represents the spinodal curve.

Fig. 3 shows the $\tilde{T}-\tilde{v_{2}}$ projection of the phase diagram.
The solid black curve denotes the isobar $\Pi/\Pi_{c}=0$. Coexistence and
spinodal curves may be seen in $\tilde{\Pi}-\tilde{v}_{2}$ projection as well.
The swollen limb of the coexistence curve where predominately swollen phase
exists is represented by the solid red curve while the shrunken limb of the
coexistence curve is represented by the solid blue curve in Fig. 3.

Finally, the $\tilde{\Pi}-\tilde{T}$ projection is shown in Fig. 4. In this
figure, the entire coexistence or the phase transition region described in the
earlier diagrams (Fig. 2 and Fig. 3) appears as the solid red line. The
shrunken and swollen phases outside the phase transition region are also shown
in the figure. Also shown is the critical point and the supercritical gel
region. The spinodal curve has also been depicted by the black dash-dotted
curve in Fig. 3 and Fig. 4.

Figure 4: $\tilde{\Pi}-\tilde{T}$ phase diagram denoting the various phases.
The red solid curve represents the coexistence curve. The black dash-dotted
curve is the spinodal curve. The black dot represents the critical point.

For further analysis of phase transition mechanism, we formulate a method for
approaching equilibrium in a system where surface effects, as characterized by
a gradient energy given by Eq. (23), play a significant role. In line with the
work of [30], the continuity equation corresponding to the free energy of the
system as in Eq. (24) may be written as,

$\dfrac{\partial v_{2}}{\partial
t}=\nabla\cdot\\{Dv_{2}(1-v_{2})\nabla\mu_{mod}\\}$ (29)

where $D$ is the diffusion coefficient, $v_{2}$ is the volume fraction of
polymer which is the order parameter of the system and $\mu_{mod}$ is the
chemical potential associated with it. In nonequilibrium systems, the chemical
potential $\mu_{mod}$ varies across distinct points, serving as the driving
force for the process of diffusion. The chemical potential $\mu_{mod}$ in
Eq.(29) is calculated as,

$\begin{split}\mu_{mod}=&\Psi_{mod}+(1-v_{2})\dfrac{\delta\Psi_{mod}}{\delta
v_{2}}\\\ =&\Psi_{mod}+(1-v_{2})\left[\dfrac{\partial\Psi_{mod}}{\partial
v_{2}}-\nabla\dfrac{\partial\Psi_{mod}}{\partial(\nabla
v_{2})}\right]\end{split}$ (30)

where,

$\dfrac{\partial\Psi_{mod}}{\partial
v_{2}}=\left[(1-\kappa)\dfrac{\partial\Psi}{\partial
n_{1}}+\kappa\dfrac{\partial\Psi_{grad}}{\partial
n_{1}}+(\Psi_{grad}-\Psi)\dfrac{\partial\kappa}{\partial
n_{1}}\right]\dfrac{\partial n_{1}}{\partial v_{2}}$ (31)

$\dfrac{\partial\Psi_{mod}}{\partial(\nabla
v_{2})}=\kappa\dfrac{\partial\Psi_{grad}}{\partial(\nabla
n_{1})}\dfrac{\partial(\nabla n_{1})}{\partial(\nabla v_{2})}$ (32)

Eq. (29) is analogous to Cahn-Hilliard equation [23] with the difference that
the term $Dv_{2}(1-v_{2})$ substitutes for the mobility $M$. Eq. (29) can be
used to gain substantial insight into the dynamical aspects of phase change
using, say, numerical schemes such as finite difference or Euler’s method,
among others.

##  IV Thermodynamic geometry of gels

The phase diagrams afforded a certain measure of understanding of phase
transition in our system. Now, using Ruppeiner’s approach, we attempt to track
the evolution of the microstructure throughout this process. Indeed, the
nonlinear nature of the manifold on which the microstructural evolution
happens should rule out the emergence of unstable regions. Writing Eq. (2) in
terms of the fluctuating variables $T$ and $n_{1}$, we obtain the following
expression for the squared length of an incremental line element,

$(\Delta
l)^{2}=-\dfrac{1}{k_{B}T}\left(\dfrac{\partial^{2}{\Psi_{mod}}}{\partial{T}^{2}}\right)(\Delta
T)^{2}+\dfrac{1}{k_{B}T}\left(\dfrac{\partial^{2}{\Psi_{mod}}}{\partial{n_{1}}^{2}}\right)(\Delta
n_{1})^{2}$ (33)

The (matrix form of the) metric tensor in Eq. (4) is given by

$\mathbf{g}=\dfrac{1}{k_{B}T}\left[\begin{array}[]{cc}-\dfrac{\partial^{2}{\Psi_{mod}}}{\partial{T}^{2}}&0\\\
0&\dfrac{\partial^{2}{\Psi_{mod}}}{\partial{n_{1}}^{2}}\\\ \end{array}\right]$
(34)

In terms of $T$ and $n_{1}$, the Ricci scalar curvature $R$ may be determined
using Eq. (9),

$R=\dfrac{1}{\sqrt{g}}\left[\dfrac{\partial}{\partial
T}\left(\dfrac{1}{\sqrt{g}}\dfrac{\partial g_{22}}{\partial
T}\right)+\dfrac{\partial}{\partial
n_{1}}\left(\dfrac{1}{\sqrt{g}}\dfrac{\partial g_{11}}{\partial
n_{1}}\right)\right]$ (35)

where
$g=-\dfrac{1}{(k_{B}T)^{2}}\dfrac{\partial^{2}{\Psi_{mod}}}{\partial{T}^{2}}\dfrac{\partial^{2}{\Psi_{mod}}}{\partial{n_{1}}^{2}}$
is the determinant of the metric tensor.

|
---|---
|
|

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Plots of normalized scalar curvature $R_{N}$ with volume fraction
($\tilde{v_{2}}$) at temperatures (a) $\tilde{T}=0.986$, (b)
$\tilde{T}=0.993$, (c) $\tilde{T}=0.996$, (d) $\tilde{T}=1$, (e)
$\tilde{T}=1.006$, (f) $\tilde{T}=1.023$

(a) Swollen phase (SW)

(b) Shrunken phase (SH)

(c) Heterogeneous state with both swollen phase and shrunken phase (Unstable
SW + SH)

(d) Gel phase at critical point

Figure 6: Schematic representation of the microstructure of different phases
of the gel Figure 7: Plot showing the behavior of normalized scalar curvature
$R_{N}$ along the swollen limb (top red line) and shrunken limb (bottom blue
line) of the coexistence curve.

---

(a)

(b)

Figure 8: Plot of $\ln|R_{N}|$ vs $\ln|1-\tilde{T}|$. Black markers denote the
numerically obtained data and the solid line denotes the fitted curve. (a) The
coexistence swollen phase curve (red line)—the slope is -2.00902. (b) The
coexistence shrunken phase curve (blue line)—the slope is -1.99123.

Using the free energy density (Eq. (24)) in Eq. (35), the expression for $R$
is obtained as,

$R=\dfrac{p}{q}$ (36)

where,

$p=-k_{B}\Theta
n_{1}\bar{v}s\begin{bmatrix}C_{v}\tilde{T}(1+n_{1_{0}}\bar{v})(r+s)+C_{v}\left\\{\tilde{T}s+\tilde{T}\bar{v}n_{1}r+\tilde{T}r+\bar{v}n_{1}s\left(\tilde{T}-\dfrac{\Theta}{T_{c}}\right)\right\\}\\\
+\tilde{T}\left\\{\tilde{T}r+\tilde{T}\bar{v}n_{1}r+\tilde{T}s+\bar{v}n_{1}s\left(\tilde{T}-\dfrac{\Theta}{T_{c}}\right)\right\\}\left\\{CT_{c}+\dfrac{0.05A(\tilde{T}-1)}{|\tilde{T}-1|^{1.95}}+\dfrac{0.55AD(\tilde{T}-1)}{|\tilde{T}-1|^{1.45}}\right\\}\end{bmatrix}$
(37)
$q=2C_{v}^{2}\begin{bmatrix}\tilde{T}r+\tilde{T}\bar{v}n_{1}r+\tilde{T}s+\bar{v}n_{1}s\left(\tilde{T}-\dfrac{\Theta}{T_{c}}\right)\end{bmatrix}^{2}$
(38)

  * $r=S\bar{v}n_{1}\left\\{-2+3\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{2/3}\right\\}$,

  * $s=6(1+n_{1_{0}}\bar{v})^{2}\left(\dfrac{1+n_{1_{0}}\bar{v}}{1+n_{1}\bar{v}}\right)^{2/3}$.

where $S$ is the constant defined in (14).

We define the normalized scalar curvature $R_{N}=\dfrac{R}{k_{B}}$ in line
with [44] and observe its variation with reduced volume $\tilde{v}_{2}$ for
reduced temperatures $\tilde{T}=0.986,0.993,0.996,1,1.006$ and $1.023$; the
results are shown in Fig. 5. The curvature assumes small to moderate values
over a majority of the parametric space except in the vicinity of the spinodal
curve and the coexistence region, where it diverges. For $\tilde{T}=1$,
divergence of $R_{N}$ is observed only at the critical point $v_{2_{c}}=1$
(see Fig. 5 (d)). For $\tilde{T}>1$, divergence of $R_{N}$ is not observed
(see Figs. 5 (e) and (f)).

In the context of black holes, divergence of $R_{N}$ is speculated to have
some relation with divergence of the correlation length $\xi$ [44]. We believe
that the same holds true for volumetric phase transition and divergence in
curvature is directly related to that in the fluctuations of the order
parameter $\tilde{v}_{2}$. Thus, density fluctuations are divergent in the
coexistence region and at the critical point. Within the coexistence region,
density fluctuations manifest as phase segregation at the macroscopic scale,
attributed to spinodal decomposition within the spinodal or unstable region
and, nucleation and growth in the metastable region. At the critical point,
fluctuations cause segregation at the microscopic scale which manifests as
critical opalescence in light scattering experiments. Photon correlation
spectroscopic studies on polyacrylamide gels [4] have established the
existence of a phase segregated microstructure on the coexistence curve which
confirms our understanding. In Fig. 6, we have shown some schematic
representations of the microstructure in swollen, shrunken, critical and
unstable regimes based on our interpretation of curvature and aided by the
work in [4]. The gel appears homogeneous in the swollen and shrunken phases
and heterogeneous in the unstable and critical regimes. Based on what is known
of the Van der Waals fluid system [44], we may expect a homogeneous
microstructure in the supercritical regime but with large microscale density
fluctuations; see Figs. 5 (e) and (f). But owing to a lack of experimental
data pertaining to $C_{v}$ at higher temperatures, further insight at this
stage is not available.

Curvature has also been exploited to reveal the nature of interactions among
micro-constituents during different stages of phase transition in van der
Waals fluids and black holes. For example, the sign of curvature in the case
of van der Waals fluid [41] indicates whether the predominant force between
fluid molecules is attractive or repulsive for a particular temperature or
volume fraction. We examine the variation in the sign of curvature along the
coexistence curve on lines similar to [44]. The curvature remains negative for
all temperatures leading to the conclusion that no radical change occurs in
the nature of interaction with change in volume fraction $v_{2}$ over the
entire range of temperatures considered.

Finally, we would like to examine the critical behaviour of $R_{N}$ near the
critical point to further establish the universality class of phase
transition. We start by examining the critical exponent with which the
curvature diverges. To evaluate the exponent, we resort to a numerical
technique described in [45]. First we assume a critical-temperature dependent
functional form for $R_{N}$ in the vicinity of the critical point and fit the
curve to discrete points lying on the swollen and shrunken limbs of the
coexistence curve.

$R_{N}\sim-(1-\tilde{T})^{-\alpha}$ (39)

or, equivalently,

$\ln{|R_{N}|}=-\alpha\ln{|1-\tilde{T}|}+\beta$ (40)

near the critical point. Fig. 8 shows the discrete points and the fitted
curves. Along the coexistence swollen phase curve, we get the following
fitting result,

$\ln{|R_{N}|}=-2.00902\ln{|1-\tilde{T}|}-5.01464$ (41)

The result along the coexistence shrunken phase curve is,

$\ln{|R_{N}|}=-1.99123\ln{|1-\tilde{T}|}-4.7542$ (42)

Comparing Eqs. (41) and (42), we can deduce that the slope of $\ln{|R_{N}|}$
vs. $\ln{|1-\tilde{T}|}$ is -2. Thus, the critical exponent associated with
scalar curvature is 2. Also, $R_{N}$ as a function of $\tilde{T}$ is given by,

$R_{N}\sim-(1-\tilde{T})^{-2}$ (43)

Based on the mean field theory, the critical exponent associated with the
correlation length $\xi$ is $0.5$ [3]. So one may define the following
relationship between the curvature and the correlation length,

$R_{N}\sim-\xi^{4}$ (44)

The critical exponent of 2 and the scaling law of Eq. (44) establish that
phase transition in gels indeed belongs to the same universality class as van
der Waals fluids and black hole systems [44, 45, 54].

##  V Conclusion

In view of various practical applications, discontinuous volumetric phase
transition in gels is of interest to both scientific and industrial
communities. A substantive characterization of phase transition based on an
understanding of the microstructure is thus crucial. In the present article,
we have studied phase separation dynamics in the light of a Landau-Ginzburg
theory in order to gain an understanding of the transformative processes
underlying phase transition. We have used Ruppeiner’s thermodynamic approach
to probe the microstructure of a polyacrylamide gel undergoing phase
transition. The divergence points in the Ricci scalar curvature indicate
diverging correlation length and hence large density fluctuations. This
implies phase segregation and a heterogeneous microstructure during
transition. Further, the uniform negative sign of the curvature indicates no
drastic change in the molecular interactions for the range of temperatures
considered. Finally, we have numerically calculated the critical exponent
associated with the scalar curvature and established a scaling law between
curvature and correlation length. We have thus theoretically confirmed that
this phase transition belongs to the same universality class as gas-liquid
phase transition in van der Waals fluids. We may also draw inferences on the
gel microstructure across transitions and across critical points. While the
microstructure is phase-segregated during transition and at critical point, no
such conclusion can be reached in the supercritical regime with no
singularities in curvature post the critical point. In the supercritical
regime, from the present approach, we may merely assume that the nature of
interactions does not change and that the microstructure is homogeneous,
pending a more detailed study where curvature needs to be analyzed at even
higher temperatures with more accurate notions of specific heat in line with
[55]. In our future work, we would also like to exploit our observations to
design efficient supercritical systems for biomedical and structural
applications.

## References

  * Koetting  _et al._ [2015] M. C. Koetting, J. T. Peters, S. D. Steichen, and N. A. Peppas, Stimulus-responsive hydrogels: Theory, modern advances, and applications, Materials Science and Engineering: R: Reports 93, 1 (2015).
  * Kamata  _et al._ [2015] H. Kamata, X. Li, U.-i. Chung, and T. Sakai, Design of hydrogels for biomedical applications, Advanced healthcare materials 4, 2360 (2015).
  * Shibayama and Tanaka [1993] M. Shibayama and T. Tanaka, Volume phase transition and related phenomena of polymer gels, Responsive gels: volume transitions I , 1 (1993).
  * Tanaka  _et al._ [1977] T. Tanaka, S. Ishiwata, and C. Ishimoto, Critical behavior of density fluctuations in gels, Physical Review Letters 38, 771 (1977).
  * Tanaka [1978] T. Tanaka, Collapse of gels and the critical endpoint, Physical review letters 40, 820 (1978).
  * Dušek and Patterson [1968] K. Dušek and D. Patterson, Transition in swollen polymer networks induced by intramolecular condensation, Journal of Polymer Science Part A-2: Polymer Physics 6, 1209 (1968).
  * Amiya and Tanaka [1987] T. Amiya and T. Tanaka, Phase transitions in crosslinked gels of natural polymers, Macromolecules 20, 1162 (1987).
  * Katayama and Ohata [1985] S. Katayama and A. Ohata, Phase transition of a cationic gel, Macromolecules 18, 2781 (1985).
  * Hirokawa and Tanaka [1984] Y. Hirokawa and T. Tanaka, Volume phase transition in a non-ionic gel, in  _AIP Conference Proceedings_ , Vol. 107 (American Institute of Physics, 1984) pp. 203–208.
  * Hirotsu  _et al._ [1987] S. Hirotsu, Y. Hirokawa, and T. Tanaka, Volume-phase transitions of ionized n-isopropylacrylamide gels, The Journal of chemical physics 87, 1392 (1987).
  * Otake  _et al._ [1988] K. Otake, T. Tsuji, M. Konno, and S. Saito, Preparation of a new hydrogel and porous glass composite membrane, Journal of chemical engineering of Japan 21, 443 (1988).
  * Inomata  _et al._ [1990] H. Inomata, S. Goto, and S. Saito, Phase transition of n-substituted acrylamide gels, Macromolecules 23, 4887 (1990).
  * Otake  _et al._ [1990] K. Otake, H. Inomata, M. Konno, and S. Saito, Thermal analysis of the volume phase transition with n-isopropylacrylamide gels, Macromolecules 23, 283 (1990).
  * Siegel  _et al._ [1988] R. A. Siegel, M. Falamarzian, B. A. Firestone, and B. C. Moxley, ph-controlled release from hydrophobic/polyelectrolyte copolymer hydrogels, Journal of Controlled Release 8, 179 (1988).
  * Hirokawa  _et al._ [1984] Y. Hirokawa, T. Tanaka, and S. Katayama, Microbial adhesion and aggregation” ed. kc marshall (1984).
  * Ricka and Tanaka [1984] J. Ricka and T. Tanaka, Swelling of ionic gels: quantitative performance of the donnan theory, Macromolecules 17, 2916 (1984).
  * Ohmine and Tanaka [1982] I. Ohmine and T. Tanaka, Salt effects on the phase transition of ionic gels, The Journal of Chemical Physics 77, 5725 (1982).
  * Siegel and Firestone [1988] R. A. Siegel and B. A. Firestone, ph-dependent equilibrium swelling properties of hydrophobic polyelectrolyte copolymer gels, Macromolecules 21, 3254 (1988).
  * Osada [2005] Y. Osada, Conversion of chemical into mechanical energy by synthetic polymers (chemomechanical systems), in  _Polymer Physics_ (Springer, 2005) pp. 1–46.
  * Tanaka  _et al._ [1982] T. Tanaka, I. Nishio, S.-T. Sun, and S. Ueno-Nishio, Collapse of gels in an electric field, Science 218, 467 (1982).
  * Giannetti  _et al._ [1988] G. Giannetti, Y. Hirose, Y. Hirokawa, and T. Tanaka, Molecular electronic devices (1988).
  * Flory [1953] P. J. Flory,  _Principles of polymer chemistry_ (Cornell university press, 1953).
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system. i. interfacial free energy, The Journal of chemical physics 28, 258 (1958).
  * Cahn and Hilliard [1959] J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system. iii. nucleation in a two-component incompressible fluid, The Journal of chemical physics 31, 688 (1959).
  * Nauman and Balsara [1989] E. Nauman and N. P. Balsara, Phase equilibria and the landau—ginzburg functional, Fluid Phase Equilibria 45, 229 (1989).
  * Binder [1987a] K. Binder, Theory of first-order phase transitions, Reports on progress in physics 50, 783 (1987a).
  * Binder [1987b] K. Binder, Dynamics of phase separation and critical phenomena in polymer mixtures, Colloid and Polymer Science 265, 273 (1987b).
  * Ariyapadi and Nauman [1990] M. Ariyapadi and E. Nauman, Gradient energy parameters for polymer–polymer–solvent systems and their application to spinodal decomposition in true ternary systems, Journal of Polymer Science Part B: Polymer Physics 28, 2395 (1990).
  * He  _et al._ [1996] D. Q. He, S. Kwak, and E. B. Nauman, On phase equilibria, interfacial tension and phase growth in ternary polymer blends, Macromolecular theory and simulations 5, 801 (1996).
  * Nauman and He [2001] E. Nauman and D. Q. He, Nonlinear diffusion and phase separation, Chemical Engineering Science 56, 1999 (2001).
  * Li and Tanaka [1992] Y. Li and T. Tanaka, Phase transitions of gels, Annual Review of Materials Science 22, 243 (1992).
  * Hochberg  _et al._ [1979] A. Hochberg, T. Tanaka, and D. Nicoli, Spinodal line and critical point of an acrylamide gel, Physical Review Letters 43, 217 (1979).
  * Hirotsu [1987] S. Hirotsu, Phase transition of a polymer gel in pure and mixed solvent media, Journal of the Physical Society of Japan 56, 233 (1987).
  * Tokita and Tanaka [1991] M. Tokita and T. Tanaka, Friction coefficient of polymer networks of gels, The Journal of chemical physics 95, 4613 (1991).
  * Li and Tanaka [1989] Y. Li and T. Tanaka, Study of the universality class of the gel network system, The Journal of chemical physics 90, 5161 (1989).
  * Ruppeiner [1995] G. Ruppeiner, Riemannian geometry in thermodynamic fluctuation theory, Reviews of Modern Physics 67, 605 (1995).
  * Das  _et al._ [2022] S. Das, A. Raza, and D. Roy, Geometric thermodynamics of strain-induced crystallization in polymers, Phys. Rev. E 106, 015005 (2022).
  * Ruppeiner [2010] G. Ruppeiner, Thermodynamic curvature measures interactions, American Journal of Physics 78, 1170 (2010).
  * Ruppeiner and Bellucci [2015] G. Ruppeiner and S. Bellucci, Thermodynamic curvature for a two-parameter spin model with frustration, Physical Review E 91, 012116 (2015).
  * Ruppeiner [2008] G. Ruppeiner, Thermodynamic curvature and phase transitions in kerr-newman black holes, Physical Review D 78, 024016 (2008).
  * Ruppeiner and Seftas [2020] G. Ruppeiner and A. Seftas, Thermodynamic curvature of the binary van der waals fluid, Entropy 22, 1208 (2020).
  * Ruppeiner  _et al._ [2021] G. Ruppeiner, P. Mausbach, and H.-O. May, Thermodynamic geometry of the gaussian core model fluid, Fluid Phase Equilibria , 113033 (2021).
  * Kumara  _et al._ [2022] A. N. Kumara, C. L. A. Rizwan, K. Hegde, M. S. Ali, and K. M. Ajith, Microstructure of five-dimensional neutral gauss–bonnet black hole in anti-de sitter spacetime via $p-v$ criticality, General Relativity and Gravitation 55, 10.1007/s10714-022-03050-y (2022).
  * Wei  _et al._ [2019a] S.-W. Wei, Y.-X. Liu, and R. B. Mann, Repulsive interactions and universal properties of charged anti–de sitter black hole microstructures, Phys. Rev. Lett. 123, 071103 (2019a).
  * Wei  _et al._ [2019b] S.-W. Wei, Y.-X. Liu, and R. B. Mann, Ruppeiner geometry, phase transitions, and the microstructure of charged ads black holes, Phys. Rev. D 100, 124033 (2019b).
  * Das and Roy [2021] S. Das and D. Roy, A poroviscoelasticity model based on effective temperature for water and temperature driven phase transition in hydrogels, International Journal of Mechanical Sciences 196, 106290 (2021).
  * Tanaka [1979] T. Tanaka, Phase transitions in gels and a single polymer, Polymer 20, 1404 (1979).
  * Pallas and Pethica [1985] N. R. Pallas and B. A. Pethica, Liquid-expanded to liquid-condensed transition in lipid monolayers at the air/water interface, Langmuir 1, 509 (1985), https://doi.org/10.1021/la00064a019 .
  * Hifeda and Rayfield [1985] Y. Hifeda and G. Rayfield, Phase transitions in fatty acid monolayers containing a single double bond in the fatty acid tail, Journal of Colloid and Interface Science 104, 209 (1985).
  * Hifeda and Rayfield [1992] Y. Hifeda and G. Rayfield, Evidence for first-order phase transitions in lipid and fatty acid monolayers, Langmuir 8, 197 – 200 (1992).
  * Crane  _et al._ [1999] J. M. Crane, G. Putz, and S. B. Hall, Persistence of phase coexistence in disaturated phosphatidylcholine monolayers at high surface pressures, Biophysical Journal 77, 3134 (1999).
  * Tanaka [1996] H. Tanaka, Coarsening mechanisms of droplet spinodal decomposition in binary fluid mixtures, The Journal of chemical physics 105, 10099 (1996).
  * Clerk-Maxwell [1875] J. Clerk-Maxwell, On the dynamical evidence of the molecular constitution of bodies, Nature 11, 357 (1875).
  * Kumar and Sarkar [2014] P. Kumar and T. Sarkar, Geometric critical exponents in classical and quantum phase transitions, Phys. Rev. E 90, 042145 (2014).
  * Bolmatov  _et al._ [2013] D. Bolmatov, V. Brazhkin, and K. Trachenko, Thermodynamic behaviour of supercritical matter, Nature communications 4, 2331 (2013).


# Dark matter admixed neutron star as a possible compact component in the
GW190814 merger event

H. C. Das1,2<EMAIL_ADDRESS>Ankit Kumar1,2 S. K. Patra1,2 1Institute of
Physics, Sachivalaya Marg, Bhubaneswar 751005, India 2Homi Bhabha National
Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India

###### Abstract

We put constraints on the secondary component of GW190814 by analyzing the
observational data of the event. The relativistic mean-field models are used
to calculate the mass-radius profile and tidal deformability of the compact
object, considering it as a massive neutron star with the presence of dark
matter particles inside it. With the increase of dark matter percentage, the
maximum mass, radius, and tidal deformability of the neutron star decreases.
We observe that the predicted properties are well consistent with GW190814
observational data, suggesting the possibility of a dark matter admixed
neutron star if the underlying nuclear equation of state is sufficiently
stiff.

##  I Introduction

The LIGO and Virgo collaboration detected the most enigmatic compact binary
coalescence event (GW190814) involving a black hole and a compact object with
a mass range of 22.2–24.3 and 2.50–2.67 $M_{\odot}$ respectively Abbott  _et
al._ (2020). Since the electromagnetic counterpart has not been found and has
no measurable signature of tidal deformations from the gravitational wave
(GW), the secondary component might be the lightest black hole or the heaviest
neutron star (NS). Therefore, it is a big challenge to explain the nature of
such compact object. There is a lot of debate to understand the mystery of the
secondary component of the GW190814 event Most  _et al._ (2020); Vattis  _et
al._ (2020); Tews  _et al._ (2020); Zhang and Li (2020); Lim  _et al._ (2020);
Godzieba  _et al._ (2020); Huang  _et al._ (2020); Tan  _et al._ (2020);
Fattoyev  _et al._ (2020); Roupas (2021); Biswas  _et al._ (2021). The
analysis of GW190814 data implies the possible nature of such a compact object
as a NS only when (i) the equation of state (EOS) is very stiff Huang  _et
al._ (2020) or (ii) it is a rapidly rotating compact object below the mass
shedding frequency Zhang and Li (2020); Biswas  _et al._ (2021). Thus, one can
consider either of the two facts, which can produce the mass of the secondary
component of GW190814 around 2.50 $M_{\odot}$.

The binary NS merger event GW170817 Abbott  _et al._ (2017) has evoked to put
an upper bound on the maximum mass of nonrotating NS. By combining the total
binary mass of GW170817 inferred from GW signal with electromagnetic
observations, an upper limit of $M_{{\rm max}}\leq 2.17\ M_{\odot}$ is
predicted by Margalit et al. Margalit and Metzger (2017). Rezzolla et al. have
put the upper bound by combining the GW observations and quasiuniversal
relations as $M_{{\rm max}}\leq 2.16_{-0.15}^{+0.17}\ M_{\odot}$ Rezzolla  _et
al._ (2018). Further analysis by employing both energy and momentum
conservation laws and the numerical-relativity simulations show that the
maximum mass of a cold NS is bound to be less than $2.3\ M_{\odot}$ Shibata
_et al._ (2019). Also, various massive pulsar discoveries constrained the EOSs
of the supra-nuclear matter inside the core of the NS Demorest  _et al._
(2010); Antoniadis  _et al._ (2013); Cromartie  _et al._ (2020). These
observational data also suggest strong constraints on the maximum mass of the
slowly rotating NS with a lower bound of $\sim 2\ M_{\odot}$, which discarded
many EOSs. Recently, the Neutron Star Interior Composition Explorer (NICER)
data also imposed stringent constraints on the mass and radius of a canonical
star from the analysis of PSR J0030+0451 data Miller  _et al._ (2019); Riley
_et al._ (2019); Raaijmakers  _et al._ (2019). However, one cannot exclude the
existence of a super-massive NS as the secondary component of GW190814.

Many energy density functionals have been formulated to study the finite
nuclei, the nuclear matter, and the NS properties, out of which only a few
EOSs can reproduce the properties consistently with empirical/experimental
data Dutra  _et al._ (2012, 2014). To achieve the mass of the secondary
component of the GW190814 event, one has to include stiff EOS, which may or
may not satisfy different constraints such as flow Danielewicz  _et al._
(2002), GW170817 Abbott  _et al._ (2017, 2018), x-ray Antoniadis  _et al._
(2013); Cromartie  _et al._ (2020), and NICER data Miller  _et al._ (2019);
Riley  _et al._ (2019). For this, we search for stiff relativistic mean-field
(RMF) EOSs and find that the linear Walecka models (L-RMF) Serot and Walecka
(1986), the standard nonlinear $\sigma$-self-coupling interactions (NL-RMF)
Boguta and Bodmer (1977) and the density-dependent RMF (DD-RMF) Huang  _et
al._ (2020); Rather  _et al._ (2021) sets generally reproduce the maximum mass
of the NS more than 2.50 $M_{\odot}$. The BigApple force is designed to
reproduce the maximum mass of the NS more than 2.50 $M_{\odot}$ Fattoyev  _et
al._ (2020); Das  _et al._ (2020a). Although the L-RMF model predicts larger
mass, it fails to reproduce the properties of the finite nuclei and nuclear
matter. The recent work of Huang et al. indicates that the possibility of the
secondary object in GW190814 event as pure hadronic NS with DD-RMF parameter
sets Huang  _et al._ (2020). In DD-RMF family, the DD-LZ1 and DD-MEX forces
reproduce the NS masses as 2.554 and 2.556 $M_{\odot}$, respectively. Thus,
the DD-RMF sets could predict only the lower limit mass of the secondary
component in the GW190814 event. Hence, it ruled out the possibility of the
hybrid star or any admixture of exotic components such as hyperons, dark
matter, etc., inside the NS. Therefore, the L-RMF and DD-RMF forces are not
very useful for the investigation.

Figure 1: The EOS are shown for NL3 and NL-RA1 with DM Fermi momenta 0–0.06
GeV. The joint constraints from the gravitational wave data GW170817 and
GW190814 in shaded region are adopted from Ref. Abbott  _et al._ (2020).

In the present study, we consider the NL-RMF parameter sets as our primary
input to achieve the mass of the secondary component. Some of these forces are
NL3 Lalazissis  _et al._ (1997), NL3* Lalazissis  _et al._ (2009), NL1
Reinhard  _et al._ (1986), NL-SH Sharma  _et al._ (1993), NL3-II Lalazissis
_et al._ (1997) and NL-RA1 Rashdan (2001), which are generally considered to
be stiff EOSs. The maximum masses of all these parametrizations fall in the
range of 2.7–2.8 $M_{\odot}$, which are slightly more massive than the
secondary component of the GW190814 event. Therefore, we have two possible
scenarios to reduce the maximum mass and radius of the NS (i) with the
addition of hyperons or (ii) with the addition of DM. With the addition of
exotic particles, the EOS becomes softer, which reduces the mass, radius, and
tidal deformability of the NS Weissenborn  _et al._ (2012); Biswal  _et al._
(2019); Biswal (2019); Das  _et al._ (2019); Quddus  _et al._ (2020); Das  _et
al._ (2020b, 2021a). However, in the case of hyperons inside the NS, the
predicted canonical radius for NL-RMF type sets found to be $\sim 15$ km,
which does not lie in the range given by the NICER (11.96–14.26 km) Miller
_et al._ (2019) and GW190814 (12.2–13.7 km) Abbott  _et al._ (2020)
predictions. Therefore, in this case, we exclude the hyperons and involve the
dark matter as an extra candidate inside the NS.

##  II Model

Dark matter (DM) accounts for $\sim$85% of the composition in the Universe.
When a compact star rotates in the Galaxy through the DM halo, it captures
some of the DM due to its enormous gravitational potential and huge baryon
density Goldman and Nussinov (1989); Kouvaris and Tinyakov (2011). We take
nonannihilating weakly interacting massive particles (WIMPs) as a DM
candidate, because WIMPS are the thermal relics and most abundant DM particle.
Different models have already been incorporated to know its actual nature and
its effects on the NS properties Sandin and Ciarcelluti (2009); De Lavallaz
and Fairbairn (2010); Kouvaris and Tinyakov (2011); Ciarcelluti and Sandin
(2011); Leung  _et al._ (2011); Li  _et al._ (2012); Panotopoulos and Lopes
(2017); Ellis  _et al._ (2018); Das  _et al._ (2019); Bhat and Paul (2020);
Ivanytskyi  _et al._ (2019); Quddus  _et al._ (2020); Das  _et al._ (2020b,
2021b, 2021a). We modeled a Lagrangian density by assuming that the DM
particles interact with nucleons via standard model (SM) Higgs Panotopoulos
and Lopes (2017); Das  _et al._ (2019, 2020b, 2021b, 2021a).

Figure 2: The speed of sound for NL3 parameter set with DM Fermi momenta
0–0.06 GeV. The QCD conformal limit ($c_{s}^{2}=1/3$) and maximum speed limit
($c_{s}^{2}=1$) are shown with dashed magenta line. Figure 3: The $M-R$
relations for NL3 and NL-RA1 with different DM Fermi momenta. The horizontal
bars represent the maximum mass constraints from the PSR J0740+6620 (light
pink) and the GW190814 event (dark cyan). The NICER data are also shown with
two boxes from two different analyses Miller  _et al._ (2019); Riley  _et al._
(2019). The double-headed green line represents the radius constraints by the
GW190814 Abbott  _et al._ (2020) for 1.4 $M_{\odot}$ NS.

$\displaystyle{\cal{L}}$ $\displaystyle=$
$\displaystyle{\cal{L}}_{NM}+\bar{\chi}\left[i\gamma^{\mu}\partial_{\mu}-M_{\chi}+yh\right]\chi+\frac{1}{2}\partial_{\mu}h\partial^{\mu}h$
(1)
$\displaystyle-\frac{1}{2}M_{h}^{2}h^{2}+f\frac{M_{nucl.}}{v}\bar{\varphi}h\varphi,$

where ${\cal{L}}_{NM}$ represents the nuclear matter Lagrangian as given in
Refs. Das  _et al._ (2021b, a). $\varphi$ and $\chi$ are the nucleon and DM
wave function respectively. The mass of the DM ($M_{\chi}$) is 200 GeV taken
in the present calculation. The $y$ and $f$ are the DM-Higgs coupling
($y=0.07$) and nucleon-Higgs form factor ($f=0.35$) respectively. $M_{nucl.}$
is the mass of the nucleon as 939 MeV. The quantities such as $h$, $M_{h}$,
and $v$ are the Higgs field, mass of the Higgs ($M_{h}=125$ GeV) and Higgs
vacuum expectation value ($v=246$ GeV) respectively. We constrained the values
of $y$ and $f$ by calculating the scattering cross sections of the DM nucleons
which can be written as Cline  _et al._ (2013)

$\sigma_{SI}=\frac{y^{2}f^{2}M_{nucl.}^{2}}{4\pi}\frac{\mu_{r}}{v^{2}M_{h}^{2}},$
(2)

where $\mu_{r}$ is the reduced mass. The calculated cross sections for
$M_{\chi}=200$ GeV is found to be $9.70\times 10^{-46}$ cm2, which is well
consistent with XENON-1T Aprile  _et al._ (2016), PandaX-II Tan  _et al._
(2016) and LUX Akerib  _et al._ (2017) with 90% confidence level. The LHC had
also given a limit on the WIMP-nucleon scattering cross section in the range
from $10^{-40}$ to $10^{-50}$ cm2 Djouadi  _et al._ (2012). From the lattice
QCD Czarnecki  _et al._ (2010) and MILC collaboration results Toussaint and
Freeman (2009), the value of $f$ is constrained to be $f=0.33_{-0.07}^{+0.30}$
Aad  _et al._ (2015). The taken value $f=0.35$ in this calculation lies in the
region. Thus, we constrain the values of $y$ and $f$ from both direct
detection experiments and LHC searches.

Inside the NS, a neutron decays to proton, electron, and antineutrino. The
reverse process is also occurring to maintain both $\beta$ equilibrium and
charge neutrality conditions. Hence, the total energy density $\cal{E}$ and
pressure $P$ for the DM admixed NS involves the contribution from nucleons,
leptons, and DM, as given in Das  _et al._ (2020b, 2021b, 2021a).

$\displaystyle{\cal{E}}$ $\displaystyle=$
$\displaystyle{\cal{E}}_{NM}+{\cal{E}}_{l}+\frac{1}{\pi^{2}}\int_{0}^{k_{f}^{DM}}dk\
k^{2}\sqrt{k^{2}+M_{\chi}^{\star 2}}$ (3) $\displaystyle+$
$\displaystyle\frac{1}{2}M_{h}^{2}h_{0}^{2},$

and

$\displaystyle P$ $\displaystyle=$ $\displaystyle
P_{NM}+P_{l}+\frac{1}{3\pi^{2}}\int_{0}^{k_{f}^{DM}}\frac{dk\
k^{4}}{\sqrt{k^{2}+M_{\chi}^{\star 2}}}$ (4) $\displaystyle-$
$\displaystyle\frac{1}{2}M_{h}^{2}h_{0}^{2},$

where the ${\cal{E}}_{NM}$ (${\cal{E}}_{l}$) and $P_{NM}$ ($P_{l}$) are the
energy density and pressure for nuclear matter (lepton) Das  _et al._ (2020b,
2021b, 2021a). $k_{f}^{DM}$ and $M_{\chi}^{*}$ are the Fermi momentum and the
effective mass of the DM respectively.

We calculate the core EOS for DM admixed NS with NL-RMF family of parameter
sets. For the crust part, we use the Barcelona-Catania-Paris-Madrid (BCPM)
equation of state Sharma  _et al._ (2015) and form an unified EOS for the NS.
The outer crust is composed of the neutron-rich isotopes embedded in the
Coulomb lattice in the degenerate electron gas. As a result, the nuclei
capture electrons and form different isotopes in the region $26\leq Z\lesssim
40$. The degenerate electrons solely determine the pressure in this region.
Hence, the EOSs for this region are relatively well known Baym  _et al._
(1971); Sharma  _et al._ (2015). With a density of more than $2.6155\times
10^{-4}$ fm-3, the nuclei are highly neutron rich, and they do not bind the
additional neutrons. Therefore, the neutrons are dripped out from the nuclei.
The region defines the boundary between the outer and inner crust. As the
densities increase further, the nucleons form clusters in Coulomb crystals of
neutron-rich nuclei embedded by uniform electrons. The minimization of the
systems’ energy produces different exotic structures known as pasta Ravenhall
_et al._ (1983); Lorenz  _et al._ (1993). Thus, the detailed knowledge of the
EOSs in this region is essential for studying NS properties.

##  III Results and Discussions

In this section, we calculate the NS properties by using NL-RMF-type EOSs. As
it is mentioned, all these sets predict the maximum NS mass more than $2.67\
M_{\odot}$. In Fig. 1, the EOSs (such as NL3 and NL-RA1) are depicted with
$k_{f}^{DM}$ in the range of 0–0.06 GeV, which are shortly represented as
DM0–DM6. The EOSs become softer with the increase of $k_{f}^{DM}$. The
GW190814 discovery provides a joint constraint for NS-black hole (NSBH) merger
Abbott  _et al._ (2020). It is assumed to be a NS only when the maximum mass
is not less than the secondary component of GW190814. The EOSs that correspond
to DM3 and DM4 almost pass through the joint constraints imposed by GW170817
and GW190814. It indicates that this amount of DM (0.03–0.04 GeV) may be
available inside the NS. To test the validity of the EOS, we calculate the
speed of sound ($c_{s}^{2}=\partial P/\partial\mathcal{E}$) for DM admixed NS.
We find that our model respects the causality, which is shown in Fig. 2. The
value of $c_{s}^{2}$ increases up to 0.4 fm-3, and it becomes constant beyond
that.

Figure 4: (colour online) The tidal deformability of the NS with the mass for
NL3 and NL-RA1 parameter set for different DM Fermi momenta. The dark cyan
(green) line implies the constraint given by GW190814 (GW170817) data for 1.4
$M_{\odot}$ NS.

The macroscopic properties such as mass ($M$), radius ($R$), and dimensionless
tidal deformability ($\Lambda$) are also calculated with different DM
fractions. By solving Tolman-Oppenheimer-Volkoff equations Tolman (1939);
Oppenheimer and Volkoff (1939); one can find the mass-radius ($M$-$R$)
relation for each EOSs by assuming the $P(0)=P_{C}$ and $P(R)=0$ as the
boundary conditions. We find that the NS’s mass and radius reduces with DM’s
addition, as shown in Fig. 3. The maximum mass of the secondary object and the
canonical radius constraints are 2.50–2.67 $M_{\odot}$ and 12.20–13.70 km as
given in GW190814 Abbott  _et al._ (2020). The $M_{{\rm max}}$ and $R_{1.4}$
range correspond to DM3 are found to be 2.65–2.67 $M_{\odot}$ and 12.38–12.75
km respectively for the considered parameter sets (see Table 1), which is
consistent with GW190814 data. Under the assumption of a stiff nuclear
equation of state, the value of maximum mass with assumed parameter set hints
towards the possibility of the secondary component of the GW190814 event to be
a dark matter admixed neutron star. The NICER Miller  _et al._ (2019); Riley
_et al._ (2019) constraints are also respected by DM3 and DM4. Thus we
conclude that given a specific model of the nuclear EOS, one can fix the
fraction of DM inside the NS from the observational data.

When the NS presents in the tidal field created by its companion star, then it
is deformed. The relation between quadruple deformation $Q_{ij}$ and tidal
field $\epsilon_{ij}$ is written as Hinderer (2008); Kumar  _et al._ (2017)

$Q_{ij}=-\lambda\epsilon_{ij},$ (5)

where the $\lambda$ is defined as tidal deformability. The dimensionless tidal
deformability can be defined as $\Lambda=\lambda/M^{5}$. We calculate the
$\Lambda$ for different DM Fermi momenta as shown in Fig. 4. The values of
$\Lambda_{1.4}$ for DM3 are 842.12, 784.59, 932.99, 804.36 and 818.48 for NL3,
NL3*, NL-SH, NL3-II and NL-RA1 respectively. The predictions of
$\Lambda_{1.4}$ from GW170817 is $190_{-120}^{+390}$ at 90% confidence level
Abbott  _et al._ (2018). The GW190814 also put a constraint on $\Lambda_{1.4}$
under NSBH scenario, $\Lambda_{1.4}=616_{-158}^{+273}$ Abbott  _et al._
(2020). The calculated values of $\Lambda_{1.4}$ for DM3 are well matched with
GW190814 data, and, for DM4 case, it is well consistent with GW170817 data.
This also gives sufficient hints to constraint the amount of DM inside the NS.

##  IV Summary

The mass gap of 2.5–5.0 $M_{\odot}$ in the secondary component of GW190814
motivates us to explore its true nature, and we find the possibility that the
compact object could be a DM admixed NS. To know the mystery of DM inside the
NS, we modeled a Lagrangian density by assuming that the DM particles interact
with nucleons via SM Higgs. To calculate neutron star EOS, we adopt the RMF
model with NL-RMF parameter sets. The DM admixed NS properties are computed
with different DM fractions inside it.

To achieve the mass of the secondary object, one should take stiff EOSs. In
RMF parametrizations, only NL-RMF and DD-RMF sets are stiff in character,
which predict the maximum masses more than 2.50 $M_{\odot}$. Huang et al.
Huang  _et al._ (2020) also suggested that the possibility of the secondary
object as a NS (with DD-RMF sets without any admixture of DM or strange
particles). However, the NL-RMF forces with neutron star EOSs predict masses
more than 2.67 $M_{\odot}$ are candidates for DM admixed NS. This scenario is
only seen for NL-RMF type sets because no DD-RMF forces predict the mass more
than 2.67 $M_{\odot}$. Hence, we assume that those stars that correspond to
DD-RMF sets are only NS without DM inside it.

Table 1: The quantities such as $M_{{\rm max}}$, $R_{{\rm max}}$, $R_{1.4}$ and $\Lambda_{1.4}$ are given with $k_{f}^{DM}=$ 0, 0.03 and 0.04 GeV for NL3, NL3*, NL1, NL-SH, NL3-II and NL-RA1 respectively. $k_{f}^{DM}$ (GeV) |  NL3 |  NL3* |  NL1 |  NL-SH |  NL3-II |  NL-RA1
---|---|---|---|---|---|---
0.00 | 0.03 | 0.04 | 0.00 | 0.03 | 0.04 | 0.00 | 0.03 | 0.04 | 0.00 | 0.03 | 0.04 | 0.00 | 0.03 | 0.04 | 0.00 | 0.03 | 0.04
| $M_{{\rm max}}$
---
$(M_{\odot})$
2.78 | 2.65 | 2.50 | 2.76 | 2.64 | 2.49 | 2.84 | 2.71 | 2.55 | 2.79 | 2.66 | 2.51 | 2.77 | 2.65 | 2.50 | 2.78 | 2.66 | 2.51
| $R_{{\rm max}}$
---
(km)
13.40 | 12.46 | 11.53 | 13.10 | 12.38 | 11.60 | 13.63 | 12.75 | 11.67 | 13.53 | 12.63 | 11.60 | 13.14 | 12.42 | 11.52 | 13.19 | 12.47 | 11.57
| $R_{1.4}$
---
(km)
14.08 | 13.16 | 11.82 | 14.03 | 13.13 | 12.13 | 14.73 | 13.46 | 11.90 | 14.35 | 13.41 | 11.90 | 14.05 | 13.16 | 11.81 | 14.11 | 13.22 | 11.85
$\Lambda_{1.4}$ | 1311.16 | 842.12 | 499.04 | 1250.45 | 784.59 | 492.40 | 1503.77 | 912.65 | 548.83 | 1528.42 | 932.99 | 563.90 | 1282.90 | 804.36 | 503.89 | 1338.98 | 818.48 | 516.96

With the addition of DM, the EOS becomes softer and $M_{{\rm max}}$, $R_{{\rm
max}}$, $R_{1.4}$ and $\Lambda_{1.4}$ reduce with increase of DM percentage as
given in Table 1. The EOS corresponds to DM3 and DM4 well passes through the
joint constraint inferred from GW170817+GW190814 by assuming that the
secondary component is a neutron star only when its maximum mass is not less
than the 2.50 $M_{\odot}$. The causality is also respected by this model for
all values of DM Fermi momenta. Other properties such as $M_{{\rm max}}$ and
$R_{1.4}$ are well consistent with the GW190814 data for the DM3 case. The
values of $R_{1.4}$ for DM3 and DM4 are well constrained with the NICER
results in almost all assumed sets. The calculated $\Lambda_{1.4}$ for DM3 and
DM4 lie in the range given by both GW190814 and GW170817 data. Thus, one can
constrain the DM percentage inside the NS from these observational data.

Therefore, the possibility of the secondary object in GW190814 as a DM admixed
NS only when the maximum mass should be more than 2.67 $M_{\odot}$. We hope
future GW and x-ray observations may answer the presence of DM inside the NS
and constrain its properties. Thus, we suggest that the LIGO/Virgo have to
include DM inside the compact objects when they inferred some properties such
as mass, radius, and tidal deformability of the binary neutron star.

## References

  * Abbott  _et al._ (2020) R. Abbott, T. D. Abbott, S. Abraham, _et al._ , The Astrophysical Journal 896, L44 (2020).
  * Most  _et al._ (2020) E. R. Most, L. J. Papenfort, L. R. Weih, and L. Rezzolla, MNRAS 499, L82 (2020).
  * Vattis  _et al._ (2020) K. Vattis, I. S. Goldstein, and S. M. Koushiappas, Phys. Rev. D 102, 061301 (2020).
  * Tews  _et al._ (2020) I. Tews, P. T. H. Pang, T. Dietrich, _et al._ , (2020), arXiv:2007.06057 [astro-ph.HE] .
  * Zhang and Li (2020) N.-B. Zhang and B.-A. Li, The Astrophysical Journal 902, 38 (2020).
  * Lim  _et al._ (2020) Y. Lim, A. Bhattacharya, J. W. Holt, and D. Pati, (2020), arXiv:2007.06526 [nucl-th] .
  * Godzieba  _et al._ (2020) D. A. Godzieba, D. Radice, and S. Bernuzzi, (2020), arXiv:2007.10999 [astro-ph.HE] .
  * Huang  _et al._ (2020) K. Huang, J. Hu, Y. Zhang, and H. Shen, The Astrophysical Journal 904, 39 (2020).
  * Tan  _et al._ (2020) H. Tan, J. Noronha-Hostler, and N. Yunes, Phys. Rev. Lett. 125, 261104 (2020).
  * Fattoyev  _et al._ (2020) F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and B. Reed, Phys. Rev. C 102, 065805 (2020).
  * Roupas (2021) Z. Roupas, Astrophysics and Space Science 366 (2021), 10.1007/s10509-021-03919-5.
  * Biswas  _et al._ (2021) B. Biswas, R. Nandi, P. Char, S. Bose, and N. Stergioulas, Monthly Notices of the Royal Astronomical Society 505, 1600 (2021).
  * Abbott  _et al._ (2017) B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (The LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * Margalit and Metzger (2017) B. Margalit and B. D. Metzger, The Astrophysical Journal 850, L19 (2017).
  * Rezzolla  _et al._ (2018) L. Rezzolla, E. R. Most, and L. R. Weih, The Astrophysical Journal 852, L25 (2018).
  * Shibata  _et al._ (2019) M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi, Phys. Rev. D 100, 023015 (2019).
  * Demorest  _et al._ (2010) P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Nature 467, 1081–1083 (2010).
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. C. Freire, _et al._ ,  340 (2013), 10.1126/science.1233232.
  * Cromartie  _et al._ (2020) H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. Demorest, _et al._ , Nature Astronomy 4, 72–76 (2020).
  * Miller  _et al._ (2019) M. C. Miller, F. K. Lamb, A. J. Dittmann, _et al._ , APJ 887, L24 (2019).
  * Riley  _et al._ (2019) T. E. Riley, A. L. Watts, S. Bogdanov, _et al._ , APJ 887, L21 (2019).
  * Raaijmakers  _et al._ (2019) G. Raaijmakers, T. E. Riley, A. L. Watts, _et al._ , Astrophys. J 887, L22 (2019).
  * Dutra  _et al._ (2012) M. Dutra, O. Lourenço, J. S. Sá Martins, _et al._ , Phys. Rev. C 85, 035201 (2012).
  * Dutra  _et al._ (2014) M. Dutra, O. Lourenço, S. S. Avancini, _et al._ , Phys. Rev. C 90, 055203 (2014).
  * Danielewicz  _et al._ (2002) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002).
  * Abbott  _et al._ (2018) B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018).
  * Serot and Walecka (1986) B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16, 1 (1986).
  * Boguta and Bodmer (1977) J. Boguta and A. Bodmer, Nucl. Phys. A 292, 413 (1977).
  * Rather  _et al._ (2021) I. A. Rather, U. Rahaman, M. Imran, _et al._ , (2021), arXiv:2102.04067 [nucl-th] .
  * Das  _et al._ (2020a) H. C. Das, A. Kumar, B. Kumar, S. K. Biswal, and S. K. Patra, (2020a), arXiv:2009.10690 [nucl-th] .
  * Lalazissis  _et al._ (1997) G. A. Lalazissis, J. König, and P. Ring, Phys. Rev. C 55, 540 (1997).
  * Lalazissis  _et al._ (2009) G. Lalazissis, S. Karatzikos, R. Fossion, D. P. Arteaga, A. Afanasjev, and P. Ring, Physics Letters B 671, 36 (2009).
  * Reinhard  _et al._ (1986) P. G. Reinhard, M. Rufa, J. Maruhn, W. Greiner, and J. Friedrich, Z. Phys. A Atomic Nuclei 323, 13 (1986).
  * Sharma  _et al._ (1993) M. Sharma, M. Nagarajan, and P. Ring, Physics Letters B 312, 377 (1993).
  * Rashdan (2001) M. Rashdan, Phys. Rev. C 63, 044303 (2001).
  * Weissenborn  _et al._ (2012) S. Weissenborn, D. Chatterjee, and J. Schaffner-Bielich, Phys. Rev. C 85, 065802 (2012).
  * Biswal  _et al._ (2019) S. K. Biswal, S. K. Patra, and S.-G. Zhou, APJ 885, 25 (2019).
  * Biswal (2019) S. K. Biswal, AIP Conf. Proceedings 2127, 020031 (2019).
  * Das  _et al._ (2019) A. Das, T. Malik, and A. C. Nayak, Phys. Rev. D 99, 043016 (2019).
  * Quddus  _et al._ (2020) A. Quddus, G. Panotopoulos, B. Kumar, S. Ahmad, and S. K. Patra, Journal of Physics G: Nuclear and Particle Physics 47, 095202 (2020).
  * Das  _et al._ (2020b) H. C. Das, A. Kumar, B. Kumar, _et al._ , MNRAS 495, 4893 (2020b).
  * Das  _et al._ (2021a) H. C. Das, A. Kumar, and S. K. Patra, (2021a), arXiv:2104.01815 [astro-ph.HE] .
  * Goldman and Nussinov (1989) I. Goldman and S. Nussinov, Phys. Rev. D 40, 3221 (1989).
  * Kouvaris and Tinyakov (2011) C. Kouvaris and P. Tinyakov, Phys. Rev. D 83, 083512 (2011).
  * Sandin and Ciarcelluti (2009) F. Sandin and P. Ciarcelluti, Astroparticle Phys. 32, 278–284 (2009).
  * De Lavallaz and Fairbairn (2010) A. De Lavallaz and M. Fairbairn, Phys. Rev. D 81, 123521 (2010).
  * Ciarcelluti and Sandin (2011) P. Ciarcelluti and F. Sandin, Phys. Lett. B 695, 19–21 (2011).
  * Leung  _et al._ (2011) S.-C. Leung, M.-C. Chu, and L.-M. Lin, Phys. Rev. D 84, 107301 (2011).
  * Li  _et al._ (2012) A. Li, F. Huang, and R.-X. Xu, Astroparticle Phys. 37, 70–74 (2012).
  * Panotopoulos and Lopes (2017) G. Panotopoulos and I. Lopes, Phys. Rev. D 96, 083004 (2017).
  * Ellis  _et al._ (2018) J. Ellis, G. Hütsi, K. Kannike, L. Marzola, M. Raidal, and V. Vaskonen, Phys. Rev. D 97, 123007 (2018).
  * Bhat and Paul (2020) S. A. Bhat and A. Paul, The European Physical Journal C 80, 544 (2020).
  * Ivanytskyi  _et al._ (2019) O. Ivanytskyi, V. Sagun, and I. Lopes, (2019), arXiv:1910.09925 [astro-ph.HE] .
  * Das  _et al._ (2021b) H. C. Das, A. Kumar, B. Kumar, S. K. Biswal, and S. K. Patra, JCAP 2021, 007 (2021b).
  * Cline  _et al._ (2013) J. M. Cline, P. Scott, K. Kainulainen, and C. Weniger, Phys. Rev. D 88, 055025 (2013).
  * Aprile  _et al._ (2016) E. Aprile, J. Aalbers, F. Agostini, _et al._ , JCAP 2016, 027–027 (2016).
  * Tan  _et al._ (2016) A. Tan, M. Xiao, X. Cui, _et al._ (PandaX-II Collaboration), Phys. Rev. Lett. 117, 121303 (2016).
  * Akerib  _et al._ (2017) D. S. Akerib, S. Alsum, H. M. Araújo, X. Bai, _et al._ (LUX Collaboration), Phys. Rev. Lett. 118, 021303 (2017).
  * Djouadi  _et al._ (2012) A. Djouadi, O. Lebedev, Y. Mambrini, and J. Quevillon, Physics Letters B 709, 65 (2012).
  * Czarnecki  _et al._ (2010) A. Czarnecki, J. G. Körner, and J. H. Piclum, Phys. Rev. D 81, 111503 (2010).
  * Toussaint and Freeman (2009) D. Toussaint and W. Freeman (MILC Collaboration), Phys. Rev. Lett. 103, 122002 (2009).
  * Aad  _et al._ (2015) G. Aad, B. Abbott, J. Abdallah, _et al._ , Journal of High Energy Physics 2015 (2015), 10.1007/jhep11(2015)206.
  * Sharma  _et al._ (2015) B. K. Sharma, M. Centelles, X. Viñas, M. Baldo, and G. F. Burgio, A&A 584, A103 (2015).
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, apj 170, 299 (1971).
  * Ravenhall  _et al._ (1983) D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys. Rev. Lett. 50, 2066 (1983).
  * Lorenz  _et al._ (1993) C. P. Lorenz, D. G. Ravenhall, and C. J. Pethick, Phys. Rev. Lett. 70, 379 (1993).
  * Tolman (1939) R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * Hinderer (2008) T. Hinderer, The Astrophysical Journal 677, 1216 (2008).
  * Kumar  _et al._ (2017) B. Kumar, S. K. Biswal, and S. K. Patra, Phys. Rev. C 95, 015801 (2017).


# Paucity phenomena for polynomial products

Victor Y. Wang  Courant Institute of Mathematical Sciences, New York
University, New York, NY, USA<EMAIL_ADDRESS>and  Max Wenqiang Xu
Department of Mathematics, Stanford University, Stanford, CA, USA
<EMAIL_ADDRESS>
###### Abstract.

Let $P(x)\in\mathbb{Z}[x]$ be a polynomial with at least two distinct complex
roots. We prove that the number of solutions
$(x_{1},\dots,x_{k},y_{1},\dots,y_{k})\in[N]^{2k}$ to the equation

$\prod_{1\leq i\leq k}P(x_{i})=\prod_{1\leq j\leq k}P(y_{j})\neq 0$

(for any $k\geq 1$) is asymptotically $k!N^{k}$ as $N\to+\infty$. This solves
a question first proposed and studied by Najnudel. The result can also be
interpreted as saying that all even moments of random partial sums
$\frac{1}{\sqrt{N}}\sum_{n\leq N}f(P(n))$ match standard complex Gaussian
moments as $N\to+\infty$, where $f$ is the Steinhaus random multiplicative
function.

###### Key words and phrases:

Multiplicative Diophantine equations, paucity, random multiplicative
functions, Gaussian behavior

###### 1991 Mathematics Subject Classification:

Primary 11K65; Secondary 11D45, 11D57, 11D79, 11N37

##  1\. Introduction

We are interested in counting solutions to Diophantine equations involving
products of polynomial values. Let $P(x)\in\mathbb{Z}[x]$ be a polynomial.
Consider the equation

$\prod_{1\leq i\leq k}P(x_{i})=\prod_{1\leq j\leq k}P(y_{j})\neq 0$ (1.1)

for a given integer $k\geq 1$. For integers $N\geq 1$, let
$[N]:=\\{1,2,\dots,N\\}$ and define

$A_{P,2k}([N]):=\\#\\{(x_{1},x_{2},\dots,x_{k},y_{1},y_{2},\dots,y_{k})\in[N]^{2k}:\text{solutions
to \eqref{eqn: S}}\\}.$

We say that (1.1) has a _paucity_ of nontrivial solutions if the number of
trivial solutions (solutions where the variables are equal in pairs), which is
of exact order $k!N^{k}$, is asymptotically the total number of solutions
$A_{P,2k}([N])$. One might expect a typical equation of the form (1.1) to have
this paucity property, which is a problem raised and studied in some special
cases of $P(x)$ by Najnudel in [14]. A recent result of Klurman, Shkredov and
Xu [13, Proposition 1.4] confirms that such phenomena exist in the $k=2$ case.

###  1.1. Main results

The main result of this paper is to prove the paucity phenomenon occurs for
all $P(x)$ satisfying mild necessary conditions. Let $P(x)\in\mathbb{Z}[x]$
have degree $d\geq 2$ and at least two distinct complex roots. (In other
words, we assume that $P$ is not of the form $c(x-r)^{d}$ for any integers
$c,r$.) Let $e_{P}$ be the maximum multiplicity of any complex root of $P$;
then $1\leq e_{P}\leq d-1$. (For example, if $P(x)=x^{2}(x+1)$, then
$e_{P}=2$.)

######  Theorem 1.1.

Let $P,d,e_{P}$ be as above. Then for integers $k,N\geq 1$, we have

$A_{P,2k}([N])=k!N^{k}+O_{P,k,\varepsilon}(N^{k-1/6e_{P}+\varepsilon}).$

A related paucity result (where one takes $P$ to be a linear polynomial with
_irrational_ coefficients) was obtained independently by Bourgain, Garaev,
Konyagin, and Shparlinski [2] and Heap, Sahay and Wooley [9]. Though it may be
possible (by considering certain norms) to use our result to obtain a weaker
version of their result (weaker in the error term), we are not aware of any
(known or easy) implication in the other direction.

###  1.2. Linear case and optimality of the results

Since (1.1) is “homogeneous” in $P$, the condition $P(x)\neq c(x-r)^{d}$ in
Theorem 1.1 basically requires $P(x)$ to not be a linear transformation, i.e.
the degree $d\neq 0,1$. There has been a lot of literature in the case $d=1$,
e.g. [4, 6, 8, 10]. It is proved that in this case, $A_{P,2k}([N])$ is larger
than the number of diagonal solutions by a power of $\log N$ with degree
depending on $k$. Thus, there is no paucity phenomenon in the linear case and
our mild condition is necessary.

However, we remark that in the companion paper [15] of the authors and Pandey,
the analogous counting problem is solved for $P(x)=x$ with short interval
support $[N,N+H]$, $H\leq N/(\log N)^{Ck^{2}}$,111Bourgain, Garaev, Konyagin,
and Shparlinski [2] handled a shorter range of the form $H\leq
N^{1-\varepsilon}$. and this is one of the key steps in making progress
towards a recent question of Harper [7]. The “paucity” in the linear case
comes from the “shortness” of the interval; in comparison, the “paucity” in
the present paper comes from the fact that $P(x)$ is “genuinely of degree
$\geq 2$” (i.e. it has at least two distinct complex roots).

###  1.3. Application to random multiplicative functions

One origin of the study of Diophantine equations of the form (1.1) (when
$P(x)=x$) comes from the moments of the Steinhaus random multiplicative
function: let $f(p)$ be independent random variables (for all primes $p$)
taking values from the complex unit circle, and define $f(n)$ completely
multiplicatively for all $n$. We can use orthogonality to derive that

$\mathbb{E}_{f}\left|\sum_{n\leq N}f(n)\right|^{2k}=A_{P(n)=n,2k}([N]).$

Computing the moments here is exactly the counting problem to (1.1). As
mentioned before, we have now a good understanding of the high moments of
$\frac{1}{\sqrt{N}}\sum_{n\leq N}f(n)$. When the support of $f(n)$ is the full
set $[N]$, these moments blow up and grow like powers of $\log N$ as $k$
grows. Our recent work [15] with Pandey proves that when $f(n)$ is supported
on short intervals $[N,N+H]$ with $H\to+\infty$ and $H\ll N/(\log N)^{Ck^{2}}$
as $N\to+\infty$, these moments are consistent with Gaussian moments for all
$k\geq 1$.

Theorem 1.1 can be interpreted as the following moment asymptotic for the
Steinhaus random multiplicative function:

$\mathbb{E}_{f}\left|\frac{1}{\sqrt{N}}\sum_{n\leq
N}f(P(n))\right|^{2k}=k!+O_{k,P,\varepsilon}(N^{-1/6e_{P}+\varepsilon}).$

The moments of random sums are strongly connected to their limiting
distribution. In particular, to prove the random sums have Gaussian limiting
distribution, it is sufficient to show that all the moments match the Gaussian
moments (although this is not a necessary condition). Thus, Theorem 1.1
potentially gives a way to establish the central limit theorem for
$\frac{1}{\sqrt{N}}\sum_{n\leq N}f(P(n))$ subject to establishing the odd
moments, which one may be able to establish by using similar arguments in the
present paper. See [15] for a reference for treating odd moments in the
special case $d=1$. In fact, Klurman, Shkredov and Xu [13] already proved that
$\frac{1}{\sqrt{N}}\sum_{n\leq N}f(P(n))$ has Gaussian limiting distribution
by exploiting the martingale structure, which successfully reduces the task of
establishing a central limit theorem to a fourth moment estimate. There is a
lot more literature on the study of the distribution when $P(n)$ is linear and
we refer readers to [3, 16, 5, 11].

###  1.4. Notation

Our notation is standard. For any two functions $f,g:\mathbb{R}\to\mathbb{R}$,
we write $f\ll g,g\gg f,g=\Omega(f)$ or $f=O(g)$ if there exists a positive
constant $C$ such that $|f|\leq Cg$. We use $O_{k}$ to indicate that the
implied constant depends on $k$. Let $\omega(n)$ be the number of distinct
prime factors of $n$, and $\tau_{k}(n)$ be the number of ways to write $n$ as
an ordered product of $k$ positive integers. We use $\lfloor x\rfloor$ to
denote the largest integer $\leq x$.

### Acknowledgement

We thank Oleksiy Klurman, Ilya Shkredov, and Igor Shparlinski for helpful
comments on earlier versions of the paper, and thank Yotam Hendel for
providing a reference for Lemma 2.1.

##  2\. Proof of the main result

The proof of Theorem 1.1 shares some ideas with the authors’ recent work [15]
with Pandey, particularly the use of “factored” congruence conditions and
induction on $k$.222The work [2, proof of Lemma 22] also uses congruence
conditions, but does not “factor” them as we do in Corollary 2.3 below (and as
in [15, §2 and §3]). But the present proof also requires a novel use of some
algebra (e.g. a result of Huxley [12]) and algebraic geometry (e.g. [1,
Theorem 4] of Bombieri and Pila), among other ingredients. More specifically,
our argument requires us to separately consider three types of solutions:
trivial solutions (for which $x_{i}$ and $y_{j}$ are equal in pairs), “almost-
trivial” nontrivial solutions (for which some
$\operatorname{gcd}(P(x_{i}),P(y_{j}))$ is very large), and “generic”
nontrivial solutions. For the two types of nontrivial solutions, see
especially Corollary 2.3 and its proof. Roughly speaking, we use [12] to show
that there are relatively few “generic” nontrivial solutions, while we use [1,
Theorem 4] to show that _on average_ there are relatively few “almost-trivial”
nontrivial solutions.

We now build up to a proof of Theorem 1.1. Let $P,d,e_{P}$ be as in Theorem
1.1. Furthermore, assume that $P$ has positive leading coefficient. (Theorem
1.1 is invariant under scaling $P$, so we may certainly make this assumption.)

Let $Q(x)\in\mathbb{Z}[x]$ be a polynomial, with no repeated complex roots,
such that

$Q(x)\mid P(x)\mid Q(x)^{e_{P}}.$ (2.1)

Let $\Delta_{Q}\neq 0$ be the discriminant of $Q$. (On a first reading, we
suggest restricting to the special case $P(x)=x(x+1)$ considered by Najnudel
in [14], where $e_{P}=1$ and one can take $Q(x)=x(x+1)$ with $\Delta_{Q}=1$.)

The following result on polynomial congruences is due to Huxley [12]. See also
the stronger result in Stewart [17, Corollary 2], and the ensuing historical
discussion there.

######  Lemma 2.1 (Huxley).

Let $P,d,e_{P},Q,\Delta_{Q}$ be as above. Let $\ell\geq 1$ be an integer. Then

$\\#\\{x\bmod{\ell}:Q(x)\equiv 0\bmod{\ell}\\}\leq
d^{\omega(\ell)}\lvert\Delta_{Q}\rvert^{1/2}.$

This result leads to the following proposition:

######  Proposition 2.2.

Let $P,d,e_{P},Q,\Delta_{Q}$ be as above. Let $z,N\geq 1$ be integers. Then

$\\#\\{x\in[N]:z\mid P(x)\\}\ll
d^{\omega(z)}\lvert\Delta_{Q}\rvert^{1/2}\cdot(1+N/z^{1/e_{P}}).$

###### Proof.

Let $\ell$ be the smallest positive integer such that $z\mid\ell^{e_{P}}$.
Clearly $\ell\mid z$ and $\ell\geq z^{1/e_{P}}$.

If $x\in[N]$ and $z\mid P(x)$, then $z\mid Q(x)^{e_{P}}$ by (2.1), so
$\ell\mid Q(x)$. This implies that

$\\#\\{x\in[N]:z\mid P(x)\\}\leq\\#\\{x\in[N]:\ell\mid Q(x)\\}.$ (2.2)

But a short calculation using Lemma 2.1 gives the bound

$\\#\\{x\in[N]:\ell\mid Q(x)\\}\ll
d^{\omega(\ell)}\lvert\Delta_{Q}\rvert^{1/2}\cdot(1+N/\ell).$ (2.3)

(For $N=\ell$, the bound (2.3) follows directly from Lemma 2.1; for $N<\ell$
it then follows by upon enlarging $N$ to $\ell$, while for $N>\ell$ it follows
by reduction modulo $\ell$.)

Since $\ell\mid z$ and $\ell\geq z^{1/e_{P}}$, the proposition follows from
(2.2) and (2.3). ∎

We next deduce the following general counting result. Here and later, we let

$G_{P,\lambda}([N],z):=\\#\\{(x,a,b)\in[N]\times[\lfloor\lambda\rfloor]^{2}:a\cdot
z=b\cdot P(x)\text{ and }a<b\\},$ (2.4)

whenever $z,N,\lambda\in\mathbb{Z}$ with $N,\lambda\geq 1$. (We allow $z\leq
0$ just to keep Lemma 2.6 clean.)

######  Corollary 2.3.

Let $P,d,e_{P},Q,\Delta_{Q}$ be as above. Let $k,z,N,\lambda\geq 1$ be
integers. Then

$\begin{split}T:=\\#\\{(x_{1},x_{2},\dots,x_{k})\in[N]^{k}:{}&z\mid
P(x_{1})P(x_{2})\cdots P(x_{k}),\\\ &\min(P(x_{1}),\dots,P(x_{k}))>0,\\\
&\max(P(x_{1}),\dots,P(x_{k}))<z\\}\end{split}$ (2.5)

is at most

$k\cdot G_{P,\lambda}([N],z)\cdot N^{k-1}+\tau_{k}(z)\cdot
O(d^{\omega(z)})^{k}\lvert\Delta_{Q}\rvert^{k/2}\cdot\left(\frac{N^{k}}{z^{1/e_{P}}}+\frac{N^{k-1}}{\lambda^{1/e_{P}}}+N^{k-2}\right).$
(2.6)

###### Proof.

For each $j\in[k]$, let $I_{j}$ be the contribution to $T$ from tuples
$(x_{1},\dots,x_{k})$ with $\operatorname{gcd}(P(x_{j}),z)\geq z/\lambda$.
Note that if $x\in[N]$ has the properties $\operatorname{gcd}(P(x),z)\geq
z/\lambda$ and $0<P(x)<z$, then $\operatorname{gcd}(P(x),z)=z/b$ for some
positive integer $b\leq\lambda$ dividing $z$, and thus $P(x)=a\cdot z/b$ for
some positive integer $a<b$. It follows that for each $j\in[k]$,

$I_{j}\leq G_{P,\lambda}([N],z)\cdot N^{k-1}.$ (2.7)

Now let $S$ be the contribution to $T$ from tuples $(x_{1},\dots,x_{k})$ with
$\operatorname{gcd}(P(x_{j}),z)<z/\lambda$ for all $j\in[k]$. Given positive
integers $u_{1},\dots,u_{k}$ with $u_{1}u_{2}\cdots u_{k}=z$, let
$S_{u_{1},\dots,u_{k}}$ be the contribution to $T$ from tuples
$(x_{1},\dots,x_{k})$ with $\operatorname{gcd}(P(x_{j}),z)<z/\lambda$ and
$u_{j}\mid P(x_{j})$ for all $j\in[k]$. Then $S_{u_{1},\dots,u_{k}}=0$ unless
$u_{1},\dots,u_{k}<z/\lambda$, in which case we must have $k\geq 2$ and (by
Proposition 2.2 and the fact that $u_{1}u_{2}\cdots u_{k}=z$)

$\begin{split}S_{u_{1},\dots,u_{k}}&\leq
O(d^{\omega(z)})^{k}\lvert\Delta_{Q}\rvert^{k/2}\cdot\left(N^{k-2}+\frac{N^{k-1}}{\min_{i\in[k]}(z/u_{i})^{1/e_{P}}}+\frac{N^{k}}{z^{1/e_{P}}}\right)\\\
&\leq
O(d^{\omega(z)})^{k}\lvert\Delta_{Q}\rvert^{k/2}\cdot\left(N^{k-2}+\frac{N^{k-1}}{\lambda^{1/e_{P}}}+\frac{N^{k}}{z^{1/e_{P}}}\right).\end{split}$
(2.8)

As $S\leq\sum_{u_{1}u_{2}\cdots u_{k}=z}S_{u_{1},\dots,u_{k}}$, we conclude
that

$S\leq\tau_{k}(z)\cdot
O(d^{\omega(z)})^{k}\lvert\Delta_{Q}\rvert^{k/2}\cdot\left(N^{k-2}+\frac{N^{k-1}}{\lambda^{1/e_{P}}}+\frac{N^{k}}{z^{1/e_{P}}}\right).$
(2.9)

Since $T\leq I_{1}+\dots+I_{k}+S$, the desired result follows from (2.7) and
(2.9). ∎

To further estimate the upper bound (2.6) for the quantity (2.5), we show that
on average over certain values of $z$, the quantity $G_{P,\lambda}([N],z)$
appearing in (2.6) is not large. We show this in Lemma 2.6, with Lemma 2.4 and
Lemma 2.5 as preparations.

The next two lemmas imply that when $a\neq b$, the polynomial $a\cdot
P(y)-b\cdot P(x)$ has no linear factor, and therefore has few integral zeros
in boxes. Lemma 2.5 is a well-known counting result due to Bombieri and Pila
[1, Theorem 4], which is also used in [13].

######  Lemma 2.4.

Let $P,d$ be as above. Let $a,b$ be distinct positive integers. Then the
polynomial $a\cdot P(y)-b\cdot P(x)$ has no factor in $\mathbb{C}[x,y]$ of
degree $1$.

###### Proof.

Suppose for contradiction that $a\cdot P(y)-b\cdot P(x)$ has a linear factor,
say $fx+gy+h$ (with $f,g,h\in\mathbb{C}$ and $(f,g)\neq(0,0)$). By symmetry,
we may assume $g\neq 0$. By scaling, we may then assume $g=-1$. The factor
theorem in $(\mathbb{C}[x])[y]$ now implies

$aP(fx+h)=bP(x).$

Thus $f^{d}=b/a$ (so $|f|=(b/a)^{1/d}\notin\\{0,1\\}$), and the (invertible)
affine map $L(x):=fx+h$ on $\mathbb{C}$ permutes the roots of $P$. Thus each
root of $P$ is fixed by $L^{d!}$. But since $|f|\neq 1$, the map
$L^{d!}=f^{d!}x+O_{d,f,h}(1)$ has at most one fixed point. Thus $P$ has at
most one root, contradicting our original requirement that $P$ have at least
two distinct roots. ∎

######  Lemma 2.5 (Bombieri–Pila [1, Theorem 4]).

Let $\mathcal{C}$ be an absolutely irreducible curve (over the rationals) of
degree $d\geq 2$ and $N\geq\exp(d^{6})$. Then the number of integral points on
$\mathcal{C}$ and inside the square $[0,N]\times[0,N]$ does not exceed
$N^{1/d}\exp(12\sqrt{d\log N\log\log N})$.

Recall the definition of $G_{P,\lambda}([N],z)$ from (2.4). The following
lemma shows that $G_{P,\lambda}([N],P(y))$ is not large when we average over
$y$.

######  Lemma 2.6.

Let $P,d$ be as above. Let $N,\lambda\geq 1$ be integers. Then

$\sum_{y\in[N]}G_{P,\lambda}([N],P(y))\ll_{d,\varepsilon}\lambda^{2}\cdot
N^{1/2+\varepsilon}.$

###### Proof.

Plugging in (2.4), we see that the lemma concerns the number of solutions
$(x,y)\in[N]^{2}$ to $a\cdot P(y)=b\cdot P(x)$ on average over positive
integers $a,b\leq\lambda$ with $a<b$. Consider any such integers $a,b$. By
Lemma 2.4, each irreducible complex component of the curve $a\cdot P(y)-b\cdot
P(x)=0$ has degree $\geq 2$, and thus by Lemma 2.5, at most
$O_{d,\varepsilon}(N^{1/2+\varepsilon})$ integral points $(x,y)\in[N]^{2}$.
The lemma follows upon summing over $a,b$. ∎

With Corollary 2.3 and Lemma 2.6, we are ready to prove Theorem 1.1. The proof
proceeds by using a “congruence condition and induction” strategy.

###### Proof of Theorem 1.1.

Call a solution to (1.1) _trivial_ if the first $k$ variables
$x_{1},\dots,x_{k}$ equal the last $k$ variables $y_{1},\dots,y_{k}$ in pairs.
Clearly the contribution to $A_{P,2k}([N])$ from trivial solutions is

$k!N^{k}+O_{k}(N^{k-1}),$ (2.10)

since it is at most $k!N^{k}$ and at least $k!N(N-1)\cdots(N-k+1)$. It remains
to bound $N_{P,2k}([N])$, the number of nontrivial solutions
$(x_{1},\dots,x_{k},y_{1},\dots,y_{k})\in[N]^{2k}$ to (1.1).

Let $R_{P,2k}([N])$ be the contribution to $N_{P,2k}([N])$ from tuples with

$y_{k}=\max(y_{1},\dots,y_{k})=\max(x_{1},\dots,x_{k})=x_{k}.$

Let $N^{\prime}_{P,2k}([N])$ be the contribution to $N_{P,2k}([N])$ from
tuples with

$y_{k}=\max(y_{1},\dots,y_{k})>\max(x_{1},\dots,x_{k}).$ (2.11)

Then by symmetry,

$N_{P,2k}([N])\leq k^{2}\cdot R_{P,2k}([N])+2k\cdot N^{\prime}_{P,2k}([N]).$
(2.12)

Since $R_{P,2k}([N])\leq N\cdot N_{P,2k-2}([N])$ and $R_{P,2}([N])=0$, it
follows from (2.12) that

$N_{P,2k}([N])\leq k\cdot\max_{1\leq j\leq k}((k^{2}N)^{k-j}\cdot 2j\cdot
N^{\prime}_{P,2j}([N])).$ (2.13)

We now bound $N^{\prime}_{P,2k}([N])$. We assume $P$ has positive leading
coefficient. Let $M(P)$ be a positive integer (depending only on $P$) such
that for all integers $n\geq M(P)$, we have

$P(n)>\max(P(0),P(1),\dots,P(n-1))\qquad\text{and}\qquad P(n)\geq n^{d}/2.$
(2.14)

Let $M\geq M(P)$ and $\lambda\geq 1$ denote integers to be chosen later. Then
by Corollary 2.3 (applied with $z:=P(y_{k})$ for $y_{k}\geq M$, using the
conditions (2.11) and (2.14)) and the divisor bound, the quantity
$N^{\prime}_{P,2k}([N])$ is at most

$M^{2k}+O_{k,d,\Delta_{Q},\varepsilon}(N^{\varepsilon})\cdot\sum_{M\leq y\leq
N}\left(N^{k-1}G_{P,\lambda}([N],P(y))+\frac{N^{k}}{y^{d/e_{P}}}+\frac{N^{k-1}}{\lambda^{1/e_{P}}}+N^{k-2}\right).$
(2.15)

By Lemma 2.6, we conclude that

$N^{\prime}_{P,2k}([N])\leq
M^{2k}+O_{k,d,\Delta_{Q},\varepsilon}(N^{\varepsilon})\cdot(\lambda^{2}N^{k-1/2}+N^{k}M^{1-d/e_{P}}+N^{k}\lambda^{-1/e_{P}}+N^{k-1}).$
(2.16)

Choosing $M:=\lfloor M(P)\cdot N^{1/4}\rfloor$ (say) and $\lambda:=\lfloor
N^{1/6}\rfloor$, we get (since $1-d/e_{P}\leq-1/e_{P}$)

$N^{\prime}_{P,2k}([N])\leq
O_{k,d,\Delta_{Q},M(P),\varepsilon}(N^{\varepsilon})\cdot(N^{k/2}+N^{k-1/4e_{P}}+N^{k-1/6e_{P}})\ll_{k,P,\varepsilon}N^{k-1/6e_{P}+\varepsilon}.$
(2.17)

Finally, plug (2.17) into (2.13) to get
$N_{P,2k}([N])\ll_{k,P,\varepsilon}N^{k-1/6e_{P}+\varepsilon}$. This, plus the
estimate (2.10) for trivial solutions, establishes the theorem. ∎

## References

  * [1] E. Bombieri and J. Pila.  The number of integral points on arcs and ovals.  Duke Math. J., 59(2):337–357, 1989.
  * [2] J. Bourgain, M. Z. Garaev, S. V. Konyagin, and I. E. Shparlinski.  Multiplicative congruences with variables from short intervals.  J. Anal. Math., 124:117–147, 2014.
  * [3] S. Chatterjee and K. Soundararajan.  Random multiplicative functions in short intervals.  Int. Math. Res. Not. IMRN, 3:479–492, 2012.
  * [4] A. Granville and K. Soundararajan.  Large character sums.  J. Amer. Math. Soc., 14(2):365–397, 2001.
  * [5] A. J. Harper.  On the limit distributions of some sums of a random multiplicative function.  J. Reine Angew. Math., 678:95–124, 2013.
  * [6] A. J. Harper.  Moments of random multiplicative functions, II: High moments.  Algebra Number Theory, 13(10):2277–2321, 2019.
  * [7] A. J. Harper.  A note on character sums over short moving intervals, 2022.  arXiv:2203.09448.
  * [8] A. J. Harper, A. Nikeghbali, and M. Radziwill.  A note on Helson’s conjecture on moments of random multiplicative functions.  In Analytic number theory, pages 145–169. Springer, Cham, 2015\.
  * [9] W. Heap, A. Sahay, and T. D. Wooley.  A paucity problem associated with a shifted integer analogue of the divisor function.  J. Number Theory, 242:660–668, 2023.
  * [10] W. P. Heap and S. Lindqvist.  Moments of random multiplicative functions and truncated characteristic polynomials.  Q. J. Math., 67(4):683–714, 2016.
  * [11] B. Hough.  Summation of a random multiplicative function on numbers having few prime factors.  Math. Proc. Cambridge Philos. Soc., 150(2):193–214, 2011.
  * [12] M. N. Huxley.  A note on polynomial congruences.  In Recent progress in analytic number theory, Vol. 1 (Durham, 1979), pages 193–196. Academic Press, London-New York, 1981.
  * [13] O. Klurman, I. D. Shkredov, and M. W. Xu.  On the random Chowla conjecture.  arXiv e-prints, page arXiv:2202.08767, Feb. 2022.
  * [14] J. Najnudel.  On consecutive values of random completely multiplicative functions.  Electron. J. Probab., 25:Paper No. 59, 28, 2020.
  * [15] M. Pandey, V. Y. Wang, and M. W. Xu.  Partial sums of typical multiplicative functions over short moving intervals, 2022.  arXiv:2207.11758.
  * [16] K. Soundararajan and M. W. Xu.  Central limit theorems for random multiplicative functions, in preparation.
  * [17] C. L. Stewart.  On the number of solutions of polynomial congruences and Thue equations.  J. Amer. Math. Soc., 4(4):793–835, 1991.