# Evolution of binary seeds in collapsing protostellar gas clouds

Tatsuya Satsuka,1,2 Toru Tsuribe,2 Suguru Tanaka,1 and Kentaro Nagamine1,3
1Theoretical Astrophysics, Department of Earth and Space Science, Graduate
School of Science, Osaka University, 1-1 Machikaneyama, Toyonaka 560-0043,
Japan
2 College of Science, Ibaraki University, 2-1-1 Bunkyo, Mito, Ibaraki
310-8512, Japan
3 Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505
S. Maryland Pkwy, Las Vegas, NV 89154-4002, USA
<EMAIL_ADDRESS>
(Accepted xxx. Received xxx; in original form xxx)

###### Abstract

We perform three dimensional smoothed particle hydrodynamics (SPH) simulations
of gas accretion onto the seeds of binary stars to investigate their short-
term evolution. Taking into account of dynamically evolving envelope with non-
uniform distribution of gas density and angular momentum of accreting flow,
our initial condition includes a seed binary and a surrounding gas envelope,
modelling the phase of core collapse of gas cloud when the fragmentation has
already occurred. We run multiple simulations with different values of initial
mass ratio $q_{0}$ (the ratio of secondary over primary mass) and gas
temperature. For our simulation setup, we find a critical value of $q_{\rm
c}=0.25$ which distinguishes the later evolution of mass ratio q as a function
of time. If $q_{0}\ga q_{\rm c}$, the secondary seed grows faster and $q$
increases monotonically towards unity. If $q_{0}\la q_{\rm c}$, on the other
hand, the primary seed grows faster and $q$ is lower than $q_{0}$ at the end
of the simulation. Based on our numerical results, we analytically calculate
the long-term evolution of the seed binary including the growth of binary by
gas accretion. We find that the seed binary with $q_{0}\ga q_{\rm c}$ evolves
towards an equal-mass binary star, and that with $q_{0}\la q_{\rm c}$ evolves
to a binary with an extreme value of $q$. Binary separation is a monotonically
increasing function of time for any $q_{0}$, suggesting that the binary growth
by accretion does not lead to the formation of close binaries.

###### keywords:

accretion, accretion discs – hydrodynamics – circumstellar matter – binaries:
general – stars: formation

††pagerange: Evolution of binary seeds in collapsing protostellar gas
clouds–LABEL:lastpage††pubyear: xxxx

##  1 Introduction

It is widely recognised that the majority of main sequence stars and pre-main
sequence stars are in binary systems and that these binaries have various
distributions in mass ratio $q$, binary separation, and binary frequency
(Duquennoy & Mayor, 1991; Ghez et al., 1993; Kouwenhoven et al., 2005;
Raghavan et al., 2010; Kraus et al., 2011; De Rosa et al., 2014). For main
sequence stars, Duquennoy & Mayor (1991) provided statistics of G-dwarf binary
systems. They found that binaries have log-normal separation distribution, and
the number of binaries with wide separations decreases with increasing $q$,
while the number of binaries with close separations is roughly constant or
increase with increasing $q$. Raghavan et al. (2010) and De Rosa et al. (2014)
found a similar distribution of $q$ for A-type binaries as Duquennoy & Mayor
(1991).

For young stellar objects (YSOs), the number of binaries with intermediate
mass YSOs decreases with increasing $q$ (Kouwenhoven et al., 2005). In
contrast, the number of binaries with low-mass YSOs increases with increasing
$q$ (Kraus et al., 2011), and have a log-normal separation distribution.
Although observers have investigated various distributions of binary
properties in detail as mentioned above, formation process of these binaries
is not well understood.

The favoured scenario of binary formation is the fragmentation during runaway
collapse of cloud core (Boss & Bodenheimer, 1979; Miyama et al., 1984; Tsuribe
& Inutsuka, 1999a, b), and the fragmentation in protostellar disc after the
runaway collapse (Williams & Tohline, 1988; Adams et al., 1989; Bonnell, 1994;
Bonnell & Bate, 1994; Woodward et al., 1994; Vorobyov, 2010). Details of the
fragmentation process during cloud core collapse is described in Tohline
(2002) and the references therein. In the present paper, we do not discuss the
fragmentation process, but focus on the evolution of a seed binary after
fragmentation.

After the fragmentation, the binary seeds start to accrete the surrounding
envelope and grow towards main sequence stars. Since the initial mass of seed
binary is lower than $1$ per cent of the stellar mass (Bonnell & Bate, 1994),
observed physical properties of a binary (i.e., properties at the end of
accretion phase) are completely different from those of the seed binaries.

Recently, large-scale hydrodynamical simulations have been performed in order
to explain the observed binary properties (Bate et al., 2002a, b, 2003;
Attwood et al., 2009; Bate, 2009a, b; Offner et al., 2009). Bate (2009a)
simulated a large-scale, homogeneous, isothermal, and turbulent cloud core
collapse ignoring radiative and magnetic effects, and succeeded in reproducing
the observed multiplicity and binary properties such as the frequency of very
low-mass binaries, the log-normal separation distribution, and the
$q$-distribution of wide and close binaries. However, they did not mention
about how and what physical processes determine the distributions of binaries.
In order to understand that, we need to investigate the gas accretion onto a
binary in more detail.

A number of simulations of gas accretion onto a binary have been performed
using SPH codes in two dimensions (2D) (Dunhill et al., 2015; Young et al.,
2015; Young & Clarke, 2015) and 3D (Artymowicz & Lubow, 1996; Bate & Bonnell,
1997), as well as using grid codes in 2D (Ochi et al., 2005; Hanawa et al.,
2010; D’Orazio et al., 2013; Farris et al., 2014). For example, Bate & Bonnell
(1997) investigated a steady, isothermal, non-self-gravitating gas accretion
onto a binary assuming constant angular momentum and density of gas at the
outer boundary. They found that, in the case of low gas angular momentum, the
$q$-value decreases, because the primary is closer to the mass centre of
binary than the secondary. On the other hand, in the case of high angular
momentum, the $q$-value increases because the infalling gas encounters the
secondary first.

Ochi et al. (2005) also investigated in a similar model to Bate & Bonnell
(1997), and found contradictory results that the $q$-value can decrease even
when the angular momentum of gas is high. They argued that this discrepancy
was caused by the lower numerical resolution in the simulation of Bate &
Bonnell (1997), which may have enhanced the accretion onto the secondary.
However, there were other differences between the simulations of Bate &
Bonnell (1997) and Ochi et al. (2005), such as the gas temperature, computing
method, and gravitational potential of the binary, therefore the ultimate
cause of the difference was unclear at that time.

Young et al. (2015) returned to this problem and investigated the dependence
on gas temperature in the same model as Bate & Bonnell (1997) and Ochi et al.
(2005), where they fixed the angular momentum of gas and changed only the gas
temperature. As a result, they concluded that the disagreement between Bate &
Bonnell (1997) and Ochi et al. (2005) was simply caused by the difference in
the gas temperature. Young et al. (2015) found that, in the case of cold gas,
the gas from circum-binary disc is easily trapped inside the secondary’s Roche
lobe, whereas in the hot case, the flow from the secondary’s Roche lobe to the
primary’s Roche lobe emerges, suppressing the growth of $q$-value.

These previous works that we described above were limited to the most simple
situations, i.e., isothermal non-self-gravitating gas, constant angular
momentum and density of gas at the outer boundary, and an isolated binary
(i.e., no growth of binary by accretion). However, if we want to compare the
simulations with real binary systems in star-forming regions, we need to
consider the unsteady gas accretion caused by the non-uniform distribution of
angular momentum and density of the envelope. Motivated by this current
situation, we perform SPH simulations to investigate the unsteady gas
accretion onto a seed binary considering the non-uniform distribution of
angular momentum and density of infalling envelope. Note that we still ignore
the self-gravity of gas and the growth of binary by accretion, with the aim of
understanding physically how these distributions of gas affect the binary
evolution. If the accreted mass exceeds the initial mass of the seed binary,
the self-gravity of gas and the growth of binary by accretion might cause non-
negligible effects on binary evolution. Therefore, we focus on the short-term
evolution until the accreted mass exceeds the initial binary mass. Sections 2
and 3 describe our model and calculation method. We present the results of our
simulations in Section 4. In Section 5, we discuss the results and give
estimates of the long-term evolution of seed binaries. We conclude in Section
6.

##  2 MODEL OF SEED BINARY AND ENVELOPE

Figure 1: Schematic diagrams of the initial condition of our simulations in
the centre-of-mass frame of the envelope and binary. The cross shows the
origin, two black circles in the centre are the binary seeds, and the gray
regions show the envelope. The top (bottom) panel shows the cross section in a
face-on (edge-on) view. The specific angular momentum vector of gas is aligned
with $z$-axis. The angle $\theta$ is determined by equation (7).

In the present work, we skip the formation process of the seed binary in order
to focus on the gas accretion onto the seed binary. The seed binary and the
surrounding envelope are set up in the initial condition as we describe below
and in Fig. 1.

For the density and specific angular momentum distribution of the envelope, we
take the following profiles:

$\displaystyle\rho(r)$ $\displaystyle=$ $\displaystyle\rho_{\rm
in}\left(\frac{R_{\rm in}}{r}\right)^{2}\quad{\rm and}$ (1) $\displaystyle
j(r)$ $\displaystyle=$ $\displaystyle j_{\rm in}\left(\frac{r}{R_{\rm
in}}\right),$ (2)

where $r$ is distance from the centre-of-mass of the envelope, the constants
$\rho_{\rm in}$ and $j_{\rm in}$ represent the density and specific angular
momentum at $r=R_{\rm in}$, and $R_{\rm in}$ is the radius of the inner edge
of the envelope. This power-law distribution is motivated by the numerical
results often seen for the self-gravitating isothermal cloud core collapse
with rotation (Narita et al., 1984; Matsumoto et al., 1997; Saigo & Hanawa,
1998).

