text:
# Hybrid stars in the quark-meson coupling model with superconducting quark matter

## I Introduction

It is widely believed that hadronic matter undergoes a phase transition to quark matter at high densities and/or high temperatures. The high temperature limit is important in heavy ion collisions and/or cosmology, whereas the high baryon density behavior is important for the study of neutron, hybrid and quark stars.

After a gravitational colapse of a massive star takes place, a neutron star with practically zero temperature can be born few seconds after deleptonization. The correct calculations of the star properties depend on the appropriate equations of state (EOS) that describe its crust and interior. The crust of the neutron star, where density is low is believed to be adequately described by hadronic matter. Its interior, however, where density is of the order of $`510`$ times nuclear saturation density remains to be properly understood. Whether the central part of the star is composed of quark matter only, of mixed matter or of paired quark matter is one of the subjects of the present work.

In the present paper we are interested in building the EOS for mixed matter of quark and hadron phases. We employ the quark-meson coupling model (QMC) guichon ; ST  including hyperons in order to describe the hadron phase. In the QMC model, baryons are described as a system of non-overlapping MIT bags which interact through the effective scalar and vector mean fields, very much in the same way as in the Walecka model (QHD) qhd . Many applications and extensions of the model have been made in the last years tsht ; pal ; recent ; temp ; phase ; jin96 ; mj98 .

While the QMC model shares many similarities with QHD-type models, it however offers new opportunities for studying nuclear matter properties. One of the most attractive aspects of the model is that different phases of hadronic matter, from very low to very high baryon densities and temperatures, can be described within the same underlying model, namely the MIT bag model. In the QMC, matter at low densities and temperatures is a system of nucleons interacting through meson fields, with quarks and gluons confined within MIT bag. For matter at very high density and/or temperature, one expects that baryons and mesons dissolve and the entire system of quarks and gluons becomes confined within a single, big MIT bag. Another important aspect of the QMC is that the internal structure of the nucleon is introduced explicitly. It is found that the EOS for infinite nuclear matter at zero temperature derived from the QMC model is much softer than the one obtained in the Walecka model qhd . Also, the QMC model nucleon effective mass lies in the range $`0.7`$ to $`0.8`$ of the free nucleon mass, which agrees with results derived from non-relativistic analysis of scattering of neutrons from lead nuclei mahaux87  and is larger in comparison with Walecka model effective mass. Consequently, we expect that contrary to the NL3 and TM1 parameterizations of the non linear Walecka model (NLWM) with hyperons, the mass of the nucleon does not become zero at densities $`\rho <10\rho _0`$, $`\rho _0`$ being the density of nuclear matter at saturation. At finite temperature, there arises yet another difference between predictions of QMC and QHD, namely the behavior of the effective nucleon mass with the temperature at fixed baryon density. While in QHD-type models the nucleon mass always decreases with temperature, in the QMC it increases. The difference arises because of the explicit treatment of the internal structure of the nucleon in the QMC. When the bag is heated up, quark-antiquark pairs are excited in the interior of the bag, increasing the internal energy of the bag phase . In what follows we consider the QMC with a constant bag constant, e. g. not dependent on the scalar meson field as in the modified quark meson coupling models (MQMC) jin96 ; mj98 .Contrary to the QMC result, the bag radius increases with density for all MQMC models with a meson dependent bag constant mj98 . For densities not much larger than nuclear matter saturation density $`\rho _0`$ the bags start to overlap which implies a breakdown of the model. In the present work we are interested in describing properties of nuclear matter for densities which go beyond the saturation density and therefore have chosen to consider the QMC model.

For the quark phase we have chosen to use both unpaired quark matter (UQM) described by the MIT bag model bag ; freedman ; farhi ; kapusta  and paired quarks described by the color-flavor locked (CFL) phase. Recently many authors shh ; ar ; bo ; srp ; arrw ; kw  have discussed the possibility that the quark matter is in a color-superconducting phase, in which quarks near the Fermi surface are paired, forming Cooper pairs which condense and break the color gauge symmetry mga . At sufficiently high density the favored phase is called CFL, in which quarks of all three colors and all three flavors are allowed to pair.

Once the mixed EOS is built by enforcing appropriate Gibbs criteria and chemical equilibrium conditions, the properties of the stars are calculated and discussed. We have restricted ourselves to the investigation of neutron stars ($`T=0`$) and hence, neutrino trapping is not an expected mechanism because it is important only during the first seconds after gravitational collapse of the core of a massive star bl86 ; kj95 .

In recent works recente1 ; recente2 , two of us have studied the properties of mixed stars whose equations of state were built with the non-linear Walecka model (NLWM) for the hadron matter qhd  and the MIT Bag bag ; freedman ; farhi ; kapusta  and the Nambu-Jona-Lasinio njl  models for the unpaired quark matter. In the first work recente1  the effects of temperature were investigated and in the second one the consequences of imposing neutrino trapping at fixed entropies were studied. In both works, only in very special cases, the interior of the stars were made of quarks only. In general, a mixed phase of hadrons and quarks was favored.

In what follows we compare the properties of neutron stars obtained within the QMC model with hyperons and with both quark models discussed above, namely, the onset of hyperons, mixed and quark phases, strangeness content, the maximum allowed masses and core composition.

The present paper is organized as follows: in section II the QMC model with hyperons is reviewed. In sections III and IV the unpaired quark matter and the CFL phase are described. In section V the mixed phase is implemented and the results are shown and discussed in section VI. Finally, in the last section, the conclusions are drawn.

## II The quark-meson coupling model for hadronic matter

In the QMC model, the nucleon in nuclear medium is assumed to be a static spherical MIT bag in which quarks interact with the scalar and vector fields, $`\sigma `$, $`\omega `$ and $`\rho `$ and these fields are treated as classical fields in the mean field approximation. The quark field, $`\psi _q(x)`$, inside the bag then satisfies the equation of motion:

$$\left[i\text{/}(m_q^0g_\sigma ^q\sigma )g_\omega ^q\omega \gamma ^0+\frac{1}{2}g_\rho ^q\tau _z\rho _{03}\right]\psi _q(x)=0,q=u,d,s,$$
(1)
where $`m_q^0`$ is the current quark mass, and $`g_\sigma ^q`$, $`g_\omega ^q`$ and $`g_\rho ^q`$ and denote the quark-meson coupling constants. The normalized ground state for a quark in the bag is given by guichon ; ST

$$\psi _q(𝐫,t)=𝒩_q\mathrm{exp}\left(iϵ_qt/R_B\right)\left(\begin{array}{c}j_0\left(x_qr/R_B\right)\\ i\beta _q\stackrel{}{\sigma }\widehat{r}j_1\left(x_qr/R_B\right)\end{array}\right)\frac{\chi _q}{\sqrt{4\pi }},$$
(2)
where

$$ϵ_q=\mathrm{\Omega }_q+R_B\left(g_\omega ^q\omega +\frac{1}{2}g_\rho ^q\tau _z\rho _{03}\right);\beta _q=\sqrt{\frac{\mathrm{\Omega }_qR_Bm_q^{}}{\mathrm{\Omega }_q+R_Bm_q^{}}},$$
(3)
with the normalization factor given by

$$𝒩_q^2=2R_B^3j_0^2(x_q)\left[\mathrm{\Omega }_q(\mathrm{\Omega }_q1)+R_Bm_q^{}/2\right]/x_q^2,$$
(4)
where $`\mathrm{\Omega }_q\sqrt{x_q^2+(R_Bm_q^{})^2}`$, $`m_q^{}=m_q^0g_\sigma ^q\sigma `$, $`R_B`$ is the bag radius of the baryon, and $`\chi _q`$ is the quark spinor. The quantities $`\psi _q,ϵ_q,\beta _q,𝒩_q,\mathrm{\Omega }_q,m_q^{}`$ all depend on the baryon considered. The bag eigenvalue, $`x_q`$, is determined by the boundary condition at the bag surface

$$j_0(x_q)=\beta _qj_1(x_q).$$
(5)
The energy of a static bag describing baryon $`B`$ consisting of three ground state quarks can be expressed as

$$E_B^{\mathrm{bag}}=\underset{q}{}n_q\frac{\mathrm{\Omega }_q}{R_B}\frac{Z_B}{R_B}+\frac{4}{3}\pi R_B^3B_B,$$
(6)
where $`Z_B`$ is a parameter which accounts for zero-point motion and $`B_B`$ is the bag constant. The set of parameters used in the present work is given in table 1. The effective mass of a nucleon bag at rest is taken to be

$$M_B^{}=E_B^{\mathrm{bag}}.$$
(7)
The equilibrium condition for the bag is obtained by minimizing the effective mass, $`M_B^{}`$ with respect to the bag radius

$$\frac{dM_B^{}}{dR_B^{}}=0.$$
(8)
For the QMC model, the equations of motion for the meson fields in uniform static matter are given by

$$m_\sigma ^2\sigma =\underset{B}{}g_{\sigma B}C_B(\sigma )\frac{2J_B+1}{2\pi ^2}_0^{k_B}\frac{M_B^{}(\sigma )}{\left[k^2+M_B^2(\sigma )\right]^{1/2}}k^2𝑑k,$$
(9)
$$m_\omega ^2\omega _0=\underset{B}{}g_{\omega B}\left(2J_B+1\right)k_B^3/(6\pi ^2),$$
(10)
$$m_\rho ^2\rho _{03}=\underset{B}{}g_{\rho B}I_{3B}\left(2J_B+1\right)k_B^3/(6\pi ^2).$$
(11)
In the above equations $`J_B`$, $`I_{3B}`$ and $`k_B`$ are respectively the spin, isospin projection and the Fermi momentum of the baryon species $`B`$. The hyperon couplings are not relevant to the ground state properties of nuclear matter,but information about them can be available from the levels in $`\mathrm{\Lambda }`$ hypernuclei chrien .

$$g_{\sigma B}=x_{\sigma B}g_{\sigma N},g_{\omega B}=x_{\omega B}g_{\omega N},g_{\rho B}=x_{\rho B}g_{\rho N}$$

and $`x_{\sigma B}`$, $`x_{\omega B}`$ and $`x_{\rho B}`$ are equal to $`1`$ for the nucleons and acquire different values in different parameterizations for the other baryons. Note that the $`s`$-quark is unaffected by the sigma and omega mesons i.e. $`g_\sigma ^s=g_\omega ^s=0.`$

In Eq. (9) we have

$$g_{\sigma B}C_B(\sigma )=\frac{M_B^{}(\sigma )}{\sigma }=\frac{E_B^{\mathrm{bag}}}{\sigma }=\underset{q=u,d}{}n_qg_\sigma ^qS_B(\sigma )$$
(12)
where

