# Hadron-quark phase transition in neutron star by combining the relativistic
Brueckner-Hartree-Fock theory and Dyson-Schwinger equation approach

Pianpian Qin  Department of Physics and Chongqing Key Laboratory for Strongly
Coupled Physics, Chongqing University, Chongqing 401331, China  Zhan Bai
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  Sibo Wang  Department of Physics and Chongqing Key Laboratory for
Strongly Coupled Physics, Chongqing University, Chongqing 401331, China
Chencan Wang  School of Physics and Astronomy, Sun Yat-Sen University, Zhuhai
519082, China  Si-xue Qin<EMAIL_ADDRESS>Department of Physics and Chongqing
Key Laboratory for Strongly Coupled Physics, Chongqing University, Chongqing
401331, China

###### Abstract

Starting from the relativistic Brueckner-Hartree-Fock theory for nuclear
matter and the Dyson-Schwinger equation approach for quark matter, the
possible hadron-quark phase transition in the interior of a neutron star is
explored. The first-order phase transition and crossover are studied by
performing the Maxwell construction and three-window construction
respectively. The mass-radius relation and the tidal deformability of the
hybrid star are calculated and compared to the joint mass-radius observation
of a neutron star and the constraints from gravitational wave detection. For
the Maxwell construction, no stable quark core is found in the interior of a
neutron star. For the three-window construction, the parameters of the smooth
interpolation function are chosen in such a way to keep the thermodynamic
stability and lead to a moderate crossover density region. To support a two-
solar-mass neutron star under the three-window construction, the effective
width of medium screening effects in quark matter should be around $0.35$ GeV.

##  I Introduction

At sufficiently high temperature or density, the quarks are likely to “escape”
from nucleons and become basic degrees of freedom. Such transition is known as
the hadron-quark phase transition. Unlike the transition at high temperature,
which is reachable in heavy ion collisions, the transition at large density is
far beyond the scope of terrestrial experiments. The corresponding density can
only be found in astronomical compact objects, such as neutron stars. The
neutron star is the remnant of a supernova explosion, and the density inside a
neutron star could reach about $5\sim 10\rho_{\text{sat}}$, where
$\rho_{\text{sat}}=0.16~{}\text{fm}^{-3}$ is the nuclear saturation density
Lattimer (2004, 2021). At such a density, it is very likely that the hadron-
quark phase transition will take place and deconfined quark matter appears
Annala  _et al._ (2020); Bauswein  _et al._ (2019). Therefore, it is possible
to use an astronomical observation of neutron stars to constrain the theories
about the dense matter and phase transition.

Currently, it is difficult to unitedly describe the hadron matter, quark
matter and hadron-quark phase transition within a single theoretical
framework. To study the hadron-quark phase transition inside the hybrid star,
it is common to describe the hadron and quark matter separately with
corresponding approaches, and then use different construction schemes to
combine them and get a complete equation of state (EOS).

For the hadron matter, there have already been many nuclear many-body methods,
among which the nonrelativistic Rikovska Stone  _et al._ (2003); Fattoyev  _et
al._ (2013) and relativistic Fattoyev  _et al._ (2018); Li  _et al._ (2020);
Fu  _et al._ (2022); Sun  _et al._ (2008); Tong  _et al._ (2020) density
functional theories (DFTs) are very important. They are based on effective
nucleon-nucleon $(NN)$ interactions, where the parameters are determined by
fitting the ground state properties of finite nuclei and infinite nuclear
matter at saturation. The predictions of the neutron star properties with
different DFTs are rather divergent due to the loose constraints of the
effective $NN$ interactions at higher densities. In comparison, one can start
from realistic $NN$ interactions where the parameters are constrained with the
$NN$ scattering data in free space, and utilize the many-body methods to deal
with the realistic $NN$ interactions, such as the relativistic Brueckner-
Hartree-Fock (RBHF) theory Shen  _et al._ (2019) as well as its
nonrelativistic counterpart BHF theory Baldo and Maieron (2007). In
particular, the RBHF theory is rooted in the relativistic framework, which
contains significant three-body force effects self-consistently Sammarruca
_et al._ (2012) and naturally avoids the problem of superluminance in
nonrelativistic methods Li  _et al._ (1992). The RBHF theory has been
successfully applied to study the neutron star with pure hadron matter Krastev
and Sammarruca (2006); Wang  _et al._ (2020a); Tong  _et al._ (2022).

In the description of quark matter, phenomenological models are widely used
(see, e.g., Refs. Hatsuda and Kunihiro (1994); Buballa and Oertel (1999);
Fukushima (2004); Benvenuto and Lugones (1995); Torres and Menezes (2013);
Chodos  _et al._ (1974); Zacchi  _et al._ (2015)). In spite of these
successful models, a study that is directly based on quantum chromodynamics
(QCD) is still needed. The most important first-principle method, the lattice
QCD Montvay and Münster (1997); Ding and et al. (2019), is powerful at zero
chemical potential but encounters the notorious “sign problem” at large
density. The perturbative QCD Kurkela  _et al._ (2010); Kurkela and Vuorinen
(2016) also loses its power in the phase transition region. Therefore, a
nonperturbative, continuum approach is required to study the phase transition
of cold-dense matter. In particular, the Dyson-Schwinger equation (DSE)
approach Roberts and Schmidt (2000); Fischer (2019) is a typical functional
method based on QCD. It can deal with the confinement and dynamical chiral
symmetry breaking simultaneously, and can be naturally used at finite
temperature and chemical potential. In recent years, there have already been
studies of neutron stars by using the DSE method Chen  _et al._ (2011, 2012,
2015); Yasutake  _et al._ (2016); Bai  _et al._ (2018); Bai and Liu (2021).
However, the relevant study is still in its early stages, and more efforts
must be paid to make this method mature.

As for the construction of phase transition, the most widely used method is
the Maxwell construction Endo  _et al._ (2006); Hempel  _et al._ (2009);
Yasutake  _et al._ (2009). It assumes that the phase transition is of first
order Glendenning (1992); Bhattacharyya  _et al._ (2010) and a stable quark
core is formed inside the neutron star. It has also been argued that with the
increase of density, the boundaries of hadrons will gradually disappear and
the transition from hadron matter to quark matter will be a smooth crossover
Baym  _et al._ (2018). Under this assumption, the three-window construction
Masuda  _et al._ (2013a, b) was proposed. In this paper, both of these schemes
will be used.

Apart from those theoretical approaches, on the experimental side, the study
of hadron-quark phase transition relies closely on the astronomical
observations of neutron stars. Several massive neutron stars have already been
detected with high-precision mass measurements Demorest  _et al._ (2010);
Fonseca and et al. (2016); Antoniadis and et al. (2013); Cromartie and et al.
(2020); Fonseca and et al. (2021). These observations provide the lower limit
for the stiffness of the neutron star EOS. The joint mass-radius observation
of the neutron star from the Neutron Star Interior Composition Explorer
(NICER) mission provides additional requirements Riley and et al. (2019);
Miller and et al. (2019); Riley and et al. (2021); Miller and et al. (2021).
The detection of gravitational wave (GW) signals from a binary neutron star
merger Abbott and et al. (2017, 2018) has also provided an important
constraint for the neutron star properties.

In previous works, numerious methods for nuclear matter and quark matter have
been implemented to study the hadron-quark phase transition in the hybrid star
Schertler  _et al._ (1999); Agrawal (2010); Logoteta  _et al._ (2013); Orsaria
_et al._ (2013); Wu and Shen (2017); Miao  _et al._ (2020); Ju  _et al._
(2021); Contrera  _et al._ (2022); Huang  _et al._ (2022); Burgio  _et al._
(2002); Klähn  _et al._ (2007); Agrawal (2010); Li  _et al._ (2015); Baym  _et
al._ (2019); Nicotra  _et al._ (2006); Bai  _et al._ (2018); Bai and Liu
(2021). In Refs. Chen  _et al._ (2011, 2012, 2015); Yasutake  _et al._ (2016),
the BHF theory and the DSE approach have been combined. Nevertheless, this
prescription is inconsistent with respect to the relativity, as the BHF theory
is nonrelativistic while the DSE approach is relativistic. This paper utilizes
the RBHF theory for nuclear matter and the DSE approach for quark matter,
which are both in the relativistic framework. To construct the EOS of the
hybrid star, the Maxwell construction and three-window construction are
employed to describe the first-order phase transition and crossover,
respectively. The mass-radius relation and tidal deformability of the hybrid
star are calculated and are compared with the astronomical observations.

This paper is organized as follows. In Sec. II, the RBHF theory for nuclear
matter, the DSE approach for quark matter, and the construction schemes for
hadron-quark phase transition are briefly introduced. The results and
discussions are presented in Sec. III. Finally, a summary is given in Sec. IV.

##  II Nuclear matter, quark matter, and hadron-quark phase transition

###  II.1 Nuclear matter

In the RBHF theory, the nuclear matter is described with the nucleons as the
basic degrees of freedom. The single-particle motion of a nucleon in nuclear
medium is described by the Dirac equation

$(\bm{\alpha}\cdot\bm{p}+\beta
M_{\tau}+\beta\hat{U}_{\tau})u_{\tau}(\bm{p},s)=E_{\bm{p},\tau}u_{\tau}(\bm{p},s)\,,$
(1)

where the subscript $\tau=n,p$ indicates neutron or proton. $M_{\tau}$ is the
mass of a free nucleon. $u_{\tau}(\bm{p},s)$ is the Dirac spinor of a nucleon
with momentum $\bm{p}$, spin $s$, and single-particle energy
$E_{\bm{p},\tau}$. The single-particle potential operator $\hat{U}_{\tau}$ is
generally divided into scalar and vector components

$\hat{U}_{\tau}(\bm{p})=U_{S,\tau}(p)+\gamma^{0}U_{0,\tau}(p)+\bm{\gamma}\cdot\hat{\bm{p}}U_{V,\tau}(p)\,.$
(2)

Here $\hat{\bm{p}}=\bm{p}/p$ is the unit vector parallel to the momentum
$\bm{p}$. The quantities $U_{S,\tau}(p)$, $U_{0,\tau}(p)$, and $U_{V,\tau}(p)$
are the scalar potential, the timelike part, and the spacelike part of the
vector potential. The momentum dependent potentials can be determined uniquely
by considering the positive- and negative-energy states simultaneously, i.e.,
the RBHF theory in the full Dirac space Wang  _et al._ (2021). However, this
method is now limited up to the density $\rho=0.57\ \text{fm}^{-3}$ Wang  _et
al._ (2022), which is not high enough to study the hadron-quark phase
transition in neutron stars. This density limitation is absent if the
momentum-independence approximation Brockmann and Machleidt (1990) is
utilized, where $U_{S,\tau}$ and $U_{0,\tau}$ are approximated to be constants
and $U_{V,\tau}$ is neglected. Within this approximation, the single-particle
potential operator $\hat{U}_{\tau}$ is written as

$\hat{U}_{\tau}=U_{S,\tau}+\gamma^{0}U_{0,\tau}\,.$ (3)

By defining the effective nucleon mass $M^{*}_{\tau}$ and effective energy
$E^{*}_{\bm{p},\tau}$

$M^{*}_{\tau}=M_{\tau}+U_{S,\tau}\,,\qquad
E^{*}_{\bm{p},\tau}=E_{\bm{p},\tau}-U_{0,\tau}\,,$ (4)

the Dirac equation (1) in nuclear medium can be rewritten as

$(\bm{\alpha}\cdot\bm{p}+\beta
M^{*}_{\tau})u_{\tau}(\bm{p},s)=E^{*}_{\bm{p},\tau}u_{\tau}(\bm{p},s).$ (5)

From the Dirac equation (5), the dispersion relation
$E^{*}_{\bm{p},\tau}=\sqrt{\bm{p}^{2}+M^{*2}_{\tau}}$ is obtained and the
Dirac spinor can be solved exactly as

$u_{\tau}(\bm{p},s)=\sqrt{\frac{E^{*}_{\bm{p},\tau}+M^{*}_{\tau}}{2M^{*}_{\tau}}}\begin{pmatrix}1\\\
\frac{\bm{\sigma}\cdot\bm{p}}{E^{*}_{\bm{p},\tau}+M^{*}_{\tau}}\end{pmatrix}\chi_{s}\chi_{\tau}\,,$
(6)

where $\chi_{s}$ and $\chi_{\tau}$ are the spin and isospin wave function
respectively. The normalization condition is
$\bar{u}_{\tau}(\bm{p},s)u_{\tau}(\bm{p},s)=1$.

In the RBHF theory, the single-particle potentials $U_{S,\tau}$ and
$U_{0,\tau}$ are self-consistently determined with the effective $NN$
interaction, $G$ matrix, which can be obtained by solving the Thompson
equation Thompson (1970); Brockmann and Machleidt (1990)

$G_{\tau\tau^{\prime}}(\bm{q^{\prime}},\bm{q}|\bm{P})=V_{\tau\tau^{\prime}}(\bm{q^{\prime}},\bm{q})+\int\frac{d^{3}k}{(2\pi)^{3}}V_{\tau\tau^{\prime}}(\bm{q^{\prime}},\bm{k})\frac{Q_{\tau\tau^{\prime}}(\bm{k},\bm{P})}{W_{\tau\tau^{\prime}}-E^{*}_{\tau\tau^{\prime}}}G_{\tau\tau^{\prime}}(\bm{k},\bm{q}|\bm{P})\,,$
(7)