Within the above gas distribution, we assume that the seed binary is formed
via fragmentation of gas, carve out the envelope within $r<R_{\rm in}$, and
convert its mass into the seed binary while conserving mass and angular
momentum. The seed binary is modelled by two point masses, and their masses
are $M_{\rm p}$ and $M_{\rm s}$ for the primary and secondary, respectively.
The total mass of binary is $M_{\rm b}\equiv M_{\rm p}+M_{\rm s}$, the initial
mass ratio is $q_{0}\equiv M_{\rm s}/M_{\rm p}$, and the initial binary
separation is $a_{0}$. We assume that these seeds move circularly around the
centre-of-mass with Keplerian velocities. By conservation laws, the mass and
specific angular momentum of the seed binary is represented by

$\displaystyle M_{\rm b}$ $\displaystyle=$ $\displaystyle\int_{0}^{R_{\rm
in}}4\pi r^{2}\rho(r){\rm d}r=4\pi{R_{\rm in}}^{3}\rho_{\rm in},$ (3)
$\displaystyle J_{\rm b}$ $\displaystyle=$ $\displaystyle\int_{0}^{R_{\rm
in}}4\pi r^{2}\rho(r)j(r){\rm d}r=\frac{1}{2}M_{\rm b}j_{\rm in}$ (4)
$\displaystyle=$ $\displaystyle\frac{q_{0}}{(1+q_{0})^{2}}j_{\rm circ}M_{\rm
b},$

where $j_{\rm circ}\equiv\sqrt{GM_{\rm b}a_{0}}$ is the reference specific
angular momentum. Once $R_{\rm in}$ and $q_{0}$ are determined, $j_{\rm in}$
and $\rho_{\rm in}$ are given by

$\displaystyle\rho_{\rm in}$ $\displaystyle=$ $\displaystyle\frac{M_{\rm
b}}{4\pi{R_{\rm in}}^{3}},$ (5) $\displaystyle j_{\rm in}$ $\displaystyle=$
$\displaystyle\frac{2q_{0}}{(1+q_{0})^{2}}\left(GM_{\rm b}a_{0}\right)^{1/2}.$
(6)

The envelope gas at $r>R_{\rm in}$ is distributed as equations (1) and (2). We
set the velocities of gas at $R_{\rm in}<r<R_{\rm out}$ as

$\frac{j(r)^{2}}{r_{\rm cyl}^{3}}<\frac{GM_{\rm b}}{r}\frac{r_{\rm cyl}}{r},$
(7)

such that the gravitational force exceeds the centrifugal force (see Fig. 1).
The initial radial velocity of gas is assumed as

$v_{r}(r)=\left(\frac{2GM_{\rm b}}{r}-\frac{{j(r)}^{2}}{{r_{\rm
cyl}}^{2}}\right)^{1/2},$ (8)

such that the kinetic energy of gas is equal to the gravitational energy,
where $r_{\rm cyl}$ is cylindrical radius in the centre-of-mass frame. For
simplicity, we assume a non-self-gravitating isothermal gas, and no magnetic
fields and radiation. We emphasize that the important difference between our
work and the previous ones (Bate & Bonnell, 1997; Ochi et al., 2005; Young et
al., 2015; Young & Clarke, 2015) is that we aim to investigate unsteady
evolution of binary by considering non-uniform distribution of gas density and
angular momentum, while we still assume isolated binary and ignore self-
gravity of gas.

##  3 NUMERICAL METHOD

We use GADGET-3 SPH code (originally described by Springel 2005) in three
dimension. In this code, smoothing length of each gas particle is determined
by the number of neighbour particles: $(4\pi/3)h^{3}\rho=N_{\rm ngb}m$, where
$h$ is the smoothing length, $m$ is the mass of SPH particle, and $N_{\rm
ngb}$ is the number of neighbour particles. We adopt $N_{\rm ngb}=50$, which
corresponds to $h\sim 2.3(m/\rho)^{1/3}$. For calculation of hydrodynamics, we
choose a polytropic index $\gamma=1$ assuming an isothermal gas, and adopt
Monaghan-Balsara form of artificial viscosity with the parameters $\alpha=1.0$
and $\beta=2.0$ (see Springel 2005 and references therein). All of our
simulations are performed with $128^{3}$ SPH particles, which is roughly two
orders of magnitude higher than that in Bate & Bonnell (1997). In some of our
simulations, a steady circum-secondary disc forms, and the smoothing length is
about one-tenth to one-fifth of the scale height at the outer edge of the
disc, satisfying the criterion of Young et al. (2015) (see Appendix).

The seed binary is treated as two sink particles with a sink radius of $R_{\rm
sink}=0.01a_{0}$, following Young et al. (2015). The SPH particles are removed
from the computational domain once they fall into the sink radius of each seed
or reach the outer boundary, without any feedback to the seeds.

In the case where the envelope has an angular momentum, circum-stellar discs
are formed (Artymowicz & Lubow, 1996; Bate, 1997; Bate & Bonnell, 1997; Ochi
et al., 2005). We define our circum-stellar discs as the gas with $J<U_{\rm
L1}$, where $J$ is the Jacobi constant of gas and $U_{\rm L1}$ is the sum of
gravitational and centrifugal potential at L1 point in the corotating frame of
the seed binary. The Jacobi constant of gas is written by

$\displaystyle J=\frac{1}{2}v^{2}-\frac{GM_{\rm p}}{r_{\rm p}}-\frac{GM_{\rm
s}}{r_{\rm s}}+c_{s}^{2}\,{\rm ln}\,\rho,$ (9)

where $r_{\rm p}$, $r_{\rm s}$, $c_{s}$ are the distance from the primary, the
distance from the secondary, and the sound speed of gas, respectively.

By our definition of the circum-stellar discs, it is possible that the gas in
the circum-stellar discs will be accreted onto each seed via viscous
evolution. Therefore we define the time-dependent mass ratio as

$q(t)=\frac{M_{\rm s}+\Delta M_{\rm s}(t)}{M_{\rm p}+\Delta M_{\rm p}(t)},$
(10)

where

$\displaystyle\Delta M_{\rm p}(t)$ $\displaystyle=$ $\displaystyle M_{\rm
acc,p}(t)+M_{\rm disc,p}(t),$ (11) $\displaystyle\Delta M_{\rm s}(t)$
$\displaystyle=$ $\displaystyle M_{\rm acc,s}(t)+M_{\rm disc,s}(t).$ (12)

Here, $M_{\rm acc,p}$ is the mass accreted onto the primary, $M_{\rm acc,s}$
is the mass accreted onto the secondary, $M_{\rm disc,p}$ is the mass of
circum-primary disc, and $M_{\rm disc,s}$ is the mass of circum-secondary
disc. We terminate our simulations when $\Delta M_{\rm b}(t)\equiv\Delta
M_{\rm p}(t)+\Delta M_{\rm s}(t)$ exceeds $M_{\rm b}$, because the self-
gravity of gas and the growth of binary cannot be ignored after $\Delta M_{\rm
b}(t)>M_{\rm b}$.

With these setup, we compute the gas accretion onto the seed binary for
several values of $q_{0}$ in the range of $0.1<q_{0}<1.0$ and the sound speed
$c_{s}/\sqrt{GM_{\rm b}/a_{0}}=0.05$ (cold) and $0.25$ (hot). Hereafter, we
use the units of $G=M_{b}=a_{0}=1$, in which the orbital period of the seed
binary corresponds to $2\pi$.

##  4 RESULTS

###  4.1 Formation of discs

Figure 2: Time evolution of gas around the seed binary. The top and middle
panels show the face-on logarithmic surface density maps in the centre-of-mass
frame for $q_{0}=0.1$. The hot case is shown in the left column and the cold
case on the right column. The top panels are at $t=2\pi$, and the middle
panels show the final snapshot at $t=18.1$ (panel ($b$)) and $t=15.1$ (panel
($e$)). The white crosses show the positions of the secondary. The bottom
panels show the time evolution of mass accreted onto the primary $M_{\rm
acc,p}$ (red dashed), mass accreted onto the secondary $M_{\rm acc,s}$ (blue
dashed), mass of the circum-primary disc $M_{\rm disc,p}$ (red dotted), mass
of the circum-secondary disc $M_{\rm disc,s}$ (blue dotted), change from the
initial mass of primary $\Delta M_{\rm p}$ (red solid), change from the
initial mass of secondary $\Delta M_{\rm s}$ (blue solid), change from the
initial binary mass $\Delta M_{\rm b}$ (black solid). Note that all values for
secondary is divided by the initial mass ratio $q_{0}$. All density maps in
this paper are produced using SPLASH visualization code (Price, 2007)

Figure 3: Same as Fig. 2, but with $q_{0}=0.3$. Panels ($a$) and ($d$) are at
$t=2\pi$, panel ($b$) is at $t=6\pi$, and panel ($e$) is at $t=4\pi$.

Figure 4: Same as Fig. 2, but with $q_{0}=0.7$. Top panels are at $t=4\pi$,
and middle panels at $t=6\pi$.

In Figure 2, we show the time evolution of gas surface density (top and middle
panels), and the time evolution of the circum-stellar discs and accreted mass
onto each seed (bottom panels) for $q_{0}=0.1$ case. For both hot (left
column) and cold (right column) cases, the circum-primary disc appears at
$t\sim 2\pi$ (Fig. 2a,d). At this time, a clear density enhancement with a
bridge-like feature exists between the secondary and the primary. Inside this
bridge, we find that the Jacobi constant of gas is dissipated and becomes
$J<U_{\rm L1}$, therefore we regard this region as a shock. Neither the
circum-secondary disc nor the circum-binary disc appears until the final state
(Fig. 2b,e). In Fig. 2c,f, it is seen that the red solid line ($\Delta M_{\rm
p}$) is greater than the blue solid line ($\Delta M_{\rm s}/q$). This means
that the mass ratio decreases, as $\Delta M_{\rm p}$ and $\Delta M_{\rm s}/q$
are compared here. In this case of $q_{0}=0.1$, we see that the growth of
$\Delta M_{\rm b}(t)$ is always dominated by that of the circum-primary disc
$M_{\rm disc,p}(t)$. Therefore, the mass ratio monotonically decreases for
both hot and cold cases.