$$S_B(\sigma )=_{bag}𝑑𝐫\overline{\psi }_q\psi _q=\frac{\mathrm{\Omega }_q/2+R_Bm_q^{}(\mathrm{\Omega }_q1)}{\mathrm{\Omega }_q(\mathrm{\Omega }_q1)+R_Bm_q^{}/2};q(u,d).$$
(13)
The total energy density and the pressure including the leptons can be obtained from the grand canonical potential and they read

$`\epsilon `$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _0^2+{\displaystyle \frac{1}{2}}m_\rho ^2\rho _{03}^2`$ (14)
$`+`$ $`{\displaystyle \underset{B}{}}{\displaystyle \frac{2J_B+1}{2\pi ^2}}{\displaystyle _0^{k_B}}k^2𝑑k\left[k^2+M_B^2(\sigma )\right]^{1/2}+{\displaystyle \underset{l}{}}{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^{k_l}}k^2𝑑k\left[k^2+m_l^2\right]^{1/2},`$
$`P`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _0^2+{\displaystyle \frac{1}{2}}m_\rho ^2\rho _{03}^2`$ (15)
$`+`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{B}{}}{\displaystyle \frac{2J_B+1}{2\pi ^2}}{\displaystyle _0^{k_B}}{\displaystyle \frac{k^4dk}{\left[k^2+M_B^2(\sigma )\right]^{1/2}}}+{\displaystyle \frac{1}{3}}{\displaystyle \underset{l}{}}{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^{k_l}}{\displaystyle \frac{k^4dk}{\left[k^2+m_l^2\right]^{1/2}}}.`$
The lepton Fermi momenta are the positive real solutions of $`(k_e^2+m_e^2)^{1/2}=\mu _e`$ and $`(k_\mu ^2+m_\mu ^2)^{1/2}=\mu _\mu =\mu _e`$. The equilibrium composition of the star is obtained by solving the set of Eqs. (9)- (11) in conjunction with the charge neutrality condition (16) at a given total baryonic density $`\rho =_B(2J_B+1)k_B^3/(6\pi ^2)`$; the baryon effective masses are obtained self-consistently in the bag model.

For stars in which the strongly interacting particles are baryons, the composition is determined by the requirements of charge neutrality and $`\beta `$-equilibrium conditions under the weak processes $`B_1B_2+l+\overline{\nu }_l`$ and $`B_2+lB_1+\nu _l`$. After deleptonization, the charge neutrality condition yields

$$q_{\mathrm{tot}}=\underset{B}{}q_B(2J_B+1)k_B^3/(6\pi ^2)+\underset{l=e,\mu }{}q_lk_l^3/(3\pi ^2)=0,$$
(16)
where $`q_B`$ corresponds to the electric charge of baryon species $`B`$ and $`q_l`$ corresponds to the electric charge of lepton species $`l`$. Since the time scale of a star is effectively infinite compared to the weak interaction time scale, weak interaction violates strangeness conservation. The strangeness quantum number is therefore not conserved in a star and the net strangeness is determined by the condition of $`\beta `$-equilibrium which for baryon $`B`$ is then given by $`\mu _B=b_B\mu _nq_B\mu _e`$, where $`\mu _B`$ is the chemical potential of baryon $`B`$ and $`b_B`$ its baryon number. Thus the chemical potential of any baryon can be obtained from the two independent chemical potentials $`\mu _n`$ and $`\mu _e`$ of neutron and electron respectively.

We start by fixing the free-space bag properties for the QMC model. For the bag radius $`R_N=0.6`$, we first fixed the two unknowns $`Z_N`$ and $`B_N`$ for nucleons. These are obtained by fitting the neuclon mass $`M=939`$ MeV and enforcing the stability condition for the bag at free space. The values obtained are $`Z_N=3.98699`$ and $`B_N^{1/4}=211.303`$ MeV for $`m_u=m_d=0`$ MeV and $`Z_N=4.00506`$ and $`B_N^{1/4}=210.854`$ MeV for $`m_u=m_d=5.5`$ MeV. We then fixed these bag values, $`B_B`$, for all baryons and the parameter $`Z_B`$ and $`R_B`$ of the other baryons are obtained by reproducing their physical masses in free space and again enforcing the stability condition for their bags. The values obtained for $`Z_B`$ and $`R_B`$ are displayed in table I for $`m_u=m_d=0`$ MeV and for $`m_u=m_d=5.5`$ MeV. Note that for a fixed bag value, the equilibrium condition in free space results in an increase of the bag radius and a decrease of the parameters $`Z_B`$ for the heavier baryons.

Next we fit the quark-meson coupling constants $`g_\sigma ^q`$, $`g_\omega =3g_\omega ^q`$ and $`g_\rho =g_\rho ^q`$ for the nucleon to obtain the correct saturation properties of the nuclear matter, $`E_Bϵ/\rho M=15.7`$ MeV at $`\rho =\rho _0=0.15`$ fm<sup>-3</sup>, $`a_{sym}=32.5`$ MeV, $`K=257`$ MeV and $`M^{}=0.774M`$. We have $`g_\sigma ^q=5.957`$, $`g_{\omega N}=8.981`$ and $`g_{\rho N}=8.651`$. We take the standard values for the meson masses, $`m_\sigma =550`$ MeV, $`m_\omega =783`$ MeV $`m_\rho =770`$ MeV.

For the meson-hyperon coupling constants we have opted for three sets discussed in the literature : set a) based on quark counting arguments we take $`x_{\sigma B}=x_{\omega B}=x_{\rho B}=\sqrt{2/3}`$ as in moszk ; set b) according to gm91 ; Glen00  we choose the hyperon coupling constants constrained by the binding of the $`\mathrm{\Lambda }`$ hyperon in nuclear matter, hypernuclear levels and neutron star masses ($`x_\sigma =0.7`$ and $`x_\omega =x_\rho =0.783`$) and assume that the couplings to the $`\mathrm{\Sigma }`$ and $`\mathrm{\Xi }`$ are equal to those of the $`\mathrm{\Lambda }`$ hyperon; set c) based on the $`SU(6)`$ symmetry for the light quarks ($`u,d`$) counting rule pal  we take $`x_{\sigma B}=x_{\omega B}=2/3`$ and $`x_{\rho \mathrm{\Lambda }}=0,x_{\rho \mathrm{\Sigma }}=2,x_{\rho \mathrm{\Xi }}=1`$.

In Fig. 1 we have plotted the EOS obtained with the above parametrization of QMC and two choices of the hyperon couplings, set a) and b) together with the corresponding EOS obtained with NLWM with cubic and quartic sigma terms ($`K=`$300 MeV and $`M^{}=0.7`$ M). Although the NLWM-EOS is harder at low densities, it becomes softer at higher energies after the onset of the hyperons. This fact has consequences on the behavior of the mixed phase and on the star properties. Moreover, one can see that different choices of the meson-hyperon parameters have a greater influence on the NLWM, for which the curves separate at lower densities and become more distant one from the other at high densities, than in the QMC.

## III Unpaired quark matter equation of state

The possible existence of quark matter in the core of neutron stars is an exciting possibility freedman . Densities of these stars are expected to be high enough to force the hadron constituents or nucleons to overlap, thereby yielding quark matter. We take the quark matter equation of state as in Ref. farhi ; kapusta  in which $`u`$,$`d`$ and $`s`$ quark degrees of freedom are included in addition to electrons. Up and down quark masses are set to 5.5 MeV and the strange quark mass is taken to be 150 MeV. In chemical equilibrium $`\mu _d=\mu _s=\mu _u+\mu _e`$. In terms of neutron and electric charge chemical potentials $`\mu _n`$ and $`\mu _e`$, one has

$$\mu _u=\frac{1}{3}\mu _n\frac{2}{3}\mu _e,\mu _d=\frac{1}{3}\mu _n+\frac{1}{3}\mu _e,\mu _s=\frac{1}{3}\mu _n+\frac{1}{3}\mu _e.$$
(17)
The pressure for quark flavor $`f`$, with $`f=u,d`$ or $`s`$ is kapusta

$`P_q`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle \underset{f}{}}\left[\mu _fk_f(\mu _f^22.5m_f^2)+1.5m_f^4ln\left({\displaystyle \frac{\mu _f+k_f}{m_f}}\right)\right],`$ (18)
where the Fermi momentum is $`k_f=(\mu _f^2m_f^2)^{1/2}`$.

For the leptons, the pressure reads

$$P_l=\frac{1}{3\pi ^2}\underset{l}{}\frac{p^4dp}{\sqrt{p^2+m_l^2}}.$$
(19)
The total pressure, including the bag constant B, which simulates confinement becomes

$$P=P_l+P_qB.$$
(20)
There are only two independent chemical potentials $`\mu _n`$ and $`\mu _e`$. $`\mu _e`$ is adjusted so that the matter is electrically neutral, i.e. $`P/\mu _e=0`$.

## IV Color-flavor locked quark phase

In this section we study the equation of state taking into consideration a CFL quark paired phase. We treat the quark matter as a Fermi sea of free quarks with an additional contribution to the pressure arising from the formation of the CFL condensates.

The CFL phase can be described with the help of the thermodynamical potential which reads ar :

$$\mathrm{\Omega }_{CFL}(\mu _q,\mu _e)=\mathrm{\Omega }_{quarks}(\mu _q)+\mathrm{\Omega }_{GB}(\mu _q,\mu _e)+\mathrm{\Omega }_l(\mu _e),$$
(21)
where $`\mu _q=\mu _n/3`$ and

$$\mathrm{\Omega }_{quarks}(\mu _q)=\frac{6}{\pi ^2}_0^\nu p^2𝑑p(p\mu _q)+\frac{3}{\pi ^2}_0^\nu p^2𝑑p(\sqrt{p^2+m_s^2}\mu _q)\frac{3\mathrm{\Delta }^2\mu _q^2}{\pi ^2}+B,$$
(22)
with $`m_u=m_d`$ set to zero,

$$\nu =2\mu _q\sqrt{\mu _q^2+\frac{m_s^2}{3}},$$
(23)
$`\mathrm{\Omega }_{GB}(\mu _q,\mu _e)`$ is a contribution from the Goldstone bosons arising due to the chiral symmetry breaking in the CFL phase son ; ar :

$$\mathrm{\Omega }_{GB}(\mu _q,\mu _e)=\frac{1}{2}f_\pi ^2\mu _e^2\left(1\frac{m_\pi ^2}{\mu _e^2}\right)^2,$$
(24)
where

