\section{Introduction}
The first neutron star was discovered by Bell and Hewish in 1967. Since then, the neutron star has spawned the topics of research in numerous studies in astronomy and nuclear physics. A neutron star can be treated as a huge nucleus composed of neutrons, protons, and leptons~\cite{glendenning97,weber99}. The properties of neutron stars, like mass and radii, can be precisely measured with the gradual development of astronomical observation technology. This development has promoted theoretical studies on the nuclear many-body theory~\cite{prakash97,heiselberg00}, especially on the properties of nuclear matter at high densities, where high density refers to densities above 3 times the nuclear saturation density in the core of a neutron star.

In the Tolman-Oppenheimer-Volkoff (TOV) equation~\cite{oppenheimer39,tolman39}, which is derived by solving the Einstein equations for a spherically symmetric and time invariant system, the properties of a neutron star are fixed by the equation of state (EOS) of neutron star matter, which is a charge neutral system in $\beta$ equilibrium. The EOSs of neutron star matter have been widely studied using various nuclear many-body theories, by the application of density functional theory (DFT) and {\it ab initio} methods in nuclear physics~\cite{dutra12,dutra14}.

The effective nucleon-nucleon ($NN$) interactions are adopted in DFT, which are determined by fitting the empirical saturation properties of nuclear matter and the ground states of finite nuclei. There are usually two schemes used in DFT, i.e., relativistic and non-relativistic approximation. The most popular nuclear many-body theory in non-relativistic DFT is the Skyrme-Hartree-Fock theory~\cite{vautherin72,bender03,stone07}, which was constructed based on the point-coupling $NN$ interaction. In the relativistic version, covariant DFT was proposed based on the one-boson-exchange picture~\cite{serot86,ring96,meng06}. The EOSs of DFT around the saturation density region can be constrained very well and are kept consistent within different effective $NN$ interactions. However, their behaviors in the high density region, especially in the cores of neutron stars, show variability, which generate widely different predictions of the properties of neutron stars. With the discoveries of massive neutron stars, whose masses are around $2M_\odot$, many DFT interactions are eliminated~\cite{demorest10,antoniadis13}.

On the other hand, in the state-of-the-art nuclear many-body methods, via {\it{ab initio}} calculations, a realistic $NN$ interaction is used that is obtained by reproducing the $NN$ scattering data. The most popular realistic $NN$ potentials are constructed by the meson-exchange picture~\cite{machleidt89} and chiral effective field theory~\cite{epelbaum09,machleidt11}. The $NN$ potentials, based on chiral low-momentum expansions, are soft cores at short distances due to pion exchange; whereas, in the meson-exchange potential, such as the Bonn potential, the $\omega$ and $\rho$ meson exchanges are essential in the core region of the $NN$ interaction, which generate very strong repulsive core. There are also other high-momentum effects at short distance, such as the strong tensor force, besides strong repulsion~\cite{jastrow55} which is related with the nucleon inner structure~\cite{toki80,oka81} in the realistic $NN$ interaction. Since 1950, many microscopic many-body methods have been proposed to take such high momentum contributions into account. Akmal {\it et al.} used the hyper-netted chain-summation technique in the non-relativistic variational method~\cite{pandharipande79} and the AV18 potential~\cite{wiringa95} to calculate the EOS of neutron rich matter~\cite{akmal98} with a three-body force. Krastev and Sammarruca adopted the relativistic Brueckner-Hartree-Fock (RBHF)~\cite{brockmann90} theory to study the properties of neutron stars with Bonn potentials~\cite{krastev06}. Gandolfi {\it et al.} also compared their results with the calculations made by Akmal using the auxiliary field diffusion Monte Carlo method~\cite{gandolfi09}. Hebeler {\it et al.} discussed the properties of neutron star with $NN$ interactions from chiral effective filed theory and constrained the neutron star mass-radius relation in the framework of renormalization group theory~\cite{hebeler10,hebeler13}. All of these values of EOSs obtained from {\it{ab initio}} calculations agree well with each other in the high density region and generate similar properties of massive neutron stars.

For neutron rich matter, the tensor effect in realistic $NN$ interaction is suppressed largely in the $T=1$ isospin channel \cite{krastev06,hu10,wang12,hu13}. Therefore, it is enough to explicitly consider only the short-range correlation in neutron rich matter. Based on this motivation, a relativistic central variational (RCV) method was proposed by including a central correlation Jastrow function~\cite{hu11}, which was inspired by the series works by Panda {\it{et al.}}~\cite{panda05,panda06,panda07}. In the RCV method, the solutions of the relativistic Hartree-Fock (RHF) method~\cite{bouyssy87} were chosen as the trial wave function of the total system. The saturation properties obtained by the RCV method are $20\%$ different from the empirical data, likely because of the lack of tensor correlation in symmetric nuclear matter. However, the EOS of pure neutron matter within the Bonn potentials can reproduce the results from RBHF theory very well. The solutions of the RHF method in nuclear matter are constructed by spinors with plane wave functions, the summation of spin and isospin for which can be easily and systematically evaluated via the Feynman trace technology. Therefore, the calculation in the RCV method for neutron-rich matter is essential and economical, when compared with other time-consuming {\it{ab initio}} methods such as non-relativistic variational methods and the RBHF theory in asymmetric nuclear matter.

Furthermore, the calculation procedures of the {\it{ab initio}} aforementioned methods, when they were applied to neutron star matter in $\beta$ equilibrium, are complicated, therefore, the equations of $\beta$ equilibrium are usually solved with various approximations instead of self-consistence treatments~\cite{akmal98,krastev06}. Hence, the intention of this work is to apply the RVC method to calculate the EOSs of neutron star matter in $\beta$ equilibrium self-consistently and use these EOSs to study the properties of neutron stars with realistic $NN$ interaction, Bonn potentials. In addition, we compare our results with the calculations of other {\it{ab initio}} methods and discuss the role of central correlation in neutron stars. In Sect. 2, we give a basic theoretical formulation of the RCV method. In Sect. 3, numerical results are presented on the properties of neutron stars and are compared with other calculations. In Sect. 4, we give a summary of the present work.