Figure 3 is the same as Fig. 2, except that it is for $q_{0}=0.3$. The circum-
primary disc and the shock between the primary and the secondary appear at
$t\sim 2\pi$ (Fig. 3a,d), similarly to $q_{0}=0.1$ case. In Fig. 3b,e, the
circum-secondary disc appears. From Fig. 3c,f, it is seen that $\Delta M_{\rm
p}<\Delta M_{\rm s}/q$ (i.e., mass ratio increases) after $t\sim 3\pi$ for
both hot and cold cases.

Fig. 4 is the same as Fig. 2 and 3, but with $q_{0}=0.7$. In this case, both
of the circum-primary disc and the circum-secondary disc are simultaneously
formed after $t\sim 2\pi$. At $t\sim 4\pi$ (Fig. 4a,d), circum-primary disc
and circum-secondary disc are clearly seen. At $t\sim 6\pi$ (Fig. 4b,e), the
circum-binary disc is seen. After that, gas is accreted onto the circum-binary
disc first, before being accreted onto the seeds. Then the gas arrives at
circum-stellar discs through L1 or L2 point.

Figure 5: Properties of circum-stellar discs at the final state of $q_{0}=0.7$
and $c_{s}=0.25$ (hot) case. In each panel, thick red and thin blue solid
lines are for the circum-primary and circum-secondary disc, respectively. The
vertical black dashed and dot-dashed line denote the edge of circum-primary
and circum-secondary disc, respectively. The abscissa is the cylindrical
radius from each seed. Panel (a): Jacobi constant, with the horizontal black
solid line showing $-U_{\rm L1}$. Panel (b): Gas surface density. Panel (c):
Magnitude of vertical gravitational force divided by the magnitude of vertical
gas pressure gradient force. The black horizontal dotted line represents the
equilibrium of vertical forces. Panel (d): Specific angular momentum around
each seed, normalised by $r_{\rm cyl}v_{\rm K}(1-c_{s}^{2}/v_{\rm K}^{2})$.

Fig. 5 shows radial profiles of Jacobi constant (Fig. 5a), gas surface density
(Fig. 5b), normalised gravitational force in vertical direction (Fig. 5c), and
normalised specific angular momentum (Fig. 5d) of circum-stellar discs at the
final state in the hot case with $q_{0}=0.7$. In this figure, each
distribution is averaged in each cylindrical shell. In Fig. 5a, the outer edge
of the circum-primary disc $R_{\rm disc,p}$ (dashed line) and the outer edge
of the circum-secondary disc $R_{\rm disc,s}$ (point-dashed line) are
identified by $J=U_{\rm L1}$. In Fig. 5b, $R_{\rm disc,p}$ (or $R_{\rm
disc,s}$) is indeed consistent with the location where the surface density of
the circum-primary disc (or the circum-secondary disc) start to decline in the
outer region. Thus, the circum-stellar discs indeed consists of gas with
$J<U_{\rm L1}$. Another decline in the surface density appears at $r_{\rm
cyl}\sim 0.02$, which is caused by the sink particle approximation. In Fig.
5c, it is seen that the vertical pressure gradient force roughly balances with
gravitational force at the region of $0.02<r_{\rm cyl}<R_{\rm disc,p}$ or
$R_{\rm disc,s}$, indicating an equilibrium in the vertical direction. In Fig.
5d, it is seen that angular momentum of gas is approximately Keplerian inside
$R_{\rm disc,p}$ or $R_{\rm disc,s}$. In summary, the circum-stellar discs
consist of gas with $J<U_{\rm L1}$ and are vertically supported by pressure
gradient, and rotation velocity of disc is almost Keplerian. This structure is
often seen in the standard disc (Shakura & Sunyaev, 1973).

###  4.2 Short-term evolution of mass ratio

Figure 6: Time evolution of the mass ratio $q(t)$ in hot (thick red lines) and
cold (thin blue lines) cases. The simulations are terminated when the accreted
mass reaches $\Delta M_{\rm b}(t)=M_{\rm b}$.

Figure 6 shows the time evolution of mass ratio $q(t)$. In the cases with
$q_{0}=0.1$, $q(t)$ decreases monotonically. With $q_{0}=0.3$, $q(t)$
decreases until $t\sim 3\pi$, but it starts to increase after $t\sim 3\pi$,
and eventually exceeds $q_{0}$. In the cases with $q_{0}>0.5$, $q(t)$
increases monotonically. Taken together, our results show that the time
evolution of the mass ratio is qualitatively determined by $q_{0}$. If
$q_{0}<0.23$ (in hot case) or $q_{0}<0.26$ (in cold case), the mass ratio at
the final state is smaller than $q_{0}$. If $q_{0}>0.23$ (in hot case) or
$q_{0}>0.26$ (in cold case), on the other hand, the mass ratio at the final
state is larger than $q_{0}$. Dependence on temperature clearly appears at
$t>4\pi$ in the cases with $q_{0}>0.7$. In these cases, the circum-stellar
discs settle in a steady state. Young et al. (2015) found that, in steady
circum-stellar discs, hot gas easily crosses the L1 point from the secondary’s
Roche lobe to the primary’s Roche lobe compared to the cold gas, and in such a
case the growth of $q(t)$ is suppressed. In our simulations, the similar
tendency as described above is seen in the cases with $q_{0}>0.7$.

Bate (2000) investigated the evolution of seed binaries with various
distribution of angular momentum and density of gas, using the protobinary
evolution (PBE) code which employs the steady state solutions of Bate &
Bonnell (1997). In order to examine the results of PBE code, Bate (2000) also
performed three dimensional SPH simulations with $N_{\rm SPH}=1\times 10^{5}$
including angular momentum and density distribution, although limited to one
particular case of $q_{0}=0.6$. Although the distributions of $\rho\propto
r^{-1}$ and $j\propto r^{2}$ adopted in Bate (2000) are different from our
equations (1) and (2), our model is similar to one of the models in Bate
(2000) because both distributions have the same relation (equation 24). Thus
we can compare our results and that of Bate (2000) except for the time-scale.
Focusing on the short-term evolution until $\Delta M_{\rm b}<M_{\rm b}$, our
results mentioned in this subsection is consistent with the results of the SPH
simulation and the PBE calculations in Bate (2000). Therefore, our numerical
results from SPH simulations confirm the semi-analytical results from PBE
calculations in Bate (2000).

##  5 DISCUSSION

###  5.1 Categorising the Accreting gas

Focusing on the short-term evolution while $\Delta M_{\rm b}(t)\lid M_{\rm
b}$, the accreting gas onto the seed binary can be categorised into four
different modes as we describe below. To characterise the properties of
accreting gas in each mode, it is useful to plot the relation between initial
specific angular momentum of gas and $q_{0}$. In Fig. 7, $j_{\rm in}$ (thick
black dotted line) and $j_{\rm out}$ (thick black dot-dashed line) denotes the
initial gas specific angular momentum at $R_{\rm in}$ and $R_{\rm out}$,
respectively, and $j_{\rm M_{b}}$ (thick black dashed line) denotes the
initial gas specific angular momentum at $r_{\rm M_{b}}$, inside which the gas
mass is equal to $M_{\rm b}$. The specific angular momentum of the secondary
and the primary are defined as $j_{s}$ (thick red dashed line) and $j_{p}$
(thick red solid line). The specific angular momentum of L1 point is defined
as $j_{\rm L1}$ (blue solid line). The specific angular momentum of the
circum-binary disc $j_{\rm cb}$ (green dashed line) is defined as

$j_{\rm cb}=\sqrt{\frac{2}{1+q_{0}}}j_{\rm circ},$ (13)

such that the centrifugal potential ${j_{\rm cb}}^{2}/2{r_{\rm cyl}}^{2}$
equals to the gravitational potential $GM_{\rm b}/r_{\rm cyl}$ at $r_{\rm
cyl}=a_{0}/(1+q_{0})$ which is the distance of secondary from the mass centre
(Ochi et al., 2005). Since the initial specific angular momentum of gas is
determined by equation (2), gas with $j_{\rm in}$ is expected to fall first
onto the seed binary. At the end of the short-term evolution, gas with $j_{\rm
M_{b}}$ is expected to fall onto the circum-stellar discs if we ignore the
complex dynamics until the gas falls. In all our simulations, at the end of
the short-term evolution, more than $80\%$ of gas in the circum-stellar discs
comes from the gas whose initial angular momentum is $j_{\rm in}<j<j_{\rm
M_{b}}$. Thus $j_{\rm M_{b}}$ adequately represents the specific angular
momentum of accreted gas at the end of the short-term evolution. Here we
define the specific angular momentum of accreted gas as $j_{\rm acc}$, which
is in the region filled by backslash, mesh, and single slash in Fig. 7.

The first mode of accreting gas is the “circum-primary disc mode” (the region
filled by backslash in Fig. 7). We can see this mode when

$j_{\rm acc}<j_{\rm L1}.$ (14)

For example, when $q_{0}=0.1$, all accreted gas satisfies equation (14),
indicating that the gas easily enters inside L1 point where the Jacobi
constant of the gas is dissipated by the shock as discussed in Subsection 4.1.
Since L1 point and mass centre of the seed binary are in Roche lobe of the
primary, the gas forms a circum-primary disc (Fig. 2b,e), the primary seed
grows, and the mass ratio decreases monotonically (Fig. 6).

The second mode of accreting gas is the “marginal mode” (the region filled by
mesh in Fig. 7). This mode is seen when