$$f_\pi ^2=\frac{(218ln2)\mu _q^2}{36\pi ^2},m_\pi ^2=\frac{3\mathrm{\Delta }^2}{\pi ^2f_\pi ^2}m_s(m_u+m_d),$$
(25)
$`\mathrm{\Omega }_l(\mu _e)`$ is the negative of expression (19), and the quark number densities are equal, i.e.,

$$\rho _u=\rho _d=\rho _s=\frac{\nu ^3+2\mathrm{\Delta }^2\mu _q}{\pi ^2}.$$
(26)
In the above expressions $`\mathrm{\Delta }`$, the gap parameter, is taken to be 100 MeV ar .

The electric charge density carried by the pion condensate is given by

$$Q_{CFL}=f_\pi ^2\mu _e\left(1\frac{m_\pi ^4}{\mu _e^4}\right).$$
(27)
In the above thermodynamic potential, we have neglected the contribution due to the kaon condensation which is an effect of order $`m_s^4`$ and thereby small compared to the$`\mathrm{\Delta }^2\mu _q^2`$ contribution to the thermodynamic potential for $`\mathrm{\Delta }100`$ MeV.

## V Mixed Phase and hybrid star properties

We now consider the scenario of a mixed phase of hadronic and quark matter. In the mixed phase charge neutrality is imposed globally i.e. the quark and hadronic phases are not neutral separately but rather, the system prefers to rearrange itself so that

$$\chi \rho _c^{QP}+(1\chi )\rho _c^{HP}+\rho _c^l=0$$
(28)
where $`\rho _c^{QP}`$ and $`\rho _c^{HP}`$ are the charge density of quark and hadron phase, $`\chi `$ is the volume fraction occupied by the quark phase, $`(1\chi )`$ is the volume fraction occupied by the hadron phase and $`\rho _c^l`$ is the lepton charge density. As usual, the phase boundary of the coexistence region between the hadron and quark phase is determined by the Gibbs criteria. The critical pressure and critical neutron and electron chemical potentials are determined by the conditions,

$$\mu _{HP,i}=\mu _{QP,i}=\mu _i,i=n,e,T_{HP}=T_{QP},P_{HP}(\mu _{HP},T)=P_{QP}(\mu _{QP},T),$$

reflecting the needs of chemical, thermal and mechanical equilibrium, respectively. The energy density and the total baryon density in the mixed phase read:

$$\epsilon =\chi \epsilon ^{QP}+(1\chi )\epsilon ^{HP}+\epsilon ^l,$$
(29)
$$\rho =\chi \rho ^{QP}+(1\chi )\rho ^{HP}.$$
(30)
Notice that in all equations above the quark phase (QP) can be either the UQM or the CFL phase. The EOS for the mixed phase are then constructed. Once they are obtained, the properties of the neutron stars can be computed. The equations for the structure of a relativistic spherical and static star composed of a perfect fluid were derived from Einstein’s equations by Oppenheimer and Volkoff tov . They are

$$\frac{dP}{dr}=\frac{G}{r}\frac{\left[\epsilon +P\right]\left[M+4\pi r^3P\right]}{(r2GM)},$$
(31)
$$\frac{dM}{dr}=4\pi r^2\epsilon ,$$
(32)
with $`G`$ as the gravitational constant and $`M(r)`$ as the enclosed gravitational mass. We have used $`c=1`$. Given an EOS, these equations can be integrated from the origin as an initial value problem for a given choice of the central energy density, $`(\epsilon _0)`$. The value of $`r(=R)`$, where the pressure vanishes defines the surface of the star. We solve the above equations to study the structural properties of the star, using the EOS derived above.

## VI Results and discussion

In all figures shown, set a) for the meson-hyperon coupling constants was used, unless stated otherwise. We have omitted all results for the parameter sets b) and c) because they are very similar to the ones obtained with set a). Actually, all differences appear only in the hadron phase at densities where the mixed phase is already the dominant one. It is worth emphasizing that, as stated in sections III and IV, the $`u`$ and $`d`$ quark masses are different in the UQM and CFL models.

In figures 2a) and 2b) the EOS obtained with the unpaired quark model and the color flavor locked phase are displayed for different values of the Bag pressure $`B`$. The onset of the mixed phase and quark pure phase occurs at lower denities for smaller values of $`B`$. This effect has already been discussed in prak97  for the UQM description. A smaller value of $`B`$ gives a softer EOS in the mixed phase because the onset of the mixed phase occurs at lower densities. However at higher densities, after the onset of the quark phase it becomes harder. The mixed phase shrinks with the decrease of the $`B`$ parameter in both quark models. This fact can be also observed in table 2, where the beginning and ending energy densities of the mixed phase are displayed in the last two columns. For the CFL, no mixed phase was found for B<sup>1/4</sup>=180 MeV with $`m_s=150`$ MeV and $`\mathrm{\Delta }=100`$ MeV, giving rise to a pure quark matter star. For the sake of comparison we have plotted the EOS for the same bag pressure for the UQM and CFL phase in figure 3. We include in the same figure the EOS obtained with NLWM plus UQM (short-dashed line). One can see from this figure and table 2 that the mixed phase in the CFL appears and ends at lower energy densities than in the UQM. This result has the same qualitatively behavior as the ones shown in ar . Comparing QMC and NLWM with UQM one can see that they behave similarly at low densities. For $`ϵ>4.5\text{fm}^4`$ the NLWM-EOS becomes softer in the mixed phase. This is due to the fact that at higher densities the NLWM with hyperons becomes softer than the QMC-EOS, Fig. 1.

In figure 4 the particle population for the baryons, leptons and quarks are shown for the UQM and different bag pressures and for the CFL with one chosen Bag value. Hyperons only appear if the bag constant is very high, in our case $`B^{1/4}=210.85`$ MeV. This value corresponds the $`B_B`$ value of the bag used in the QMC model for the hadronic phase. A different behavior was obtained with NLWM in recente1  where the hyperons are present both for $`B^{1/4}=180`$ and 190 MeV. As already discussed, both the mixed phases and the quark phase appear at lower densities for lower $`B`$ values. For the CFL phase and $`B^{1/4}`$ =200 MeV, the quarks appear at $`2.67`$ nuclear saturation density. In ar , for a bag pressure of $`B^{1/4}`$ =190 MeV, the quarks appear at 2.149 $`\rho _0`$. We can see that our results are compatible with the ones shown there, obtained within a NLWM formalism with $`K=240`$ MeV and $`M^{}=0.78M`$. While in figures 4c), $`u,d`$ and $`s`$ quark populations come out different, in figure 4d), because of the imposition of equal quark densities in the CFL model, they are forced to be equal. Also, the $`e^{}`$ population, which disappears at the onset of the quark phase in all the EOS studied, in the CFL model, goes to zero at a much lower density than in the UQM model, because the onset of a pure quark phase occurs at lower densities.

As already discussed in recente1 ; recente2 , the presence of strangeness in the core and crust of neutron and proto-neutron stars has important consequences in understanding some of their properties. In figure 5 we show the strangeness fraction defined as

$$r_s=\chi r_s^{QP}+(1\chi )r_s^{HP}$$
(33)
with

$$r_s^{QP}=\frac{\rho _s}{3\rho },r_s^{HP}=\frac{_B|q_s^B|\rho _B}{3\rho },$$

where $`q_s^B`$ is the strange charge of baryon $`B`$, for the different models discussed in this work. In all cases the strangeness fraction rises steadily. If the UQM is used, at the onset of the pure quark phase it has reached 30 per cent of the baryonic matter. Although the amount of strangeness varies in the mixed phase, it is the same in the pure quark phase independently of the model used to describe the hadron phase. Nevertheless, if the CFL model is used for the quark phase, as a result of the equal quark densities imposition, the strangeness content reaches 1/3 of the total baryonic matter. Comparing QMC and NLWM with UQM, we conclude the strangeness content increases faster when the NLWM is used. This is due to the fact that in this model the hyperons also contribute.

In table 2 we show the values obtained for the maximum mass of a neutron star as function of the central density for some of the EOS studied in this work with UQM and CFL. Different bag parameters are used. We can see that the maximum mass of the star increases and its central energy decreases with increasing $`B`$. This result agrees with the fact that a larger $`B`$ value corresponds to a harder EOS at high densities as shown in Fig. 2. For $`B^{1/4}=180`$ MeV (QMC+UQM) and $`B^{1/4}=190`$ MeV (QMC+CFL), the central density of the star lies outside of the range of the mixed phase, an indication of a star with a quark core. In the last case we get a very high central density due to the very soft EOS this parametrization gives rise to. In fact the predicted maximum mass for a hybrid star within this parametrization is too low as discussed below. In all the other cases discussed the central density of the star is always within the mixed phase. We have also added some results obtained with the non-linear Walecka model recente1  instead of the QMC for the hadron phase for the sake of comparison. For the same $`B`$, the maximum mass is about the same, the central core is also made up of quarks and the mixed phase starts and finishes at lower densities. In ar , the authors have found maximum masses around 1.6 $`M_{}`$ for maximum $`B`$ values of $`B^{1/4}=185`$ MeV. For these $`B`$ values, our results come out at the same order.

The radius of the maximum mass star is sensitive to the low density EOS. In order to calculate the radius and to plot it versus the star mass, we have used the results of Baym, Pethick and Sutherland bps  for low baryonic densities.

In Fig. 6 the mass of the family of stars obtained with QMC is plotted in terms of their radii for both quark models used and $`B^{1/4}=190`$ and 200 MeV. We also include the family of stars obtained with NLWM plus UQM with $`B^{1/4}=190`$ MeV. The radius of the stars within QMC is 11.95 Km ($`B^{1/4}=190`$ MeV) and 12.47 Km ($`B^{1/4}=200`$ MeV) for UQM and 8.79 Km ($`B^{1/4}=190`$ MeV) and 13.31 Km ($`B^{1/4}=200`$ MeV) for CFL. For the NLWM family of stars we get $`R=12.53`$ Km ($`B^{1/4}=190`$ MeV). Some conclusions can be drawn. Comparing QMC and NLWM for $`B^{1/4}=190`$ MeV, the maximum mass of a stable star is similar, for both models. The fact that NLWM- $`M_{}`$ is larger shows that the main contribution to the star comes from the less dense regions, where the NLWM-EOS is harder than the QMC-EOS (see Fig. 3). We consider now the families of stars obtained within QMC for both $`B`$ values and using UQM and CFL. The quark contribution becomes more important for the smaller $`B`$ values mainly for the CFL results. In particular, for $`B^{1/4}=190`$ MeV we get a quite small maximum $`M_{}`$ and a corresponding small radius. A maximum mass of 1.32 $`M_{}`$ is too small for accounting for the presently known radio-pulsar masses tc , even after corrections due to rotation. Similar numbers where obtained with the NLWM plus CFL in ar  for slightly different values of $`B`$ and $`m_s`$. In fact, since these maximum mass stars have a small hadron exterior region the properties of the star is mainly determined by the quark model used. It is also clear from Fig. 6 that the star properties are sensitive to the $`B`$ value, in particular within the CFL formalism. Taking $`B^{1/4}=200`$ MeV the contribution of the mixed phase to stable stars is much smaller, restricted to the core of the star. In Fig. 6b) the properties of the stars which only contain hadronic matter do not coincide because we have taken $`m_u=m_d=5.5`$ MeV for the QMC parametrization for the results QMC plus UQM and zero otherwise.