where $\tau\tau^{\prime}=$ $nn$, $pp$, or $np$. In Eq. (7),
$\bm{P}=(\bm{k}_{1}+\bm{k}_{2})/2$ is the center-of-mass momentum, and
$\bm{k}=(\bm{k}_{1}-\bm{k}_{2})/2$ is the relative momentum of two interacting
nucleons with momenta $\bm{k}_{1}$ and $\bm{k}_{2}$ in the rest frame of
nuclear matter. The quantities $\bm{q}$, $\bm{q^{\prime}}$, and $\bm{k}$ are
the initial, final, and intermediate relative momenta of the two nucleons
scattering in nuclear matter, respectively.
$W_{\tau\tau^{\prime}}=E^{*}_{\bm{P}+\bm{q},\tau}+E^{*}_{\bm{P}-\bm{q},\tau^{\prime}}$
is the starting energy, and
$E^{*}_{\tau\tau^{\prime}}=E^{*}_{\bm{P}+\bm{k},\tau}+E^{*}_{\bm{P}-\bm{k},\tau^{\prime}}$
is the total single-particle energy of intermediate two-nucleon states. The
Pauli operator $Q_{\tau\tau^{\prime}}$ avoids the $NN$ scattering to occupied
states in the Fermi sea and is defined as

$Q_{\tau\tau^{\prime}}(\bm{k},\bm{P})=\left\\{\begin{aligned}
&1&\qquad&\text{if}\quad|\bm{P}+\bm{k}|>k_{F}^{\tau}\quad\text{and}\quad|\bm{P}-\bm{k}|>k_{F}^{\tau^{\prime}}\,,\\\
&0&\qquad&\text{otherwise}\,,\end{aligned}\right.$ (8)

where $k_{F}^{\tau}$ represents the Fermi momentum for nucleon $\tau$. For
nuclear matter with total nucleon density $\rho=\rho_{n}+\rho_{p}$ and isospin
asymmetry $\delta=(\rho_{n}-\rho_{p})/(\rho_{n}+\rho_{p})$, the Fermi momentum
is calculated as $k_{F}^{\tau}=\left[3\pi^{2}(1\pm\delta)\rho/2\right]^{1/3}$.

With the $G$ matrix, one can calculate the single-particle potential energy

$U_{\tau}(p)=\sum_{s^{\prime},\tau^{\prime}}\int_{0}^{k^{\tau^{\prime}}_{F}}\frac{d^{3}p^{\prime}}{(2\pi)^{3}}\frac{M^{*}_{\tau^{\prime}}}{E^{*}_{\bm{p}^{\prime},\tau^{\prime}}}\langle\bar{u}_{\tau}(\bm{p},1/2)\bar{u}_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})|\bar{G}(W)|u_{\tau}(\bm{p},1/2)u_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})\rangle\,,$
(9)

where $\bar{G}$ is the antisymmetried $G$ matrix and the factor
$M^{*}_{\tau^{\prime}}/E^{*}_{\bm{p}^{\prime},\tau^{\prime}}$ comes from the
normalization condition above. Alternatively, the single-particle potential
energy can be obtained by sandwiching the single-particle potential operator
between the Dirac spinors

$U_{\tau}(p)=\frac{M^{*}_{\tau}}{E^{*}_{\bm{p},\tau}}\langle\bar{u}_{\tau}(\bm{p},1/2)|\hat{U}_{\tau}|u_{\tau}(\bm{p},1/2)\rangle=\frac{M^{*}_{\tau}}{E^{*}_{\bm{p},\tau}}U_{S,\tau}+U_{0,\tau}\,.$
(10)

By combining Eqs. (9) and (10), one can extract the two constants $U_{S,\tau}$
and $U_{0,\tau}$ with two momenta, e.g., $0.5k^{\tau}_{F}$ and $k^{\tau}_{F}$.

Equations (1), (7), (9), and (10) constitute a coupled set of equations that
needs to be solved self-consistently. Starting from initial values of
$U^{(0)}_{S,\tau},U^{(0)}_{0,\tau}$ in vacuum, the Dirac spinors are obtained
by solving the Dirac equation (1). Then one solves the Thompson equation (7)
to get the $G$ matrix and obtain the single-particle potential energy by using
the integrals in Eq. (9). From Eq. (10) a new set of values for
$U^{(1)}_{S,\tau},U^{(1)}_{0,\tau}$ are found, which are to be used in the
next iteration. This iterative procedure is repeated until the satisfactory
convergence is reached.

Once the solution is converged, the binding energy per nucleon in nuclear
matter can be calculated as

$\begin{split}\frac{E}{A}=&\frac{1}{\rho}\sum_{s,\tau}\int_{0}^{k^{\tau}_{F}}\frac{d^{3}p}{(2\pi)^{3}}\frac{M^{*}_{\tau}}{E^{*}_{\bm{p},\tau}}\langle\bar{u}_{\tau}(\bm{p},s)|\gamma\cdot\bm{p}+M_{\tau}|u_{\tau}(\bm{p},s)\rangle\\\
&\
-\frac{1-\delta}{2}M_{p}-\frac{1+\delta}{2}M_{n}+\frac{1}{2\rho}\sum_{s,s^{\prime},\tau,\tau^{\prime}}\int_{0}^{k^{\tau}_{F}}\frac{d^{3}p}{(2\pi)^{3}}\int_{0}^{k^{\tau^{\prime}}_{F}}\frac{d^{3}p^{\prime}}{(2\pi)^{3}}\frac{M^{*}_{\tau}}{E^{*}_{\bm{p},\tau}}\frac{M^{*}_{\tau^{\prime}}}{E^{*}_{\bm{p}^{\prime},\tau^{\prime}}}\\\
&\
\times\langle\bar{u}_{\tau}(\bm{p},s)\bar{u}_{\tau^{\prime}}(\bm{p^{\prime}},s^{\prime})|\bar{G}(W)|u_{\tau}(\bm{p},s)u_{\tau^{\prime}}(\bm{p^{\prime}},s^{\prime})\rangle\,.\end{split}$
(11)

The neutron star matter is regarded as the beta-equilibrium nuclear matter
consisting of protons, neutrons, electrons, and muons. The energy density
$\varepsilon_{H}$ of the neutron star matter or pure hadron star is calculated
as $\varepsilon_{H}=\rho e_{\text{tot}}$, where the total energy
$e_{\text{tot}}$ is defined as

$e_{\text{tot}}(\rho,Y_{n},Y_{p},Y_{e},Y_{\mu})=E/A(\rho,Y_{p})+Y_{p}M_{p}+Y_{n}M_{n}+E_{e}/A+E_{\mu}/A.$
(12)

The quantities $E_{e}/A$ and $E_{\mu}/A$ are the contributions from electrons
and muons, which are treated as gas of relativistic noninteracting fermions.
The equilibrium particle concentrations $Y_{i}=\rho_{i}/\rho\ (i=n,p,e,\mu)$
can be calculated via the $\beta$-stability condition and charge neutrality
condition

$\displaystyle\mu_{n}-\mu_{p}=$ $\displaystyle\ \mu_{e},$ (13a)
$\displaystyle\mu_{n}-\mu_{p}=$ $\displaystyle\ \mu_{\mu},$ (13b)
$\displaystyle\rho_{e}+\rho_{\mu}=$ $\displaystyle\ \rho_{p},$ (13c)

where $\mu_{i}\ (i=n,p,e,\mu)$ is the chemical potential for particle $i$. For
electrons and muons, the chemical potential is obtained via

$\mu_{i}=\frac{\partial e_{\text{tot}}}{\partial Y_{i}}.$ (14)

For protons and neutrons, the chemical potential is calculated as the single-
particle energy at the Fermi surface. Once the total energy $e_{\text{tot}}$
is calculated, the pressure of the neutron star matter can be obtained as

$P_{H}=\rho^{2}\frac{\partial e_{\text{tot}}}{\partial\rho}.$ (15)

Owing to the cluster effects in nuclear matter with density $\rho<0.08\
\text{fm}^{-3}$, the RBHF theory is not applicable and the EOS introduced with
the Baym-Bethe-Pethick Baym  _et al._ (1971) and Baym-Pethick-Sutherland
models Baym  _et al._ (1971) is used. It should also be noticed that the
results of the RBHF calculation are not so smooth that numerical derivatives
can be performed. Parametrization or regression strategies Krastev and
Sammarruca (2006); Katayama and Saito (2013) are often employed to the neutron
star EOS. In this work, Gaussian process regression Huang  _et al._ (2022) is
implemented to process our RBHF results.

###  II.2 Quark matter

The DSEs are the equations of motion for fields. They can be derived by
differentiating the action of QCD and contain all the information of QCD
Lagrangian. In this paper, the DSE of the quark propagator is concerned, the
quark-gluon vertex is truncated by a symmetry preserving scheme, and the gluon
propagator is approximated by the interaction model. Therefore, the
complicated DSE set can be reduced to a solvable gap equation. The Feynman
diagram for the gap equation is shown in Fig. 1.

Figure 1: The Feynman diagram for the gap equation. The solid (wavy) line with
black thick dot is for the dressed quark (gluon) propagator, the thick gray
dot represents the dressed quark-gluon interaction vertex, and the thin black
dot represents the bare quark-gluon interaction vertex.

At zero temperature and finite chemical potential, the gap equation for the
quark propagator $S(p;\mu)$ reads as

$S^{-1}(p;\mu)=Z_{2}(i\bm{\gamma}\cdot\bm{p}+i\gamma_{4}\tilde{p}_{4}+m_{q})+\Sigma(p;\mu)\,,$
(16)

where $p=(\bm{p},p_{4})$ is the four-momentum and $\tilde{p}_{4}=p_{4}+i\mu$
with $\mu$ is the quark chemical potential. In Eq. (16), the Euclidean metric
is used, which is different from the Minkowski metric used in Eq. (3). $m_{q}$
is the current mass of the quark $q$ with $q=u,d,s$. In this paper, the
current-quark masses $m_{u/d}=0$ are chosen for simplicity, while
$m_{s}=115~{}\text{MeV}$ is obtained by fitting the $K$ meson mass in vacuum
Alkofer  _et al._ (2002). $\Sigma(p;\mu)$ is the renormalized self-energy of
the quark

$\Sigma(p;\mu)=Z_{1}\int\frac{d^{4}q}{(2\pi)^{4}}g^{2}(\mu)D_{\rho\sigma}(k;\mu)\frac{\lambda^{a}}{2}\gamma_{\rho}S(q,\mu)\frac{\lambda^{a}}{2}\Gamma_{\sigma}(q,p;\mu)\,,$
(17)

where $Z_{1}$, $Z_{2}$ are the renormalization constants. $\lambda^{a}/2$ is
the fundamental representation of $SU(3)$ color symmetry.
$D_{\rho\sigma}(k;\mu)$ with $k=p-q$ is the renormalized dressed gluon
propagator, and $\Gamma_{\sigma}(q,p;\mu)$ is the renormalized dressed quark-
gluon vertex. $g(\mu)$ is the density-dependent coupling constant.

To solve the gap equation (16), we have to know the quark-gluon vertex
$\Gamma_{\sigma}(q,p;\mu)$ and gluon propagator $D_{\rho\sigma}(k;\mu)$. In
this paper, the dressed quark-gluon vertex $\Gamma_{\sigma}(q,p;\mu)$ is
truncated by using the rainbow approximation Maris and Roberts (1997), and the
nonperturbative dressing effect of the dressed gluon propagator
$D_{\rho\sigma}(k;\mu)$ can be absorbed in the effective interaction function
$\mathcal{G}(k^{2};\mu)$ Alkofer  _et al._ (2002). Therefore, the interaction
kernel can be expressed as

$Z_{1}g^{2}(\mu)D_{\rho\sigma}(k;\mu)\Gamma_{\sigma}(q,p;\mu)=\mathcal{G}(k^{2};\mu)D_{\rho\sigma}^{\text{free}}(k)\gamma_{\sigma}\,,$
(18)

where
$D_{\rho\sigma}^{\text{free}}(k)=\frac{1}{k^{2}}\left(\delta_{\rho\sigma}-k_{\rho}k_{\sigma}/k^{2}\right)$
is the free gluon propagator in the Landau gauge. The interaction function
$\mathcal{G}(k^{2};\mu)$ is usually divided into two parts: the infrared part
and ultraviolet perturbative part. At zero temperature, the quark properties
are mainly determined by its infrared behavior, so we omit the ultraviolet
part as in Ref. Chen  _et al._ (2011) for better numerical behavior. The gluon
interaction function we applied is Chen  _et al._ (2011)

$\frac{\mathcal{G}(k^{2};\mu)}{k^{2}}=\frac{4\pi^{2}D}{\omega^{6}}k^{2}e^{-k^{2}/\omega^{2}}e^{-\mu^{2}/\omega_{\text{eff}}^{2}}\,.$
(19)

Hence, the integration is finite at the ultraviolet limit, and the
renormalization procedure can be omitted by simply setting all the
renormalization factors to unit. The parameters $D$ and $\omega$ control the
strength and width of the interaction in vacuum respectively. It is found that
the observables of vector and pseudoscalar mesons are insensitive to
variations of $\omega\in[0.4,0.6]$ as long as $D\omega=\text{constant}$ Maris
and Tandy (1999); Qin  _et al._ (2011). As in Refs. Chen  _et al._ (2011); Bai
_et al._ (2018), we choose $\omega=0.5$ GeV and $D=1.0~{}\text{GeV}^{2}$.