$j_{\rm L1}<j_{\rm acc}<j_{\rm s}.$ (15)

When $q_{0}=0.2$, for example, all accreted gas satisfies equation (15). In
this case, most of gas is trapped by the primary similarly to the circum-
primary disc mode. In the end, the mass ratio decreases. However, the gas that
satisfies equation (15) enters inside secondary’s Roche lobe more easily than
in the circum-primary disc mode. Inside secondary’s Roche lobe, the Jacobi
constant of gas is dissipated by the shock. As a result, $M_{\rm acc,s}$
becomes non-negligible in the end.

The third mode of accreting gas is “circum-stellar discs mode” (the region
filled by single slash in Fig. 7). We can see this mode when

$j_{\rm s}<j_{\rm acc}<j_{\rm cb}.$ (16)

In this case, a circum-secondary disc is formed. Once the circum-secondary
disc is formed, $\Delta M_{\rm s}/q_{0}$ dominates, and the mass ratio
increases monotonically. We can see this mode when $q_{0}=0.7$, for example
(see Fig. 4).

The fourth mode is the “circum-binary disc mode” (the region filled by double
slash in Fig. 7). We can see this mode when the specific angular momentum of
gas is larger than $j_{\rm cb}$ (equation 13),

$j>j_{\rm cb}.$ (17)

In this mode, the majority of gas cannot enter inside each Roche lobe because
of the centrifugal barrier, and gas settles down to the circum-binary disc
first. Then, the gas enters inside each Roche lobe through L2 or L3 point, and
falls onto the circum-stellar discs. This behaviour is seen at $t>6\pi$ for
$q_{0}=0.7$ (see Fig. 4b,e). Since $j_{\rm M_{b}}$ is lower than $j_{\rm cb}$
for any $q_{0}$ in our simulations (Fig. 7), the gas in circum-binary disc
mode is not accreted by the end of short-term evolution. Therefore, the
circum-binary disc mode is irrelevant for the $q$-evolution in the short term.
To investigate $q$-evolution in this mode, we need to simulate the long-term
evolution.

In our simulations, the time evolution of mass ratio qualitatively changes at
$q_{\rm c,hot}=0.23$ (hot case) or $q_{\rm c,cold}=0.26$ (cold case). The
values of $q_{\rm c,hot}$ and $q_{\rm c,cold}$ roughly correspond to the
intersection point of $j_{\rm s}$ and $j_{\rm{M_{b}}}$ in Fig. 7. Therefore,
we define a critical initial mass ratio $q_{\rm c}$ at this intersection
point, and we find $q_{\rm c}=0.25$ from Fig. 7. The value of $q_{\rm c,cold}$
is somewhat closer to $q_{\rm c}$ than $q_{\rm c,hot}$. This is because the
gas flow is closer to a ballistic motion in the cold limit than in the hot
case. With a finite gas temperature, pressure gradient force pushes out the
gas in radial direction. Therefore, even if $j_{\rm M_{b}}<j_{\rm L1}$, the
rotation radius of gas with $j_{\rm M_{b}}$ can reach ${j_{\rm
L1}}^{2}/GM_{\rm b}$. Since $j_{\rm M_{b}}$ is monotonically increasing
function of $q_{0}$, $q_{\rm c,hot}$ is somewhat lower than $q_{\rm c}$. The
difference between $q_{\rm c,hot}$ and $q_{\rm c,cold}$ is small since this
push-out effect is expected to be weak when $c_{s}/v_{\rm K}<1$. Here we
emphasize that the critical value $q_{\rm c}=0.25$ was derived only for a
particular distribution of angular momentum and density (equation 24), and
that it was evaluated when $\Delta M_{\rm b}(t)=M_{\rm b}$.

In summary, gas accretion onto the primary dominates in the circum-primary
disc mode and the marginal mode. While in the circum-stellar discs mode, a
circum-secondary disc is formed and accretion onto the secondary becomes
significant enough to increase the mass ratio. The gas in circum-binary disc
mode forms a circum-binary disc.

Figure 7: Relation between the initial mass ratio $q_{0}$ and the specific
angular momentum $j$ of the envelope. Each line shows the specific angular
momentum of secondary seed (thick red dashed), primary seed (thick red solid),
L1 point (thin blue solid), circum-binary disc (thin green dashed), initial
gas specific angular momentum at $R_{\rm in}$ (thick black dotted), at $R_{\rm
out}$ (thick-black dot-dashed), and at $r_{\rm Mb}$ (thick black dashed). The
black open circle at the intersection of $j_{\rm p}$ and $j_{\rm M_{b}}$
indicates the critical value $q_{\rm c}=0.25$. Each shaded region indicates a
different mode of gas accretion: circum-primary disc mode (backslash),
marginal mode (mesh), circum-stellar discs mode (single slash), and circum-
binary disc mode (double slash).

###  5.2 Analytic Estimate of Long-term Evolution

In our numerical simulations, we focus on the short-term evolution until
$\Delta M_{\rm b}(t)=M_{\rm b}$ assuming an isolated binary with no self-
gravity. In this subsection, we discuss the long-term evolution of binary
separation analytically including binary growth by accretion.

There are two effects which change the binary separation by accretion. One is
the increase of binary mass. When the binary mass becomes larger and if the
angular momentum is conserved, then the binary separation becomes smaller
because of stronger gravitational force. The other is the increase of binary
angular momentum, which increases the binary separation. The evolution of
binary separation is determined by the competition between above two effects.
These effects become especially important when $\Delta M_{\rm b}(t)>M_{\rm
b}$. First, we formulate the time evolution of binary in our model. Then, we
discuss one possibility in which the long-term evolution can be predicted
based on our numerical results of short-term evolution.

As for the binary, we define the time-dependent binary mass $M_{\rm b}(t)$,
binary separation $a(t)$, mass ratio $q(t)$. The reference specific orbital
angular momentum can be written as

$j_{\rm circ}(t)=\sqrt{GM_{\rm b}(t)a(t)}.$ (18)

Then the time-dependent orbital angular momentum of binary $J_{\rm b}(t)$ is
written by

$J_{\rm b}(t)=\frac{2q(t)}{(1+q(t))^{2}}M_{\rm b}(t)j_{\rm circ}(t).$ (19)

We introduce following non-dimensional variables:

$\displaystyle\tilde{M}(t)$ $\displaystyle=$ $\displaystyle\frac{M_{\rm
b}(t)}{M_{\rm b}},$ (20) $\displaystyle\tilde{J}(t)$ $\displaystyle=$
$\displaystyle\frac{J_{\rm b}(t)}{J_{\rm b}},$ (21)
$\displaystyle\tilde{a}(t)$ $\displaystyle=$
$\displaystyle\frac{a(t)}{a_{0}},$ (22)

and $j_{\rm circ}(t)$ is represented as

$j_{\rm
circ}(t)=\frac{(1+q(t))^{2}}{q(t)}\frac{q_{0}}{(1+q_{0})^{2}}\frac{{\tilde{J}}(t)}{{\tilde{M}}(t)}j_{\rm
circ}.$ (23)

Note that we stop our simulations when it becomes $\tilde{M}=2$.

As for the envelope, in our model (equations 1 and 2), the specific angular
momentum of gas $j$ and the gas mass inside the radius $r$, $M_{\rm gas}$, has
a relationship

$j\propto M_{\rm gas}\propto r.$ (24)

From equations (20) and (24), $j_{\rm in}$ as a function of time is given by

$j_{\rm in}(t)=\tilde{M}j_{\rm in}.$ (25)

From equations (6), (18) and (25), we have

$\displaystyle j_{\rm in}$ $\displaystyle=$
$\displaystyle\frac{2q_{0}}{(1+q_{0})^{2}}j_{\rm circ},$ (26) $\displaystyle
j_{\rm in}(t)$ $\displaystyle=$
$\displaystyle\frac{2q(t)}{(1+q(t))^{2}}\frac{{\tilde{M}}^{2}(t)}{{\tilde{J}(t)}}j_{\rm
circ}(t).$ (27)

Equations (26) and (27) represent the specific angular momentum at the inner
edge of the envelope. The power indices of ${\tilde{M}}$ and ${\tilde{J}}$ in
equation (27) reflect the spatial distribution of density and angular momentum
in the envelope. If the relation

$\frac{{\tilde{M}}^{2}(t)}{{\tilde{J}}(t)}=1$ (28)

holds and if $q(t)=q_{0}$, equations (26) and (27) are the same in units of
$M_{\rm b}(t)=a(t)=1$ and $M_{\rm b}=a_{0}=1$. This indicates that the
evolution of binary system is self-similar when equation (28) holds and
$q(t)=q_{0}$. Note that, in equations (25) and (27), it is implicitly assumed
that all angular momentum and mass of the envelope is converted to the orbital
angular momentum and mass of the binary.

After the above preparation, we can now discuss the time evolution of binary
separation. From equations (4) and (19), we have

$\tilde{a}(t)=\left(\frac{q(t)}{q_{0}}\right)^{-2}\left(\frac{1+q(t)}{1+q_{0}}\right)^{4}\frac{{\tilde{J}}^{2}(t)}{{\tilde{M}}^{3}(t)}.$
(29)

From equation (29), we can see that the separation becomes larger with
increasing orbital angular momentum of the binary, and that it becomes smaller
with increasing mass. Moreover, the separation also depends on $q(t)$, and
this dependence originates from equation (19). For given $J_{\rm b}(t)$ and
$M_{\rm b}(t)$, one can see from equation (19) that $a(t)$ inside $j_{\rm
circ}$ depends on $q(t)$. If equation (28) and $q(t)=q_{0}$ hold, binary
separation is proportional to accreted mass in our model:

$\tilde{a}(t)={\tilde{M}}(t).$ (30)