A straightforward method of determining neutron stars properties is by measuring the gravitational redshift of spectral lines produced in neutron star photosphere which provides a direct constraint on the mass-to-radius ratio. Recently a redshift of 0.35 from three different transitions of the spectra of the X-ray binary EXO0748-676 was obtained in cottam . This redshift corresponds to $`M_{}/R(km)=0.15`$. In Fig. 6 we have added the line corresponding to this constraint and only the EOS for QMC plus CFL with B<sup>1/4</sup>=190 MeV barely satisfies this constraint. In fact, the above constraint excludes all the EOS with hyperons, quarks or obtained within a relativistic mean-field approach, as compiled in lat01 .

## VII Conclusions

In the present paper we have studied the EOS for neutron stars using both the unpaired quark matter based on the MIT Bag model and the CFL phase for describing the quark phase and a relativistic mean-field quark-meson coupling description in which quarks interact via the exchange of $`\sigma `$-, $`\omega `$\- and $`\rho `$\- mesons for the hadron phase.

We have compared the results obtained within the QMC model with the ones obtained within the NLWM. For similar properties at nuclear saturation we conclude that contrary to NLWM, with the QMC the hyperons only appear for a very high value of the bag constant. Also hyperons make the NLWM-EOS much softer than the QMC-EOS. For the hyperon meson couplings we have used three choices and verified that they did not have any effect on the onset of hyperons. The none appearance of hyperons affects the variation of the strangeness content of the star with density: except for the quark phase we get higher fractions with the NLWM.

For the bag pressure parameter $`B`$ we have used three different values, which produce different EOS and consequently the properties of the stars are dependent on them. Small $`B`$ values give stable stars with a quark core. However, for a given gap constant $`\mathrm{\Delta }`$ we obtain a phase transition to a deconfined CFL phase only if $`B`$ is greater than a critical value. For lower values we get a EOS of pure quark matter, and pure quark matter stars. For $`\mathrm{\Delta }=100`$ MeV we should have $`B^{1/4}185`$ MeV. For $`B^{1/4}=190`$ MeV and $`m_s=150`$ MeV, with QMC plus CFL we are not able to obtain stable stars with masses equal to most of radio-pulsar masses known. We have also shown that just one of the EOS we have studied satisfies the constraint imposed by the recently measured redshift of 0.35 from three different transitions of the spectra of the X-ray binary EXO0748-676cottam . A more systematic study has to be done in order to verify whether a set of acceptable parameters for the QMC+ MIT bag model gives a family of stars which contains the measured value.

In pal , strange meson fields, namely the scalar meson field $`f_0(975)`$ and the vector meson field $`\varphi (1020)`$, were also considered in order to reproduce the observed strongly attractive $`\mathrm{\Lambda }\mathrm{\Lambda }`$ interaction. They have shown that the introduction of these strange mesons makes the EOS harder due to the repulsive effect of the $`\varphi (1020)`$, meson. A harder EOS for the hadronic matter gives rise to a onset of the mixed phase at lower densities and a smaller mixed phase. The inclusion of these mesons and their influence on the properties of the stars are under investigation.

It has been shown that the effect of temperature on the maximum mass of stable stars is small compared to the effect of neutrino trapping prak97 ; recente1 ; recente2 . Therefore, it would be interesting to include neutrino trapping even at $`T=0`$ MeV and impose leptonic number conservation in the models used here and check the properties of the arising stars.

## ACKNOWLEDGMENTS

This work was partially supported by CNPq (Brazil), CAPES (visiting researcher project), CAPES(Brazil)/GRICES (Portugal) under project 003/100 and FEDER/FCT (Portugal) under the project POCTI/35308/FIS/2000.


format: text

text:
# Damping factors for the gap-tooth scheme

## 1 Introduction

For an important class of multi-scale problems, a separation of scales exists between the (microscopic, detailed) level of description of the available model, and the (macroscopic, continuum) level at which one would like to observe the system. Consider, for example, a kinetic Monte Carlo model of bacterial chemotaxis \[SGK03\]. A stochastic biased random walk model describes the probability of an individual bacterium to run or “tumble”, based on the rotation of its flagellae. Technically, it would be possible to run the detailed model for all space and time, and observe the macroscopic variables of interest, but this would be prohibitively expensive. It is known, however, that, under certain conditions, one can write a closed deterministic model for the evolution of the *concentration* of the bacteria as a function of space and time.

The recently proposed *gap-tooth scheme* \[KGK02\] can then be used instead of performing stochastic time integration in the whole domain. A number of small intervals (boxes, *teeth*), separated by large gaps, are introduced; they qualitatively correspond to mesh points for a traditional, continuum solution of the unavailable chemotaxis equation. The scheme works as follows: We first choose a number of macroscopic grid points. We then choose a small interval around each grid point; initialize the fine scale, microsopic solver within each interval consistently with the macroscopic initial conditions; and provide each box with appropriate (as we will see, to some extent artificial) boundary conditions. Subsequently, we use the microscopic model in each interval to simulate evolution until time $`\mathrm{\Delta }t`$, and obtain macroscopic information (e.g. by computing the average density in each box) at time $`\mathrm{\Delta }t`$. This amounts to a coarse time-$`\mathrm{\Delta }t`$ map; this procedure is then repeated.

The generalized Godunov scheme of E and Engquist \[EE03\] also solves an unavailable macroscopic equation by repeated calls to a microscopic code; however, the assumption is made that the unavailable equation can be written in conservation form. In the gap-tooth scheme discussed here, the microscopic computations are performed without assuming such a form for the “right-hand-side” of the unavailable macroscopic equation; we evolve the detailed model in a subset of the domain, and try to recover macroscopic evolution through interpolation in space and extrapolation in time.

We have showed analytically, in the context of numerical homogenization, that the gap-tooth scheme is close to a finite difference scheme for the homogenized equation \[SKR03\]. However, that analysis employed simulations using an algebraic constraint, ensuring that the initial macroscopic gradient is preserved at the boundary of each box over the time-step $`\mathrm{\Delta }t`$. This requires altering an existing microscopic code, so as to impose this macroscopically-inspired constraint. This may be impractical (e.g. if the macroscopic gradient has to be estimated), undesirable (e.g. if the development of the code is expensive and time-consuming) or even impossible (e.g. if the microscopic code is a *legacy code*). Generally, a given microscopic code allows us to run with a set of pre-defined boundary conditions. It is highly non-trivial to impose macroscopically inspired boundary conditions on such microscopic codes, see e.g. \[LLY98\] for a control-based strategy. This can be circumvented by introducing buffer regions at the boundary of each small box, which shield the short-time dynamics within the computational domain of interest from boundary effects. One then uses the microscopic code with its *built-in* boundary conditions.

Here, we show we can study the gap-tooth scheme (with buffers) through its numerically obtained damping factors, by estimating its eigenvalues. Integration with nearby coarse initial conditions is used to estimate matrix-vector products of the linearization of the coarse time-$`\mathrm{\Delta }t`$ map with known perturbation vectors; these are integrated in matrix-free iterative methods such as Arnoldi eigensolvers. For a standard diffusion problem, we show that the eigenvalues of the gap-tooth scheme are approximately the same as those of the finite difference scheme. When we impose Dirichlet boundary conditions at the boundary of the buffers, we show that the scheme converges to the standard gap-tooth scheme for increasing buffer size.

## 2 Physical Problem/ Governing Equations

We consider a general reaction-convection-diffusion equation with a dependence on a small parameter $`ϵ`$,

$$\frac{}{t}u(x,t)=f(u(x,t),\frac{}{x}u(x,t),\frac{^2}{x^2}u(x,t),x,\frac{x}{ϵ}),$$
(1)
with initial condition $`u(x,0)=u_0(x)`$ and Dirichlet boundary conditions $`u(0,t)=v_o`$ and $`u(1,t)=v_1`$. We further assume that $`f`$ is 1-periodic in $`y=\frac{x}{ϵ}`$.

Since we are only interested in the macroscopic (averaged) behavior, let us define an averaging operator for $`u(x,t)`$ as follows

$$U(x,t):=𝒮_h(u)(x,t)=_{x\frac{h}{2}}^{x+\frac{h}{2}}u(\xi ,t)d\xi .$$
(2)
This operator replaces the unknown function with its local average in a small box of size $`h>>ϵ`$ around each point. If $`h`$ is sufficiently small, this amounts to the removal of the microscopic oscillations of the solution, retaining its macroscopically varying components.

The averaged solution $`U(x,t)`$ satisfies an (unknown) macroscopic partial differential equation,

$$\frac{}{t}U(x,t)=F(U(x,t),\frac{}{x}U(x,t),\frac{^2}{x^2}U(x,t),x;h),$$
(3)
which does not depend on the small scale, but instead has a dependence on the box width $`h`$.

The goal of the gap-tooth scheme is to approximate the solution $`U(x,t)`$, while only making use of the detailed model (1). For illustration purposes, consider as a microscopic model the constant coefficient diffusion equation,

$$\frac{}{t}u(x,t)=a^{}\frac{^2}{x^2}u(x,t),$$
(4)
In this case both $`U(x,t)`$ and $`u(x,t)`$ satisfy (4). The microscopic and macroscopic models are the same, which allows us to focus completely on the method and its implementation.

## 3 Multiscale/Multiresolution Method

### 3.1 The gap-tooth scheme

Suppose we want to obtain the solution of the *unknown* equation (3) on the interval $`[0,1]`$, using an equidistant, macroscopic mesh $`\mathrm{\Pi }(\mathrm{\Delta }x):=\{0=x_0<x_1=x_0+\mathrm{\Delta }x<\mathrm{}<x_N=1\}`$. To this end, consider a small interval (*tooth*, box) of length $`h<<\mathrm{\Delta }x`$ centered around each mesh point, and let us perform a time integration using the microscopic model (1) in each box. We provide each box with boundary conditions and initial condition as follows.