In the interaction model Eq. (19), the factor
$e^{-\mu^{2}/\omega_{\text{eff}}^{2}}$ depicts the effects of the medium
screening on the interaction at finite chemical potential $\mu$, where
$\omega_{\text{eff}}$ is the effective width of medium screening effects. The
larger the $\omega_{\text{eff}}$, the weaker the medium screening effects and
the stronger the gluon interaction. In Refs. Chen  _et al._ (2011, 2016); Bai
_et al._ (2018); Bai and Liu (2021), the effective width of medium screening
effects is primarily adopted as $\omega_{\text{eff}}\lesssim 0.5$ GeV to study
the hadron-quark phase transition. In this paper, we will also study the
dependence of our results on the effective width $\omega_{\text{eff}}$.

With the interaction model in Eqs. (18) and (19), the self-energy in Eq. (17)
can be reduced as

$\Sigma(p;\mu)=\int\frac{d^{4}q}{(2\pi)^{4}}\frac{4\pi^{2}D}{\omega^{6}}k^{2}e^{-k^{2}/\omega^{2}}e^{-\mu^{2}/\omega_{\text{eff}}^{2}}\left(\delta_{\rho\sigma}-k_{\rho}k_{\sigma}/k^{2}\right)\frac{\lambda^{a}}{2}\gamma_{\rho}S(q,\mu)\frac{\lambda^{a}}{2}\gamma_{\sigma}\,.$
(20)

Therefore, the gap equation (16) is simplified as an solvable equation in
terms of the quark propagator $S(q,\mu)$, which can be decomposed according to
the Lorentz structure

$S(p;\mu)^{-1}=i\bm{\gamma}\cdot\bm{p}A(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})+B(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})+i\gamma_{4}\tilde{p}_{4}C(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})\,.$
(21)

The scalar functions $A(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})$,
$B(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})$, and
$C(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})$ can then be solved with the gap
equation. It is known that the gap equation has multiple solutions. In vacuum,
massless quarks may have the solution $B(p)\equiv 0$ or $B(p)\neq 0$, which
are called the Wigner solution and Nambu solution, respectively. The Wigner
solution corresponds to the dynamical chiral symmetry (DCS) phase, where the
quarks are bare and have no dynamical mass. The Nambu solution corresponds to
the dynamical chiral symmetry breaking (DCSB) phase, where the massless quarks
are dressed and the mass functions $M(p)=B(p)/A(p)$ acquire nonzero values.
Although there are arguments that there might exist quarkyonic matter that is
DCS but still confined, it is usually believed that DCSB and confinement
appear simultaneously. Therefore, the Wigner solution corresponds to the
deconfined quark phase, and the Nambu solution corresponds to the confined
hadron phase. To describe the quark core in the interior of a neutron star at
high density, we only consider the Wigner solution in this paper. In the left
and right panel in Fig. 2, the Wigner solution of the quark propagator for the
massless quark is shown as the functions of momentum and chemical potential,
respectively. At zero chemical potential, the functions $A$ and $C$ are
identical as the Lorentz $O(4)$ symmetry is satisfied. At zero momentum, the
function $A$ shows different chemical potential dependence in comparison to
the function $C$ at intermediate chemical potential, while both functions $A$
and $C$ approach unit as the quark propagator approaches asymptotic freedom at
high chemical potential.

Figure 2: The Wigner solution of the quark propagator for massless quark with
varying momentum (left) or varying chemical potential (right).

With the Wigner solution, the distribution function $f_{1}(\bm{p};\mu)$ can be
obtained from the quark propagator

$f_{1}(\bm{p};\mu)=\frac{1}{4\pi}\int_{-\infty}^{\infty}dp_{4}\
\text{tr}_{D}[-\gamma_{4}S(p;\mu)]$ (22)

where the trace is for the spinor indices. With the quark propagator, the
integration in $f_{1}(\bm{p};\mu)$ can be converted to a contour integral on
the complex plane of $\tilde{p_{4}}$,

$f_{1}(\bm{p};\mu)=\frac{1}{\pi}\int_{-\infty}^{\infty}dp_{4}\
\frac{i(p_{4}+i\mu)C\left(\left|\bm{p}\right|^{2},(p_{4}+i\mu)^{2}\right)}{\mathcal{M}},$
(23)

where

$\mathcal{M}=\bm{p}^{2}A^{2}(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})+\tilde{p}_{4}^{2}C^{2}(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2})+B(\left|\bm{p}\right|^{2},\tilde{p}_{4}^{2}).$
(24)

Then the number density $n_{q}$ with quark flavor $q=u,d,s$ can be obtained as
a function of chemical potential Chen  _et al._ (2008)

$n_{q}(\mu)=2N_{c}N_{f}\int\frac{d^{3}p}{(2\pi)^{3}}f_{1}(\bm{p};\mu)\,,$ (25)

where $N_{f}=3$ and $N_{c}=3$ denote the number of flavor and color,
respectively.

The quark star matter is composed of quarks and leptons under the $\beta$
equilibrium and electric charge neutral condition,

$\displaystyle\mu_{d}=$ $\displaystyle\ \mu_{u}+\mu_{e}=\mu_{s}\,,$ (26a)
$\displaystyle\mu_{d}=$ $\displaystyle\ \mu_{u}+\mu_{\mu}\,,$ (26b)
$\displaystyle\frac{2n_{u}-n_{d}-n_{s}}{3}=$ $\displaystyle\ n_{e}+n_{\mu}\,,$
(26c)

where the chemical potential for quark is denoted by $\mu_{q}$ with $q=u,d,s$.
In Eq. (26c), the number density of lepton $n_{l}$ with $l=e,\mu$ is
calculated as $n_{l}=k^{3}_{Fl}/3\pi^{2}$ where the Fermi momentum $k_{Fl}$ is
related to the chemical potential $\mu_{l}$ by
$k^{2}_{Fl}=\mu^{2}_{l}-m^{2}_{l}$. In the description of the pure quark star,
we take $m_{e}=0.511~{}\text{MeV}$ and $m_{\mu}=105~{}\text{MeV}$. In
practice, for given $\mu_{n}$ with

$\mu_{n}=\mu_{u}+2\mu_{d},$ (27)

one can obtain the chemical potential $\mu_{i}$ with $i=u,d,s,e,\mu$ by
solving Eqs. (26) and (27).

Once the quark chemical potential $\mu_{q}$ is obtained, the number density
$n_{q}(\mu_{q})$ can be calculated with Eq. (25) and the baryon number density
$\rho$ is found as $\rho=\frac{1}{3}(n_{u}+n_{d}+n_{s})$. Furthermore, one can
calculate the pressure by integrating the number density

$P_{q}(\mu_{q})=P_{q}(\mu_{q,0})+\int_{\mu_{q,0}}^{\mu_{q}}d\mu\
n_{q}(\mu)\,.$ (28)

Theoretically, the starting point of the integral $\mu_{q,0}$ can be any
value. In this paper we take $\mu_{u,0}=\mu_{d,0}=0$ and choose $\mu_{s,0}$ as
the value of the starting point of the Wigner phase. Similar to Ref. Chen  _et
al._ (2011), we take the vacuum pressure
$P_{u}(\mu_{u,0})=P_{d}(\mu_{d,0})=-45~{}\text{MeV}\cdot\text{fm}^{-3}$ and
$P_{s}(\mu_{s,0})=0$. Since the leptons are treated as free Fermi gas, their
pressure is

$P_{l}=\frac{1}{24\pi^{2}}\left[k_{Fl}\mu_{l}\left(2k_{Fl}^{2}-3m_{l}^{2}\right)+3m_{l}^{4}\ln\left(\left|\frac{k_{Fl}+\mu_{l}}{m_{l}}\right|\right)\right]\,,\quad
l=e,\mu\,.$ (29)

Finally, with the chemical potential $\mu_{i}$ and number density $n_{i}$ with
$i=u,d,s,e,\mu$, the total pressure and total energy density of the pure quark
star can be obtained as

$\displaystyle P_{Q}=$ $\displaystyle\sum_{i=u,d,s,e,\mu}P_{i}(\mu_{i})\,,$
(30a) $\displaystyle\varepsilon_{Q}=$
$\displaystyle\sum_{i=u,d,s,e,\mu}\mu_{i}n_{i}-P_{Q}\,.$ (30b)

###  II.3 Hadron-quark phase transition

In this paper, we will study the hadron-quark phase transition as a first-
order transition as well as a crossover. The first-order phase transition is
described with the Maxwell construction, and the crossover is described with
the three-window construction. For the Maxwell construction, the phase
transition occurs when the baryon chemical potential and pressure of the two
phases are equal

$P_{H}(\mu_{n,c})=P_{Q}(\mu_{n,c})\,,$ (31)

where $\mu_{n,c}$ is the critical baryon chemical potential of the first-order
phase transition. Therefore, the pressure of the hybrid star under the Maxwell
construction is