The analytic result of equation (30) is consistent with the numerical work by
Bate (2000). Here, we discuss one possibility in which the long-term evolution
can be predicted by reusing the result of the short-term evolution. From
equations (26) and (27), we see that the difference between $j_{\rm in}/j_{\rm
circ}$ and $j_{\rm in}(t)/j_{\rm circ}(t)$ is caused only by the mass ratio,
if equation (28) always holds. According to our simulations, in the hot case
with $q_{0}=0.5$, $q\approx 0.7$ when $\tilde{M}=2$ from Fig. 6. Under the
above assumptions, we can reuse the former result to predict that the mass
ratio would be $q\approx 0.9$ when it reaches $\tilde{M}=3$. Repeating this
procedure, we can predict the long-term evolution of a seed binary.

Figure 8: Binary separation at the end of the short-term evolution
$\tilde{a}(t_{\rm end})$ in the cases of hot (thick red line) and cold (thin
blue line). The black horizontal line denotes $\tilde{a}(t_{\rm end})=2$.

We saw in Fig. 6 that, in the short-term evolution, $q(t)$ increases
monotonically if $q_{0}>q_{\rm c}$, and vice versa. Based on this result and
the argument in the previous paragraph, we argue that the long-term evolution
of $q(t)$ is qualitatively determined by $q_{0}$. Fig. 8 plots equation (29)
at the end of the short-term evolution (i.e., binary separation at
$\tilde{M}=2$) using our numerical results of $q(t)$ and equation (28). Fig. 8
shows that the separation reaches $\tilde{a}(t_{\rm end})=2$ in the cases with
$q_{0}\rightarrow 1.0$ and $q_{0}\simeq q_{\rm c}$, indicating that the time
evolution of a binary is self-similar in these cases (equation 30). Fig. 8
also shows that $\tilde{a}(t_{\rm end})>1$ for any $q_{0}$, which suggests
that the binary separation is a monotonically increasing function of time and
therefore close binaries are difficult to form.

Here, we note again that these analytic results are based on the assumption
that all angular momentum of the envelope is converted to the orbital angular
momentum of the binary. In other words, we are disregarding the division of
gas angular momentum into orbital angular momentum of binary and that of
circum-stellar discs. In order to investigate the growth of separation more
properly, a direct calculation of the binary orbit is needed.

##  6 CONCLUSION AND FUTURE WORK

In the present work, we investigate the short-term evolution of a seed binary
using the SPH code GADGET-3 in three dimensions. Our simulation setup includes
non-uniform distribution of gas and angular momentum with $\rho\propto r^{-2}$
and $j\propto r$, respectively. In the initial condition, the seed binary is
assumed to have formed around the mass centre of the binary by fragmentation,
conserving angular momentum and mass. The seed binary is isolated, and self-
gravity of gas is ignored. With this setup, we compute the accretion of
gaseous envelope onto the seed binary until the binary mass growth exceeds its
initial mass, surveying the parameter ranges of $0.1<q_{0}<1.0$ and the sound
speeds $c_{s}/\sqrt{GM_{\rm b}/a_{0}}=0.05$ (cold) and $0.25$ (hot).

As a result, we categorise the gas accretion into four different modes as
follows:

  1. 1.

“Circum-primary disc mode” is seen when the specific angular momentum of
accreting gas is lower than that of L1 point, i.e., $j_{\rm acc}<j_{\rm L1}$.
Most of the gas falls onto the primary and the circum-primary disc, and hence
$q(t)$ monotonically decreases. This is because the specific angular momentum
is small enough, and the gas with $j_{\rm acc}<j_{\rm L1}$ enters the
primary’s Roche lobe.

  2. 2.

“Marginal mode” is seen when $j_{\rm L1}<j_{\rm acc}<j_{\rm s}$. In this case,
although most of the gas is trapped by the primary similarly to the ”circum-
primary disc mode”, the gas is able to enter the secondary’s Roche lobe, and
the secondary starts to accrete gas. As a result $q(t)$ becomes smaller than
$q_{0}$ after the short-term evolution.

  3. 3.

“Circum-stellar discs mode” is seen when $j_{\rm s}<j_{\rm acc}<j_{\rm cb}$.
If the specific angular momentum of gas exceeds that of the secondary, gas
starts to rotate around the secondary, and a circum-secondary disc is also
formed. Once the circum-secondary disc is formed, $q(t)$ monotonically
increases.

  4. 4.

“Circum-binary disc mode” is seen when $j_{\rm cb}<j$. In this case, gas
cannot fall onto the circum-stellar discs directly because of its large
angular momentum. Therefore, the gas falls onto the circum-binary disc first,
and then later enter the Roche lobes through L2 or L3 point.

We find that the short-term evolution of $q$-value is qualitatively different
according to its initial value $q_{0}$. If $q_{0}>q_{\rm c}=0.25$, the final
mass ratio exceeds $q_{0}$. This critical value $q_{\rm c}$ is determined by
the condition $j_{\rm s}=j_{\rm M_{\rm b}}$ in Fig. 7. The critical value
$q_{\rm c}=0.25$ was derived only for a particular distribution of angular
momentum and density (equation 24), and that it was evaluated when $\Delta
M_{\rm b}=M_{\rm b}$. In circum-primary disc mode or marginal mode, the
dominant accretion onto the primary decreases the $q$-value. However, once the
circum-secondary disc is formed, the accretion onto the secondary becomes
significant enough to increase the mass ratio. The value of $q_{\rm c}$ does
not differ dramatically depending on gas temperature as long as $c_{s}/v_{\rm
K}<1$.

We also estimate the long-term evolution of a seed binary analytically.
Assuming that equation (28) holds, we argue that the evolution of binary
system would be self-similar, and the short-term evolution of $q(t)$ from our
simulations can be reused just by updating the initial mass ratio. As a
result, we find that the binary separation is a monotonically increasing
function of time for any $q_{0}$. This result suggests that close binaries are
difficult to form. In the future, we will include direct computations of
binary orbit in our simulations in order to investigate the effect of binary
growth by accretion.

## Acknowledgments

This work is partially based on the Master’s thesis of Suguru Tanaka (Tanaka,
2010). We are grateful to Volker Springel for providing us with the original
version of GADGET-3 code, and to Kengo Tomida and Fumio Takahara for useful
discussions and continuous encouragement. KN acknowledges the partial support
by JSPS KAKENHI Grant Number 26247022. Numerical simulations were in part
carried out on XC30 at the Centre for Computational Astrophysics, National
Astronomical Observatory of Japan.

## Appendix

### Resolution dependence of mass ratio evolution

Figure 9: Mean rate of mass ratio change ${\dot{q}}$ versus $h(t_{\rm end})/H$
which is the ratio of SPH smoothing length to the scale height at the outer
edge of the circum-secondary disc. The open symbols are for hot case, and the
filled symbols are for the cold case with $N_{\rm SPH}=64^{3}$ (triangle),
1283(circle), $2\times 128^{3}$ (square), $4\times 128^{3}$ (pentagon), 2563
(diamond). The value of $h(t_{\rm end})/H$ is smaller for better resolution.
All of our simulations presented in the main text is performed with $N_{\rm
SPH}=128^{3}$.

In order to investigate the resolution dependence of mass ratio evolution, we
rerun several simulations in the cases with $q_{0}=0.7$ (hot and cold) and
varying number of SPH particles $N_{\rm SPH}=64^{3},\,128^{3},\,2\times
128^{3},\,4\times 128^{3},$ and 2563. In a steady circum-stellar disc,
pressure gradient and gravitational force are in equilibrium in the vertical
direction as mentioned in Subsection 4.1. Therefore, the scale height of this
disc $H$ is represented by

$H(r_{\rm cyl})=\frac{c_{s}}{v_{\rm K}(r_{\rm cyl})}r_{\rm cyl}.$ (31)

We introduce $h/H$, the ratio of SPH smoothing length to the scale height at
the outer edge of the circum-secondary disc. Here, we define the edge of the
circum-secondary disc where $J=U_{\rm L1}$. In Fig. 9, we plot the mean rate
of change of the mass ratio, ${\dot{q}}=(q(t_{\rm end})-q_{0})/t_{\rm end}$,
as a function of $h(t_{\rm end})/H$. In both hot and cold cases, one can see
that the variation of ${\dot{q}}$ is within $25\%$ for different resolution.
Young et al. (2015) suggested that ${\dot{q}}$ is independent of resolution if
the SPH smoothing length at the outer edge of circum-secondary disc satisfies
$h/H<1$. Indeed, we also confirm this criterion with our three dimensional
simulations.