Boundary conditions. Since the microscopic model (1) is diffusive, it makes sense to impose a fixed gradient at the boundary of each small box for a time $`\mathrm{\Delta }t`$ for the macroscopic function $`U(x,t)`$. The value of this gradient is determined by an approximation of the concentration profile by a polynomial, based on the (given) box averages $`U_i^n`$, $`i=1,\mathrm{},N`$.

$$u(x,t_n)p_i^k(x;t_n),x[x_i\frac{h}{2},x_i+\frac{h}{2}],$$

where $`p_i^k(x;t_n)`$ denotes a polynomial of (even) degree $`k`$. We require that the approximating polynomial has the same box averages in box $`i`$ and in $`\frac{k}{2}`$ boxes to the left and to the right. This gives us

$$\frac{1}{h}_{x_{i+j}\frac{h}{2}}^{x_{i+j}+\frac{h}{2}}p_i^k(\xi ;t_n)d\xi =U_{i+j}^n,j=\frac{k}{2},\mathrm{},\frac{k}{2}.$$
(5)
One can easily check that

$$𝒮_h(p_i^k)(x,t_n)=\underset{j=\frac{k}{2}}{\overset{\frac{k}{2}}{}}U_{i+j}^nL_{i,j}^k(x),L_{i,j}^k(x)=\underset{\stackrel{l={\scriptscriptstyle \frac{k}{2}}}{lj}}{\overset{\frac{k}{2}}{}}\frac{(xx_{i+l})}{(x_{i+j}x_{i+l})}$$
(6)
where $`L_{i,j}^k(x)`$ denotes a Lagrange polynomial of degree $`k`$. The derivative of $`p_i^k(x,t_n)`$ at $`x_i\pm \frac{h}{2}`$ is subsequently used as a Neumann boundary condition,

$$\frac{\mathrm{d}}{\mathrm{d}x}p_i^k(x;t_n)|_{x_i\pm \frac{h}{2}}$$
(7)
In \[SKR03\], it is shown how to enforce the macroscopic gradient to be constant when the microscopic model exhibits fast oscillations.

Initial condition. For the time integration, we must impose an initial condition $`\stackrel{~}{u}^i(x,t_n)`$ in each box $`[x_i\frac{h}{2},x_i+\frac{h}{2}]`$, at time $`t_n`$. We require $`\stackrel{~}{u}^i(x,t_n)`$ to satisfy the boundary condition and the given box average. We choose a quadratic polynomial, centered around the coarse mesh point,

$$\stackrel{~}{u}^i(x,t_n)a(xx_i)^2+b(xx_i)+c.$$
(8)
Using the constraints (7) and requiring $`\frac{1}{h}_{x_i\frac{h}{2}}^{x_i+\frac{h}{2}}\stackrel{~}{u}^i(\xi ,t_n)d\xi =U_i^n`$, we obtain

$$a=\frac{s_i^+s_i^{}}{2h},b=\frac{s_i^++s_i^{}}{2},c=U_i^n\frac{h}{24}(s_i^+s_i^{}).$$
(9)
The algorithm. The complete *gap-tooth* algorithm to proceed from $`t_n`$ to $`t_{n+1}=t_n+\mathrm{\Delta }t`$ is given below:

1. At time $`t_n`$, construct the initial condition $`\stackrel{~}{u}^i(x,t_n)`$, $`i=0,\mathrm{},N`$, using the box averages $`U_j^n`$ ($`j=0,\mathrm{},N`$) as defined in (9).
2. Compute $`\stackrel{~}{u}^i(x,t)`$ by solving the equation (1) until time $`t_{n+1}=t+\mathrm{\Delta }t`$ with Neumann boundary conditions (7).
3. Compute the box average $`U_i^{n+1}`$ at time $`t_{n+1}`$.

It is clear that this amounts to a “coarse to coarse” time-$`\mathrm{\Delta }t`$ map. We write this map as follows,

$$U^{n+1}=S_k(U^n;t_n+\mathrm{\Delta }t),$$
(10)
where $`S`$ represents the numerical time-stepping scheme for the macroscopic (coarse) variables and $`k`$ denotes the degree of interpolation.

We emphasize that the scheme is also applicable if the microscopic model is not a partial differential equation. In this case, we replace step 2 with a coarse time-stepper, based on the *lift-run-restrict* procedure that was outlined in \[GKT02\]. Numerical experiments using this algorithm are presented in \[GLK03, GK02\]. Figure 1(a) gives a schematic representation of the algorithm.

### 3.2 The gap-tooth scheme with buffers

We already mentioned that, in many cases, it is not possible or convenient to constrain the macroscopic gradient. However, the only crucial issue is that the detailed system in each box should evolve *as if it were embedded in a larger domain*. This can be effectively accomplished by introducing a larger box of size $`H>>h`$ around each macroscopic mesh point, but still only use (for macro-purposes) the evolution over the smaller, “inner” box. This is illustrated in figure 1(b). Lifting and evolution (using *arbitrary* outer boundary conditions) are performed in the larger box; yet the restriction is done by taking the average of the solution over the inner, small box. The goal of the additional computational domains, the *buffers*, is to buffer the solution inside the small box from outer boundary effects. This can be accomplished over *short enough* times, provided the buffers are *large enough*; analyzing the method is tantamount to making these statements quantitative.

The idea of using a buffer region was also used in the multi-scale finite element method (oversampling) of Hou \[HW97\] to eliminate the boundary layer effect; also Hadjiconstantinou makes use of overlap regions to couple a particle simulator with a continuum code \[Had99\]. If the microscopic code allows a choice of different types of “outer” microscopic boundary conditions, selecting the size of the buffer may also depend on this choice.

## 4 Results

We first show analytically and numerically that the standard gap-tooth scheme converges for equation (4). We then analyze convergence of the scheme with buffers and Dirichlet boundary conditions through its damping factors.

### 4.1 Convergence of the gap-tooth scheme

###### Theorem 1.

The gap-tooth scheme, applied to equation (4) with exact, analytical integration within the boxes, and boundary conditions defined through interpolating polynomials of (even) order $`k`$, is equivalent to a finite difference discretization of equation (4) of order $`k`$ central differences in space and an explicit Euler time step.

###### Proof.

When using exact (analytic) integration in each box, we can find an explicit formula for the gap-tooth time-stepper. The initial profile is given by (8), $`\stackrel{~}{u}_i(x,t_n)=a(xx_i)^2+b(xx_i)+c`$. Due to the Neumann boundary conditions, time integration can be done analytically, using

$$\stackrel{~}{u}(x,t_n+\mathrm{\Delta }t)=a(xx_i)^2+b(xx_i)+c+2aa^{}\mathrm{\Delta }t$$

Averaging this profile over the box gives the following time-stepper for the box averages,

$$U_i^n=U_i^n+a^{}\frac{s_i^+s_i^{}}{h}\mathrm{\Delta }t.$$

We know that $`s_i^\pm =\frac{\mathrm{d}}{\mathrm{d}x}p_i^k(x,t_n)|_{x_i\pm \frac{h}{2}},`$ where $`p_i^k(x,t_n)`$ is determined by (5). One can easily verify that

$$\frac{\mathrm{d}^2}{\mathrm{d}x^2}𝒮_h(p_i^k)(x,t_n)=\frac{s_i^+s_i^{}}{h},$$

is a $`k`$-th order approximation of $`\frac{^2u}{x^2}`$, which concludes the proof.∎∎

As an example, we apply the gap-tooth scheme to the diffusion equation (4) with $`a^{}=1`$. We choose an initial condition $`U(x,t)=1|2x1|`$, with Dirichlet boundary conditions, and show the result of a fourth-order gap-tooth simulation with $`\mathrm{\Delta }x=0.05`$, $`\mathrm{\Delta }t=510^3`$ and $`h=0.01`$. Inside each box, we used a second order finite difference scheme with microscopic spatial mesh size $`\delta x=110^3`$ and $`\delta t=510^7`$. The results are shown in figure 2.

### 4.2 Damping factors

Convergence results are typically established by proving consistency and stability. If one can prove that the error in each time step can be made arbitrarily small by refining the spatial and temporal mesh size, and that an error made at time $`t_n`$ does not get amplified in future time-steps, one has proved convergence. This requires the solution operator to be stable as well.

In the abscense of explicit formulas, one can examine the damping factors of the time-stepper. If, for decreasing mesh sizes, all (finitely many) eigenvalues and eigenfunctions of the time-stepper converge to the dominant eigenvalues and eigenfunctions of the time evolution operator, one expects the solution of the scheme to converge to the true solution of the evolution problem.

Consider equation (4) with Dirichlet boundary conditions $`u(0,t)=0`$ and $`u(1,t)`$, and denote its solution at time $`t`$ by the time evolution operator

$$u(x,t)=s(u_0(x);t),$$
(11)
We know that

$$s(\mathrm{sin}(m\pi x);t)=\mathrm{e}^{(m\pi )^2t}\mathrm{sin}(m\pi x),m.$$

Therefore, if we consider the time evolution operator over a fixed time $`\overline{t}`$, $`s(,\overline{t})`$, then this operator has eigenfunctions $`\mathrm{sin}(m\pi x)`$, with resp. eigenvalues

$$\lambda _m=\mathrm{e}^{(m\pi )^2\overline{t}}.$$
(12)
A good (finite difference) scheme approximates well all eigenvalues whose eigenfunctions can be represented on the given mesh. We note that it is possible to decouple the time horizon $`\overline{t}`$ from the gap-tooth (or finite-difference) time-step $`\mathrm{\Delta }t`$, in order to study the effect of different discretizations on the same reporting horizon.

Since the operator defined in (11) is linear, the numerical time integration is equivalent to a matrix-vector product. Therefore, we can compute the eigenvalues using matrix-free linear algebra techniques, even for the gap-tooth scheme, for which it might not even be possible to obtain a closed expression for the matrix. We note that this analysis gives us an indication about the quality of the scheme, but it is by no means a proof of convergence.

We illustrate this with the computation of the eigenvalues of the gap-tooth scheme with Neumann box boundary conditions. In this case, we know from theorem 1 that these eigenvalues should correspond to the eigenvalues of a finite difference scheme on the same mesh. We compare the eigenvalues of the gap-tooth scheme of order $`k=2`$ for equation (4) with diffusion coefficient $`a^{}=0.45825686`$. As method parameters, we choose $`\mathrm{\Delta }x=0.05`$, $`h=510^3`$, $`\mathrm{\Delta }t=2.510^4`$ for a time horizon $`\overline{t}=410^3`$, which corresponds to 16 gap-tooth steps. Inside each box, we use a finite difference scheme of order $`2`$ with $`\delta x=110^4`$ and an implicit Euler time-step of $`510^5`$. We compare these eigenvalues to those the finite difference scheme with $`\mathrm{\Delta }x=0.05`$ and $`\mathrm{\Delta }t=2.510^4`$, and with the dominant eigenvalues of the “exact” solution (a finite difference approximation with $`\mathrm{\Delta }x=110^3`$ and $`\mathrm{\Delta }t=110^7`$). The result is shown in figure 3. The difference between the finite difference approximation and the gap-tooth scheme in the higher modes, which should be zero according to theorem 1, is due to the numerical solution inside each box and the use of numerical quadrature for the average.