\section{Relativistic central variational method}
We firstly show the Hamiltonian of the Bonn potential~\cite{machleidt89}, which is constructed based on the one-boson-exchange potential. It is defined as a sum of one-particle amplitudes of six bosons, consisting of $\pi$ and $\eta$ pseudovectors, $\sigma$ and $\delta$ scalars, and $\rho$ and $\omega$ vector mesons,
\beq
H&=&\sum^A_{i}T_i+\frac{1}{2}\sum^A_{i,j}V_{i,j}\nn
&=&\int d^3\mathbf{x}\overline \psi(\mathbf x)(-i\bm{\gamma}\cdot\bm{\nabla}+M_N)\psi(\mathbf x)\nn
&&+\frac{1}{2}\sum_{i=\sigma,\delta,\atop\eta,\pi,\omega,\rho}\int d^3\mathbf{x}'d^3\mathbf{x}\overline \psi(\mathbf x')\overline \psi(\mathbf x)\frac{\Gamma_i(1,2)}{m^2_i+\mathbf q^2}\psi(\mathbf x)\psi(\mathbf x'),
\eeq
where $\psi$ corresponds the nucleon field. $\mathbf q$ is the transfer momentum between two nucleons and $M_N,~m_i$ are the masses of nucleons and mesons, respectively. The $\Gamma_i$ matrices represent the vertex between nucleons and mesons, which are depicted below:
\beq
&&\Gamma_\sigma(1,2)=-g^2_\sigma,\nn
&&\Gamma_\delta(1,2)=-g^2_\delta\bm{\tau}_1\cdot\bm{\tau}_2,\nn
&&\Gamma_\eta(1,2)=-\left(\frac{f_\eta}{m_\eta}\right)^2(\rlap/q\gamma_5)_1(\rlap/q\gamma_5)_2,\nn
&&\Gamma_\pi(1,2)=-\left(\frac{f_\pi}{m_\pi}\right)^2(\rlap/q\gamma_5)_1(\rlap/q\gamma_5)_2\bm{\tau}_1\cdot\bm{\tau}_2,\nn
&&\Gamma_\omega(1,2)=g^2_\omega\gamma_\mu(1)\gamma^\mu(2),\nn
&&\Gamma^V_\rho(1,2)=g^2_\rho\gamma_\mu(1)\gamma^\mu(2)\bm{\tau}_1\cdot\bm{\tau}_2,\nn
&&\Gamma^T_\rho(1,2)=\left(\frac{f_\rho}{2M_N}\right)^2q_\nu\sigma^{\mu\nu}(1)q^\lambda\sigma_{\mu\lambda}(2)\bm{\tau}_1\cdot\bm{\tau}_2,\nn
&&\Gamma^{VT}_\rho(1,2)=i\left(\frac{g_\rho f_\rho}{M_N}\right)\gamma_\mu(2)\sigma^{\mu\nu}q_\nu(1)\bm{\tau}_1\cdot\bm{\tau}_2,
\eeq
where $\sigma_{\mu\nu}=\frac{i}{2}[\gamma_\mu,\gamma_\nu]$ is an antisymmetric tensor gamma matrix, and the tensor coupling part between $\omega$ meson and nucleon has been neglected because the value of $f_\omega/g_\omega$ is negligible. Meanwhile, a monopole form factor should be considered,
\beq
F_i(q^2)=\frac{\Lambda^2_i-m^2_i}{\Lambda^2_i+q^2},
\eeq
for each meson-nucleon vertex denoted by $i$. All coupling constants and cut-off momenta $\Lambda_i$ were determined by fitting the $NN$ scattering data.

In the RCV method, one introduces a central correlation function of the wave function of the RHF theory as the trial wave function of the nuclear matter system \cite{hu11}, to treat the strong short-range repulsion,
\beq
|\Psi\ket=F|\Phi\ket,
\eeq
where $|\Phi\ket$ is the RHF wave function, and the correlation factor $F$ is chosen to be a product of two-body correlation functions $f(r_{ij})$,
\beq
F=\prod^A_{i<j}f(r_{ij}).
\eeq
Here, $f(r_{ij})$ is Jastrow correlation function~\cite{jastrow55}. The total energy density with the correlation function is obtained as,
\beq
\mathcal{E}_c&=&\frac{E_c}{\Omega}=\frac{1}{\Omega}\bra\Phi|\widetilde{H}|\Phi\ket\nn
&=&\bra T\ket+\bra T_c\ket+\bra V\ket.
\eeq
where the explicit form of the correlated Hamiltonian $\widetilde{H}$ appears as
\beq
\widetilde{H}&=&\sum^A_iT_i+\frac{1}{2}\sum^A_{i,j}\widetilde {V}_{ij}\nn
&=&\sum^A_iT_i+\frac{1}{2}\sum^A_{i,j}\{f^\dag(r_{ij})[T_i+T_j+V_{ij}]f(r_{ij})-(T_i+T_j)\}.
\eeq
Furthermore, the kinetic energy density part is not only related with the original one-body kinetic operator $\bra T\ket$,
\beq
\bra T\ket=\frac{\lambda}{\pi^2}\int^{k_F}_0p^2dp[p\hat P+M_N\hat M],
\eeq
but also the two-body operator with the Jastrow correlation function, $\bra T_c\ket$,
\beq
\bra T_c\ket&=&C\rho_B\frac{\lambda}{\pi^2}\int^{k_F}_0p^2dp[p\hat P+M_N\hat M]\nn
&-&\frac{2\lambda}{(2\pi)^4}\int^{k_F}_0p^2dpp'^2dp'\{[p\hat P(p)+2M_N \hat M(p)]I(p,p')+p'\hat P(p)J(p',p)\},
\eeq
where the first and second terms correspond to the direct (Hartree) and exchange (Fock) contributions of $\bra T_c\ket$, respectively and $C=\int d^3\mathbf r[f^2(r)-1]$. The isospin degeneracy is $\lambda=2$ for symmetric nuclear matter and $\lambda=1$ for pure neutron matter.  $\hat P$ and $\hat M$ are defined as
\beq
\hat P(p)=\frac{p^*(p)}{E^*(p)},~~~~~~~~~~~~~~\hat M(p)=\frac{M^*_N(p)}{E^*(p)}.
\eeq
The potential energy density $\bra V\ket$ will be split into the direct term $\bra V_D\ket$ and exchange term $\bra V_E\ket$. The pseudo-vector mesons do not provide their contributions in the Hartree approximation. We have the result of $\bra V_D\ket$ as,
\beq
\bra V_D\ket&=&-\frac{\widetilde {F}_\sigma(0,0)}{2}(\rho^p_S+\rho^p_S)^2-\frac{\widetilde {F}_\delta(0,0)}{2}(\rho^p_S-\rho^n_S)^2\nn
&&+\frac{\widetilde {F}_\omega(0,0)}{2}(\rho^p_B+\rho^n_B)^2+\frac{\widetilde {F}_\rho(0,0)}{2}(\rho^p_B-\rho^n_B)^2,
\eeq
where $\rho_S$ and $\rho_B$ are the scalar and nucleon vector densities, respectively, given as,
\beq
\rho^{p(n)}_S&=&\frac{1}{\pi^2}\int^{k^{p(n)}_F}_0p^2dp\hat M(p),\nn
\rho^{p(n)}_B&=&\frac{1}{\pi^2}\int^{k^{p(n)}_F}_0p^2dp.
\eeq
The proton and neutron cases are distinguished by the superscript $p/n$. The exchange contribution of $\sigma$ meson $\bra V^\sigma_E\ket$ as an example is given by,
\beq
\bra V^\sigma_E\ket=\frac{\lambda}{2(2\pi)^4}\int^{k_F}_0pdpp'dp\left[A_\sigma(p,p')+\hat M(p)\hat M(p')B_\sigma(p,p')+\hat P(p)\hat P(p')C_\sigma(p,p')\right].
\eeq
Here, the moment dependent functions $A_i, B_i, C_i, I, J$ and $F_i$ are various angular integrals related with the meson-exchange potentials, which are listed in the appendix of Ref.~\cite{hu11}. The contributions of other mesons can be expressed in similar forms. Finally, the total energy density is written as,
\beq
\mathcal{E}_c=\bra T\ket+\bra T_c\ket+\bra V_D\ket+\sum_i\bra V^i_E\ket.
\eeq
Usually, several free parameters, $c_1, c_2,\dots,c_i$, will appear in the Jastrow function $f(r)$. In the present work, the correlation function is chosen as,
\beq
f(r)=1-(c_0+c_1r+c_2r^2+c_3r^3)e^{-c_4r}.
\eeq
where the exponential term makes $f(r)$ unity at large distance. A natural choice from the unitary property of the correlation function is a normalization constraint on $f(r_{ij})$,
\beq\label{nc}
\int d^3 r_{ij}[f^2(r_{ij})-1]=0,
\eeq
and should go to zero for small $r_{ij}$ because of the repulsive core of $NN$ interaction, which lead to $c_0=1$. Furthermore, we should also ensure it is a monotonously increasing property at short distance,
\beq\label{gcon}
f'(0)\geq 0.
\eeq
We can then calculate the binding energy of nuclear matter after determining the remaining parameters with the variational principle. The minimal value of the total energy should appear at $f'(0)=0$ and $f''(0)=0$ with the constraint (\ref{gcon}) to make the Jastrow correlation function increase at short distance. Therefore, we can obtain the relations between $c_1, c_2$, and $c_4$,
\beq
c_1=c_4,~~~~c_2=\frac{c^2_4}{2}.
\eeq
Now, there is only one parameter, $c_4$, in the actual calculation, because $c_3$ is fixed by the normalization condition of the Jastrow correlation function, Eq. (\ref{nc}). We would like to determine $c_4$ by the variational principle with the energy density,
\beq
\frac{\partial\mathcal{E}_c}{\partial c_4}=0.
\eeq
More detailed formulas of the RCV method can be found in Ref.~\cite{hu11}.

For neutron star matter, there are not only the nucleons, but also the leptons, like electrons and muons. All of them exist in the neutron star with the equilibrium conditions of the chemical potentials for the $\beta$ decay,
\beq\label{bteq}
\mu_n&=&\mu_p+\mu_e,\nn
\mu_\mu&=&\mu_e,
\eeq
where the chemical potentials $\mu_n,~\mu_p,~\mu_\mu$, and $\mu_e$ are determined by the relativistic energy-momentum relation at the Fermi momentum $p=k_F$,
\beq
\mu_i&=&\Sigma^i_0(k_F)+E^\ast_i(k_F),\nn
\mu_\lambda&=&\sqrt{k^2_{F,\lambda}+m^2_\lambda},
\eeq
where, $i=n,~p$ and $\lambda=e,~\mu$. $\Sigma^i_0(k_F)$ is the zero component of the self-energy of proton or neutron. Furthermore, the nucleon density conservation and charge neutrality are imposed in neutron star matter as,
\beq
&&\rho=\rho_{n}+\rho_{p},\nn
&&\rho_e+\rho_\mu=\rho_{p}.
\eeq
The pressure of the neutron star system can be obtained with thermodynamics relation, as
\beq
P(\rho)=\rho^2\frac{\partial}{\partial\rho}\frac{\varepsilon}{\rho}=\sum_{i=n,p,e,\mu}\rho_{i}\mu_i-\varepsilon.
\eeq

The stable configurations of a neutron star then can be obtained from the well known hydrostatic equilibrium equations, by Tolman, Oppenheimer and Volkoff~\cite{tolman39,oppenheimer39} for the pressure $P$ and the enclosed mass $m$,
\beq\label{tov}
\frac{dP(r)}{dr}&=&-\frac{Gm(r)\varepsilon(r)}{r^{2}}\frac{\Big[1+\frac{P(r)}{\varepsilon(r)}\Big]\Big[1+\frac{4\pi r^{3}P(r)}{m(r)}\Big]}
{1-\frac{2Gm(r)}{r}},\nn
\frac{dm(r)}{dr}&=&4\pi r^{2}\varepsilon(r),
\eeq
where, $P(r)$ is the pressure of neutron star at radius, $r$, and $m(r)$ is the total star mass inside a sphere of radius $r$. Once the EOS $P(\varepsilon)$ is specified, $\varepsilon$ being the total energy density ($G$ is the gravitational constant), for a chosen central value of the energy density, the numerical integration of Eq.(\ref{tov}) provides the mass-radius relation of neutron star.

\section{Results and discussions}
\subsection{The equations of state of neutron star matter}
After solving the equilibrium conditions of chemical potentials for various particles in the neutron star matter, Eq.~(\ref{bteq}), the binding energies per nucleon as functions of density are shown in the left panel of Fig.~\ref{re} with three Bonn potentials, Bonn A, Bonn B, and Bonn C, which were fitted by the scattering data of $NN$ system. The binding energy of the Bonn A potential is the smallest of the three Bonn potentials, which is in accordance with the conclusion of RBHF theory for symmetric nuclear matter~\cite{brockmann90}. The tensor component of Bonn A among these three potentials is the weakest, whose $D$ state probability, $P_D$, is the smallest in deuteron. Furthermore, the difference of binding energies of these three potentials in neutron star matter is obviously less than that in the symmetric nuclear matter. The effect of the tensor force becomes weaker with increasing neutron fraction and does not play any role in pure neutron matter. Similarly, the pressure as a function of energy density for neutron star matter is required as input data, when we calculate the properties of neutron stars in the TOV equation. The pressure-energy relations with Bonn potentials are given in the right panel of Fig.~\ref{re}. Their behaviors are very similar to those of the binding energy and are almost identical with each other. It demonstrates that these EOSs will generate similar properties of the neutron stars.
\begin{figure}[!hbt]
	\centering
	\subfigure[]{\includegraphics[width=8cm]{re.eps}}
	\subfigure[]{\includegraphics[width=8cm]{ep.eps}}
	\caption{The energy per nucleon as a function of density (a panel) and the pressure as a function of energy density (b panel) of neutron star matter with Bonn A, Bonn B, and Bonn C potentials.}
	\label{re}
\end{figure}

The central correlation on kinetic energy per nucleon is presented for the neutron star matter in Fig.~\ref{tc}. The variational method can be achieved based on the competition between the correlation on kinetic energy and potential. The central correlation on kinetic energy provides a repulsive effect to prevent two-nucleon approach at high densities, whereas the one on potential gives an attractive contribution to remove the repulsion of realistic $NN$ interactions in short-range region. Finally, the correlations on kinetic and potential energies determine the minimum total energy and explicitly confirm the variational parameters in the Jastrow function. It can be found that the correlation on kinetic energy contributes almost half of the binding energy per nucleon for neutron star matter and plays a very essential role in the RCV method.
\begin{figure}[!hbt]
	\centering
	\includegraphics[width=9cm]{tc.eps}
	\caption{The central correlation on kinetic energy for neutron star matter with Bonn A, Bonn B, and Bonn C potentials.}
	\label{tc}
\end{figure}

In the RCV method there is only one independent variational parameter, $c_4$~\cite{hu11}, which is shown in Fig.~\ref{c4} as a function of density with Bonn A, Bonn B, and Bonn C potentials. From this figure, we can see that the central correlation strength increases slowly at low density, reaches a maximum value around the normal nuclear saturation density, and starts to decrease with increasing density thereafter. It demonstrates that the central correlation will have to consider more variables to generate the saturation density and that it becomes weaker at high density, since the distance between two nucleons is already sufficiently compressed.
\begin{figure}[!hbt]
	\centering
	\includegraphics[width=9cm]{c4.eps}
	\caption{The variational parameter $c_4$ as a function of density for neutron star matter with Bonn A, Bonn B, and Bonn C potentials.}
	\label{c4}
\end{figure}
\subsection{The properties of neutron stars}
With the EOSs from the RCV method, we obtain the properties of neutron stars by solving the TOV equation with Bonn potentials. The mass-radius and mass-density relations of neutron stars within our present framework are illustrated in Fig.~\ref{rm}. In the left panel, we plot the mass-radius relations of neutron stars with Bonn A, Bonn B, and Bonn C potentials. The maximum masses and corresponding radii predicted are almost the same, around $2.18 M_\odot$ and $11$ km, respectively. These results are in good agreement with the previous {\it{ab initio}} calculations. The maximum masses of neutron stars are $2.2 M_\odot$ in non-relativistic variational method by Akmal {\it{et al.}}~\cite{akmal98}. Similarly, the  maximum masses and corresponding radii of neutron stars were given as $2.24 M_\odot$ and $10.8$ km in the RBHF theory with Bonn potentials~\cite{krastev06}. It indicates that our RCV method can economically describe the neutron star matter as well as both the non-relativistic full variational method and the RBHF theory. The maximum masses and the mass-radius relations of neutron stars with three Bonn potentials have obvious distinctions, since the neutron star matter includes not only neutrons but also protons. The tensor force should contribute to some effects at the $T=0$ channel here, especially at low density, which is shown clearly in the mass-radius curves at large $R$. In the right panel of Fig.~\ref{rm}, we plot the mass-density relations of neutron stars with Bonn potentials. The central densities at maximum neutron star mass are $1.078$ fm$^{-3}$ for all of three Bonn potentials, which are similar to the calculations of Akmal {\it{et al.}}~\cite{akmal98}, $1.14$ fm$^{-3}$ and RBHF theory, $1.003-1.013$ fm$^{-3}$~\cite{krastev06}. In the low density region, the masses of neutron stars display distinguishable behaviors. They are almost identical with increasing densities.
\begin{figure}[!hbt]
	\centering
	  \subfigure[]{\includegraphics[width=8cm]{rm.eps}}
	  \subfigure[]{\includegraphics[width=8cm]{rhom.eps}}
	\caption{The mass-radius relations (panel a) and mass-density relations (panel b) of neutron stars with Bonn A, Bonn B, and Bonn C potentials.}
	\label{rm}
\end{figure}

In Table~\ref{pro}, the properties of neutron stars within the present framework are listed for Bonn A, Bonn B, and Bonn C potentials. They are compared with the ones taken from the RBHF theory~\cite{krastev06}. No matter the maximum mass or the corresponding radii and central density, our results are only $3\%$ different from the ones calculated via the RBHF theory. Furthermore, the maximum  masses of neutron stars in our calculation satisfy the requirements of recent observation on massive neutron stars~\cite{demorest10,antoniadis13}, around $2M_\odot$. The corresponding radius is also located within the constraint region worked out by Hebeler {\it{et al.}}~\cite{hebeler10,hebeler13}.
\begin{table}[!hbt]
	\centering
	\begin{tabular}{c c c c c c c}
		\hline
		\hline
		Potentials&                                    &     RCV      &                        &                                    &RBHF            &                     \\
		          &~~$M_{\text{max}}$ ($M_{\odot}$) &     $R$(km)      &  $\rho_c$ (fm$^{-3}$) &~~$M_{\text{max}}$ ($M_{\odot}$)&     $R$(km)    &  $\rho_c$ (fm$^{-3}$)\\
		\hline

		Bonn A    &~~2.184                               &10.99         &1.078                  &~~2.240                               &10.74           &1.013 \\

		Bonn B    &~~2.181                               &11.08         &1.078                  &~~2.240                               &10.79           &1.008 \\

		Bonn C    &~~2.179                               &11.22         &1.078                  &~~2.238                               &10.83           &1.003 \\

		\hline
		\hline
	\end{tabular}
	\caption{The properties of neutron stars (maximum mass, corresponding radii and central density) within present framework and compared with the ones calculated by RBHF theory for Bonn A, Bonn B, and Bonn C potentials.}\label{pro}
\end{table}

The fractions of various particles appearing in neutron stars, which are neutrons, protons, electrons, and muons are plotted in Fig.~\ref{ypabc} with the three Bonn potentials. At the beginning, the muon is absent for the $\beta$ equilibrium conditions. When the electron chemical potential is larger than the muon mass, the muon will appear in the neutron star matter at densities less than $0.2$ fm$^{-3}$, which is smaller than the density of muons in the RBHF theory. The earliest appearance of the muon in the Bonn A potential around the normal saturation density. At high density, the fraction of muons will approach that of electrons. The proton fraction in Bonn A has the largest magnitude compared to the other two Bonn potentials, which should be due to its smaller tensor component.
\begin{figure}[!hbt]
	\centering
	\includegraphics[bb=0 90 290 240, width=0.7\textwidth]{ypabc.eps}
	\caption{The particle fractions in neutron stars with Bonn A, Bonn B, and Bonn C potentials.}
	\label{ypabc}
\end{figure}

We give the proton fractions with Bonn potentials alone in Fig.~\ref{urca}. From this figure, the difference of proton fractions obtained from the three Bonn potentials is clearly revealed, which becomes very obviously in the intermediate density region and decreases in the high density region. The proton fraction in our calculation continues to increase with density and is in contrast to the case in RBHF theory, where it decreases at high density and the largest fraction is about $0.13$~\cite{krastev06}. In RBHF theory, to simplify the calculation, the leptons were treated with the non-relativistic approximation. At high densities, the relativistic effect becomes more important, where the Fermi momentum is very high and can be comparable with the light speed. Therefore, the direct URCA processes related with the cooling mechanism of neutron star could not occur easily with Bonn B and Bonn C potentials in RBHF theory, which should be satisfied by a proton fraction larger than approximately $1/9$. However, in our work, the direct URCA processes can be produced in all of the Bonn potentials and the densities appearing in direct URCA processes are located between $0.25-0.35$ fm$^{-3}$. These densities will lead the neutron star to cool very rapidly.
\begin{figure}[!hbt]
	\centering
	\includegraphics[width=9cm]{urca.eps}
	\caption{The proton fraction in neutron stars with Bonn A, Bonn B, and Bonn C potentials.}
	\label{urca}
\end{figure}
\subsection{The effective nucleon masses in neutron star}
We also give the Dirac effective nucleon masses in neutron star matter in Fig.~\ref{emn}. For proton or neutron Dirac effective masses reflecting the nucleon media effect, their magnitudes with different potential are almost equal, which is very similar situation to the energy and pressure cases. If we compare the proton effective mass with the neutron effective mass, it can be found that the proton effective mass is larger than the neutron one at low density. When the density is $\rho>0.4$ fm$^{-3}$, the splitting of proton and neutron effective masses is reversed. This behavior is quite different from the effective mass in RHF theory, where the proton effective mass should be larger than the neutron one in all density regions~\cite{li16}. For RBHF theory, the splitting of proton and neutron effective masses for asymmetric nuclear matter is strongly dependent on the treatment of the $G-$matrix. Brockmann and Machleidt used the single-particle potential to extract the effective nucleon masses, where the neutron effective mass is larger than the proton one in neutron-rich matter~\cite{sammarruca10}; whereas Dalen {\it{et al.}} adopted the projection method to distinguish the spin components in the $G-$matrix and calculated the effective nucleon mass with these interaction components in the RHF model~\cite{dalen07,dalen10}. In this scheme, the proton effective mass is larger than the neutron one.
\begin{figure}[!hbt]
	\centering
	\includegraphics[width=9cm]{emn.eps}
	\caption{The effective proton and neutron masses in neutron stars with Bonn A, Bonn B, and Bonn C potential.}
	\label{emn}
\end{figure}

To discuss the splitting of proton-neutron effective masses in our framework, we show the kinetic and potential contributions of the nucleon effective mass in a neutron star with Bonn A potential in Fig.~\ref{emc}.  Without the central correlation function, the kinetic energy is a one-body operator, which does not provide any contribution to the nucleon effective mass. Once the central correlation function is included, the central correlation on kinetic energy becomes a two-body operator and contributes to the nucleon effective mass. From Fig.~\ref{emc}, we can find that the potential generates the negative contribution to the effective nucleon mass. In neutron-rich matter, the neutron effective mass will obtain more negative components compared to the proton one. Therefore, in the RHF model, the proton effective mass is larger than the neutron one in neutron-rich matter. Meanwhile, the central correlation on kinetic energy has the positive contributions to effective nucleon mass. Furthermore, its effect on neutrons is much higher than on protons at high densities. With the competition between the potential and kinetic energy, the proton effective mass is larger than the neutron one at low density and with increasing density, the neutron effective mass is larger than the proton one. Therefore, the central correlation on kinetic energy plays a very important role in the splitting of proton-neutron effective masses. For Dalen {\it{et al.}}, they obtained the effective nucleon mass from the potential part; whereas Brockmann and Machleidt considered the splitting of effective masses from the single particle potential, which is related to kinetic energy. This may be the reason why there are opposite conclusions between the two groups, in RBHF theory on effective nucleon mass.
\begin{figure}[!hbt]
	\centering
	\includegraphics[width=9cm]{emc.eps}
	\caption{The kinetic and potential contributions on the nucleon effective masses in neutron star with Bonn A potential.}
	\label{emc}
\end{figure}

\section{Conclusion}
The RCV method based on the framework of RHF theory was applied to study the properties of neutron stars with one-boson-exchange potentials, i.e., Bonn potentials, which were determined by fitting the nucleon-nucleon ($NN$) scattering data. In neutron-rich matter, the tensor force has a very small effect of its isospin feature. Therefore, it is essential to take the central correlation on the strong repulsion of $NN$ interaction at short distances for the description of neutron-rich matter. The equation of state (EOS) of pure neutron matter obtained by a novel {\it{ab initio}} calculation of RBHF theory was completely reproduced by the present RCV method.

The EOSs of neutron star matter in $\beta$ equilibrium with nucleons and leptons, were self-consistently solved in the RCV method with Bonn A, Bonn B, and Bonn C potentials. Their behaviors were almost identical for the weak effect of tensor force in neutron-rich matter, since the difference among three Bonn potentials only appeared in their tensor components. The relativistic central correlation on kinetic energy played a very important role in the process of minimizing the total binding energy, with the variational principle, and gave half of the contribution to total binding energy. Furthermore, we found that the strength of the central correlation function was the strongest at the saturation density through the variational parameters, which correspond to the saturation mechanism of symmetric nuclear matter.

The properties of neutron stars were studied with the EOSs of neutron star matter by solving the TOV equation. The maximum neutron star masses and corresponding radii were around $2.18 M_\odot$ and $11$ km, respectively, using the RCV method with Bonn potentials. The central densities of neutron stars were about $1.078$ fm$^{-3}$. These results are in good agreement with the calculations from RBHF theory and the non-relativistic variational method. It demonstrated that the RCV method can describe neutron-rich matter reasonably and economically, compared with the conventional {\it{ab initio}} calculation.

The proton fractions in neutron star matter with the three Bonn potentials showed some differences. The proton fraction in Bonn A potential was the largest and Bonn C the smallest. The direct URCA processes would be generated in all of these potentials above the densities, $0.25-0.35$ fm$^{-3}$ with a proton fraction larger than approximately $1/9$. Furthermore, the splitting of proton-neutron effective masses was reversed with increasing density of neutron star matter. In the low density region, the proton effective mass was larger than the neutron one; whereas these behaviors are opposite at high densities, since the central correlation on kinetic energy played a more important role at high densities compared with the correlation on potential. This may explain the conflicting predictions about the splitting of proton-neutron masses in the two different treatments in RBHF theory.

Although, we can describe the properties of neutron stars very well, which are comparable to the results obtained by the other {\it{ab initio}} methods, it is necessary to take the tensor correlation into account to study symmetric nuclear matter and its saturation properties to reproduce the empirical data in future.

\section*{Acknowledgments}
This work was supported in part by the National Natural Science Foundation of China (Grant No. 11375089 and Grant No. 11405090).


\section{Coalition Formation}
While the coalitional game in partition form in \eqref{eq:coalitiona_game} characterizes the conflict between the players through their utility functions, a \emph{coalition formation game} specifies the dynamics that change the coalition structures to reach a solution of \eqref{eq:coalitiona_game}. A coalition formation game is formally represented with the pair \cite{Diamantoudi2007}
\begin{equation}
(\sset{P},\gg),
\end{equation}
where $\sset{P}$ is the set of all coalition structures (Definition~\ref{def:coalition_structure}) and $\gg$ is a binary relation on $\sset{P}$, called the \emph{dominance relation} \cite{Diamantoudi2007}. The dominance relation $\gg$ compares two coalition structures in $\sset{P}$ and will be defined as the intersection of two other binary relations in the following.

\begin{definition}[Dominance relation]\label{def:dominance}
Coalition structure $\sset{C}'$ dominates $\sset{C}$, written as $\sset{C}' \gg \sset{C}$, if and only if (iff) $\sset{C}' \succ \sset{C}$, indicating that $\sset{C}'$ is preferred to $\sset{C}$, and $\sset{C}' \xleftarrow[]{} \sset{C}$, indicating that $\sset{C}'$ is reachable from $\sset{C}$.
\end{definition}
The subsequent two sections will deal with the modeling of the dominance relation defined above. Here, we briefly reveal the purpose of the two binary relations which define it: The preference relation $\succ$ on $\sset{P}$ will be used to compare two different coalition structures based on the players' utility functions, i.e., the partition function $\mat{u}$ in \eqref{eq:coalitiona_game}. The reachability binary relation $\xleftarrow[]{}$ on $\sset{P}$ will determine whether it is possible to \emph{change} a coalition structure to another coalition structure by a feasible \emph{deviation} of some players, as illustrated in~\figurename~\ref{fig:group_deviation}. The players' deviation models can be distinguished between:
\begin{itemize}
\item \emph{Group-based deviation}: A \emph{group} of players, possibly members of different coalitions, can leave their coalitions and form a new coalition.
\item \emph{Individual-based deviation}: A \emph{single} player leaves its coalition to join another.
\end{itemize}

The reachability binary relation will affect the number of possible ways to change a given coalition structure. In multiuser networks, this change would typically require communication between the deviators, cf. \figurename~\ref{fig:group_deviation}. Therefore, the feasibility of communication between the devices will directly affect the modeling of the ``reachability'' binary relation. We will explore different parameterized deviation models later in the next two sections.

By the assumption of player rationality, a coalition structure $\sset{C}$ will change to $\sset{C}'$ if $\sset{C}'$ dominates $\sset{C}$. The following notion of coalition structure stability is called \emph{deviation-proof} and is motivated by \cite{Apt2009,Diamantoudi2007}.

\begin{definition}[Stability]\label{def:stability}
Coalition structure $\sset{C} \in \sset{P}$ is stable iff there exists no $\sset{C}' \in \sset{P}$ such that $\sset{C}' \gg \sset{C}$.
\end{definition}

A stable coalition structure indicates that there exists no other coalition structure which is reachable from the stable coalition structure and is preferred by the players. While reachability depends on the feasibility of communication between the devices in multiuser networks, using the preferences of the players as a criterion for stability accounts for efficiency in the outcome. Stable coalition structures are in this regard suitable clustering solutions for multiuser networks.

On modeling the dominance relation (Definition~\ref{def:dominance}), coalition formation is realized through Algorithm~\ref{alg:coalition_formation}. The algorithm is initialized at a suitable coalition structure $\sset{C}$. First, in Step~\ref{alg_line:deviation} a coalition structure $\sset{C}'$ is found which is reachable from $\sset{C}$. Then, Step~\ref{alg_line:preference} checks whether $\sset{C}'$ is preferred by the players compared to $\sset{C}$. If so, the coalition structure changes to $\sset{C}'$. This process is repeated until no dominating coalition structures are found, i.e., stability is reached. Note that the convergence of Algorithm~\ref{alg:coalition_formation} must be studied for the specific application which defines the dominance relation.

\begin{algorithm}[t]
\caption{\label{alg:coalition_formation} Generic coalition formation algorithm.}
\begin{algorithmic}[1]
\State \textbf{Input}: $(\sset{P},\gg)$, coalition structure $\sset{C} \in \sset{P}$
\Repeat
\State Find $\sset{C}'$ such that $\sset{C}' \xleftarrow[]{} \sset{C}$ \label{alg_line:deviation}
\State If $\sset{C}' \succ \sset{C}$ then $\sset{C}$ is updated to $\sset{C}'$ \label{alg_line:preference}
\Until {stability}
\State \textbf{Output}: {$\sset{C}$}
\end{algorithmic}
\end{algorithm}

\section{Discussion}

Coalition formation games essentially provide structured procedures which are useful for solving clustering problems. The significance of coalition formation games are their ability to exploit the users' performance measures, provided these are formulated to depend on the clustering, in order to obtain efficient outcomes. Moreover, owing to the players' rationality assumption, coalition formation games particularly support the design of distributed clustering mechanisms.

The complexity of coalition formation mainly depends on the used deviation model, which ultimately impacts the overhead in communication between the players and the overhead in cooperation within a coalition. While conceptually different, nested group-based deviation ($q$-merge and $q$-split) possess higher complexity than individual-based deviation. Despite this fact, it is difficult to know in advance which of the deviation models leads to better performance in different applications. Note, that further deviation models, other than the ones mentioned in this lecture note, can be developed, as for example in~\cite{Brandt2015a}.

\begin{figure*}[t]
\centering
\begin{minipage}[c]{7.4cm}
\centering
\includegraphics[width=5.3cm,clip]{figures/topology}%
\subcaption{\label{fig:topology} Deployment of transmitters. Line intensity between transmitters is the frequency of cooperation.}
\end{minipage}
\hspace{.4cm}
\begin{minipage}[c]{9.5cm}
\centering
\includegraphics[height=5.21cm,clip]{figures/evaluation
\subcaption{\label{fig:evaluation} Performance of different coalition formation models depending on parameter $q$ for the setup in \figurename~\ref{fig:topology}. The number at the curves corresponds to the average number of MISO links which cooperate.}
\end{minipage}
\vspace{-0.1cm}
\caption[]{\label{fig:simulation} Coalition formation in the example MISO IC.\footnotemark} \vspace{-0.6cm
\end{figure*}

\footnotetext{The Matlab simulation code is available online at https://github.com/rami-mochaourab/coalition-formation}

We simulate coalition formation in the considered MISO IC example with $17$ links. The deployment of the transmitters, each with $t = 8$ antennas, is shown in \figurename~\ref{fig:topology}. The receivers' locations are uniformly distributed within $200$ meters from their associated transmitters. We assume Rayleigh small-scale fading, $46$ dBm transmission power at the transmitters, and the following parameters for the setup \cite{Brandt2015a}: \smallskip

\noindent
\begin{center}
\footnotesize
\begin{tabular}{|l|l||l|l||l|l|}
\hline
Path loss & $15.3 + 3.76$ $\log_{10}$ (distance [meters]) & Receiver noise figure & $9$ dB & Noise PSD & $-174$ dBm/Hz \\
Shadow fading & Log-normal i.i.d. with standard deviation $8$ dB & Carrier frequency & $2$ GHz & Bandwidth & $10$ MHz\\
\hline
\end{tabular}
\end{center}
\smallskip

For this setup, the average performance (over $10^3$ realizations) of stable coalition structures reached through coalition formation (Algorithm~\ref{alg:coalition_formation}) with different deviation models and associated preference relations is shown in \figurename~\ref{fig:evaluation}. Here, we always initialize Algorithm~\ref{alg:coalition_formation} with singleton coalitions

Coalition formation with ``$q$-merge \& $2$-split'' alternates between $q$-merge and $2$-split operations. In order to ensure convergence of this algorithm, as in the individual deviation model, we do not allow a set of coalitions to merge to a coalition which has previously existed. It can be observed that by combining $2$-split with $q$-merge, we allow further dynamics in coalition formation which leads to better performance compared to only $q$-merge deviation. For larger values of $q$, the average coalition sizes of both group-based deviation models increase due to the possibility of merging larger number of coalitions. This, in most cases, leads to improvement in average performance, however at the cost of higher complexity (\figurename~\ref{fig:complexity}). Compared to group-based deviation, individual deviation leads to smaller average coalition sizes and inferior performance, as is also observed in \cite{Mochaourab2014} for a similar MISO IC setting.

In \figurename~\ref{fig:topology}, the line intensity between the transmitters shows the frequency (in percent) of cooperation between any two MISO links when using ``$4$-merge \& $2$-split''. Clearly, it is beneficial to form coalitions between links which are close to each other in order to benefit most from interference management.

\section{Group-based Deviation}
Group-based deviation describes the mechanism in which a set of players, possibly in different coalitions, build a single coalition, while the remaining players do not change their coalitions. An example is given in \figurename~\ref{fig:group_deviation}, where the coalition structure $\sset{C}$ in \figurename~\ref{fig:coalition_structure} changes (in a single deviation step) because the set of deviators form a single coalition. This deviation model is considered in \cite{Diamantoudi2007} with the assumption of farsighted players, meaning that the players are able to anticipate the effects of a sequence of deviations. In this lecture note, we rather assume \emph{nearsighted} players, i.e., the players can only know the effects of a single deviation.

A special type of group-based deviation is called \emph{nested deviation} \cite{Apt2009,Diamantoudi2007} which restricts the change in the coalition structure as follows: 1) Several existing coalitions can only merge to form a single coalition. 2) A single existing coalition can split to form smaller coalitions. We will describe parameterized versions of such nested deviations in the following two subsections.

\subsection{Merge Deviation Model}
The merge deviation allows a set of at most $q$ coalitions in some coalition structure $\sset{C}$ to \emph{merge} and form a single coalition \cite{Mochaourab2014}. Accordingly, the coalition structure $\sset{C}$ changes to $\sset{C}'$. This mechanism is illustrated in~\figurename~\ref{fig:merge}, where two coalitions comprising $\sset{T}$ form a single coalition $\bigcup \sset{T}$. This deviation model defines the reachability binary relation $\sset{C}' \xleftarrow[\textrm{merge}]{} \sset{C}$ as follows.
\begin{definition}[$q$-merge]\label{def:coalition_merging}
For all $\sset{C}',\sset{C} \in \sset{P}$, $\sset{C}' \xleftarrow[\textrm{merge}]{} \sset{C}$ is true iff $\sset{C}' = (\bigcup \sset{T}) \cup (\sset{C} \setminus \sset{T})$ with ${\sset{T} \subset \sset{C}}, {\abs{\sset{T}} \leq q}$.
\end{definition}

The parameter $q$ in this model restricts the number of reachable coalition structures and thus affects the complexity of coalition formation. Specifically, the complexity of $q$-merge corresponds to the worst case number of coalition structures $\sset{C}' \in \sset{P}$ for which $\sset{C}' \xleftarrow[\textrm{merge}]{} \sset{C}$ is true. This number directly determines the number of searches in worst case required in Step~\ref{alg_line:deviation} in Algorithm~\ref{alg:coalition_formation}. Given a coalition structure $\sset{C}$, the number of possible ways to merge a set of at least two and at most $q$ coalitions is given by \cite{Mochaourab2014}
\begin{equation}\label{eq:comp_merging}
{D}(\abs{\sset{C}},q) = \sum\nolimits_{j=2}^{\min\{q,\abs{\sset{C}}\}} \binom{\abs{\sset{C}}}{j}.
\end{equation}
Since the number of coalitions in $\sset{C}$ is largest when $\sset{C}$ includes only single-player coalitions, then the worst case complexity is ${D}(\abs{\set{N}},q)$. From Lemma 1 in \cite{Mochaourab2014} we have that for fixed $q < \abs{\set{N}}$, the growth rate of ${D}(\abs{\set{N}},q)$ in $\abs{\set{N}}$ is bounded by a polynomial of degree $q$. While it is desired to set large values for $q$ in order to encompass a large set of reachable coalition structures, it is however important to account for the associated coalition formation complexity which increases exponentially in $q$.

The convergence of Algorithm~\ref{alg:coalition_formation} using the $q$-merge reachability model in Definition~\ref{def:coalition_merging} is guaranteed since only merging of coalitions is allowed and the number of possible merging steps is finite.

The application of $q$-merge in the MISO IC requires synchronization between the links in order to determine the set of coalitions which deviate. This can be done by consistent enumeration of all \emph{reachable} coalition structures at all links and considering each one in a sequential manner. The consistent enumeration of all reachable coalition structures can be done, e.g., through lexicographic ordering~\cite{Mochaourab2014}. Then, for each reachable coalition structure, the members of the coalitions involved in the deviation send each other a cooperation request. The links' replies to the requests depend on the preference relation, which we discuss later below.

\subsection{Split Deviation Model}
In the split deviation model, a single coalition in some coalition structure $\sset{C}$ can \emph{split} into at most $q$ coalitions. In doing so, the coalition structure $\sset{C}$ changes to $\sset{C}'$ as is illustrated in~\figurename~\ref{fig:split}.
\begin{definition}[$q$-split]\label{def:coalition_splitting}
For all $\sset{C}',\sset{C}\in\sset{P}$, $\sset{C}' \xleftarrow[\textrm{split}]{} \sset{C}$ is true iff $\sset{C}' = (\sset{C} \setminus \set{T}) \cup \br{\set{S}_1,\ldots,\set{S}_k}$, where $\set{T} \in \sset{C}$, $\set{T} = \bigcup_{j=1}^k \set{S}_j$ and $k \leq q$.
\end{definition}

Analogous to $q$-merge, the convergence of Algorithm~\ref{alg:coalition_formation} with $q$-split is guaranteed since only splitting of coalitions is allowed and the number of splitting steps is finite.

\begin{figure*}[t]

  \centering
  \includegraphics[width=10cm,clip]{figures/complexity_all
  \caption{\label{fig:complexity} Complexity of $q$-merge and $q$-split deviations, with $D$ and $T$ defined in \eqref{eq:comp_merging} and \eqref{eq:comp_spliting}, respectively.}
  \vspace{-.6cm}
\end{figure*}
The complexity of the $q$-split deviation model is the worst case number of coalition structures that have to be considered in Step~\ref{alg_line:deviation} in Algorithm~\ref{alg:coalition_formation}. The number of ways to split a set $\set{T}$ into at least two and at most $q$ subsets is \cite{Mochaourab2014}
\begin{equation}\label{eq:comp_spliting}
T(\abs{\set{T}},q) = \sum\nolimits_{\ell=2}^{q} S(\abs{\set{T}},\ell),
\end{equation}
where $S(\abs{\set{T}},\ell) = \frac{1}{\ell!} \sum\nolimits_{k=0}^{\ell}(-1)^{k}\binom{\ell}{k} (\ell - k)^{\abs{\set{T}}}$ are the stirling numbers of the second kind \cite[Ch. 8]{Brualdi2004}. Then, the complexity of Step~\ref{alg_line:deviation} in Algorithm~\ref{alg:coalition_formation} is quantified by $\sum_{\set{T}\in\sset{C}} T(\abs{\set{T}},q)$. For illustration purposes, we consider the complexity of $q$-split for $\sset{C}$ being the grand coalition $\set{N}$ and compare its complexity to $q$-merge for increasing number of players $|\set{N}|$ and different values of $q$ in \figurename~\ref{fig:complexity}. It can be observed that $q$-split has significantly larger complexity compared to $q$-merge; The least complex choice of $q=2$ in the $q$-split deviation model gives an exponential growth ($T(\abs{\set{N}},2)$) which occurs for the most complex $q$-merge deviation model ($D(\abs{\set{N}},\abs{\set{N}})$).

It is possible to define coalition deviation to include a combination of merging and splitting steps~\cite{Apt2009}. The advantages with this deviation model, e.g., allowing more adaptive clustering and guaranteeing efficient clustering outcomes, are discussed in Section ``Coalition Formation Games'' in~\cite{Bacci2016}. Note however that the convergence of coalition formation with sequential merging and splitting is not implied but must be studied depending on the characteristics of the application in multiuser networks.

For the application of $q$-split in the MISO IC example, similar synchronization requirements are needed at the links as in the $q$-merge case, along with a consistent enumeration of all reachable coalition structures.

\subsection{Preference Relations}
A preference relation determines whether a coalition structure $\sset{C}'$ is preferred to another coalition structure $\sset{C}$ based on the utility functions of the players in~\eqref{eq:coalitiona_game}. One preference relation which is appropriate for distributed implementation in multiuser networks, e.g., the MISO IC, is the \emph{Pareto order relation}~$\succ_\set{T}$:
\begin{equation}\label{eq:preference_group}
\sset{C}' \succ_\set{T} \sset{C} \Leftrightarrow u_i(\sset{C}') \geq u_i(\sset{C}) \text{ for all } i \in \set{T},
\text{ and } u_j(\sset{C}') > u_j(\sset{C}) \text{ for some } j \in \set{T}.
\end{equation}
\noindent Note that only a specific set of players $\set{T}$ are considered in $\sset{C}' \succ_\set{T} \sset{C}$, which typically will correspond to the players which deviate. The modeling of the preference relation in~\eqref{eq:preference_group} is suitable for distributed algorithms in multiuser networks since each device in $\set{T}$ knows its utility function and can do the comparison locally. Then, the devices only need to exchange their \emph{binary decisions} to implement the preference relation.

\section{Individual-based deviation}

In individual-based deviation, a coalition structure can change only if a single player leaves its coalition and joins another \cite{Bogomolnaia2002}. This mechanism is illustrated in \figurename~\ref{fig:deviation}
\begin{definition}[Individual Deviation]\label{def:deviation}
For all $\sset{C}',\sset{C} \in \sset{P}$, $\sset{C}' \xleftarrow[\textrm{individual}]{} \sset{C}$ is true iff
\begin{equation}
\sset{C}' = (\sset{C} \setminus \{\set{C}_i,\set{S}\}) \cup \{\set{C}_i\setminus \{i\}, \set{S} \cup \{i\}\},\text{ for some } i \in \set{N} \text{ and } \set{S} \in (\sset{C} \setminus \set{C}_i) \cup \{\}.
\end{equation}
\end{definition}

The convergence of Algorithm~\ref{alg:coalition_formation} with the individual-based deviation model will generally depend on the used preference relation and the utilities of the players~\cite{Bogomolnaia2002}. It is however possible to enforce convergence through suitable restrictions on the deviation model. One widely used approach prevents a player from joining a coalition which it has been a member of before~\cite{Saad2012,Brandt2015a}. According to Theorem 1 in~\cite{Saad2012}, this restriction ensures convergence of the algorithm irrespective of the used preference relation.

In the MISO IC example, individual deviation requires selecting a single link at a time, according to a random or fixed procedure, which asks the members of an existing coalition to join their coalition. If the preference relation allows the deviation to take place, then the link informs the members of its coalition.

While there exist several preference relations which are suitable for individual-based deviation~\cite{Bogomolnaia2002}, we mention in this lecture note one preference relation which considers both the utilities of the \emph{deviating} player $i$ and the members of the coalition $\set{S}$ it requests to join:
\begin{equation}\label{eq:preference_individual}
\sset{C}' \succ_{i,\set{S}} \sset{C} \Leftrightarrow u_i(\sset{C}') > u_i(\sset{C}), \textrm{ and } u_j(\sset{C}') \geq u_j(\sset{C}) \text{ for all } j \in \set{S}.
\end{equation}
The above preference relation is suitable for application in multiuser networks, because the members of $\set{S}$ are required to cooperate with player $i$, such as by updating their beamforming strategies in the example MISO IC. Then, it is reasonable to demand that the utilities of the players in $\set{S}$ do not decrease

\section{Relevance}

Clustering mechanisms are essential in certain multiuser networks for achieving efficient resource utilization. The main purpose of this lecture is to present the theory of coalition formation which is useful for distributed clustering problems. We reveal the generality of the theory, which is currently missing in the literature, and study complexity aspects which must be considered in multiuser networks.
\section{Problem Formulation}

Let $\set{N}$ be the set\footnote{Throughout, we will use calligraphic font for sets and boldface calligraphic font for sets of sets.} of devices in the network, e.g., the set of links in the MISO IC. A clustering of the network, as is illustrated in \figurename~\ref{fig:coalition_structure}, will be called a \emph{coalition structure} \cite{Thrall1963}.
\begin{definition}[Coalition structure]\label{def:coalition_structure}
A {coalition structure} $\sset{C}$ is a partition of $\set{N}$ into a set of pairwise disjoint coalitions. Let $\sset{P}$ be the set containing all coalition structures.
\end{definition}%
Each element of a coalition structure $\sset{C}$ is called a \emph{coalition} and comprises a set of cooperating devices. Let $\set{C}_i \in \sset{C}$ denote the coalition which includes device $i\in \set{N}$ in coalition structure~$\sset{C}$. For suitable application in communication networks, the devices' cooperation strategies within a coalition must be designed taking into account the available local information along with the communication constraints between the devices. For the MISO IC, we assume that each transmitter has perfect local channel state information (CSI), i.e., each transmitter $i$ knows $\mat{h}_{ij}$ for all $j$ perfectly. With this information and assuming a total power constraint $p_i$, one conventional cooperative beamforming strategy for transmitter $i$ nulls the interference at the receivers in its coalition $\set{C}_i$ while maximizing the desired signal power:
\begin{equation}\label{eq:zeroforcing}
\mat{w}_{i}(\set{C}_i) = \argmax_{||\mat{w}_i||^2 \leq p_i} ~ |\mat{h}_{ii}^\dag \mat{w}_i| \text{ s.t. } \mat{h}_{ij}^\dag \mat{w}_i = 0, \text{ for } j\in\set{C}_i\backslash \{i\}.
\end{equation}
The cooperative beamforming strategy in \eqref{eq:zeroforcing} is called \emph{zero forcing beamforming} and admits a simple closed form expression. Observe that \eqref{eq:zeroforcing} also includes the \emph{noncooperation} strategy of a transmitter $i$, namely when $\set{C}_i = \{i\}$, which corresponds to maximum ratio transmission $\mat{w}_{i}(\{i\}) = \sqrt{p_i} |\mat{h}_{ii}|/||\mat{h}_{ii}||$.

\begin{figure*}[t]

  \centering
  \begin{minipage}[t]{5.3cm}
  \centering
  \includegraphics[width=5.3cm,clip]{illustrations/coalition_structure_light}%
  \subcaption{\label{fig:coalition_structure} Coalition structure $\sset{C}$}
  \end{minipage}
  \hspace{1cm}
  \begin{minipage}[t]{5.3cm}
  \centering
  \includegraphics[width=5.3cm,clip]{illustrations/group_deviation_light}%
  \subcaption{\label{fig:group_deviation} General deviation}
  \end{minipage}

  \begin{minipage}[t]{5.3cm}
  \includegraphics[width=5cm,clip]{illustrations/merge_light}
  \subcaption{\label{fig:merge} Merge}
  \end{minipage}%
  \hfill

  \begin{minipage}[t]{5.3cm}
  \includegraphics[width=5cm,clip]{illustrations/split_light}
  \subcaption{\label{fig:split} Split}
  \end{minipage}%
  \hfill
  \begin{minipage}[t]{5.3cm}
  \includegraphics[width=5cm,clip]{illustrations/individual_deviation_light}
  \subcaption{\label{fig:deviation} Individual deviation}
  \end{minipage}
  \vspace{-.2cm}
  \caption{Illustrations for the deviation models.}
  \vspace{-.6cm}
\end{figure*}

The devices' performance measures, called \emph{utility functions}, are functions of the devices' strategies. In the MISO IC example, we model the utility function of link $i$ as its achievable rate
\begin{equation}\label{eq:utility}
u_i(\sset{C}) = \log_2 \Bigg( 1 + \frac{ |\mat{h}_{ii}^\dag \mat{w}_i(\set{C}_i)|^2}{\sum\nolimits_{\substack{j \in \set{N}\setminus\set{C}_i}} |\mat{h}_{ji}^\dag \mat{w}_j(\set{C}_j)|^2 + \underbrace{\sum\nolimits_{\substack{k \in \set{C}_i\setminus\{i\}}} |\mat{h}_{ki}^\dag \mat{w}_k(\set{C}_k)|^2}_{=0} + \sigma^2}\Bigg),
\end{equation}
\noindent which depends on the beamforming strategy in~\eqref{eq:zeroforcing}. Here we have assumed a block flat-fading channel model and single-user decoding at the receivers. Observe, that the utility function in \eqref{eq:utility} depends on the overall coalition structure $\sset{C}$ through the selected cooperation strategy in \eqref{eq:zeroforcing}. Then, the effect of cooperation within coalition $\set{C}_i$ is evident through interference-nulling by the transmitters in the same~coalition.

We will assume that the devices are able to make local decisions and adapt their transmission strategies in a way to increase their associated utility functions. Relating to the MISO IC example, a transmitter~$i$ would seek to cooperate with a set of links in $\set{C}_i$ using~\eqref{eq:zeroforcing} only if this improves its achievable rate  in~\eqref{eq:utility}. Such modeled behavior for the devices conforms with the game theoretic assumption of players' rationality, and is a main characteristic for supporting distributed mechanisms in multiuser networks \cite{Bacci2016}. To this end, we will refer to the devices as \emph{players}.

On characterizing the players in the network and modeling their utility functions depending on the coalition structure, the clustering problem can be formulated as a \emph{coalitional game in partition form} \cite{Thrall1963}
\begin{equation}\label{eq:coalitiona_game}
(\set{N}, \mat{u}).
\end{equation}
\noindent In \eqref{eq:coalitiona_game}, $\set{N}$ is the set of players and $\mat{u}:\sset{P} \rightarrow \mathbb{R}^{\abs{\set{N}}}$ is called the \emph{partition function} which assigns to each player a utility depending on the coalition structure, e.g., the utility function in~\eqref{eq:utility}. The solution of \eqref{eq:coalitiona_game} is a set of coalition structures which are stable according to a certain stability criterion. The problem of finding stable coalition structures is the study of coalition formation games, which we consider next.