## References

  * Adams et al. (1989) Adams F. C., Ruden S. P., Shu F. H., 1989, ApJ, 347, 959
  * Artymowicz & Lubow (1996) Artymowicz P., Lubow S. H., 1996, ApJ, 467, L77
  * Attwood et al. (2009) Attwood R. E., Goodwin S. P., Stamatellos D., Whitworth A. P., 2009, A&A, 495, 201
  * Bate (1997) Bate M. R., 1997, MNRAS, 285, 16
  * Bate (2000) Bate M. R., 2000, MNRAS, 314, 33
  * Bate (2009a) Bate M. R., 2009a, MNRAS, 392, 590
  * Bate (2009b) Bate M. R., 2009b, MNRAS, 392, 1363
  * Bate & Bonnell (1997) Bate M. R., Bonnell I. A., 1997, MNRAS, 285, 33
  * Bate et al. (2002a) Bate M. R., Bonnell I. A., Bromm V., 2002a, MNRAS, 332, L65
  * Bate et al. (2002b) Bate M. R., Bonnell I. A., Bromm V., 2002b, MNRAS, 336, 705
  * Bate et al. (2003) Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS, 339, 577
  * Bonnell (1994) Bonnell I. A., 1994, MNRAS, 269
  * Bonnell & Bate (1994) Bonnell I. A., Bate M. R., 1994, MNRAS, 271
  * Boss & Bodenheimer (1979) Boss A. P., Bodenheimer P., 1979, ApJ, 234, 289
  * De Rosa et al. (2014) De Rosa R. J., Patience J., Wilson P. A., Schneider A., Wiktorowicz S. J., Vigan A., Marois C., Song I., Macintosh B., Graham J. R., Doyon R., Bessell M. S., Thomas S., Lai O., 2014, MNRAS, 437, 1216
  * D’Orazio et al. (2013) D’Orazio D. J., Haiman Z., MacFadyen A., 2013, MNRAS, 436, 2997
  * Dunhill et al. (2015) Dunhill A. C., Cuadra J., Dougados C., 2015, MNRAS, 448, 3545
  * Duquennoy & Mayor (1991) Duquennoy A., Mayor M., 1991, A&A, 248, 485
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Ghez et al. (1993) Ghez A. M., Neugebauer G., Matthews K., 1993, ApJ, 106, 2005
  * Hanawa et al. (2010) Hanawa T., Ochi Y., Ando K., 2010, ApJ, 708, 485
  * Kouwenhoven et al. (2005) Kouwenhoven M. B. N., Brown A. G. A., Zinnecker H., Kaper L., Portegies Zwart S. F., 2005, A&A, 430, 137
  * Kraus et al. (2011) Kraus A. L., Ireland M. J., Martinache F., Hillenbrand L. A., 2011, ApJ, 731, 8
  * Matsumoto et al. (1997) Matsumoto T., Hanawa T., Nakamura F., 1997, ApJ, 478, 569
  * Miyama et al. (1984) Miyama S. M., Hayashi C., Narita S., 1984, ApJ, 279, 621
  * Narita et al. (1984) Narita S., Hayashi C., Miyama S. M., 1984, Progress of Theoretical Physics, 72, 1118
  * Ochi et al. (2005) Ochi Y., Sugimoto K., Hanawa T., 2005, ApJ, 623, 922
  * Offner et al. (2009) Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ, 703, 131
  * Price (2007) Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159
  * Raghavan et al. (2010) Raghavan D., McAlister H. A., Henry T. J., Latham D. W., Marcy G. W., Mason B. D., Gies D. R., White R. J., ten Brummelaar T. A., 2010, ApJS, 190, 1
  * Saigo & Hanawa (1998) Saigo K., Hanawa T., 1998, ApJ, 493, 342
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Tanaka (2010) Tanaka S., 2010, Master’s thesis, Osaka University
  * Tohline (2002) Tohline J. E., 2002, ARA&A, 40, 349
  * Tsuribe & Inutsuka (1999a) Tsuribe T., Inutsuka S., 1999a, ApJ, 526, 307
  * Tsuribe & Inutsuka (1999b) Tsuribe T., Inutsuka S., 1999b, ApJ, 523, L155
  * Vorobyov (2010) Vorobyov E. I., 2010, ApJ, 723, 1294
  * Williams & Tohline (1988) Williams H. A., Tohline J. E., 1988, ApJ, 334, 449
  * Woodward et al. (1994) Woodward J. W., Tohline J. E., Hachisu I., 1994, ApJ, 420, 247
  * Young et al. (2015) Young M. D., Baird J. T., Clarke C. J., 2015, MNRAS, 447, 2907
  * Young & Clarke (2015) Young M. D., Clarke C. J., 2015, MNRAS, 452, 3085


# Effect of Confinement of Gluons and One Pion Exchange in Nucleon-Nucleon
Interaction

V K Nilakanthan111Corresponding author<EMAIL_ADDRESS>1,2, V C
Shastry1, S Raghavendra1,3, and K B Vijaya Kumar1

###### Abstract

The Nucleon-Nucleon interaction in the singlet (${}^{1}S_{0}$) and triplet
(${}^{3}S_{1}$) channels have been studied in the framework of Relativistic
Harmonic Model using the Resonating Group Method in Born-Oppenheimer
approximation. The full Hamiltonian consists of the kinetic energy, two body
confinement potential, confined one gluon exchange potential and one pion
exchange potential. Contribution of Confined One Gluon Exchange Potential and
One Pion Exchange Potential to the adiabatic nucleon-nucleon interaction
potential is discussed.

1Department of Studies in Physics, Mangalore University, Mangalore, India -
574199.

2Department of Post Graduate Studies and Research in Physics, St Aloysius
College (Autonomous), Mangalore, India - 575003.

3Department of Post Graduate Studies in Physics, S.D.M. College (Autonomous),
Ujire, India - 574240.

Keywords: Nucleon-Nucleon Interaction, Relativistic Harmonic Model, Confined
One Gluon Exchange, One Pion Exchange, Resonating Group Method.

PACS: 13.75.Cs, 12.39.Ki, 21.45.Bc

##  1 Introduction

The nucleon-nucleon (NN) interaction has always been an important and never
ending challenge. It is known that nucleus is treated as a system of nucleons
and that nucleons are made up of quarks, bound by the strong interaction.
Quantum Chromodynamics (QCD) is the widely accepted theory for strong
interaction. The two key features in the QCD are color confinement and
asymptotic freedom. Since the exact form of the confinement is not known from
QCD, phenomenological models have to be employed to study NN interaction.
These phenomenological models include either relativistic models [1, 2, 3, 4,
5, 6] where the interaction is treated perturbatively or non relativistic
quark models (NRQM) [7, 8, 9, 10, 11, 12, 13, 14] where a suitable potential
is used.

In this work we have investigated the singlet and triplet NN potential in the
frame work of relativistic harmonic model (RHM). In an earlier work, the NN
scattering potentials and phase shifts were obtained without taking into
account the contribution of one pion exchange potential (OPEP) [15, 16]. The
aim was to test the effect of confined gluons on NN scattering phase shifts.
Since it is well established that the OPEP provides state-independent
repulsion and contributes substantially to the NN adiabatic potential [17,
18], we feel any reasonable quark model should include OPEP to obtain the
partially conserved axial current, consistent with the chiral symmetry, which
specifies the coupling of pions with quarks. Hence, a picture of the nucleon
of core radii of about 0.6 fm and a pion of small size coupling to quarks
would help to clarify why the long-range NN potential could be accurately
described by the OPEP, which is of vital importance for nuclear physics.

The pion is treated as an elementary field which couples to the quarks in the
framework of NRQM. The strength of the coupling is same as that of the
experimental $\pi$-N coupling strength at zero momentum transfer [17]. In an
alternative approach, a suitable form factor with a cutoff mass $\Lambda$ is
introduced to take care of the effect of the inner structure of pions at short
range [19].

In conventional quark models, the Hamiltonian consists of kinetic energy term,
one gluon exchange potential (OGEP) term and confinement potential. In all
these models, the exchange part of the color magnetic interaction is
considered responsible for the short range repulsion [15, 16]. The OGEP is
obtained from the QCD Lagrangian in the non relativistic limit by retaining
terms to the order of $\frac{1}{c^{2}}$. The gluon propagators used to derive
the OGEP are similar to the free photon propagators in QED, which is used to
obtain the Fermi-Breit interaction. Since confinement of color implies
confinement of quarks and gluons, to determine the nature of NN interaction, a
decisive role should be played by the confined dynamics of the gluons. For the
confinement of quarks we make use of the relativistic harmonic model (RHM) and
for the confinement of gluons we make use of the current confinement model
(CCM) [15]. The confined one gluon exchange potential (COGEP) used in the
present work is derived in the CCM using confined gluon propagators (CGP) [15,
16].

Recently, nuclear potentials were calculated by lattice QCD by utilizing the
long distance behavior of Nambu-Bethe-Salpeter wave functions [20]. Further
insights in this work can be seen in References [21, 22]. In a work by Vinh
Mau et al., it was shown that the behavior of NN interaction at short
distances is dependent on additional terms added for the intermediate and long
range forces [23]. With the experimental data provided by the J-PARC [24],
PANDA [25], NICA [26], and HIAF [27] projects, there is good progress in this
regime.

The paper has been divided into 4 sections. The Relativistic Harmonic Model
and Resonating group method are discussed in the section 2. The results and
discussions are presented in the section 3 and the conclusions in section 4.

##  2 Model

###  2.1 Relativistic Harmonic Model

In the RHM [29] the quarks in a nucleon are considered to be confined by a
Lorentz scalar plus vector harmonic oscillator potential,

$\frac{1}{2}(1+\gamma_{0})\alpha^{2}r^{2}+M$ (1)

where $\gamma_{0}=\left[\begin{array}[]{cc}1&0\\\ 0&-1\\\ \end{array}\right]$
is the Dirac’s matrix, M is a constant mass and $\alpha^{2}$ is the
confinement strength parameter. In RHM the quark wave function $\psi$ is given
by,

$\psi=N\left[\begin{array}[]{c}\phi\\\
\frac{{\bm{\sigma}}.{\bm{p}}}{E+M}\phi\\\ \end{array}\right]$ (2)

where

$N=\sqrt{\frac{2(E+M)}{3E+M}}$

Here $E$ is an eigen value of the single particle Dirac equation. The lower
component of $\psi$ is eliminated to obtain a harmonic oscillator wave
equation in $\phi$

$(\frac{p^{2}}{E+M}+\alpha^{2}r^{2})\phi=(E-M)\phi$ (3)

The full Hamiltonian used in this work is

$H=K+V_{int}+V_{conf}-K_{CM}$ (4)

where $K$ is the Kinetic energy, $V_{int}$ is the interaction potential,
$V_{conf}$ is the confinement potential and $K_{CM}$ is the kinetic energy of
the center of mass.

$K=\sum_{i=1}^{6}\frac{p_{i}^{2}}{(E+M)}$ $K_{CM}=\frac{P^{2}}{6(E+M)}$

where $\frac{(E+M)}{2}$ is the dynamic effective mass of the quarks, $p_{i}$
is momentum of the $i^{th}$ quark and $P$ is the momentum of the center of
mass. The interaction potential is given by,