We now examine the effect of introducing a buffer region, as described in section 3.2. We consider again equation (4) with $`a^{}=0.45825686`$, and we take the gap-tooth scheme with parameters $`\mathrm{\Delta }x=0.05`$, $`h=510^3`$, $`\mathrm{\Delta }t=2.510^4`$ for a time horizon $`\overline{t}=410^3`$, and an internal time-stepper as above. We introduce a buffer region of size $`H`$, and we impose Dirichlet boundary conditions at the outer boundary of the buffer region. Lifting is done in identically the same way as for the gap-tooth scheme without buffers; we only use (9) as the initial condition in the larger box $`[x_i\frac{H}{2},x_i+\frac{H}{2}]`$. We compare the eigenvalues again with the equivalent finite difference scheme and the exact solution, for increasing sizes of the buffer box $`H`$. Figure 4 shows that, as $`H`$ increases, the eigenvalues of the scheme converge to those of the original gap-tooth scheme. We see that, in this case, we would need a buffer of size $`H=410^2`$, i.e. 80% of the original domain, for a good approximation of the damping factors. Note that it is conceptually possible, though inefficient, to let the buffer boxes overlap.

## 5 Summary/Conclusions

We described the gap-tooth scheme for the numerical simulation of multi-scale problems. This scheme simulates the macroscopic behaviour over a macroscopic domain when only a microscopic model is explicitly available. In the case of diffusion, we showed equivalence of our scheme to standard finite differences of arbitrary (even) order, both theoretically and numerically.

We showed that it is possible, even without analytic formulas, to study the properties of the gap-tooth scheme and generalizations through the damping factors of the resulting coarse time-$`\mathrm{\Delta }t`$ map. We illustrated this for the original gap-tooth scheme and for an implementation using Dirichlet boundary conditions in a buffer box. We showed that, as long as the buffer region is “large enough” to shield the internal region from the boundary effects over a time $`\mathrm{\Delta }t`$, we get a convergent scheme. Therefore, we are able to use microscopic codes in the gap-tooth scheme *without modification*.

In a forthcoming paper, we will explore, using these damping factors for many different types of boundary conditions, the relation between the quality of the boundary conditions, the size of the buffer region and the time-step before reinitialization. We will investigate the trade-off between the effort required to impose a particular type of boundary conditions (and the eventual macroscopically inspired control-based strategy) and the efficiency gain due to smaller buffer sizes and/or longer possible time-steps before reinitialization. Here, we showed that this investigation is made possible by studying the damping factors of the resulting coarse time-$`\mathrm{\Delta }t`$ map.

Acknowledgements GS is a Research Assistant of the Fund of Scientific Research - Flanders. This work has been partially supported by an IUAP grant and by the Fund of Scientific Research through Research Project G.0130.03 (GS, DR), and by the AFOSR and the NSF (IGK). The authors thank Olof Runborg for discussions that improved this text and the organizers of the Summer School in Multi-scale Modeling and Simulation in Lugano.

## Index


text:
# A possible route to grow a (Mn:Si(1-x)Gex)-based diluted magnetic semiconductor

## Abstract

We systematically investigate, using ab initio density-functional theory calculations, the properties of interstitial and substitutional Mn in both Si and Ge, as well as in the Si1-xGex alloy. We show that volume effects are not the main reason Mn prefers to be a subsitutional impurity in pure Ge, and chemical effects, therefore, play an important role. Using realistic models of Si1-xGex, we show that for $`x0.16`$ substitutional Mn in Ge-rich neighborhoods become more stable than interstitial Mn, which may allow the growth of Si-based diluted magnetic semiconductors.

Diluted magnetic semiconductor (DMS) materials dietl ; science-dietl ; ohno1 ; ohno2  have been know for a long time book , however, it was the relatively recent growth of III-V ohno1 ; III-V  Mn-based ferromagnetic DMS that have brought these materials to the spotlight, motivated by the possibility of fabricating useful spintronic devices. Although these accomplishments are rather significant, if one considers the possible technological impact there is no doubt that the synthesis of a Si-based DMS material would be of great relevance. Ferromagnetism in the MnxGe1-x compound has been reported by more than one group mnge , which makes one ask the obvious question: would it be possible to grow a similar MnxSi1-x ferromagnetic material?

Considering this question in a recent paper nosso , we have argued that: (i) as Si and Ge have similar electronic structures, it is reasonable to expect that for large enough $`x`$, MnxSi1-x would become ferromagnetic below some critical temperature, in the same way MnxGe1-x does it. This hypothesis is supported by a recent calculation freeman  where both MnxGe1-x and MnxSi1-x were studied. The question then is if the required values of $`x`$ can be achieved; (ii) in Si, a Mn impurity favors an interstitial site, whereas in Ge it prefers a substitutional site (see below); (iii) as a consequence, a Mn substitutional impurity in Ge cannot diffuse as easily as an interstitial Mn in Si woodbury ; difmnsi ; fastdiff , allowing the introduction of a large enough number of impurities without their diffusion and subsequent clustering. Therefore, it seems that a crucial issue to obtain a MnxSi1-x material is to find a way to have the Mn occupying substitutional sites. We have previously found nosso  that there are substitutional sites at the Si(100) surface with formation energies comparable to interstitial sites, opening up the possibility of placing MnSi through some non-equlibrium growth process, akin to the low-temperature molecular beam epitaxy procedures that are crucial to obtain III-V based DMSs III-V .

In the present paper we address the following questions: (i) Is it simply a volume effect that would make Mn prefer a tetrahedral interstitial site in Si whereas it prefers a substitutional site in Ge? (ii) Will a Ge rich environment in the SiGe alloy also favor a substitutional Mn site? Regarding question (i), we find that the cause Mn is substitutional in Ge is not simply a volume effect, and therefore finding ways to increase the Si lattice parameter would not suffice to stabilize the substitutional Mn impurity. However, since specific chemical interactions are the possible explanation for the differences between Mn in Si and Ge, this implies that the answer to question (ii) above is yes. We actually find that the substitutional Mn site shoud be the predominant one for Ge concentrations in the Si1-xGex alloy for $`x0.16`$. These findings have the important consequence that a Mn doped SiGe alloy may present similar ferromagnetic properties as the MnxGe1-x material, opening up the possibility of Si-based spintronics.

All our results are based on ab initio calculations, based on spin-polarized density functional theory within the generalized gradient approximation (GGA) gga . We have used ultrasoft pseudopotentials vanderbilt  and a plane wave expansion up to 230 eV, as implemented in the VASP code vasp . We have used a FCC-based supercell containing 128 sites, and the $`L`$-points for the Brillouin zone sampling. In all calculations the positions of all atoms in the supercell were relaxed until all the forces components were smaller than 0.02 eV/Å.

We initially considered pure Si and Ge crystals. We placed a Mn atom both at an interstitial as well as at a substitutional site. The formation energy of a neutral interstitial impurity neutral , $`E_f^I`$, is calculated as

$$E_f^I=E_{def}E_{bulk}\mu _{Mn}$$
(1)
where $`E_{def}`$ is the total energy of the supercell with the defect, $`E_{bulk}`$ is the total energy of the supercell without the defect (either pure Si or pure Ge), and $`\mu _{Mn}`$ is the Mn chemical potential chempot . For neutral substitutional impurities neutral  the formation energy ($`E_f^S`$) is given by

$$E_f^S=(E_{def}+\mu _X)E_{bulk}\mu _{Mn}$$
(2)
where $`\mu _X`$ is the chemical potential of either Si or Ge chempot .

In Table 1 we present the results for Mn atoms in Si and Ge at their equilibrium lattice constants (we find $`a_{Si}=5.445`$ Å for Si and $`a_{Ge}=5.750`$ Å for Ge). As expected difmnsi ; fastdiff , for silicon the Mn atoms prefer to be at interstitial sites, with an energy difference $`\mathrm{\Delta }E^{SI}`$ between the substitutional and interstitial sites of $`\mathrm{\Delta }E^{SI}=0.3`$ eV. For Mn in Ge, on the other hand, we obtain that the Mn substitutional impurity has a lower formation energy by $`\mathrm{\Delta }E^{SI}=0.6`$ eV. In order to investigate how much of this difference between Si and Ge is a result of their lattice parameter difference, i.e., how much is due to a volume effect, we repeated the same calculations for Si (Ge) in the lattice parameter of Ge (Si) chempot . The results are also presented in Table 1. In silicon, no significant changes were observed, with a small reduction in the substitutional Mn formation energy, which led to $`\mathrm{\Delta }E^{SI}=0.2`$ eV instead of 0.3 eV. For Mn in Ge, the interstitial site becomes highly unfavorable, with the formation energy changing by more than 1 eV. The formation energy for the substitutional site, on the other hand, changed very little. This leads to a large increase in the energy difference between the substitutional and interstitial sites: $`\mathrm{\Delta }E^{SI}=1.7`$ eV. Therefore, even though increasing the Si lattice parameter has a small tendency towards favoring a Mn atom in a substitutional site, the difference in lattice parameter between Si and Ge cannot account, by itself, for the distinct behavior of Mn in these materials. Thus, specific chemical interactions between the Mn atom and the host atoms are crucial in determining the distinct behavior in Si and Ge.

This suggests that in the Si1-xGex alloy the Mn atoms may prefer to be at a substitutional site with a Ge rich environment, instead of at an interstitial site. In order to confirm this possibility, we performed calculations for both interstitial as well as substitutional Mn in Si1-xGex, for x=0.25, 0.5 and 0.75. Since the Si1-xGex is a substitutionally disordered alloy, the vicinity of a Mn impurity is not uniquely determined, and for either the substitutional or the interstitial Mn, there are five different types of sites, if only the nearest-neighborhood is considered, i.e., a Mn surrounded by a configuration of Si and Ge atoms that can be labeled as SiνGe4-ν, for $`\nu `$ varying from 0 to 4.

To simulate the Si1-xGex disordered alloys, we used supercells with 128 atoms where the atoms were distributed as Special Quasirandom Structures (SQS) Wei . Details of the preparation of the supercells were given elsewhere liga1 , and a similar procedure has already been used to study vacancies in this alloy liga2 ; liga3 . It should be mentioned that due to the SQS approach, the disorder of the alloy is taken into account in an explicit way. As we have shown before liga1 , the alloy lattice parameter follows very closely the Vegard’s law, and we therefore use $`a_{Si_{1x}Ge_x}=(1x)a_{Si}+xa_{Ge}`$.