$P(\mu_{n})=\left\\{\begin{aligned}
P_{H}\,,\qquad\text{if}\quad\mu_{n}\leq\mu_{n,c}\,,\\\
P_{Q}\,,\qquad\text{if}\quad\mu_{n}>\mu_{n,c}\,.\end{aligned}\right.$ (32)

Correspondingly, the energy density of the hybrid star under the Maxwell
construction is

$\varepsilon(\mu_{n})=\left\\{\begin{aligned}
\varepsilon_{H}\,,\qquad\text{if}\quad\mu_{n}\leq\mu_{n,c}\,,\\\
\varepsilon_{Q}\,,\qquad\text{if}\quad\mu_{n}>\mu_{n,c}\,.\end{aligned}\right.$
(33)

For the three-window construction, a smooth interpolation of the energy
density between the hadron and quark phases is performed

$\varepsilon(\rho)=f_{-}(\rho)\varepsilon_{H}(\rho)+f_{+}(\rho)\varepsilon_{Q}(\rho)\,.$
(34)

The interpolation function $f_{\pm}$ is chosen as Masuda  _et al._ (2013a, b)

$f_{\pm}=\frac{1}{2}\left(1\pm\tanh\left(\frac{\rho-\bar{\rho}}{\Gamma}\right)\right)\,,$
(35)

where the parameters $\bar{\rho}$ and $\Gamma$ describe the center density and
the width of the transition region.

Taking parameters ($\bar{\rho},\Gamma$)=($3.5,1.5$) as an example, the
variation behavior of the function $f_{\pm}$ in terms of the baryon density is
shown in Fig. 3. Once the energy density of the hybrid star has been obtained,
the pressure can be determined with the thermodynamic relation,

$\displaystyle
P(\rho)=\rho^{2}\frac{\partial(\varepsilon/\rho)}{\partial\rho}.$ (36)

Figure 3: The interpolation functions $f_{\pm}$ with respect to the baryon
density $\rho/\rho_{\text{sat}}$. The center density and the width of the
transition region are $(\bar{\rho},\Gamma$)=($3.5,1.5$). The shaded region
denotes the crossover region $\bar{\rho}-\Gamma\leq\rho\leq\bar{\rho}+\Gamma$.

With the EOS of the neutron star matter, the mass-radius relation of the
neutron star can be obtained by solving the Tolman–Oppenheimer–Volkov equation
Oppenheimer and Volkoff (1939); Tolman (1939), and the tidal deformability can
be calculated as in Ref. Hinderer  _et al._ (2010).

##  III Results and discussion

Figure 4: Left panel: energy density and pressure as functions of density for
the pure hadron star described by the RBHF theory with potentials pvCD-Bonn A,
B, and C. Right panel: similar to the left panel, but for the pure quark star
described by the DSE approach with different effective width
$\omega_{\text{eff}}$ of medium screening in the interaction model (19).

In the left panel of Fig. 4, the energy density $\varepsilon$ and pressure $P$
of the pure hadron star are calculated by the RBHF theory with $NN$
interactions pvCD-Bonn A, B, and C Wang  _et al._ (2019). It is clear that the
difference of the results from pvCD-Bonn A, B, and C is negligible. This is
reasonable since the main difference between the three parametrizations is in
the tensor force strength, which is mostly reflected in the ($T=0$)
${}^{3}S_{1}$-${}^{3}D_{1}$ states with $T$ the total isospin. This partial
wave does not contribute to the ($T=1$) neutron-neutron state Krastev and
Sammarruca (2006), which is dominant in the neutron star. In the following
discussion, the potential pvCD-Bonn A is used, since the empirical nuclear
saturation properties can be described satisfactorily by the RBHF theory with
pvCD-Bonn A Wang  _et al._ (2020b). In the right panel, the energy density
$\varepsilon$ and pressure $P$ of the pure quark star are calculated by the
DSE approach with different effective width $\omega_{\text{eff}}$ of medium
screening in the interaction model (19). It is found that with the effective
width $\omega_{\text{eff}}$ increasing, i.e., the increasing of the gluon
interaction, the energy density and pressure become larger. Besides, the
differences of the energy density and pressure become more evident at higher
density.

Figure 5: The $P\mbox{-}\mu_{n}$ relation (left) and $P\mbox{-}\rho$ relation (right) of the hybrid star under Maxwell construction with nuclear matter described by the RBHF theory with potential pvCD-Bonn A and quark matter described by the DSE approach with different effective width $\omega_{\text{eff}}$ of medium screening in the interaction model (19). In comparison, the results of the hadron star described by the RBHF theory with potential pvCD-Bonn A and the $P\mbox{-}\mu$ relation obtained by the BHF theory with potential Bonn B from Ref. Chen  _et al._ (2011) are also given. Table 1: Critical chemical potential $\mu_{n,c}$, critical pressure $P(\mu_{n,c})$, and the phase transition density region $[\rho_{H},\rho_{Q}]$ obtained with the RBHF theory and DSE approach under the Maxwell construction. The results obtained with the BHF theory and DSE approach from Ref. Chen  _et al._ (2011) are also shown. Model | $\mu_{n,c}$ | $P(\mu_{n,c})$ | $\rho_{H}$ | $\rho_{Q}$
---|---|---|---|---
[GeV] | [GeV$\cdot\text{fm}^{-3}$] | [$\rho_{\text{sat}}$] | [$\rho_{\text{sat}}$]
A-0.5 | 1.995 | 0.916 | 8.105 | 16.689
A-0.35 | 1.640 | 0.471 | 6.078 | 10.940
A-0.25 | 1.366 | 0.225 | 4.379 | 7.189
BHF-0.35 | 1.416 | 0.193 | 3.556 | 5.744

Figure 5 shows the $P\mbox{-}\mu_{n}$ relation and $P\mbox{-}\rho$ relation of
the hybrid star under Maxwell construction. The legend A-$\omega_{\text{eff}}$
denotes that the nuclear matter is described by the RBHF theory with potential
pvCD-Bonn A and the quark matter is described by the DSE approach with
effective width $\omega_{\text{eff}}$ of medium screening. For comparison, the
results of the pure hadron star described by the RBHF theory are also given.
In the left panel, the intersection points of the EOSs in $P\mbox{-}\mu_{n}$
plane are the critical points of the first-order phase transition in Eq. (31).
In the right panel, the pressure is constant at the critical point of the
first-order phase transition, while the density jumps from $\rho_{H}$ to
$\rho_{Q}$. The values of $\mu_{n,c}$, $P(\mu_{n,c})$, $\rho_{H}$, and
$\rho_{Q}$ are listed in Table 1. It is found that for a larger
$\omega_{\text{eff}}$, the critical baryon chemical potential $\mu_{n,c}$ and
density region $[\rho_{H},\rho_{Q}]$ are higher, and the corresponding
pressure $P(\mu_{n,c})$ is larger. For $\omega_{\text{eff}}=0.5$ GeV, the
corresponding density of the critical point is about
$8.1~{}\rho_{\text{sat}}$. Considering that the nucleon degrees of freedom
becomes less available at higher density, the parameter $\omega_{\text{eff}}$
should be smaller than $0.5$ GeV to obtain a reasonable EOS of the hybrid
star. For comparison, in the left panel of Fig. 5, we also show the
$P\mbox{-}\mu$ relation calculated by the BHF theory with the potential Bonn B
and corresponding three-body forces from Ref. Chen  _et al._ (2011). The
critical properties obtained with the Maxwell construction are also shown in
the last row in Table 1. It is found that the nonrelativistic calculations
lead to a softer EOS and the first-order phase transition happens much earlier
than the relativistic ones.

Figure 6: The mass-radius relation of the hybrid star under Maxwell
construction with nuclear matter described by the RBHF theory with potential
pvCD-Bonn A and quark matter described by the DSE approach with different
effective width $\omega_{\text{eff}}$ of medium screening in the interaction
model (19). In comparison, the results of the pure hadron star described by
the RBHF theory are also given.

The mass-radius relation of the hybrid star under the Maxwell construction are
plotted in Fig. 6. It can be found that, once the first-order transition
happens, the mass and radius decrease simultaneously. It is known that, if the
mass decreases with respect to the increase of central density, the neutron
star is unstable against oscillation. Therefore, this sharpness reflects that
the appearance of quark matter leads to an unstable hybrid star. In other
words, under the Maxwell construction with the present models, there is no
such a quark core in the interior of a neutron star. We note that similar
results are also found by combining the BHF theory with the DSE approach in
Ref. Chen  _et al._ (2011).

Figure 7: The pressure of the hybrid star as a function of density under the
three-window construction with four representative parameter sets
$(\bar{\rho},\Gamma)$. The results for the pure hadron star are also shown.
The effective width $\omega_{\text{eff}}=0.35$ GeV is fixed for the
description of quark matter with the DSE approach.

Under the three-window construction, a phenomenological interpolation of the
energy density between the hadron and quark phases is performed, where the
interpolation function is given in Eq. (35). The parameters $\bar{\rho}$ and
$\Gamma$ characterize the center and effective width of the crossover region
in baryon density. Inspired by and basing them on Ref. Masuda  _et al._
(2013a), we consider two constraints on the choice of the two parameters: (1)
the system is always thermodynamically stable, i.e., $dP/d\rho>0$, and (2) the
crossover density region should be moderate to avoid the failure of the
nuclear matter method and quark matter method, i.e.,
$\bar{\rho}-\Gamma>\rho_{\text{sat}}$ and
$\bar{\rho}+\Gamma<6\rho_{\text{sat}}$. In addition, a stiffer EOS is favored
to satisfy the observation of massive neutron stars.

Figure 7 depicts the $P\mbox{-}\rho$ relation of the hybrid star under the
three-window construction. The quark matter effective width
$\omega_{\text{eff}}=0.35$ GeV is chosen as an example and four representative
parameter sets ($\bar{\rho},\Gamma$) are shown. The parameter sets $(3.5,1.0)$
and $(4.0,1.5)$ lead to nonmonotonic pressure and should be discarded. The
maximum density of the parameter set $(4.0,2.0)$ is $6\rho_{\text{sat}}$,
which is somehow too large for the RBHF theory. The parameter set ($3.5,1.5$)
satisfies the conditions (1) and (2) and is chosen in the following
calculations.

Figure 8: The $\varepsilon\mbox{-}\rho$ relation (left) and $P\mbox{-}\rho$
relation (right) of the hybrid star under the three-window construction with
nuclear matter described by the RBHF theory with potential pvCD-Bonn A and
quark matter described by the DSE approach with different effective width
$\omega_{\text{eff}}$ of medium screening in the interaction model (19). In
comparison, the results of the hadron star described by the RBHF theory are
also given.

In Fig. 8, the $\varepsilon\mbox{-}\rho$ relation and $P\mbox{-}\rho$ relation
of the hybrid star under the three-window construction are shown. In the left
panel, the $\varepsilon\mbox{-}\rho$ relations are interpolated in the
crossover region. Outside the crossover region, pure nuclear matter or pure
quark matter dominates. In the right panel, the pressure in the crossover
region increases monotonically with an increasing baryon density for the
effective width $\omega_{\text{eff}}$ of different quarks. From Fig. 8, it is
clear that with a lager $\omega_{\text{eff}}$, the energy density
$\varepsilon$ and pressure $P$ are larger in the crossover region.

Figure 9: The mass-radius relation of the hybrid star under the three-window construction with nuclear matter described by the RBHF theory with potential pvCD-Bonn A and quark matter described by the DSE approach with different effective width $\omega_{\text{eff}}$ of medium screening in the interaction model (19). In comparison, the results of the hadron star described by the RBHF theory are also given. Table 2: The maximum mass $M_{\text{max}}$ of the neutron star, the corresponding radius $R_{M_{\text{max}}}$, and the radius for $1.4M_{\odot}$ neutron star $R_{1.4M_{\odot}}$ under the three-window construction with different effective width $\omega_{\text{eff}}$. Model | $M_{\text{max}}$ | $R_{M_{\text{max}}}$ | $R_{1.4M_{\odot}}$
---|---|---|---
[$M_{\odot}$] | [km] | [km]
pvCD-Bonn A | 2.198 | 11.15 | 12.47
A-0.5 | 2.530 | 13.54 | 13.93
A-0.35 | 2.102 | 12.91 | 13.20
A-0.25 | 1.539 | 12.39 | 12.57

The mass-radius relation of the hybrid star under the three-window
construction is plotted in Fig. 9. The maximum mass with different models, the
corresponding radius, and the radius for $1.4M_{\odot}$ neutron star are
listed in Table 2. For $\omega_{\text{eff}}=0.25$ GeV, the maximum mass of the
hybrid star is about $1.5~{}M_{\odot}$, which cannot support the astrophysical
observation of massive neutron stars. For $\omega_{\text{eff}}=0.5$ GeV, the
maximum mass of the hybrid star is $2.53~{}M_{\odot}$. This is in
contradiction to the results from binary neutron star mergers, which require
that the EOS cannot be too stiff, and provide an upper bound for the maximum
mass Annala  _et al._ (2018). For $\omega_{\text{eff}}=0.35$ GeV, the maximum
mass of the hybrid star is $2.1~{}M_{\odot}$, which is consistent with the
current constraints from astrophysical observation. Therefore, to obtain a
hybrid star supporting $2~{}M_{\odot}$ with three-window construction, the
effective width $\omega_{\text{eff}}$ of medium screening effects should be
close to $0.35$ GeV.

Figure 10: The $P\mbox{-}\varepsilon$ relation of the pure hadron star, the
pure quark star and the hybrid star under the Maxwell construction and three-
window construction. Nuclear matter is described by the RBHF theory with
potential pvCD-Bonn A and quark matter is described by the DSE approach with
effective width $\omega_{\text{eff}}=0.35$ GeV.

In Fig. 10, the EOSs of the pure hadron star, the pure quark star, and the
hybrid star under the Maxwell construction and three-window construction are
compared. For the Maxwell construction, there is a clear plateau of pressure
as a function of energy density. This corresponds to the latent heat of the
first-order phase transition. For the three-window construction, the central
density of the phase transition region is $\bar{\rho}=3.5\rho_{\text{sat}}$
with an effective width $\Gamma=1.5\rho_{\text{sat}}$. Outside the crossover
region, the EOS of the hybrid star asymptotically approaches that of the pure
hadron star or the pure quark star. We also find the crossover region is lower
than the first-order phase transition region in terms of energy density.

Figure 11: The mass-radius relation of the pure hadron star, the pure quark
star, and the hybrid star under the Maxwell construction and three-window
construction. Nuclear matter is described by the RBHF theory with potential
pvCD-Bonn A and quark matter is described by the DSE approach with the
effective width $\omega_{\text{eff}}=0.35$ GeV. The dark (light) green and
purple regions indicate the $68\%(95\%)$ confidence intervals constrained by
the NICER analysis of PSR J0030+0451 Miller and et al. (2019) and PSR
J0740+6620 Miller and et al. (2021).

In Fig. 11, the mass-radius relation of the pure hadron star, the pure quark
star, and the hybrid star under the Maxwell construction and three-window
construction is depicted. Except for the pure quark star modeled by the DSE
approach, both the pure hadron star and the hybrid star can support
$2M_{\odot}$. In comparison, the joint constraints of the mass and radius of
neutron stars are also shown. The $68\%$ and $95\%$ contours of the joint
probability density distribution of the mass and radius of PSR J0030+0451
Miller and et al. (2019) and PSR J0740+6620 Miller and et al. (2021) from the
NICER analysis are also shown. It can be found that the mass radius of the
pure hadron star, the hybrid stars with the Maxwell, and three-window
constructions are all consistent with the recent constraints by NICER. The
radii of a $1.4M_{\odot}$ hybrid star $R_{1.4M_{\odot}}$ under the Maxwell
construction and three-window construction are $12.47$ km and $13.20$ km,
respectively.

Figure 12: Tidal deformability of the pure hadron star, the pure quark star,
and the hybrid star under the Maxwell construction and three-window
construction. Nuclear matter is described by the RBHF theory with potential
pvCD-Bonn A and quark matter is described by the DSE approach with the
effective width $\omega_{\text{eff}}=0.35$ GeV. The constraints from GW170817
$\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ Abbott and et al. (2018) are also
shown.

The tidal deformability of the pure hadron star, the pure quark star, and the
hybrid star under the Maxwell construction and three-window construction are
shown in Fig. 12. The tidal deformability of the pure hadron star and that of
the hybrid star under the Maxwell construction share the same value, $473$,
which is consistent with the constraints
$\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ from GW170817 Abbott and et al.
(2018). The tidal deformability of the hybrid star under the three-window
construction is $715$, which is slightly higher than that from the
astronomical observation. We notice that the RBHF theory used in the present
work can be improved by considering the negative-energy states in the full
Dirac space Wang  _et al._ (2021); Tong  _et al._ (2022), where the tidal
deformability can be reduced to a value which is much closer to the center
value from GW170817. Therefore, the three-window construction with the nuclear
matter described by the RBHF theory in the full Dirac space might lead to a
tidal deformability which is consistent with the GW constraints.

##  IV Summary and prospects

The possible hadron-quark phase transition is explored by combining the RBHF
theory for nuclear matter and the DSE approach for quark matter. The Maxwell
construction and three-window construction are implemented to study the first-
order phase transition and crossover, respectively. For the Maxwell
construction, the phase transition occurs where the baryon chemical potential
and pressure of the two phases are equal. With the RBHF theory and DSE
approach, there is no stable quark core in the interior of a neutron star,
which confirms previous studies with nonrelativistic Brueckner-Hartree-Fock
theory and the DSE approach. For the three-window construction, a smooth
interpolation of the energy density between the hadron and quark phases is
performed, where the parameters in the interpolation function are chosen in
such a way as to keep the thermodynamic stability and lead to a moderate
crossover density region. To support a two-solar-mass neutron star under the
three-window construction, the effective width of medium screening effects in
quark matter should be around $\omega_{\text{eff}}=0.35$ GeV. The mass-radius
relation of the hybrid star is consistent with the joint mass-radius
observation, while the tidal deformability of a $1.4$ solar mass is found
slightly higher than the constraints from gravitational wave detection.

In the future, this work can be extended by improving the theoretical methods
for hadron matter, quark matter, and construction schemes. The RBHF theory can
be improved by considering the negative energy states in the full Dirac space.
For the DSE framework, the parameters cannot be changed arbitrarily as they
are determined by the hadron properties as well as the phase transition at
finite temperature. To obtain a stiffer EOS for the quark matter, it is
necessary to improve the vertex and gluon truncation in the DSE approach, for
example, by using the Ball-Chiu or Chang-Liu-Roberts vertex Ball and Chiu
(1980); Qin  _et al._ (2013) and including more interaction channels Eichmann
_et al._ (2016) as well as higher order corrections Gao  _et al._ (2021). As
for the construction schemes, considering the studies by combing BHF and DSE
approaches Chen  _et al._ (2011), the Gibbs construction is likely to support
a stable quark core. After the corrections and improvements are realized, it
is hopeful to achieve a better understanding of the hadron-quark phase
transition in neutron stars.

###### Acknowledgements.

This work was partly supported by the China Postdoctoral Science Foundation
under Grants No. 2021M700610 and No. 2022M723230; the CAS Project for Young
Scientists in Basice Research (YSBR060); the National Natural Science
Foundation of China under Grants No. 12147102, No. 12005060, No. 12205353, and
No. 12205030; and the Fundamental Research Funds for the Central Universities
under Grants No. 2020CDJQY-Z003 and No. 2021CDJZYJH-003.

## References

  * Lattimer (2004) J. M. Lattimer, Science 304, 536 (2004).
  * Lattimer (2021) J. Lattimer, Annu. Rev. Nucl. Part. Sci. 71, 433 (2021).
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nat. Phys. 16, 907 (2020).
  * Bauswein  _et al._ (2019) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019).
  * Rikovska Stone  _et al._ (2003) J. Rikovska Stone, J. C. Miller, R. Koncewicz, P. D. Stevenson, and M. R. Strayer, Phys. Rev. C 68, 034324 (2003).
  * Fattoyev  _et al._ (2013) F. J. Fattoyev, J. Carvajal, W. G. Newton, and B.-A. Li, Phys. Rev. C 87, 015806 (2013).
  * Fattoyev  _et al._ (2018) F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Phys. Rev. Lett. 120, 172702 (2018).
  * Li  _et al._ (2020) J. J. Li, A. Sedrakian, and F. Weber, Phys. Lett. B 810, 135812 (2020).
  * Fu  _et al._ (2022) H. R. Fu, J. J. Li, A. Sedrakian, and F. Weber, Phys. Lett. B 834, 137470 (2022).
  * Sun  _et al._ (2008) B. Y. Sun, W. H. Long, J. Meng, and U. Lombardo, Phys. Rev. C 78, 065805 (2008).
  * Tong  _et al._ (2020) H. Tong, P. Zhao, and J. Meng, Phys. Rev. C 101, 035802 (2020).
  * Shen  _et al._ (2019) S. Shen, H. Liang, W. Long, J. Meng, and P. Ring, Prog. Part. Nucl. Phys. 109, 103713 (2019).
  * Baldo and Maieron (2007) M. Baldo and C. Maieron, J. Phys. G 34, R243 (2007).
  * Sammarruca  _et al._ (2012) F. Sammarruca, B. Chen, L. Coraggio, N. Itaco, and R. Machleidt, Phys. Rev. C 86, 054317 (2012).
  * Li  _et al._ (1992) G. Q. Li, R. Machleidt, and R. Brockmann, Phys. Rev. C 45, 2782 (1992).
  * Krastev and Sammarruca (2006) P. G. Krastev and F. Sammarruca, Phys. Rev. C 74, 025808 (2006).
  * Wang  _et al._ (2020a) C. Wang, J. Hu, Y. Zhang, and H. Shen, Astrophys. J. 897, 96 (2020a).
  * Tong  _et al._ (2022) H. Tong, C. Wang, and S. Wang, Astrophys. J. 930, 137 (2022).
  * Hatsuda and Kunihiro (1994) T. Hatsuda and T. Kunihiro, Phys. Rep. 247, 221 (1994).
  * Buballa and Oertel (1999) M. Buballa and M. Oertel, Phys. Lett. B 457, 261 (1999).
  * Fukushima (2004) K. Fukushima, Phys. Lett. B 591, 277 (2004).
  * Benvenuto and Lugones (1995) O. G. Benvenuto and G. Lugones, Phys. Rev. D 51, 1989 (1995).
  * Torres and Menezes (2013) J. R. Torres and D. P. Menezes, Europhys. Lett. 101, 42003 (2013).
  * Chodos  _et al._ (1974) A. Chodos, R. L. Jaffe, K. Johnson, C. B. Thorn, and V. F. Weisskopf, Phys. Rev. D 9, 3471 (1974).
  * Zacchi  _et al._ (2015) A. Zacchi, R. Stiele, and J. Schaffner-Bielich, Phys. Rev. D 92, 045022 (2015).
  * Montvay and Münster (1997) I. Montvay and G. Münster,  _Quantum Fields on a Lattice_ (Cambridge University Press, Cambridge, England, 1997).
  * Ding and et al. (2019) H. T. Ding and et al. (HotQCD Collaboration), Phys. Rev. Lett. 123, 062002 (2019).
  * Kurkela  _et al._ (2010) A. Kurkela, P. Romatschke, and A. Vuorinen, Phys. Rev. D 81, 105021 (2010).
  * Kurkela and Vuorinen (2016) A. Kurkela and A. Vuorinen, Phys. Rev. Lett. 117, 042501 (2016).
  * Roberts and Schmidt (2000) C. D. Roberts and S. M. Schmidt, Prog. Part. Nucl. Phys. 45, S1 (2000).
  * Fischer (2019) C. S. Fischer, Prog. Part. Nucl. Phys. 105, 1 (2019).
  * Chen  _et al._ (2011) H. Chen, M. Baldo, G. F. Burgio, and H.-J. Schulze, Phys. Rev. D 84, 105023 (2011).
  * Chen  _et al._ (2012) H. Chen, M. Baldo, G. F. Burgio, and H.-J. Schulze, Phys. Rev. D 86, 045006 (2012).
  * Chen  _et al._ (2015) H. Chen, J.-B. Wei, M. Baldo, G. F. Burgio, and H.-J. Schulze, Phys. Rev. D 91, 105002 (2015).
  * Yasutake  _et al._ (2016) N. Yasutake, H. Chen, T. Maruyama, and T. Tatsumi, J. Phys. Conf. Ser. 665, 012068 (2016).
  * Bai  _et al._ (2018) Z. Bai, H. Chen, and Y.-x. Liu, Phys. Rev. D 97, 023018 (2018).
  * Bai and Liu (2021) Z. Bai and Y.-X. Liu, Eur. Phys. J. C 81, 612 (2021).
  * Endo  _et al._ (2006) T. Endo, T. Maruyama, S. Chiba, and T. Tatsumi, Prog. Theor. Phys. 115, 337 (2006).
  * Hempel  _et al._ (2009) M. Hempel, G. Pagliara, and J. Schaffner-Bielich, Phys. Rev. D 80, 125014 (2009).
  * Yasutake  _et al._ (2009) N. Yasutake, T. Maruyama, and T. Tatsumi, Phys. Rev. D 80, 123009 (2009).
  * Glendenning (1992) N. K. Glendenning, Phys. Rev. D 46, 1274 (1992).
  * Bhattacharyya  _et al._ (2010) A. Bhattacharyya, I. N. Mishustin, and W. Greiner, J. Phys. G 37, 025201 (2010).
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rep. Prog. Phys. 81, 056902 (2018).
  * Masuda  _et al._ (2013a) K. Masuda, T. Hatsuda, and T. Takatsuka, Astrophys. J. 764, 12 (2013a).
  * Masuda  _et al._ (2013b) K. Masuda, T. Hatsuda, and T. Takatsuka, Prog. Theor. Exp. Phys. 2013, 073D01 (2013b).
  * Demorest  _et al._ (2010) P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Nature 467, 1081 (2010).
  * Fonseca and et al. (2016) E. Fonseca and et al., Astrophys. J. 832, 167 (2016).
  * Antoniadis and et al. (2013) J. Antoniadis and et al., Science 340, 1233232 (2013).
  * Cromartie and et al. (2020) H. T. Cromartie and et al., Nat. Astron. 4, 72 (2020).
  * Fonseca and et al. (2021) E. Fonseca and et al., Astrophys. J. Lett. 915, L12 (2021).
  * Riley and et al. (2019) T. E. Riley and et al., Astrophys. J. 887, L21 (2019).
  * Miller and et al. (2019) M. C. Miller and et al., Astrophys. J. 887, L24 (2019).
  * Riley and et al. (2021) T. E. Riley and et al., Astrophys. J. Lett. 918, L27 (2021).
  * Miller and et al. (2021) M. C. Miller and et al., Astrophys. J. Lett. 918, L28 (2021).
  * Abbott and et al. (2017) B. P. Abbott and et al. (LIGO Scientific Collaboration and VIRGO Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott and et al. (2018) B. P. Abbott and et al. (The LIGO Scientific and the Virgo Collaborations), Phys. Rev. Lett. 121, 161101 (2018).
  * Schertler  _et al._ (1999) K. Schertler, S. Leupold, and J. Schaffner-Bielich, Phys. Rev. C 60, 025801 (1999).
  * Agrawal (2010) B. K. Agrawal, Phys. Rev. D 81, 023009 (2010).
  * Logoteta  _et al._ (2013) D. Logoteta, C. m. c. Providência, and I. Vidaña, Phys. Rev. C 88, 055802 (2013).
  * Orsaria  _et al._ (2013) M. Orsaria, H. Rodrigues, F. Weber, and G. A. Contrera, Phys. Rev. D 87, 023001 (2013).
  * Wu and Shen (2017) X. H. Wu and H. Shen, Phys. Rev. C 96, 025802 (2017).
  * Miao  _et al._ (2020) Z. Miao, A. Li, Z. Zhu, and S. Han, Astrophys. J. 904, 103 (2020).
  * Ju  _et al._ (2021) M. Ju, X. Wu, F. Ji, J. Hu, and H. Shen, Phys. Rev. C 103, 025809 (2021).
  * Contrera  _et al._ (2022) G. A. Contrera, D. Blaschke, J. P. Carlomagno, A. G. Grunfeld, and S. Liebing, Phys. Rev. C 105, 045808 (2022).
  * Huang  _et al._ (2022) K. Huang, J. Hu, Y. Zhang, and H. Shen, Astrophys. J. 935, 88 (2022).
  * Burgio  _et al._ (2002) G. F. Burgio, M. Baldo, P. K. Sahu, and H.-J. Schulze, Phys. Rev. C 66, 025802 (2002).
  * Klähn  _et al._ (2007) T. Klähn, D. Blaschke, F. Sandin, C. Fuchs, A. Faessler, H. Grigorian, G. Röpke, and J. Trümper, Phys. Lett. B 654, 170 (2007).
  * Li  _et al._ (2015) A. Li, W. Zuo, and G. X. Peng, Phys. Rev. C 91, 035803 (2015).
  * Baym  _et al._ (2019) G. Baym, S. Furusawa, T. Hatsuda, T. Kojo, and H. Togashi, Astrophys. J. 885, 42 (2019).
  * Nicotra  _et al._ (2006) O. E. Nicotra, M. Baldo, G. F. Burgio, and H.-J. Schulze, Phys. Rev. D 74, 123001 (2006).
  * Wang  _et al._ (2021) S. Wang, Q. Zhao, P. Ring, and J. Meng, Phys. Rev. C 103, 054319 (2021).
  * Wang  _et al._ (2022) S. Wang, H. Tong, Q. Zhao, C. Wang, P. Ring, and J. Meng, Phys. Rev. C 106, L021305 (2022).
  * Brockmann and Machleidt (1990) R. Brockmann and R. Machleidt, Phys. Rev. C 42, 1965 (1990).
  * Thompson (1970) R. H. Thompson, Phys. Rev. D 1, 110 (1970).
  * Baym  _et al._ (1971) G. Baym, H. A. Bethe, and C. J. Pethick, Nucl. Phys. A 175, 225 (1971).
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Katayama and Saito (2013) T. Katayama and K. Saito, Phys. Rev. C 88, 035805 (2013).
  * Alkofer  _et al._ (2002) R. Alkofer, P. Watson, and H. Weigel, Phys. Rev. D 65, 094026 (2002).
  * Maris and Roberts (1997) P. Maris and C. D. Roberts, Phys. Rev. C 56, 3369 (1997).
  * Maris and Tandy (1999) P. Maris and P. C. Tandy, Phys. Rev. C 60, 055214 (1999).
  * Qin  _et al._ (2011) S.-X. Qin, L. Chang, Y.-X. Liu, C. D. Roberts, and D. J. Wilson, Phys. Rev. C 84, 042202 (2011).
  * Chen  _et al._ (2016) H. Chen, J. B. Wei, and H. J. Schulze, Eur. Phys. J. A 52, 291 (2016).
  * Chen  _et al._ (2008) H. Chen, W. Yuan, L. Chang, Y.-X. Liu, T. Klähn, and C. D. Roberts, Phys. Rev. D 78, 116015 (2008).
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * Tolman (1939) R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Hinderer  _et al._ (2010) T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Phys. Rev. D 81, 123016 (2010).
  * Wang  _et al._ (2019) C. Wang, J. Hu, Y. Zhang, and H. Shen, Chin. Phys. C 43, 114107 (2019).
  * Wang  _et al._ (2020b) C. Wang, J. Hu, Y. Zhang, and H. Shen, J. Phys. G 47, 105108 (2020b).
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018).
  * Ball and Chiu (1980) J. S. Ball and T.-W. Chiu, Phys. Rev. D 22, 2542 (1980).
  * Qin  _et al._ (2013) S.-X. Qin, L. Chang, Y.-X. Liu, C. D. Roberts, and S. M. Schmidt, Phys. Lett. B 722, 384 (2013).
  * Eichmann  _et al._ (2016) G. Eichmann, C. S. Fischer, and C. A. Welzbacher, Phys. Rev. D 93, 034013 (2016).
  * Gao  _et al._ (2021) F. Gao, J. Papavassiliou, and J. M. Pawlowski, Phys. Rev. D 103, 094013 (2021).


# An Estimator for Matching Size in Low Arboricity Graphs with Two
Applications 111This work is supported by the Iranian Institute for Research
in Fundamental Sciences (IPM), Project Number 98050014.

Hossein Jowhari 222 Department of Computer Science and Statistics, Faculty of
Mathematics, K. N. Toosi University of Technology. Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a new simple degree-based estimator for the size of
maximum matching in bounded arboricity graphs. When the arboricity of the
graph is bounded by $\alpha$, the estimator gives a $\alpha+2$ factor
approximation of the matching size. For planar graphs, we show the estimator
does better and returns a $3.5$ approximation of the matching size.

Using this estimator, we get new results for approximating the matching size
of planar graphs in the streaming and distributed models of computation. In
particular, in the vertex-arrival streams, we get a randomized
$O(\frac{\sqrt{n}}{\varepsilon^{2}}\log n)$ space algorithm for approximating
the matching size within $(3.5+\varepsilon)$ factor in a planar graph on $n$
vertices. Similarly, we get a simultaneous protocol in the vertex-partition
model for approximating the matching size within $(3.5+\varepsilon)$ factor
using $O(\frac{n^{2/3}}{\varepsilon^{2}}\log n)$ communication from each
player.

In comparison with the previous estimators, the estimator in this paper does
not need to know the arboricity of the input graph and improves the
approximation factor for the case of planar graphs.

##  1 Introduction

A matching in a graph $G=(V,E)$ is a subset of edges $M\subseteq E$ where no
two edges in $M$ share an endpoint. A maximum matching of $G$ has the maximum
number of edges among all possible matchings. Here we let $m(G)$ denote the
matching size of $G$, i.e. the size of a maximum matching in $G$. In this
paper, we present algorithms for approximating $m(G)$ in the sublinear models
of computation. In particular, our results fit the vertex-arrival stream model
(also known as the adjacency list streams). In the vertex-arrival model, in
contrast with the edge-arrival version where the input stream is an arbitrary
ordering of the edges, here each item in the stream is a vertex of the graph
followed by a list of its neighbors.

We also focus on graphs with bounded arboricity. A graph $G=(V,E)$ has
arboricity bounded by $\alpha$ if the edge set $E$ can be partitioned into at
most $\alpha$ forests. A well-known fact (known as the Nash-William theorem
[NW64]) states that a graph has arboricity $\alpha$, if and only if every
induced subgraph on $t$ vertices has at most $\alpha(t-1)$ number of edges.
Graphs with low arboricity cover a wide range of graphs such as constant
degree graphs, planar graphs, and graphs with small tree-widths. In particular
planar graphs have arboricity bounded by $3$.

A simple reduction from counting distinct elements implies that computing
$m(G)$ exactly requires $\Omega(n)$ space complexity even for trees and
randomized algorithms (see [AMS99] for the lower bound on distinct elements
problem.) This has initiated the study of finding computationally-light
estimators for $m(G)$ that take small space to compute. With this focus,
following the work by Esfandiari et al. [EHL+15], there has been a series of
papers [MV16, CJMM17, MV18, BGM+19] that have designed estimators for the
matching size based on the degrees of vertices, edges and the arboricity of
the input graph. In this paper, we design another degree-based estimator for
$m(G)$ in low arboricity graphs that has certain advantages in comparison with
the previous works and leads to new algorithmic results. Before describing our
estimator we briefly review some of the previous ideas. In the discussions
below, we assume $G$ has arboricity bounded by $\alpha$.

#### Shallow edges, high degree vertices

Esfandiari et al. [EHL+15] were first to observe that one can approximately
characterize the matching size of low arboricity graphs based on the degree
information of the vertices and the local neighborhood od the edges. Let $H$
denote the set of vertices with degree more than $h=2\alpha+3$ and let $F$
denote the set of edges with both endpoints having degree at most $h$.
Esfandiari et al. have shown that $m(G)\leq|H|+|F|\leq(5\alpha+9)m(G).$ Based
on this estimator, the authors in [EHL+15] have designed a
$\tilde{O}(\varepsilon^{-2}\alpha n^{2/3})$ space algorithm for approximating
$m(G)$ within $5\alpha+9+\varepsilon$ factor in the edge-arrival model.

#### Fractional matchings

By establishing an interesting connection with fractional matchings and the
Edmonds Polytope theorem, Mcgregor and Vorotnikova [MV16] have shown the
following quantity approximates $m(G)$ within $(\alpha+2)$ factor.

$(\alpha+1)\sum_{(u,v)\in
E}\min\\{\frac{1}{\deg(u)},\frac{1}{deg(v)},\frac{1}{\alpha+1}\\}.$

Based on this estimator, the authors in [MV16] have given a
$\tilde{O}(\varepsilon^{-2}n^{2/3})$ space streaming algorithm (in the edge-
arrival model) that approximate $m(G)$ within $\alpha+2+\varepsilon$ factor.
Also in the same work, another degree-based estimator is given that returns a
$\frac{(\alpha+2)^{2}}{2}$ factor approximation of $m(G)$. A notable property
of this estimator is that it can be implemented in the vertex-arrival stream
model in $O(\log n)$ bits of space.

####  $\alpha$-Last edges

Cormode et al. [CJMM17] (later revised by Mcgregor and Vorotnikova [MV18])
have designed an estimator that depends on a given ordering of the edges.
Given a stream of edges $S=e_{1},\ldots,e_{m}$, let $E_{\alpha}(S)$ denote a
subset of edges where $(u,v)\in E_{\alpha}(S)$ iff the vertices $u$ and $v$
both appear at most $\alpha$ times in $S$ after the edge $(u,v)$. It is shown
that $m(G)\leq|E_{\alpha}(S)|\leq(\alpha+2)m(G).$ Moreover a
$O(\frac{1}{\varepsilon^{2}}\log^{2}n)$ space streaming algorithm is given
that approximates $|E_{\alpha}(S)|$ within $1+\varepsilon$ factor in the edge-
arrival model.

###  1.1 The estimator in this paper

The new estimator is purely based on the degree of the vertices in the graph
without any dependence on $\alpha$. To estimate the matching size, we count
the number of what we call locally superior vertices in the graph. Namely,

###### Definition 1

In graph $G=(V,E)$, we call $u\in V$ a locally superior vertex if $u$ has a
neighbor $v$ such that $deg(u)\geq deg(v)$. We let $\ell(G)$ denote the number
of locally superior vertices in $G$.

We show if the arboricity of $G$ is bounded by $\alpha$, then $\ell(G)$
approximates $m(G)$ within $(\alpha+2)$ factor (Lemma 2.) This repeats the
same bound obtained by the estimators in [MV16] and [CJMM17], however for
planar graphs, we show that the approximation factor is at most $3.5$ which
beats the previous bounds (Lemma 5) 333We do not have tight examples for our
analysis. In fact, we conjecture that $\ell(G)$ approximates $m(G)$ within $3$
factor when $G$ is planar.. As an evidence, consider the 4-regular planar
graph on $9$ vertices. Both of the estimators in [MV16] and [CJMM17], report
$18$ as the estimation for $m(G)$ while the exact answer is $4$. It follows
their approximation factor is at least $4.5$.

Unfortunately, the new estimator, in spite of its simplicity, does not seem to
be applicable in the edge-arrival model without an extra pass over the stream.
To decide if a vertex is locally superior, we need to know its neighbors and
learn their degrees which becomes burdonsome in one pass. However in the
vertex-partition model, we can obtain this information in one pass and
consequently can achieve sublinear space bounds. More formally, we get a
randomized $O(\frac{\sqrt{n}}{\varepsilon^{2}}\log n)$ space algorithm for
approximating $m(G)$ within $(3.5+\varepsilon)$ factor in this model. In terms
of approximation factor, this improves over existing sublinear algorithms
[MV16, MV18].

As another application of our estimator, we get a sublinear simultaneous
protocol in the vertex-partition model for approximating $m(G)$ when $G$ is
planar. In this model, vertex set $V$ is partitioned into $t$ subsets
$V_{1},\ldots,V_{t}$ where each subset is given to a player. The $i$-th player
knows the edges on $V_{i}$. The players do not communicate with each other.
They only send one message to a referee whom at the end computes an
approximation of the matching size. (The referee does not get any part of the
input.) We assume the referee and the players have a shared source of
randomness. Within this setting, we design a protocol that approximates $m(G)$
within $3.5+\varepsilon$ factor using $O(\frac{n^{2/3}}{\varepsilon^{2}}\log
n)$ communication from each player. Note that for $t>3$ and $t=o(n^{1/3})$,
this result is non-trivial. The best previous result implicit in the works of
[CCE+16, MV16] computes a $5+\varepsilon$ factor approximation using
$\tilde{O}(n^{4/5})$ communication from each player. We should also mention
that, based on the estimator in [MV16], there is a simultaneous protocol that
reports a $12.5$ factor approximation of $m(G)$ using $O(\log n)$
communication.

##  2 Graph properties

In the following proofs, we let $M\subseteq E$ denote a maximum matching in
graph $G$. When the underlying graph is clear from the context, for the vertex
set $S$, we use $N(S)$ to denote the neighbors of the vertices in $S$
excluding $S$ itself. For vertex $u$, we simply use $N(x)$ to denote the
neighbors of $u$. The vertex $v$ is a neighbor of the edge $(x,y)$ if $v$ is
adjacent with $x$ or $y$. When $x$ is paired with $y$ in the matching $M$,
abusing the notation, we define $M(x)=y$.

###### Lemma 2

Let $G=(V,E)$ be a graph with arboricity $\alpha$. We have

$m(G)\leq\ell(G)\leq(\alpha+2)m(G).$

Proof: The left hand side of the inequality is easy to show. For every edge in
$E$, at least one of the endpoints is locally superior. Since edges in $M$ are
disjoint, at least $|M|$ number of endpoints must be locally superior. This
proves $m(G)\leq\ell(G)$.

To show the right hand side, we use a charging argument. Let $L$ denote the
locally superior vertices in $G$. Our goal is to show an upper bound on $|L|$
in terms of $|M|$ and $\alpha$. Let $X\subseteq L$ be the set of locally
superior vertices that are NOT endpoints of a matching edge. The challenge is
to prove an upper on $|X|$.

The vertices in $X$ do not contribute to the maximum matching. However all the
vertices in $N(X)$ must be endpoints of matching edges (otherwise $M$ would
not be a maximal matching.) For the same reason, there cannot be an edge
between the vertices in $X$. To prove an upper bound on $|X|$, in the first
step, using an assignment procedure, we assign a subset of vertices in $X$ to
edges in $M$ in a way any target edge gets at most $\alpha-1$ locally superior
vertices. We do the assignments in the following way.

#### The Assignment Procedure

If we find a $y\in N(X)$ with at most $\alpha-1$ neighbors in $X$, we assign
all the neighbors of $y$ in $X$ to the matching edge $(y,M(y))$. We repeat
this process, every time picking a vertex in $N(X)$ with less than $\alpha$
neighbors in $X$ and do the assignment that we just described, until we cannot
find such a vertex in $N(X)$. Note that when we assign a locally superior
vertex $x$, we remove the edges on $x$ before continuing the procedure.

Here we emphasize the fact that if $y$ has a neighbor $x\in X$, then $M(y)$
cannot have neighbors in $X\setminus\\{x\\}$ (otherwise it would create an
augmenting path and contradict the optimality of $M$.)

Let $X_{1}\subseteq X$ be the assigned locally superior vertices and
$M_{1}\subseteq M$ be the used matching edges in the assignment procedure. We
have

$|X_{1}|\leq(\alpha-1)|M_{1}|.$ (1)

Let $X_{2}=X\setminus X_{1}$ be the unassigned vertices in $X$. Now we try to
prove an upper bound on $|X_{2}|$. For this, we need to make a few
observations.

###### Observation 3

Let $Y_{2}=N(X_{2})$. The pair $y$ and $M(y)$ cannot be both in $Y_{2}$.

Proof: Suppose $y$ and $M(y)$ are both in $N(X_{2})$. Let $B$ and $C$ be the
neighbors of $y$ and $M(y)$ in $X_{2}$ respectively. If $|B\cup C|>1$, then
one can find an augmenting path of length $3$ (with respect to $M$.) A
contradiction.

On the other hand, if $|B\cup C|=1$, then $y$ and $M(y)$ have only a shared
neighbor $x\in X_{2}$ which means the edge $e=(y,M(y))$ should have been used
by the assignment procedure and as result $x\in X_{1}$. Another contradiction.
$\Box$

###### Observation 4

Every vertex $x\in X_{2}$ has degree at least $\alpha+1$.

Proof: Consider $x\in X_{2}$. Suppose, for the sake of contradiction,
$\deg(x)$ is $k$ where $k\leq\alpha$. Since $x$ is a locally superior vertex,
there must be a $y\in N(x)$ with degree at most $k$ in $G$. We know that $y$
is an endpoint of a matching edge. In the assignments procedure, whenever we
used an edge $e\in M$ all the neighbors of its endpoints (in $X$) were
assigned. Since $x$ is not assigned yet, it means the edge $(y,M(y))$ has not
been used. Consequently $y$ must have at least $\alpha$ neighbors in $X_{2}$.
Counting the edge $(y,M(y))$, we should have $\deg(y)\geq\alpha+1$. A
contradiction. $\Box$

Let $G^{\prime}=(X_{2}\cup Y_{2},E^{\prime})$ be a bipartite graph where
$E^{\prime}$ is the set of edges between $X_{2}$ and $Y_{2}$. From Observation
4, we have

$(\alpha+1)|X_{2}|\leq|E^{\prime}|.$ (2)

Since $G^{\prime}$ is a subgraph of $G$, its arboricity is bounded by
$\alpha$. As result,

$|E^{\prime}|\leq\alpha(|X_{2}|+|Y_{2}|).$ (3)

Recall that $Y_{2}$ are endpoints of matching edges. Let $M_{2}$ be those
matching edges. Observation 3 implies that $|Y_{2}|=|M_{2}|$. As result,
combining (2) and (3), we get the following.

$|X_{2}|\leq\alpha|Y_{2}|=\alpha|M_{2}|.$ (4)

To prove an upper bound on $|L|$, we also need to count the locally superior
vertices that are endpoints of matching edges. Let $Z=L\setminus X$. We have
$|Z|\leq 2|M|.$ Summing up, we get

$\displaystyle|L|$ $\displaystyle=|X_{1}|+|X_{2}|+|Z|$
$\displaystyle\leq(\alpha-1)|M_{1}|+\alpha|M_{2}|+2|M|$
$\displaystyle=\alpha(|M_{1}|+|M_{2}|)+2|M|-|M_{1}|$
$\displaystyle\leq(\alpha+2)|M|-|M_{1}|$ $\displaystyle\leq(\alpha+2)|M|$

This proves the lemma. $\Box$

###### Lemma 5

Let $G=(V,E)$ be a planar graph. We have $\ell(G)\leq 3.5m(G)$.

Proof: For planar graphs, similar to what we did in the proof of Lemma 2, we
first try to assign some of the vertices in $X$ to the matching edges using a
simple assignment procedure. (Recall that $X$ is the set of vertices in $L$
that are not endpoints of edges in $M$.)

#### The Assignment Procedure

Let $Y_{1}=\emptyset$. If we find a $y\in N(X)$ with only $1$ neighbor $x\in
X$, we assign $x$ to the matching edge $(y,M(y))$. Also we add $y$ to $Y_{1}$.
We continue the procedure until we cannot find such a vertex in $N(X)$. Note
that when we assign a locally superior vertex $x$, we remove the edges on $x$.

Let $X_{1}\subseteq X$ be the assigned locally superior vertices and
$M_{1}\subseteq M$ be the used matching edges in the assignment procedure.
Note that $|Y_{1}|=|M_{1}|$. We have

$|X_{1}|\leq|M_{1}|.$ (5)

Let $X_{2}=X\setminus X_{1}$. Using a similar argument that we used for
proving Observation 4, we can show every vertex in $X_{2}$ has degree at least
$3$. Also letting $Y_{2}=N(X_{2})$, we observe that $y\in Y_{2}$ and $M(y)$
cannot be both in $Y_{2}$ as we noticed in the Observation 3. Let
$M_{2}\subseteq M$ be the matching edges with one endpoint in $Y_{2}$. We have
$|Y_{2}|=|M_{2}|$.

Now consider the bipartite graph $G^{\prime}=(X_{2}\cup Y_{2},E^{\prime})$
where $E^{\prime}$ is the set of edges between $X_{2}$ and $Y_{2}$. Every
planar bipartite graph with $n$ vertices has at most $2n-4$ edges 444For a
short proof of this, combine the Euler’s formula $|V|-|E|+|F|=2$ with the
inequality $2|E|\geq 4|F|$ caused by each face having at least 4 sides (since
there are no odd cycles) and we get $|E|\leq 2|V|-4.$. Since $G^{\prime}$ is a
bipartite planar graph, it follows,

$3|X_{2}|\leq|E^{\prime}|<2(|X_{2}|+|Y_{2}|)=2(|X_{2}|+|M_{2}|).$ (6)

This shows $|X_{2}|<2|M_{2}|$. Letting $Z=L\setminus X$ and
$M_{3}=M\setminus(M_{1}\cup M_{2})$, we get

$|L|=|X_{1}|+|X_{2}|+|Z|\leq|M_{1}|+2|M_{2}|+2|M|\leq 3|M|+|M_{2}|-|M_{3}|.$
(7)

This already proves $|L|$ is bounded by $4|M|$. To prove the bound claimed in
the lemma, we also show that $|L|\leq 3|M|+|M_{1}|+|M_{3}|$. Combined with the
inequality (7), this proves the lemma.

Let $Y=Y_{1}\cup Y_{2}$. Note that $Y$ are one side of the matching edges in
$M_{1}\cup M_{2}$. Let $Y^{\prime}=\\{M(y)\;|\;y\in Y\\}.$ We use a special
subset of $Y^{\prime}$, named $Y^{\prime\prime}$ which is defined as follows.
We let $Y^{\prime\prime}$ denote the locally superior vertices in $Y^{\prime}$
that have degree $2$ or they are adjacent with both endpoints of an edge in
$M_{3}$. We make the following observation regarding the vertices in
$Y^{\prime\prime}$.

###### Observation 6

We can assign each vertex $y^{\prime}\in Y^{\prime\prime}$ to a distinct $e\in
Y_{1}\cup M_{3}$ where $e$ has no neighbor in
$Y^{\prime}\setminus\\{y^{\prime}\\}$.

Proof: Consider $y^{\prime}\in Y^{\prime\prime}$. If $y^{\prime}$ is adjacent
with both endpoints of an edge $e=(z,z^{\prime})\in M_{3}$, we assign
$y^{\prime}$ to $e$ (when there are multiple edges with this condition we pick
one of them arbitrarily.) Note that $z$ and $z^{\prime}$ cannot have neighbors
in $Y^{\prime}$ other than $y^{\prime}$ because otherwise it would create an
augmenting path.

Now suppose $y^{\prime}$ has degree $2$. Since $y^{\prime}$ is a locally
superior vertex, it must have a neighbor $z$ of degree at most $2$. The
neighbor $z$ cannot be in $Y_{2}\cup X_{2}$ because the vertices in $Y_{2}\cup
X_{2}$ have degree at least $3$. We distinguish between two cases.

  * •

$M(y^{\prime})\in Y_{2}$. In this case, $z$ cannot be in $Y_{1}$ either
because the vertices in $Y_{1}$ are already of degree $2$ without
$y^{\prime}$. Also $z\notin X_{1}$ because otherwise it would create an
augmenting path. The only possibility is that $z$ is an endpoint of a matching
edge in $M_{3}$. We assign $y^{\prime}$ to the matching edge
$(z,z^{\prime})\in M_{3}$. Note that $z^{\prime}$ cannot have a neighbor in
$Y^{\prime}\setminus\\{y^{\prime}\\}$ because it would create an augmenting
path.

  * •

$M(y^{\prime})\in Y_{1}$. Here $z$ could be in $X_{1}$. If this is the case,
then $M(y^{\prime})$ cannot have a neighbor in
$Y^{\prime}\setminus\\{y^{\prime}\\}$ because it would create an augmenting
path. In this case, we assign $y^{\prime}$ to $M(y^{\prime})$. If
$z=M(y^{\prime})$, then again we assign $y^{\prime}$ to $M(y^{\prime})$. The
only remaining possibility is that $z$ an endpoint of a matching edge in
$M_{3}$ which we handle it similar to the previous case.

$\Box$

$Y_{2}$$Y_{1}$$Y^{\prime}$$X_{2}$$X_{1}$$M_{3}$ Figure 1:  A demonstration of
the construction in the proof of lemmas 2 and 5. Thick edges represent
matching edges. The unfilled vertices belong to the set $Y^{\prime\prime}$.

Now, assume we assign the vertices in $Y^{\prime\prime}$ to the elements in
$Y_{1}\cup M_{3}$ according to the above observation. Let
$Y_{1}^{\prime}\subseteq Y_{1}$ and $M_{3}^{\prime}\subseteq M_{3}$ be the
vertices and edges that were used in the assignment. Let
$Y^{\prime\prime\prime}$ be the remaining locally superior vertices in
$Y^{\prime}$. Namely, $Y^{\prime\prime\prime}=(L\cap Y^{\prime})\setminus
Y^{\prime\prime}$. Before making the final point, we observe that only one
endpoint of the edges in $M_{3}$ are adjacent with vertices in
$Y^{\prime\prime\prime}$. Let $Y_{3}$ be the endpoint of edges in
$M_{3}\setminus M_{3}^{\prime}$ that have neighbors in
$Y^{\prime\prime\prime}$. Consider the bipartite graph
$G^{\prime\prime}(V^{\prime\prime},E^{\prime\prime})$ where

$V^{\prime\prime}=(X_{2}\cup
Y^{\prime\prime\prime})\cup\big{(}Y_{2}\cup(Y_{1}\setminus Y_{1}^{\prime})\cup
Y_{3}\big{)}$

and $E^{\prime\prime}$ is the set of edges between $X_{2}$ and $Y_{2}$, and
the edges between $Y^{\prime\prime\prime}$ and $Y_{2}\cup(Y_{1}\setminus
Y_{1}^{\prime})\cup Y_{3}$.

Relying on the facts that $G^{\prime\prime}$ is a planar bipartite graph,
$Y^{\prime\prime\prime}$ is composed of vertices with degree at least $3$, and
the edges on $Y^{\prime\prime\prime}$ are all in $E^{\prime\prime}$, we have

$3|X_{2}|+3|Y^{\prime\prime\prime}|\leq|E^{\prime\prime}|\leq
2(|X_{2}|+|Y_{2}|+|Y_{1}\setminus
Y_{1}^{\prime}|+|Y^{\prime\prime\prime}|+|Y_{3}|).$

It follows,

$\displaystyle|X_{2}|+|Y^{\prime\prime\prime}|$ $\displaystyle\leq
2(|Y_{2}|+|Y_{1}\setminus Y_{1}^{\prime}|+|Y_{3}|)$ $\displaystyle\leq
2(|M_{2}|+|M_{1}|-|Y_{1}^{\prime}|+|M_{3}|-|M_{3}^{\prime}|)$
$\displaystyle=2(|M|-|Y_{1}^{\prime}|-|M_{3}^{\prime}|)$

Since $|Y^{\prime\prime}|=|Y^{\prime}_{1}|+|M^{\prime}_{3}|$, we get

$|X_{2}|+|Y^{\prime\prime\prime}|\leq 2|M|-2|Y^{\prime\prime}|$ (8)

Let $Z_{1}$, $Z_{2}$ and $Z_{3}$ denote the locally superior vertices that are
endpoints of matching edges in $M_{1}$, $M_{2}$ and $M_{3}$ respectively. From
the definition of $Y^{\prime\prime}$ and $Y^{\prime\prime\prime}$, we have

$|Z_{1}|+|Z_{2}|\leq|M_{1}|+|M_{2}|+|Y^{\prime\prime}|+|Y^{\prime\prime\prime}|$
(9)

From (8) and (9), we get

$\displaystyle|L|$ $\displaystyle=|X_{1}|+|X_{2}|+|Z_{1}|+|Z_{2}|+|Z_{3}|$
$\displaystyle\leq|M_{1}|+|X_{2}|+(|M_{1}|+|M_{2}|+|Y^{\prime\prime}|+|Y^{\prime\prime\prime}|)+2|M_{3}|$
$\displaystyle=2|M_{1}|+(|X_{2}|+|Y^{\prime\prime\prime}|)+|M_{2}|+|Y^{\prime\prime}|+2|M_{3}|$
$\displaystyle\leq 2|M_{1}|+|M_{2}|+2|M|-|Y^{\prime\prime}|+2|M_{3}|$
$\displaystyle=3|M|+|M_{1}|+|M_{3}|-|Y^{\prime\prime}|$ $\displaystyle\leq
3|M|+|M_{1}|+|M_{3}|$

This finishes the proof of the lemma. $\Box$

##  3 Algorithms

We first present a high-level sampling-based estimator for $\ell(G)$. Then we
show how this estimator can be implemented in the streaming and distributed
settings using small space and communication. For our streaming result, we use
a combination of the estimator for $\ell(G)$ and the greedy maximal matching
algorithm. For the simultaneous protocol, we use the estimator for $\ell(G)$
in combination with the edge-sampling primitive in [CCE+16] and an estimator
in [MV16].

The high-level estimator (described in Algorithm 1) samples a subset of
vertices $S\subseteq V$ and computes the locally superior vertices in $S$. The
quantity $\ell(G)$ is estimated from the scaled ratio of the locally superior
vertices in the sample set.

Run the following estimator $r=\lceil\frac{8}{\epsilon^{2}}\rceil$ number of
times in parallel. In the end, report the average of the outcomes.

1\. Sample $s$ vertices (uniformly at random) from $V$ without replacement.

2\. Let $S$ be the set of sampled vertices.

3\. Compute $S^{\prime}$ where $S^{\prime}$ is the set of locally superior
vertices in $S$.

4\. Return $\frac{n}{s}|S^{\prime}|$ as an estimation for $\ell(G)$.

Algorithm 1 The high-level description of the estimator for $\ell(G)$

###### Lemma 7

Assuming $s\geq\frac{n}{\ell(G)}$, the high-level estimator in Algoirthm 1
returns a $1+\varepsilon$ factor approximation of $\ell(G)$ with probability
at least $7/8$.

Proof: Fix a parallel repetition of the algorithm and let $X$ denote the
outcome of the associated estimator. Assuming an arbitrary ordering on the
locally superior vertices, let $X_{i}$ denote the random variable associated
with $i$-th locally superior vertex. We define $X_{i}=1$ if the $i$-th locally
superior vertex has been sampled, otherwise $X_{i}=0$. We have
$X=\frac{n}{s}\sum_{i=1}^{\ell(G)}X_{i}$. Since $Pr(X_{i}=1)=\frac{s}{n}$, we
get $E[X]=\ell(G)$. Further we have

$\displaystyle
E[X^{2}]=\frac{n^{2}}{s^{2}}E\Big{[}\sum_{i,j}^{\ell(G)}X_{i}X_{j}\Big{]}$
$\displaystyle=\frac{n^{2}}{s^{2}}\Big{[}\sum_{i}^{\ell(G)}E[X_{i}^{2}]+\sum_{i\neq
j}^{\ell(G)}E[X_{i}X_{j}]\Big{]}$
$\displaystyle=\frac{n^{2}}{s^{2}}\Big{[}\frac{s}{n}\ell(G)+{\ell(G)\choose
2}\frac{s(s-1)}{n(n-1)}\Big{]}$
$\displaystyle=\frac{n}{s}\ell(G)+{\ell(G)\choose 2}\frac{n(s-1)}{s(n-1)}$
$\displaystyle<\frac{n}{s}\ell(G)+\ell^{2}(G)$

Consequently, $Var[X]=E[X^{2}]-E^{2}[X]<\frac{n}{s}\ell(G).$

Let $Y$ be the average of the outcomes of $r$ parallel and independent
repetitions of the basic estimator. We have $E[Y]=\ell(G)$ and
$Var[Y]<\frac{n}{sr}\ell(G)$. Using the Chebyshev’s inequality,

$Pr(|Y-E[Y]|\geq\varepsilon
E[Y])\leq\frac{Var[Y]}{\varepsilon^{2}E^{2}[X]}<\frac{n/s}{r\varepsilon^{2}\ell(G)}.$

Setting $r=\frac{8}{\varepsilon^{2}}$ and $s\geq\frac{n}{\ell(G)}$, the above
probability will be less than $1/8$. $\Box$

###  3.1 The streaming algorithm

We first note that we can implement the high-level estimator of Algorithm 1 in
the vertex-arrival stream model using $O(\frac{s}{\varepsilon^{2}}\log n)$
space. Consider a single repetition of the estimator. The sampled set $S$ is
selected in the beginning of the algorithm (before the stream.) This can be
done using a reservoir sampling strategy [Vit85] in $O(|S|\log n)$ space. To
decide if $u\in S$ is locally superior or not, we just need to store $\deg(u)$
and the minimum degree of the neighbors that are visited so far. This takes
$O(\log n)$ bits of space. As result, the whole space needed to implement a
single repetition is $O(s\log n)$ bits.

The streaming algorithm runs two threads in parallel. In one thread it runs
the streaming implementation of Algorithm 1 after setting
$s=\lceil\sqrt{n}\>\rceil$. In the other thread, it runs a greedy algorithm to
find a maximal matching in the input graph. We stop the greedy algorithm
whenever the size of the discovered matching $F$ exceeds $\sqrt{n}$. In the
end, if $|F|<\sqrt{n}$, we output $|F|$ as an approximation for $m(G)$,
otherwise we report the outcome of the first thread.

Note that if $|F|<\sqrt{n}$, $F$ is a maximal matching in $G$. Hence
$|F|\geq\frac{1}{2}m(G)$. Assume $|F|\geq\sqrt{n}$. In this case the algorithm
outputs the result of first thread. In this case, by Lemma 2, we know
$\ell(G)\geq\sqrt{n}$. Consequently, it follows from Lemma 7, the first thread
returns a $1+O\varepsilon$ approximation of $\ell(G)$ and consequently it
returns a $3.5+O(\varepsilon)$ approximation of $m(G)$. Since the greedy
algorithm takes at most $O(\sqrt{n})$ space, the space complexity of the
algorithm is dominated by the space usage of the first thread. We get the
following result.

###### Theorem 8

Let $G$ be a planar graph. There is a streaming algorithm (in the vertex-
arrival model) that returns a $3.5+\epsilon$ factor approximation of $m(G)$
using $O(\frac{\sqrt{n}}{\epsilon^{2}})$ space.

###  3.2 A simultaneous communication protocol

To describe the simultaneous protocol, we consider two cases separately: (a)
when the matching size is low; to be precise, when it is smaller than some
fixed value $k=n^{1/3}$, and (b) when the matching size is high, i.e. at least
$\Omega(k)$. For each case, we describe a separate solution. The overall
protocol will be these solutions (run in parallel) combined with a sub-
protocol (in parallel) to distinguish between the cases.

#### Graphs with large matching size

In the case when matching size is large, similar to what was done in the
streaming model, we run an implementation of Algorithm 1 in the simultaneous
model. To see the implementation, in the simultaneous model all the players
(including the referee) know the sampled set $S$. This results from access to
the shared randomness. For each $u\in S$, the players send the minimum degree
of the neighbors of $u$ in his input to the referee. The player that owns $u$,
also sends $\deg(u)$ to the referee. Having received this information, the
referee can decide if $u$ is a locally superior vertex or not. As result, we
can implement Algorithm 1 in the simultaneous model using a protocol with
$O(\frac{s}{\varepsilon^{2}}\log n)$ message size.

#### Graphs with small matching size

In the case where the matching size is small, we use the edge-sampling method
of [CCE+16]. We review their basic sampling primitive in its general form.
Given a graph $G(V,E)$, let $c:V\rightarrow[b]$ be a totally random function
that assigns each vertex in $V$ a random number (color) in
$[b]=\\{1,\ldots,b\\}$. The set $\textrm{Sample}_{b,d,1}$ is a random subset
of $E$ picked in the following way. Given a subset $K\subseteq[b]$ of size
$d\in\\{1,2\\}$, let $E_{K}$ be the edges of $G$ where the color of their
endpoints matches $K$. For example when $K=\\{3,4\\}$, the set $E_{\\{3,4\\}}$
contains all edges $(u,v)$ such that $\\{c(u),c(v)\\}=\\{3,4\\}$. For all
$K\subseteq[b]$ of size $d$, the set $\textrm{Sample}_{b,d,1}$ picks a random
edge from $E_{K}$. Finally, the random set $\textrm{Sample}_{b,d,r}$ is the
union of $r$ independent instances of $\textrm{Sample}_{b,d,1}$. We have the
following lemma from [CCE+16] (see Theorems 4 in the reference.)

###### Lemma 9

Let $G=(V,E)$ be a graph. When $m(G)\leq k$, with probability
$1-1/{\mathrm{poly}}(k)$, the random set $\textrm{Sample}_{100k,2,O(\log k)}$
contains a matching of size $m(G)$.

Note that, in the simultaneous vertex-partition model, the referee can obtain
an instance of $\textrm{Sample}_{b,d,1}$ via a protocol with $O(b^{d}\log n)$
message size. To see this, using the shared randomness, the players pick the
random function $c:V\rightarrow[b]$. Let $E^{(i)}$ be the subset of edges
owned by the $i$-th player. We have $E=\bigcup_{i=1}^{t}E^{(i)}$. To pick a
random edge from $E_{K}$ for a given $K\subseteq[b]$, the $i$-th player
randomly picks an edge $e\in E_{K}\cap E^{(i)}$ and sends it along with
$|E_{K}\cap E^{(i)}|$ to the referee. After receiving this information from
all the players, the referee can generate a random element of $E_{K}$. Since
there are $O(b^{d})$ different $d$-subsets of $[b]$, the size of the message
from a player to the referee is bounded by $O(b^{d}\log n)$ bits.
Consequently, the referee can produce a rightful instance of
$\textrm{Sample}_{b,d,r}$ using $O(rb^{d}\log n)$ communication from each
player.

#### How to distinguish between the cases?

To accomplish this, here we use a degree-based estimator by Mcgregor and
Vorotnikova [MV16] described in the following lemma.

###### Lemma 10

Let $G$ be a planar graph. We have

$m(G)\leq A^{\prime}(G)=\sum_{u\in V}\min\\{\deg(u)/2,4-\deg(u)/2\\}\leq
12.5\>m(G).$

It is easy to see that, in the simultaneous vertex-partition model, we can
implement this estimator with $O(\log n)$ bits communication from each player.

#### The final protocol

Let $k=\lceil n^{1/3}\rceil$. We run the following threads in parallel.

  1. (1)

A protocol that implements the high-level estimator (Algorithm 1) with
$s=\lceil 12.5n/k\rceil$ as its input parameter according to the discussions
above. Let $z_{1}$ be the output of this protocol.

  2. (2)

A protocol to compute an instance of $\textrm{Sample}_{b,d,r}$ for $b=100k$
and $d=2$ and $r=O(\log k)$. Let $z_{2}$ be the size of maximum matching in
the sampled set.

  3. (3)

A protocol to compute $A^{\prime}(G)$. Let $z_{3}$ be the output of this
thread.

In the end, if $z_{3}\geq\frac{k}{12.5}$, the referee outputs $z_{1}$ as an
approximation for $m(G)$, otherwise the referee reports $z_{2}$ as the final
answer.

###### Theorem 11

Let $G$ be a planar graph on $n$ vertices. The above simultaneous protocol
with probability $3/4$ returns a $3.5+O(\varepsilon)$ approximation of $m(G)$
where each player sends $O(\frac{n^{2/3}}{\varepsilon^{2}})$ bits to the
referee.

Proof: First we note that by choosing the constants large enough, we can
assume the thread (2) errs with probability at most $1/8$. If
$z_{3}\geq\frac{k}{12.5}$, then we know $m(G)\geq\frac{k}{12.5}$. This follows
from Lemma 10. Consequently by Lemma 2, we have $\ell(G)\geq\frac{k}{12.5}$.
Therefore from Lemma 7, we have $|z_{1}-\ell(G)|\leq\varepsilon\ell(G)$ with
probability at least $7/8$. It follows from Lemma 5 that
$(1-\varepsilon)m(G)\leq z_{1}\leq(3.5+3.5\varepsilon)m(G).$

On the other hand, if $z_{3}<\frac{k}{12.5}$, by Lemma 10 we know that $m(G)$
must be less than $k$. Having this, from Lemma 9, with probability at least
$7/8$, we get $z_{2}=m(G)$. In this case the protocol computes the exact
matching size of the graph.

The communication complexity each player is dominated by the cost of the first
thread which is $O(n^{2/3}\varepsilon^{-2}\log n)$. The total error
probability is bounded by $1/4$. This finishes the proof. $\Box$

## References

  * [AMS99] Noga Alon, Yossi Matias, and Mario Szegedy, _The space complexity of approximating the frequency moments_ , Journal of Computer and system sciences 58 (1999), no. 1, 137–147.
  * [BGM+19] Marc Bury, Elena Grigorescu, Andrew McGregor, Morteza Monemizadeh, Chris Schwiegelshohn, Sofya Vorotnikova, and Samson Zhou, _Structural results on matching estimation with applications to streaming_ , Algorithmica 81 (2019), no. 1, 367–392.
  * [CCE+16] R. Chitnis, G. Cormode, H. Esfandiari, M.T. Hajiaghayi, A. McGregor, M. Monemizadeh, and S. Vorotnikova, _Kernelization via sampling with applications to finding matchings and related problems in dynamic graph streams_ , Proceedings of the 27th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 2016, pp. 1326–1344.
  * [CJMM17] Graham Cormode, Hossein Jowhari, Morteza Monemizadeh, and S. Muthukrishnan, _The sparse awakens: Streaming algorithms for matching size estimation in sparse graphs_ , 25th Annual European Symposium on Algorithms, ESA 2017, September 4-6, 2017, Vienna, Austria, 2017, pp. 29:1–29:15.
  * [EHL+15] H. Esfandiari, M.T. Hajiaghyi, V. Liaghat, M. Monemizadeh, and K. Onak, _Streaming algorithms for estimating the matching size in planar graphs and beyond_ , Proceedings of the 26th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 2015.
  * [MV16] A. McGregor and S. Vorotnikova, _Planar matching in streams revisited_ , Proceedings of the 19th International Workshop on Approximation Algorithms for Combinatorial Optimization Problems (APPROX), 2016.
  * [MV18] Andrew McGregor and Sofya Vorotnikova, _A simple, space-efficient, streaming algorithm for matchings in low arboricity graphs_ , 1st Symposium on Simplicity in Algorithms, SOSA 2018, January 7-10, 2018, New Orleans, LA, USA, 2018, pp. 14:1–14:4.
  * [NW64] C. St. J. A. Nash-Williams, _Decomposition of finite graphs into forests_ , Journal of the London Mathematical Society 39 (1964), no. 1, 12.
  * [Vit85] Jeffrey Scott Vitter, _Random sampling with a reservoir_ , ACM Trans. Math. Softw. 11 (1985), no. 1, 37–57.