$V_{int}=V_{COGEP}+V_{OPEP}$

where

$V_{COGEP}=\frac{\alpha_{s}}{4}N^{4}[D_{0}({\bm{r}})+\frac{1}{(E+M)^{2}}(4\pi\delta^{3}({\bm{r}})-c^{2}D_{0}({\bm{r}}))(1-\frac{2}{3}{\bm{\sigma}_{i}}.{\bm{\sigma}_{j}})]{\bm{\lambda}_{i}}.{\bm{\lambda}_{j}}$
(5)

where

$D_{0}(\bm{r})=exp(-\frac{r^{2}c_{0}^{2}}{2})(\alpha_{1}\frac{1}{r}+\alpha_{2})$

In the above equation, $\bm{\lambda}_{i}$ and $\bm{\lambda}_{i}$ are the
generators of the color $SU(3)$ group for the $i^{th}$ and the $j^{th}$
quarks, $\bm{\sigma}_{i}$ and $\bm{\sigma}_{j}$ are the Pauli spin operators
of the $i^{th}$ and the $j^{th}$ quarks and $\alpha_{s}$ is the strong
coupling constant.

The OPEP [17] is given by

$V_{OPEP}=\frac{f_{q}^{2}}{3}\sum_{i<j}\frac{e^{-m_{\pi}r_{ij}}}{r_{ij}}(\bm{\sigma}_{i}.\bm{\sigma}_{j})(\bm{\tau}_{i}.\bm{\tau}_{j})$
(6)

where $\bm{\tau}_{i}$ and $\bm{\tau}_{j}$ are the isospins of the $i^{th}$ and
the $j^{th}$ quarks and $f_{q}$ is the OPEP strength parameter.

###  2.2 Resonating Group Method

The NN interaction exists only when there is exchange of quarks between the
nucleons. If the quarks are not exchanged between the nucleons, there can be
no NN interaction arising from the quark-quark (qq) interaction. This can be
seen from the fact that the matrix element of
$\bm{\lambda}_{i}.\bm{\lambda}_{j}$ vanishes when $i^{th}$ quark is in one
nucleon and $j^{th}$ is in the other nucleon in accordance with the Wigner-
Eckart theorem. Hence the only way NN interaction can arise from a qq
interaction is by constructing a totally antisymmetric wave function for the
six quark system in which case the exchange terms arising solely out of anti-
symmetrization. The RGM employs a totally antisymmetric wave function and
treats the motion of the center of mass correctly. Using the RGM technique we
solve the equation [28],

$\langle\psi|(H-E){\bm{A}}|\psi\rangle=0$ (7)

to get the energy ($E$) of the interacting nucleons. Here, $H$ is the
Hamiltonian, $\psi$ is the wave function of the nucleons and ${\bm{A}}$ is the
anti-symmetrization operator given by,

${\bm{A}}=\frac{1}{10}(1-9P^{OSTC}_{36})$ (8)

where $P^{OSTC}_{36}$ is the permutation operator for the quarks $3$ and $6$
and $OSTC$ stands for orbital, spin, isospin and color. Thus $P^{OSTC}_{36}$
operator exchanges the orbital, spin, isospin and color quantum numbers of the
quarks $3$ and $6$.

The anti-symmetrization operator splits each term in the Hamiltonian into two
parts: direct part and the exchange part. The direct part corresponds to the
self energy of the nucleons and the exchange part gives the interaction
between the nucleons. At asymptotic distances, the exchange part of the
interaction vanishes since the overlap of wave functions is absent.

The harmonic oscillator wave function used here is,

$\phi({\bm{r}}_{i})=\frac{1}{(\pi
b^{2})^{3/4}}exp(-\frac{1}{2b^{2}}({\bm{r}}_{i}-\frac{{\bm{s}}_{I}}{2})^{2})$
(9)

where $b$ is the oscillator size parameter and $s_{I}$ is the generator
coordinate.

The energy is then given by,

$E=\frac{\langle\psi|H{\bm{A}}|\psi\rangle_{l}}{\langle\psi|{\bm{A}}|\psi\rangle_{l}}$
(10)

where the subscript $l$ indicates that the quantities have been projected to
the angular momentum value $l$.

Here, each nucleon is considered as a cluster of three quarks and the two
nucleon system is considered as cluster $A$ and cluster $B$. In RGM, the total
wave function is expressed as an anti-symmetric product of the single particle
wave functions. The total wave function of the six quark system is,

$\psi_{TOT}(\xi_{A},\xi_{B},{\bm{R}}_{AB})={\bm{A}}[\phi_{A}(\xi_{A})\phi_{B}(\xi_{B})\chi({\bm{R}}_{AB})]$
(11)

where $\phi_{A}$ and $\phi_{B}$ are the internal wave functions of the
individual clusters $A$ and $B$ respectively and $\chi$ is the relative wave
function between the two clusters and $\bm{A}$ is the total anti-symmetric
operator of the six quark system. To separate the total wave function given in
equation (11), the following choice of coordinate is made,

$\displaystyle\xi_{1}={\bm{r}}_{1}-{\bm{r}}_{2},\xi_{2}={\bm{r}}_{3}-\frac{{\bm{r}}_{1}+{\bm{r}}_{2}}{2},{\bm{R}}_{A}=\frac{1}{3}({\bm{r}}_{1}+{\bm{r}}_{2}+{\bm{r}}_{3}),$
$\displaystyle\xi_{3}={\bm{r}}_{4}-{\bm{r}}_{5},\xi_{4}={\bm{r}}_{6}-\frac{{\bm{r}}_{4}+{\bm{r}}_{5}}{2},{\bm{R}}_{B}=\frac{1}{3}({\bm{r}}_{4}+{\bm{r}}_{5}+{\bm{r}}_{6}),$
$\displaystyle{\bm{R}}_{AB}={\bm{R}}_{A}-{\bm{R}}_{B},{\bm{R}}_{G}=\frac{1}{2}({\bm{R}}_{A}+{\bm{R}}_{B}).$

Here ${\bm{r}}_{i}$ is the coordinate of the $i^{th}$ quark, the coordinates
$\xi_{A}=(\xi_{1},\xi_{2})$ and $\xi_{A}=(\xi_{3},\xi_{4})$ are the internal
coordinates of the two clusters $A$ and $B$ respectively, ${\bm{R}}_{AB}$ is
the relative coordinate between the two clusters and ${\bm{R}}_{G}$ is the
center of mass coordinate of the total system.

Since the Hamiltonian is translationally invariant, $L^{l}_{ij}$ can be
written as,