To study the interstitial Mn defects, we randomly selected five sites with distinct first neighborhoods. As these sites will also have different second, third, etc., neighborhoods, an averaging procedure should be performed. From previous studies of vacancies in the alloy liga2 ; liga3 , we estimate an overall variation in the formation energies of the order of $`\pm `$ 0.1 eV due the different vicinities, and hence this averaging will not alter our conclusions in any significative way. The formation energy $`E_f^I(\nu ,x)`$ of a neutral interstitial neutral  Mn in the Si1-xGex alloy, in a given neighborhood SiνGe4-ν, labeled by $`\nu `$, is given by

$$E_f^I(\nu ,x)=E_I(\nu ,x)E_{bulk}(x)\mu _{Mn}(x).$$
(3)
Here $`E_I(\nu ,x)`$ is the total energy of the SQS structure with the Mn interstitital, $`E_{bulk}(x)`$ is the total energy of the SQS alloy without any defects, and $`\mu _{Mn}(x)`$ is the Mn chemical potential in the Si1-xGex alloy, which we determine chempot2  as

$$\mu _{Mn}(x)=(1x)\mu _{MnSi}+x\mu _{MnGe}\mu _{Si_{(1x)}Ge_x}.$$
(4)
To study substitutional Mn defects a similar procedure as described above was employed. Since an atom from the original SQS structure must now be removed, there is the extra possibility of having the Mn replacing either a Si or a Ge atom. The formation energy $`E_f^S(X,\nu ,x)`$ of a neutral neutral  Mn substituting an atom $`X`$ ($`X`$=Si,Ge), in the Si1-xGex alloy, with a neighborhood SiνGe4-ν, labeled by $`\nu `$, is given by

$$E_f^S(X,\nu ,x)=E_S(X,\nu ,x)E_{bulk}(x)+\mu _X(x)\mu _{Mn}(x).$$
(5)
Here $`E_S(X,\nu ,x)`$ is the total energy of the SQS structure with the Mn in place of an atom $`X`$ ($`X`$=Si,Ge), which has a chemical potential $`\mu _X(x)`$ in the Si1-xGex alloy chempot3 .

All our results are presented in Table 2 and Fig. 1. For the substitutional Mn calculations, we always removed a Si atom to place the Mn impurity, except for $`x=0.5`$ where a Ge atom was also removed. These latter results are shown in parethesis in Table 2. As can be seen, the formation energies are always very similar, indicating that our conclusions do not depend on Mn replacing either a Si or a Ge atom. For interstitial Mn, the differences between the formation energies for the local configurations Si0Ge4 and Si4Ge0 are always 0.4 eV, for all $`x`$, with the lowest values occuring in Si-rich vicinities (Si4Ge0). For substitutional Mn this formation-energy spread is also independent of $`x`$, and has a value of 0.3 eV. Ge-rich neighborhoods, however, have lower formation energies. An important result is the fact that the lowest formation energy for a substitutional Mn is always smaller than the lowest formation energy for an interstitial Mn. The difference between these two lowest formation energies is 0.1 eV for $`x=0.25`$, 0.4 eV for $`x=0.5`$, and 0.5 eV for $`x=0.75`$. This indicates that already for the Si0.75Ge0.25 alloy, there are local configurations that make substitutional Mn the lowest energy structure. Another significant result is the fact that for $`x0.5`$, all substitutional configurations have formation energies that are smaller than the smallest formation energy for interstitial Mn. The overall picture of our results can be appreciated in Fig. 1. Curve fittings to both the average values (continuous curves) of the formation energies as well as to the lowest values (long-dashed curves) of the formation energies are also presented in Fig. 1. The average values curves cross at $`x0.14`$, whereas the lowest values curves cross at $`x0.18`$. These results indicate that for $`x0.16\pm 0.02`$ there should be Ge-rich neighborhoods in Si(1-x)Gex where substitutional Mn atoms becomes more stable than interstitial Mn.

At thermodynamic equilibrium, and assuming that the defects are independent of each other, the population of Mn interstitial defects in Si1-xGex is given by liga2

$$N^I(\nu )=N_SP^\nu (x)exp(G_f^I(\nu )/k_BT),$$
(6)
where $`N_S`$ is the total number of sites in the lattice, $`P^\nu (x)`$ is the probability prob  for the $`\nu `$ configuration to exist in the $`x`$-concentration alloy, $`G_f^I(\nu )`$ is the Gibbs free formation energy of interstitial Mn at configuration $`\nu `$, and $`T`$ is the temperature. Note that $`N_SP^\nu (x)`$ is the effective number of interstitial sites of type $`\nu `$. Under the same assumptions, the population of Mn substitutional defects in Si1-xGex are given by

$$N^S(\nu )=N_SP^\nu (x)f(X,x)exp(G_f^S(X,\nu )/k_BT),$$
(7)
where $`f(X,x)=(1x)`$ if $`X`$=Si and $`f(X,x)=x`$ if $`X`$=Ge, and $`G_f^S(X,\nu )`$ is the Gibbs free formation energy of substitutional Mn at configuration $`\nu `$, with Mn replacing an atom $`X`$ ($`X`$=Si,Ge). Assuming that, for a given temperature and composition, the vibrational entropies of all defects are similar and independent on their vicinities, we can calculate the relative population of the Mn defects as

$$n^D(\nu ,x,T)=\frac{P^\nu (x)exp[E_f^D(\nu )/k_BT]}{_{\nu ,D}P^\nu (x)exp[E_f^D(\nu )/k_BT]}.$$
(8)
where $`D=I,S`$ represents both the interstitial and substitutional defects. For substitutional defects, since the formation energies for Mn replacing either a Ge or a Si atom are quite similar, and given that $`_Xf(X,x)=1`$, the above expression for $`n^S(\nu ,x,T)`$ is obtained after a summation over $`X`$ is performed.

In Table 3 we present the relative populations for the interstitial and substitutional Mn impurities in the Si(1-x)Gex alloy, for $`x=0.25,0.5`$ and 0.75, for two temperatures, T=300 K and T=600 K. In all cases the overall population of substitutional Mn is larger than the interstitial one. At 300 K, even for $`x=0.25`$ we already have approximately 78 % of Mn at substitutional sites. As the temperature is increased to 600 K, this percentage decreases to 56 %. The significance of this result can be appreciated by noting that in bulk Si, at T=300 K, basically 100 % of Mn are at interstitial sites, whereas at T=600 K there are approximately 99.7 % of interstitial Mn. This means that by alloying with Ge one can revert the overall population of defects from interstitial to substitutional Mn.

In conclusion, we have shown that volume effects cannot be the main reason Mn prefers to be a subsitutional impurity in pure Ge. Chemical effects, therefore, must clearly play an important role. Through a systematic study of interstitial and substitutional Mn atoms in realistic models of the Si(1-x)Gex alloy, we have shown that for $`x0.16`$ the substitutional Mn in Ge-rich neighborhoods becomes more stable than the interstitial Mn. By playing with the temperature and the concentration $`x`$, and maybe also using non-equilibrium growth conditions nosso , it should be possible to obtain Si(1-x)Gex alloys with a large enough concentration of substitutional Mn atoms. Considering that in all calculations in the alloy we obtained a Mn local moment with the same value as in the pure crystals, and given that recent studies freeman  have shown that both MnxSi1-x and MnxGe1-x should have similar magnetic properties, all these facts indicate that Mn:Si(1-x)Gex is potentially a magnetic material like MnxGe1-x mnge , opening up in this way the road towards Si-based DMS.

This research was supported by the agencies FAPESP and CNPq. We thank CENAPAD-SP for computer time. We also acknowledge G. M. Dalpian for useful discussions.


text:
# Monitoring The Evolution of SNR 1987A With Chandra

## 1. Introduction

We continue our monitoring program of SNR 1987A with the Advanced CCD Imaging Spectrometer (ACIS) on board the Chandra X-Ray Observatory. As of 2002 December, we have performed a total of seven observations (Table 1). Results from the first six observations have been reported in the literature (Burrows et al. 2000; Park et al. 2002, 2003; Michael et al. 2002). The X-ray morphology was ring-like and asymmetric (brighter in the east) with the emergence of X-ray-bright spots. The X-ray spectrum was described with a plane-parallel shock in a non-equilibrium ionization (NEI) state with an electron temperature of $`kT`$ $``$ 2.5 keV. A blast wave shock velocity of $`v`$ $``$ 3500 km s-1 was derived from the detected X-ray line profiles, which also provided direct evidence of an electron-ion non-equilibrium behind the shock. The X-ray flux was non-linearly increasing as the blast wave approaches the dense inner ring. We here report on the results from the lastest Chandra observations of SNR 1987A.

## 2. Morphology

The ACIS images (Fig. 1) have been generated by the data reduction process described in Burrows et al. (2000) and Park et al. (2002). Fig. 1 shows the overall brightening and development of new X-ray bright spots, particularly in the western side of the SNR. SNR 1987A is expected to be a complete ring in X-rays. In 2000, the strong correlation between the soft X-ray and the optical images was interpreted as emission from the slow shock encountering protrusions of the dense inner ring. The correlation between the hard X-ray and the radio images was consistent with emission from the fast shock propagating into an HII region interior to the ring. As of 2002-12, correlations between the X-ray and the optical/radio images are more complex than such a simple picture, which is perhaps expected as the blast wave is reaching the main body of the inner ring.

## 3. Spectrum and Flux

The X-ray spectrum from each observation can be fitted with a single NEI plane shock model with an electron temperature of $`kT`$$``$2$``$3 keV and low abundances (0.1$``$0.5 solar). Based on the data with the most significant photon statistics, we find that the electron temperature has decreased slightly, while the volume emission measure (EM) has substantially increased for the last two years (Fig. 2). These results are consistent with our physical picture of the blast wave encountering a dense circumstellar medium (CSM). A growing contribution to the observed spectrum from the slow shock is then expected (Michael et al. 2002).