$\displaystyle
L^{l}_{ij}=\int[\phi^{+SM}_{A}({\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{3},\frac{{\bm{s}}_{I}}{2})\phi^{+SM}_{B}({\bm{r}}_{4},{\bm{r}}_{5},{\bm{r}}_{6},\frac{-{\bm{s}}_{I}}{2})Y^{*}_{lm}({\hat{\bm{s}}}_{I})][H-E]A$
$\displaystyle[\phi^{SM}_{A}({\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{3},\frac{{\bm{s}}_{J}}{2})\phi^{SM}_{B}({\bm{r}}_{4},{\bm{r}}_{5},{\bm{r}}_{6},\frac{-{\bm{s}}_{J}}{2})Y_{lm}({\hat{\bm{s}}}_{J})]\prod_{k=1}^{6}d^{3}{\bm{r}}_{k}d{\hat{\bm{s}}}_{I}d{\hat{\bm{s}}}_{J}$

To take into account all possible interactions between the quarks, we have to
consider seven different types of operators for the potential $V_{ij}$ in the
Hamiltonian. They are $V_{12DR}$, $V_{36DR}$, $V_{12EX}$, $V_{13EX}$,
$V_{16EX}$, $V_{14EX}$ and $V_{36EX}$, where $DR$ stands for direct part of
the quark interaction between quarks $i$ and $j$ and $EX$ stands for the
corresponding exchange part.

##  3 Results and Discussions

To analyze the contributions of the various components of the Hamiltonian, we
have plotted the diagonal elements of the various kernels of the singlet and
triplet NN potentials as a function of the relative distance between the
nucleons $(s_{I})$. The adiabatic potential is calculated using the Born-
Oppenheimer approximation given by,

$V^{Ad}_{12}=\langle\psi_{1}(s_{I})|H|\psi_{2}(s_{I})\rangle-\langle\psi_{1}(\infty)|H|\psi_{2}(\infty)\rangle$

M (MeV) | 160.6
---|---
E (MeV) | 40.0
$\alpha^{2}$ (MeV fm2) | 200.01
b (fm) | 0.6
$\alpha_{s}$ | 6.5
c (fm-1) | 0.3
$\alpha_{1}$ | 1.035994
$\alpha_{2}$ (MeV) | 2.016150
c0 (fm-1) | 1.7324
$m_{\pi}$ (MeV) | 140.0
Table 1: List of Parameters

There are ten parameters in our model - the masses of the quarks ($M$), the
confinement strength($\alpha^{2}$), the harmonic oscillator size parameter
($b$), the quark-gluon coupling constant ($\alpha_{s}$), the current
confinement model parameter ($c$) and mass of the pion ($m_{\pi}$). The
coupling constant $\alpha_{s}$ is fixed by the N-$\Delta$ mass splitting which
comes from the color magnetic term of COGEP. We have chosen the value of
oscillator size parameter to be 0.6 fm which is consistent with the
experimental results of the charge distribution of the nucleons and the axial
charge distribution [17]. The fixing of the other parameters are discussed in
[15, 16]. The parameters used are listed in Table 1.

Figure 1: Adiabatic NN potential without OPEP

(a)  (b)

Figure 2: Adiabatic NN potential with OPEP; $f_{q}^{2}=12.6$ (left),
$f_{q}^{2}=20.6$ (right)

Figure 1 gives the plot of adiabatic NN potential without OPEP and Figure 2
gives the plot of adiabatic NN potential with OPEP. The intermediate range
attraction for ${}^{1}S_{0}$ state is reduced in the presence of OPEP. The
change in the potential with change in quark pion coupling constant value
$f_{q}^{2}$ can also be seen from Figure 2. The short range repulsion for
${}^{1}S_{0}$ state is larger than that of ${}^{3}S_{1}$ state. The color
magnetic part of the COGEP can be considered as the reason for this
difference.

Figure 3: Kinetic energy

Figure 3 gives the plot of adiabatic kinetic energy vs $s_{I}$ (fm). It is
clear from the figure that there is substantial contribution to the short
range repulsion from the kinetic energy. Since direct part of kinetic energy
do not contribute to the NN interaction the whole of the repulsion arises from
the exchange part of the kinetic energy.

At short distances, the exchange kernels of $\delta^{3}(\bm{r})$ is dominant
over the exchange kernels of $c^{2}D_{0}(\bm{r})$ and hence provide the short
range repulsion. The short range repulsion is larger for ${}^{1}S_{0}$ state
than the ${}^{3}S_{1}$ state. In the intermediate and long ranges, the
exchange kernels of $c^{2}D_{0}(\bm{r})$ dominates over the exchange kernels
of $\delta^{3}(\bm{r})$ and hence provide the intermediate and long range
attraction. The intermediate range attraction is larger for the ${}^{3}S_{1}$
state than the ${}^{1}S_{0}$.

(a)  (b)

Figure 4: Direct and exchange components of NN potential (left) and color
magnetic exchange part of COGEP (right).

Figure 4 is a plot of the direct and exchange parts of the Hamiltonian and the
exchange part of the color magnetic interaction in the adiabatic limit. The
exchange part of the potentials of ${}^{1}S_{0}$ and ${}^{3}S_{1}$ states show
repulsion in the short range. The ${}^{1}S_{0}$ state exchange potential is
completely repulsive and that of the ${}^{3}S_{1}$ state shows a small
attraction in the intermediate range. There is consistency with established
results in the case of repulsive contribution to the adiabatic potential at
short range, both to ${}^{1}S_{0}$ and ${}^{3}S_{1}$ states [17, 30, 31].
There is no contribution from the color electric term to the NN interaction.
Since the energy difference between $2(0s)^{3}$ and $(0s)^{6}$ configuration
must come from the expectation value of
${\bm{\lambda}}_{i}.{\bm{\lambda}}_{j}$, the radial matrix elements are same
for $2(0s)^{3}$ configuration and $(0s)^{6}$ configuration. The expectation
value of the ${\bm{\lambda}}_{i}.{\bm{\lambda}}_{j}$ depends only on the
number of quarks and hence the color electric elements of the COGEP and the
confinement potential do not contribute NN adiabatic potential. The
expectation value of ${\bm{\lambda}}_{i}.{\bm{\lambda}}_{j}$
${\bm{\sigma}}_{i}.{\bm{\sigma}}_{j}$ in the color magnetic part of the COGEP
for the $2(0s)^{3}$ and $(0s)^{6}$ configuration does not vanish and thus, the
color magnetic part provides short range repulsion [17].

##  4 Conclusions

In this work, we have investigated the role played by the COGEP and OPEP on
${}^{1}S_{0}$ and ${}^{3}S_{1}$ NN adiabatic potentials using the Born-
Oppenheimer approximation without evoking $\sigma$ meson in the framework of
the RHM and the CCM by employing the RGM formalism.

The aim was to understand the role played by the confined gluon exchange and
one pion exchange on NN interaction. The COGEP was derived using the CGP. The
major conclusions drawn from the present work are:

  1. 1.

The term arising out of the confinement of gluons is responsible for
intermediate range attraction and hence plays the role of $\sigma$ mesons. The
color magnetic terms of the COGEP yield both short range repulsion and
intermediate range attraction for ${}^{1}S_{0}$ and ${}^{3}S_{1}$ states. The
calculation clearly demonstrates the importance of the confinement of gluons
on ${}^{1}S_{0}$ and ${}^{3}S_{1}$ state NN potentials.

  2. 2.

The OPEP and exchange terms of Kinetic energy provides state independent
repulsion.

To conclude, we have obtained ${}^{1}S_{0}$ and ${}^{3}S_{1}$ NN adiabatic
potentials starting from the dynamics of quarks and gluons including
confinement of color and OPEP as a phenomenological input. Further study is in
progress.

## References

  * [1] P L Ferreira and N Zagury, Gluon Exchange Corrections for Baryons in a Model with Confining Potentials, Letters in Nuovo Cimento 20, 511 (1977).
  * [2] P L Ferreira, J A Helayel, and N Zagury, A Linear Potential Model for Quark Confinement, Nuovo Cimento A 55, 215 (1980).
  * [3] M Bander, D Silverman, B Klima, and U Maor, Relativistic Treatment of Quark - Anti-quark Spectra in QCD, Phys. Rev. D 29, 2038 (1984).
  * [4] S Godfrey and N Isgur, Mesons in a relativized quark model with chromodynamics, Phys. Rev. D 32, 189 (1985).
  * [5] S Capstick and N Isgur, Baryons in a Relativized Quark Model with Chromodynamics, Phys. Rev. D 34, (1986) 2809.
  * [6] S B Khadkikar and S K Gupta, Magnetic Moments of Light Baryons in the Harmonic Model, Phys. Lett. B 124, (1983) 523.
  * [7] O W Greenberg, Quarks, Ann. Rev. Nucl. Part. Phys. 28, 327 (1978).
  * [8] R K Bhaduri, L E Cohler, and Y Nogami, A Unified Potential for Mesons and Baryons, Nuovo Cimento A 65, 376 (1981).
  * [9] D B Lichtenberg, Energy Levels of Quarkonia in Potential Models, Int. J. Mod. Phys. A 2, 1669 (1987).
  * [10] M Narodetskii, R Ceuleneer, and C Semay, Hyperfine interaction in the nonrelativistic quark model and the convergence of harmonic oscillator variational method, J. Phys. G 18, 1901 (1992).
  * [11] N Isgur, D Scora, B Grinstein, and M B Wise, Semileptonic B and D Decays in the Quark Model, Phys. Rev. D 39, 799 (1989).
  * [12] D Scora and N lsgur, Semileptonic meson decays in the quark model: An update, Phys. Rev. D 52, 2783 (1995).
  * [13] M G Olsson, S Veseli, and K Williams,Observations on the potential confinement of a light fermion, Phys. Rev. D 51, 5079 (1995).
  * [14] R F Labed, Meson mass splittings in the nonrelativistic model, Phys. Rev. D 47, 1134 (1993).
  * [15] S. B. Khadkikar and K. B. Vijaya Kumar, N N scattering with exchange of confined gluons, Phys. Lett. B 254, 320 (1991).
  * [16] K. B. Vijaya Kumar and S. B. Khadkikar, N N interaction in a relativistic harmonic model with confined gluons, Nucl. Phys. A 556, 396 (1993).
  * [17] K. Shimizu, Study of Baryon Baryon Interactions and Nuclear Properties in the Quark Cluster Model, Phys. Lett. B 148, 418 (1984).
  * [18] C. S. Vanamali and K. B. Vijaya Kumar, Nucleon-nucleon interaction with one-pion exchange and instanton-induced interactions, Phys. Rev. C. 94 054002 (2016).
  * [19] F. Fernandez and E. Oset, A Model of the Double Spin Flip $NN$ Amplitude Based on a One Pion Exchange Potential With Quark Exchange, Nucl. Phys. A 455, 720 (1986).
  * [20] N. Ishii, NN interation from the lattice, AIP Conf. Proc. 1388, 623 (2011).
  * [21] S. Aoki, T. Hatsuda, and N. Ishii, Nuclear Force from Monte Carlo Simulations of Lattice Quantum Chromodynamics, Comput. Sci. Dis. 1, 015009 (2008).
  * [22] K. Murano, N. Ishii, S. Aoki,and T Hatsuda, Energy dependence of nucleon-nucleon potentials in lattice QCD, PoS LAT2009, 126 (2009).
  * [23] R. Vinh Mau, C. Semay, B. Loiseau, and M. Lacombe, Nuclear forces and quark degrees of freedom, Phys. Rev. Lett. 67, 1392 (1991).
  * [24] S. Kumano, Hadron physics at J-PARC, Nucl. Phys. A 782, 442 (2007).
  * [25] W. Erni et al. (PANDA Collaboration), Physics Performance Report for PANDA: Strong Interaction Studies with Antiprotons, arXiv:0903.3905 [hep-ex].
  * [26] A. N. Sissakian et al. (NICA Collaboration), The nuclotron-based ion collider facility (NICA) at JINR: New prospects for heavy ion collisions and spin physics, J. Phys. G 36, 064069 (2009).
  * [27] J. C. Yang et al., High Intensity heavy ion Accelerator Facility (HIAF) in China, Nucl. Instrum. Methods Phys. Res. Sec. B 317, 263 (2013).
  * [28] M. Oka and K. Yazaki, Short Range Part of Baryon Baryon Interaction in a Quark Model, Prog. Theor. Phys. 66 (1981) 556.
  * [29] P. C. Vinodkumar, K. B. Vijaya Kumar, and S. B. Khadkikar, Effect of the confined gluons in quark-quark interaction, Pramana J. Phys. 39, 47 (1992).
  * [30] F. Myhrer and J. Wroldnen, The nucleon-nucleon force and the quark degrees of freedom, Rev. Mod. Phys. 60, 629 (1988).
  * [31] A. Valcarce, A. Buchmann, F. Fernandez, and A. Faessler, Can one simultaneously describe the deuteron properties and the nucleon-nucleon phase shifts in the quark cluster model?, Phys. Rev. C 50, 2246 (1994).