A two-temperature model indeed statistically improves the fit, particularly with the latest data (Fig. 3). The soft component ($`kT`$$``$0.25 keV) implies a slow shock velocity ($`v`$$``$400 km s-1) with a highly advanced ionization state ($`n_et`$$``$1013 cm-3 s). This is consistent with the decelerated shock front in the dense CSM. The 0.5$``$5 keV band flux from the soft component has tripled since 2000-12, while the flux from the hard component has doubled. We note that the flux increase from the slow shock in the western half (by a factor of 6) is much larger than in the east ($``$30%). The flux increase from the fast shock is constant between the east and the west. These results suggest that the blast wave has now begun to encounter the dense CSM in the western side of the SNR a few years after it did in the eastern side. The X-ray lightcurve shows a non-linear flux increase rate (Fig. 4). We have thus attempted non-linear fits to the data: e.g., in Fig. 4, we present a fit by assuming an exponential density distribution along the radius of the inner ring. The best-fit model implies an $``$20 times higher density in the inner ring than in the HII region. The actual density contrast, which might be higher than this, could be measured as the shock-CSM interaction proceeds.

## References

Burrows, D. N., Michael, E., Hwang, U., McCray, R., Chevalier, R. A., Petre, R., Garmire, G. P., Holt, S. S., & Nousek, J. A. 2000, ApJ, 543, L149

Hasinger, G., Aschenbach, B., & Trümper, J. 1996, A&A, 312, L9

Park, S., Burrows, D. N., Garmire, G. P., Nousek, J. A., McCray, R., Michael, E., Zhekov, S. 2002, ApJ, 567, 314

Park, S., Burrows, D. N., Zhekov, S. A., Michael, E., McCray, R., & Hasinger, G. 2003, To appear in Adv. in Space Res, the proceedings of the 34th COSPAR Synposium E1.4 “High Energy Studies of Supernova Remnants and Neutron Stars”, ed. W. Becker and W. Hermsen, in press

Michael, E., Zhekov, S., McCray, R., Hwang, U., Burrows, D. N., Park, S., Garmire, G. P., Holt, S. S., Hasinger, G. 2002, ApJ, 574, 166


text:
# Observation of 𝐾_𝑠⁰𝐾_𝑠⁰ resonances in deep inelastic scattering at HERA

## 1 Introduction

The $`K_s^0K_s^0`$ system is expected to couple to scalar and tensor glueballs. This has motivated intense experimental and theoretical study during the past few years rmp:v71:1411 ; hep-ex/0101031 . Lattice QCD calculations glueball1 ; glueball2  predict the existence of a scalar glueball with a mass of $`1730\pm 100`$ MeV and a tensor glueball at $`2400\pm 120`$ MeV. The scalar glueball can mix with $`q\overline{q}`$ states with $`I=0`$ from the scalar meson nonet, leading to three $`J^{PC}=0^{++}`$ states whereas only two can fit into the nonet. Experimentally, four states with $`J^{PC}=0^{++}`$ and $`I=0`$ have been established pdg : $`f_0(980)`$, $`f_0(1370)`$, $`f_0(1500)`$ and $`f_0(1710)`$.

The state most frequently considered to be a glueball candidate is $`f_0(1710)`$ pdg  , but its gluon content has not yet been established. This state was first observed in radiative $`J/\psi `$ decays BES96  and its angular momentum $`J=0`$ was established by the WA102 experiment using a partial-wave analysis in the $`K^+K^{}`$ and $`K_s^0K_s^0`$ final statesWAf1710  . Observation of $`f_0(1710)`$ in $`\gamma \gamma `$ collisions may indicate a large quark content. A recent publication from L3 L3  reports the observation of two states in $`\gamma \gamma `$ collisions above $`1500`$ MeV, the well-established $`f_2^{}(1525)`$ pdg  and a broad resonance at 1760 MeV. It is not clear if the latter state is the $`f_0(1710)`$.

The $`ep`$ collisions at HERA provide an opportunity to study resonance production in a new environment. Production of $`K_s^0`$ particles has been studied previously at HERA epj:c2:77 ; zfp:c68:29 ; np:b480:3 . In this contribution, the first observation of resonances in the $`K_s^0K_s^0`$ final state in inclusive deep inelastic $`ep`$ scattering (DIS) is reported hep-0308006 .

## 2 Event selection and $`K_s^0`$-pair candidates

A detailed description of the ZEUS detector can be found elsewhere zeus:1993:bluebook .

The data used for this study correspond to a total integrated luminosity of 120 pb<sup>-1</sup> collected in ZEUS during the 1996-2000 running period.

The inclusive neutral current DIS process $`e(k)+p(P)e(k^{})+X`$ can be described in terms of the following variables: the negative of the invariant-mass squared of the exchanged virtual photon, $`Q^2=q^2=`$ $`(kk^{})^2`$; the fraction of the lepton energy transferred to the proton in the proton rest frame, $`y=(qP)/(kP)`$; and the Bjorken scaling variable, $`x=Q^2/(2Pq)`$.

A three-level trigger system was used to select events online zeus:1993:bluebook . The inclusive DIS selection was defined by requiring an electron found in the Uranium Calorimeter, and further requirements were applied to ensure a well defined data sample hep-0308006 .

Oppositely charged track pairs reconstructed by the ZEUS central tracking detector (CTD) and assigned to a secondary vertex were selected and combined to form $`K_s^0`$ candidates. Both tracks were assigned the mass of a charged pion and the invariant-mass $`M(\pi ^+\pi ^{})`$ was calculated. Only events with at least one pair of $`K_s^0`$ candidates were selected. The invariant mass of the $`K_s^0`$ pair candidate $`M(K_s^0,K_s^0)`$ was reconstructed in the range $`0.995<M(K_s^0K_s^0)<2.795`$ GeV. A detailed description of the $`K_s^0`$ pair candidate selection can be found in hep-0308006 .

Figure 1 shows the distribution in $`x`$ and $`Q^2`$ of selected events containing at least one pair of $`K_s^0`$ candidates.

Figure 2 shows the $`M(\pi ^+\pi ^{})`$ distribution in the range $`0.45<M(\pi ^+\pi ^{})<0.55`$ GeV after the $`K_s^0`$ pair candidate selection.

## 3 Results

The $`K_s^0K_s^0`$ spectrum may have a strong enhancement near the $`K_s^0K_s^0`$ threshold due to the $`f_0(980)`$/$`a_0(980)`$ state prd:v32:189 ; pl:b489:24 ; hep-ph/0202157 . Since the high $`K_s^0K_s^0`$ mass is the region of interest for this analysis, the complication due to the threshold region is avoided by imposing the cut $`cos\theta _{K_s^0K_s^0}<0.92`$, where $`\theta _{K_s^0K_s^0}`$ is the opening angle between the two $`K_s^0`$ candidates in the laboratory frame.

After applying all selections, 2553 $`K_s^0`$-pair candidates were found in the range $`0.995<M(K_s^0K_s^0)<2.795`$ GeV, where $`M(K_s^0K_s^0)`$ was calculated using the $`K_s^0`$ mass of $`497.672`$MeV pdg . The momentum resolution of the CTD leads to an average $`M(K_s^0K_s^0)`$ resolution which ranges from $`7`$ MeV in the $`1300`$ MeV mass region to 10 MeV in the $`1700`$ MeV region. Figure 3 shows the measured $`K_s^0K_s^0`$ invariant-mass spectrum. Two clear peaks are seen, one around 1500 MeV and the other around 1700 MeV, along with an enhancement around 1300 MeV. The data for $`\mathrm{cos}\theta _{K_s^0K_s^0}>0.92`$ are also shown.

The distribution of Fig. 3 was fitted using three modified relativistic Breit-Wigner (MRBW) distributions and a background function $`U(M)`$;

$$F(M)=\underset{i=1}{\overset{3}{}}(\frac{m_{,i}\mathrm{\Gamma }_{d,i}}{(m_{,i}^2M^2)^2+m_{,i}^2\mathrm{\Gamma }_{d,i}^2})+U(M),$$
(1)
where $`\mathrm{\Gamma }_{d,i}`$ is the effective resonance width, which takes into account spin and large width effects bwigner , $`m_{,i}`$ is the resonance mass, and $`M`$ is $`K_s^0K_s^0`$ invariant mass. The background function is

$$U(M)=A(M2m_{K_s^0})^Be^{C\sqrt{M2m_{K_s^0}}},$$
(2)
where $`A`$, $`B`$ and $`C`$ are free parameters and $`m_{K_s^0}`$ is the $`K_s^0`$ mass defined by Hagiwara et al pdg . Monte Carlo studies showed that effects of the track-momentum resolution on the mass reconstruction were small compared to the measured widths of the states. Therefore the resolution effects were ignored in the fit.

Below $`1500`$ MeV, a region strongly affected by the $`\mathrm{cos}\theta _{K_s^0K_s^0}`$ cut, a peak is seen around $`1300`$ MeV where a contribution from $`f_2(1270)/a_2^0(1320)`$ is expected. This mass region was fitted with a single Breit-Wigner.

Above $`1500`$ MeV, the lower-mass state has a fitted mass of $`1537_8^{+9}`$ MeV and a width of $`50_{22}^{+34}`$ MeV, in good agreement with the well established $`f_2^{}(1525)`$. The higher-mass state has a fitted mass of $`1726\pm 7`$ MeV and a width of $`38_{14}^{+20}`$ MeV. The widths reported here were stable, within statistical errors, to a wide variation of fitting methods including those using 30 MeV bins rather than the default 15 MeV bins. The width is narrower than the PDG value ($`125\pm 10`$ MeV) pdg  reported for $`f_0(1710)`$, but when it is fixed to this value, the fit is still acceptable.

It was found that 93% of the $`K_s^0`$-pair candidates selected within the detector and trigger acceptance are in the target region of the Breit frame bk:1972:1 ; zerw:zfp:237 , the hemisphere containing the proton remnant. Of the $`K_s^0`$-pair candidates in the target region, 78% are in the region $`x_p=2p_B/Q>1`$, where $`p_B`$ is the absolute momentum of the $`K_s^0K_s^0`$ in the Breit frame. High $`x_p`$ corresponds to production of the $`K_s^0`$-pair in a region where sizeable initial state gluon radiation may be expected. This is in contrast to the situation at $`e^+e^{}`$ colliders where the particles entering the hard scattering are colourless.

## 4 Conclusions

The first observation in $`ep`$ deep inelastic scattering of a state at $`1537`$ MeV, consistent with $`f_2^{}(1525)`$, and another at 1726 MeV, close to $`f_0(1710)`$, is reported. There is also an enhancement near 1300 MeV which may arise from the production of $`f_2(1270)`$ and/or $`a_2^0(1320)`$ states. The width of the state at $`1537`$ MeV is consistent with the PDG value for the $`f_2^{}(1525)`$. The state at 1726 MeV has a mass consistent with the glueball candidate $`f_0(1710)`$, and is found in a gluon-rich region of phase space. However, it’s width of $`38_{14}^{+20}`$ MeV is narrower than the PDG value of 125$`\pm `$10 MeV for the $`f_0(1710)`$.