text:
# Borel–Padé vs Borel–Weniger method: a QED and a QCD example

## I Introduction

In this letter we want to present some numerical results which allow us to compare the efficiency of the Borel–Padé method with that of the Borel–Weniger method for resummation of truncated perturbation series (TPS) in some physically significant scenarios. The scenarios we are referring to are those when the function, which we want to find through a resummation, is known to have certain singularity structure in the Borel plane. If there are singularities on the positive axis of the Borel plane, then we implicitly assume that in such cases we either know the correct prescription for integration in the Laplace–Borel integral, or we simply adhere to a certain adopted prescription.

(1) We will first illustrate the efficiency of the two methods on the QED example of the Euler–Heisenberg Lagrangian density, i.e., the one–loop fermion–induced effective action density in a strong uniform electromagnetic field . In this case, the solution is known, and its real part can be written in the following form:

$`\mathrm{Re}\delta \stackrel{~}{}(\stackrel{~}{a};p)`$ $`=`$ $`\mathrm{Re}{\displaystyle _0^{\mathrm{}}}𝑑w\mathrm{exp}\left({\displaystyle \frac{w}{\stackrel{~}{a}}}\right){\displaystyle \frac{(1)}{w}}\left[{\displaystyle \frac{p\mathrm{cos}(w)}{\mathrm{sin}(w+\mathrm{i}ϵ^{})}}\mathrm{coth}(pw)+{\displaystyle \frac{1}{3}}(1p^2){\displaystyle \frac{1}{w^2}}\right],`$ (1)
where we use notations

$$\stackrel{~}{a}\frac{ga}{m^2},\stackrel{~}{b}\frac{gb}{m^2},p\frac{b}{a}\frac{\stackrel{~}{b}}{\stackrel{~}{a}},\delta \stackrel{~}{}\delta /\left(\frac{m^4\stackrel{~}{a}^2}{8\pi ^2}\right).$$
(2)
Here, $`\delta `$ is the actual Lagrangian density induced by the one–loop fluctuations of the fermions in the field; $`g`$ is the field–to–fermion coupling parameter (in QED it is the positron charge $`e_0`$); $`m`$ is the mass of the fermion (electron); $`a`$ and $`b`$ are Lorentz–invariant expressions characterizing the electric and the magnetic fields $`\stackrel{}{}`$ and $`\stackrel{}{}`$, respectively

$$\left(\genfrac{}{}{0pt}{}{a}{b}\right)=\left[\pm \stackrel{}{}^2\stackrel{}{}^2+\sqrt{\left(\stackrel{}{}^2\stackrel{}{}^2\right)^2+4\left(\stackrel{}{}\stackrel{}{}\right)^2}\right]^{1/2}/\sqrt{2}.$$
(3)
Expression (1) can be obtained, for example, directly by integrating out the fermionic degreees of freedom in the path integral expression of the full effective action, then employing the proper–time integral representation for the difference of logarithms, evaluating the traces in the integrand, and subsequently performing Wick rotation by $`\pi /4`$ in the plane of the proper–time $`s`$: $`ags\mathrm{i}w+ϵ^{}`$. We refer to  for more details on the latter point. The perturbative expansion of the full solution (1), in powers of $`\stackrel{~}{a}`$, is

$$\delta \stackrel{~}{}^{\mathrm{pert}.}(\stackrel{~}{a};p)=\left[c_1(p)1!\stackrel{~}{a}^2+c_3(p)3!\stackrel{~}{a}^4+c_5(p)5!\stackrel{~}{a}^6+\mathrm{}\right],$$
(4)
with coefficients

$`c_1(p)`$ $`=`$ $`{\displaystyle \frac{1}{45}}\left[(1p^2)^2+7p^2\right],c_3(p)={\displaystyle \frac{1}{945}}\left[2(1p^2)^3+13p^2(1p^2)\right],\mathrm{etc}.`$ (5)
In the case of the pure magnetic field (p.m.f.), the corresponding expressions are simpler

$`\delta \stackrel{~}{}(\stackrel{~}{b})_{\stackrel{~}{a}=0}{\displaystyle \frac{8\pi ^2\delta _{a=0}}{m^4\stackrel{~}{b}^2}}={\displaystyle _0^{\mathrm{}}}𝑑w\mathrm{exp}\left({\displaystyle \frac{w}{\stackrel{~}{b}}}\right){\displaystyle \frac{(1)}{w}}\left[{\displaystyle \frac{\mathrm{coth}(w)}{w}}{\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{w^2}}\right],`$ (6)
$$\delta \stackrel{~}{}^{\mathrm{pert}.}(\stackrel{~}{b})_{\stackrel{~}{a}=0}=\left[\stackrel{~}{c}_11!\stackrel{~}{b}^2+\stackrel{~}{c}_33!\stackrel{~}{b}^4+\mathrm{}\right],\stackrel{~}{c}_1=\frac{1}{45},\stackrel{~}{c}_3=\frac{2}{945},\mathrm{}$$
(7)
We can now use (4)–(5), and (7), as a laboratory for resummation methods, since the full (resummed) solutions (1) and (6) are known. Since (1) and (6) are Laplace–Borel integrals, it is natural to use these examples for testing combined resummation techniques which involve Borel transformation. Borel transform $`B_L`$ of series (4) is

$$B_L(w;p)=c_1(p)w+c_3(p)w^3+c_5(p)w^5+\mathrm{},$$
(8)
and analogously for (7). In Ref. , we used Borel–Padé technique for resummation, i.e., we applied various Padé approximants $`[N/M]_B(w;p)`$ to (8)333 $`[N/M]_B(w;p)`$, being ratio of polynomials in $`w`$ of powers $`N`$ and $`M`$, respectively , is based solely on the truncated perturbation series (TPS) of (8) involving only terms with $`c_n`$: $`nN+M`$. and then employed the Laplace–Borel integral to obtain the resummed value

$$BP^{[\mathrm{N}/\mathrm{M}]}\left[\delta \stackrel{~}{}^{\mathrm{pert}.}\right](\stackrel{~}{a};p)=_0^{\mathrm{}}𝑑w\mathrm{exp}\left(\frac{w}{\stackrel{~}{a}}\right)[N/M]_\mathrm{B}(w;p).$$
(9)
The integration over poles in (9) was carried out according to the Cauchy principal value prescription, since the full solution (1) requires it.444 Various QCD and QED applications of the Borel–Padé approach with the principal value prescription can be found in . The novel method of Ref.  is, in addition, well suited for obtaining the imaginary part of $`\delta `$.

Recently, the authors of  proposed the use of Weniger (delta sequence) transformations as an alternative to the use of Padé approximants, for direct resummation of truncated perturbation series. For a truncated perturbation series (TPS) of the form

$$F_{[n+1]}(z)=\underset{0}{\overset{n+1}{}}\gamma _jz^j$$
(10)
it is defined as

$$\delta _n^{(0)}(\zeta ;\gamma _0,\mathrm{},\gamma _{n+1})=\frac{{\displaystyle \underset{j}{\overset{n}{}}}(1)^j\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right){\displaystyle \frac{(\zeta +j)_{n1}}{(\zeta +n)_{n1}}}{\displaystyle \frac{z^{nj}F_{[j]}(z)}{\gamma _{j+1}}}}{{\displaystyle \underset{j}{\overset{n}{}}}(1)^j\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right){\displaystyle \frac{(\zeta +j)_{n1}}{(\zeta +n)_{n1}}}{\displaystyle \frac{z^{nj}}{\gamma _{j+1}}}},$$
(11)
where $`(\zeta +j)_{n1}\mathrm{\Gamma }(\zeta +j+n1)/\mathrm{\Gamma }(\zeta +j)`$ are the Pochhammer symbols and $`\zeta =1`$ is usually taken. The approximant (11) is a ratio of two polynomials in $`z`$ of power $`n`$ each, and when expanded back in powers of $`z`$ it reproduces all the terms of $`F_{[n+1]}`$.

The authors  applied (11) directly to the TPS’s of $`\delta \stackrel{~}{}^{\mathrm{pert}.}(\stackrel{~}{a};p)/\stackrel{~}{a}^2`$ of (4), and when re-expanding the approximant in powers of $`\stackrel{~}{a}`$ they were able to predict the next coefficient in the series with a better precision than the one provided by the corresponding diagonal (or almost diagonal) Padé approximant. Further, in the case of the pure magnetic field they showed that the method (11), when applied directly to the TPS’s in $`\stackrel{~}{b}`$ of the induced Lagrangian density,555 The approximants (11) applied to the TPS’s of the series (7) divided by $`\stackrel{~}{b}^2`$. gave better results of resummation than the corresponding Padé approximants.

We now combine the method (11) with the Borel transformation (4) $``$ (8), and compare the results of resummation obtained in this way with the results of the corresponding Borel–Padé approximants of Ref. . Formula (11) is applied to the Borel transform (8) divided by $`w`$. We identify $`zw^2`$ (we thank the authors of  for pointing out that this clarification was missing in the original version of the preprint). In the ensuing Borel–Weniger approximant, we integrate in the Laplace–Borel integral over the poles of the integrand with the Cauchy principal value prescription, just as in Borel–Padé approximant (9), in accordance with the full known solution (1).

The results of these calculations are presented in Figs. 1(a)–(d), as functions of the electric field strength parameter $`\stackrel{~}{a}`$, for various values of $`p\stackrel{~}{b}/\stackrel{~}{a}=0.,0.5,1.5,5.0`$. In Fig. 2 we present the analogous results for the case of the pure magnetic field (p.m.f.), as function of the magnetic field parameter $`\stackrel{~}{b}`$. N3 and $`[3/4]`$ denote the Borel–Weniger and the Borel–Padé resummations based on the truncated Borel transform (8) with the first four nonzero terms (i.e., three terms beyond the leading order); N5 and $`[5/6]`$ are based on the first six terms in (8). Comparison with the exact solutions, also present in the Figures, shows that Borel–Padé is better than the corresponding Borel–Weniger, except in the case of $`p=5.0`$ (electric field combined with a much stronger magnetic field). Fig. 2 suggests that Borel–Padé is better than Borel–Weniger for resummation of functions whose Borel transforms have singularities only outside the positive axis. Further, comparison of Fig. 2 with the results of Table I of Ref.  suggests that Borel–Padé and Borel–Weniger methods are much more efficient than Weniger method in resumming series with singularities in the Borel plane. Weniger method in the p.m.f. case is better than Padé method .

We can also do analogous calculations for the induced energy densities $`\delta U`$

$`\delta 𝒰`$ $`=`$ $`a{\displaystyle \frac{\mathrm{Re}\delta }{a}}|_b\mathrm{Re}\delta ,`$ (12)
$`\delta \stackrel{~}{𝒰}(\stackrel{~}{a};p)`$ $`=`$ $`\mathrm{Re}{\displaystyle _0^{\mathrm{}}}𝑑w\mathrm{exp}\left({\displaystyle \frac{w}{\stackrel{~}{a}}}\right){\displaystyle \frac{(1)}{w}}\left[{\displaystyle \frac{pw}{\mathrm{sin}^2(w+\mathrm{i}ϵ^{})}}\mathrm{coth}(pw)+{\displaystyle \frac{1}{3}}(1+p^2)+{\displaystyle \frac{1}{w^2}}\right],`$ (13)
$`\delta \stackrel{~}{𝒰}^{\mathrm{pert}.}`$ $`=`$ $`\left[d_1(p)1!\stackrel{~}{a}^2+d_3(p)3!\stackrel{~}{a}^4+\mathrm{}\right],`$ (14)
$`d_1(p)`$ $`=`$ $`{\displaystyle \frac{1}{45}}\left[3+5p^2p^4\right],d_3(p)={\displaystyle \frac{1}{945}}\left[10+21p^27p^4+2p^6\right],\mathrm{etc}.`$ (15)
where $`\delta \stackrel{~}{𝒰}8\pi ^2\delta 𝒰/(m^4\stackrel{~}{a}^2)`$. In that case, the simple Borel transform has a double–pole structure on the positive real axis, and the Padé and Weniger approximants have trouble simulating such multiple poles adequately. Therefore, we employ a slightly modified Borel transform in the case of the induced energy densities

$$MB_U(w;p)=d_1(p)\frac{w^2}{2}+d_3(p)\frac{w^4}{4}+d_5(p)\frac{w^6}{6}+\mathrm{},$$
(16)
which has no multiple–pole structure – all the poles are simple. The (modified) Laplace–Borel integral in this case is

$$\delta \stackrel{~}{𝒰}(\stackrel{~}{a};p)=\frac{1}{\stackrel{~}{a}}_0^{\mathrm{}}𝑑w\mathrm{exp}\left(\frac{w}{\stackrel{~}{a}}\right)MB_U(w;p),$$
(17)
where again the Cauchy principal value has to be taken, once $`MB_U(w;p)`$ is replaced in (17) by its Padé or Weniger approximants. For details, we refer to Ref.  where Borel–Padé was employed also for the induced energy densities. Weniger formula (11) is now applied to the modified Borel transform (16) divided by $`w^2`$. The results are presented in Figs. 3(a)–(d), as functions of $`\stackrel{~}{a}`$ at fixed $`p=0.,0.5,1.5,5.0`$, respectively.666 In the case of the pure magnetic field, the energy density is the same as the Lagrangian density, except for the sign change. We present the solutions of Borel–Weniger and Borel–Padé based on the first four (N3, \[4/4\]) and six (N5, \[6/6\]) nonzero terms of the modified Borel transform of the energy density. We see that for the induced energy density the situation is less clear. In the cases $`p=0`$, $`0.5`$ and $`5.0`$ the Borel–Padé and Borel–Weniger resummations are apparently of comparable quality, while at $`p=1.5`$ the Borel–Padé appears to work better.

We can see these trends also if we compare the perturbation coefficients predicted by these two methods with the exact ones. These results are written in Table I for the case of the Lagrangian density (predicted $`c_9`$ and $`c_{13}`$) and in Table II for the case of the energy density (predicted $`d_9`$ and $`d_{13}`$) . Predictions of Borel–Padé and Borel–Weniger are of comparable quality in the cases of $`p=0.,0.5,5.0`$ for energy density and in the case of $`p=5.0`$ for Lagrangian density. In other cases, predictions of Borel–Padé are better.

In fact, in the case $`p=5.0`$ of the energy density, the modified Borel–Weniger is slightly, but discernibly, better than the modified Borel–Padé. Comparing predictions of Table I (for $`p=0.0`$) with those of Tables II and III of Ref.  suggests strongly that the discussed Borel–Padé and Borel–Weniger methods are better than Weniger method in predicting the coefficients $`c_n`$. Weniger method is better than Padé method in predicting $`c_n`$’s .

(2) The second example to compare the efficiency of the Borel–Padé and Borel–Weniger methods will be taken from QCD, and it will have to do with the “fixing” of a pole of a Borel transform rather than with a resummation. We look at the Bjorken polarized sum rule (BjPSR), which involves the isotriplet combination of the first moments over $`x_{\mathrm{Bj}}`$ of proton and neutron polarized structure functions

$$_0^1𝑑x_{\mathrm{Bj}}\left[g_1^{(p)}(x_{\mathrm{Bj}};Q_{\mathrm{ph}}^2)g_1^{(n)}(x_{\mathrm{Bj}};Q_{\mathrm{ph}}^2)\right]=\frac{1}{6}|g_A|\left[1S(Q_{\mathrm{ph}}^2)\right].$$
(18)
Here, $`p^2=Q_{\mathrm{ph}}^2`$$`<0`$ is $`\gamma ^{}`$ momentum transfer. At $`Q_{\mathrm{ph}}^2=3\mathrm{G}\mathrm{e}\mathrm{V}^2`$ where three quarks are assumed active ($`n_f=3`$), and if taking $`\overline{\mathrm{MS}}`$ scheme and renormalization scale (RScl) $`Q_0^2=Q_{\mathrm{ph}}^2`$, we have the following TPS of the BjPSR observable $`S(Q_{\mathrm{ph}}^2)`$ available :

$`S_{[2]}(Q_{\mathrm{ph}}^2;Q_0^2=Q_{\mathrm{ph}}^2;c_2^{\mathrm{MS}},c_3^{\mathrm{MS}})`$ $`=`$ $`a_0(1+3.583a_0+20.215a_0^2),`$ (19)
$`\mathrm{with}:a_0=a(\mathrm{ln}Q_0^2;c_2^{\mathrm{MS}},c_3^{\mathrm{MS}},\mathrm{}),n_f`$ $`=`$ $`3,c_2^{\mathrm{MS}}=4.471,c_3^{\mathrm{MS}}=20.99.`$ (20)
Here we denoted by $`a`$ the strong coupling parameter $`a\alpha _s/\pi `$.

It is known from  that the Borel transform $`B_S(z)`$ of $`S`$ has the lowest positive pole at $`z_{\mathrm{pole}}=1/\beta _0=4/9`$ (leading infrared renormalon) and that this pole has a much stronger residuum than the highest negative pole at $`z_{\mathrm{pole}}=1/\beta _0`$ (leading ultraviolet renormalon). The question we raise here is: How well can Padé and Weniger approximants to the Borel transform $`B_S(z)`$ determine the next coefficient $`r_3`$ of the term $`r_3a_0^4`$ in the TPS (19), via the requirement that $`z_{\mathrm{pole}}=4/9`$? For that, we have to know well the actual $`r_3`$. That term can be determined reasonably well on the basis of two approximants discussed in  – the effective charge approximant (ECH) $`𝒜_S^{(\mathrm{ECH})}(c_3)`$ with $`c_320.`$, and another, also RScl– and scheme–independent approximant $`𝒜_{S^2}^{1/2}(c_3)`$ with $`c_315.5`$. These two approximants give the correct location of the leading infrared renormalon pole, and when we expand them back in powers of $`a_0`$ we obtain $`r_3129.4`$ and $`r_3130.8`$, respectively. Therefore, we can estimate with high confidence the actual $`r_3`$: $`r_3=130.\pm 1`$.

It is important to consider the RScl– and scheme–invariant Borel transform when we want to apply Padé or Weniger approximants to it, so that the predicted values of $`r_3`$ will be independent of the RScl– and scheme in which we work at the intermediate stage. Such a Borel transform has been used in , and we use its variant $`\stackrel{~}{B}_S(z)`$ as specified in  \[cf. Eqs. (18)–(20) there\]. Such a Borel transform reduces (up to a $`z`$–dependent nonsingular factor) to the usual Borel transform in the approximation of the one–loop evolution. The resulting power expansion of $`\stackrel{~}{B}_S(z)`$ up to $``$$`z^3`$ will depend on the coefficient $`r_3`$

$$\stackrel{~}{B}_S(z)=1+\frac{32}{81}(\gamma 1)y+0.02078\mathrm{}y^2+\frac{8}{729}(21.88\mathrm{}+\frac{1}{6}r_3)y^3+𝒪(y^4),$$
(21)
where $`\gamma =0.577\mathrm{}`$ is Euler constant, and $`y2\beta _0z`$. If we apply $`[2/1]`$ and $`[1/2]`$ Padé approximants to the TPS (21) and demand $`z_{\mathrm{pole}}=1/\beta _0`$ ($`y_{\mathrm{pole}}=2`$), we obtain predictions $`r_3=137.0`$ and $`r_3=128.0`$, respectively. The prediction of $`[1/2]`$ is significantly better, and this could possibly be explained with the more involved denominator structure of $`[1/2]`$ in comparison to $`[2/1]`$. When applying to (21) Weniger formula (11) ($`\delta _2^{(0)}`$ with $`\zeta =1`$), we obtain $`r_3=135.3`$. This is further away from the actual value of $`130.\pm 1.`$ than the prediction of $`[1/2]`$. In both $`[1/2]`$ and $`\delta _2^{(0)}`$, the denominators are polynomials of quadratic degree in $`z`$.

To summarize this QCD example: We applied Padé and Weniger approximants to a (TPS of a) Borel transform of the Bjorken polarized sum rule and demanded that the leading infrared renormalon pole be reproduced correctly. Weniger approximant $`\delta _2^{(0)}`$ then apparently gives a somewhat worse prediction for the next coefficient than the corresponding Padé approximant $`[1/2]`$.

The work of G.C. was supported by the Korean Science and Engineering Foundation (KOSEF). The work of J.-Y.Y. was supported by the German Federal Ministry of Science (BMBF).


format: text

text:
# Stability properties of trapped Bose-Fermi gases mixture

## Abstract

The stability of Bose-Fermi gases trapped in an isotropic potentials at ultracold temperature is strongly influenced by the interaction between the fermions and the bosons. At zero temperature, the stability criterion is given in this paper using variation method, the results show that whether a fermion-boson mixture is stable depends mainly on the interaction between the fermions and the bosons. For finite temperature, however, the stability is not only related to the coupling constants, but also to the temperature. The stability conditions for finite temperature are also derived and discuss in details in this paper.
PACS number(s):03.75.Fi, 05.30.Fk,05.30.Jp

Since the realization of dilute alkali atomic vapor condensates(Bose-Einstein condensation or BEC) in 1995, large efforts have been make to study many-body effects and macroscopic properties of the gases, which may be more transparently demonstrated in BEC than in other many-body systems. For fermionic atomic vapor, however, it is difficult to achieve a degenerate gas. Since the evaporative cooling of a pure fermionic gas is ineffective at temperature sufficiently low due to the suppression of $`s`$ wave scattering between identical fermions. As theory and experiment advanced, a new rich phenomenology has appeared in which new conditions arise, which are not accessible in other BEC systems. One of the most stunning of these is the recent experimental demonstration of a condensate mixture composed of two spin states of $`{}_{}{}^{87}Rb`$. The realization of two condensates mixture is related to the sympathetic cooling mechanism, i.e., the exchange of energy due to elastic collisions between atoms of cooled and thermal samples. Most recently, B.DeMarco and D.S.Jin  report their observation of degenerate Fermi gas using an evaporative cooling strategy. Although the strategy uses a two-component Fermi gas, the mixture of Bose and Fermi gas attracts a lot of attention from the viewpoint of both experiment and theoretical study.

The mixed system of Bose and Fermi particles is itself an interesting subject for investigation. The hydrogen deuterium system has been studied at the early stage of these investigations, and there is now a lot of literature devoted to the properties of pure degenerate trapped atomic Fermi gases\[5-9\].

In a recent paper, M$`\varphi `$lmer has used a simple mean field models to study the spatial distribution of a Bose-Fermi gas mixture at $`T=0K`$ within Thomas-Fermi approximation. The results show that the distributions depend strongly on the relative sign and magnitude of the boson-boson and boson-fermion scattering lengths. Here, we shall study the Bose-Fermi gas mixture using a variation method at zero temperature, this method was first introduced in to study the BEC ground state in a harmonic trap of a Bose system, and later generalized by H.Shi and W.M.Zheng to study BEC with attractive interactions. In addition, we study the stability of the Bose-Fermi gas mixture at finite temperature. The results show that there is a region of temperature in which the phase separation of the mixture happens. And the span of the region depends on the coupling constants.

To begin, we consider a second-quantized grand canonical Hamiltonian of interacting Bose and Fermi gases

$`H`$ $`=`$ $`H_b+H_f+V_{bf},`$ (1)
$`H_b`$ $`=`$ $`{\displaystyle 𝑑r\varphi ^+(r)(\frac{p^2}{2m_b}\mu _b+\frac{1}{2}m_b\omega _br^2)\varphi (r)}+{\displaystyle \frac{g_{bb}}{2}}{\displaystyle 𝑑r𝑑r^{^{}}\varphi ^+(r)\varphi ^+(r^{^{}})\varphi (r^{^{}})\varphi (r)},`$ (2)
$`H_f`$ $`=`$ $`{\displaystyle 𝑑r\psi ^+(r)(\frac{p^2}{2m_f}\mu _f+\frac{1}{2}m_f\omega _fr^2)\psi (r)},`$ (3)
$`V_{bf}`$ $`=`$ $`g_{bf}{\displaystyle 𝑑r𝑑r^{^{}}\varphi ^+(r)\psi ^+(r^{^{}})\delta (rr^{^{}})\psi (r^{^{}})\varphi (r)},`$ (4)
where $`\varphi (r)`$ and $`\psi (r)`$ denote boson and fermion field operators with masses $`m_b`$ and $`m_f`$, respectively. For weakly interacting dilute gases, the interactions between the bosonic atoms are modeled by $`\delta `$ potentials and the interactions among the fermionic atoms are neglected, since the interactions between atoms at very low temperature is suppressed for polarized systems. $`g_{bb}`$ and $`g_{ba}`$ stand for boson-boson and boson-fermion coupling constant, respectively.

$$g_{bb}=\frac{4\pi \mathrm{}^2}{m_b}a_{bb},g_{bf}=\frac{2\pi \mathrm{}^2}{m_{bf}}a_{bf},$$

$`a_{bb}`$($`a_{bf}`$) are $`s`$-wave scattering length between boson and boson (boson and fermion), and $`m_{bf}`$ is a reduced mass of the boson and the fermion. The chemical potentials $`\mu _b`$ and $`\mu _f`$ are determined through the conditions

$$N_b=𝑑r\varphi ^+(r)\varphi (r),N_f=𝑑r\psi ^+(r)\psi (r).$$
(5)
At $`T=0`$, self-consistent mean field theory, assuming that all $`N`$ bosonic particles in a gas populated the same state denoted by single particle wave function $`\mathrm{\Phi }(r)`$, lead to a nonlinear Schrödinger equation (or the Gross-Pitaevskii equation) for $`\mathrm{\Phi }(r)=\varphi (r)`$

$$[\frac{\mathrm{}^2}{2m_b}^2+\frac{1}{2}m_b\omega _b^2r^2+g_{bb}n_b(r)]\mathrm{\Phi }(r)=E_b\mathrm{\Phi }(r),$$
(6)
we here omit quantities $`g_{bf}n_f(r)`$, which is smaller than $`g_{bb}n_b(r)`$ in the case of $`N_b>>N_f`$. In order to get a degenerate fermionic gas, the boson particles appear in the system only as a coolant, so the number of bosons is always much larger than the number of fermions. In the same approximation, the fermionic wave function is given by a Slater determinant

$$\mathrm{\Psi }(r_1,r_2,\mathrm{},r_{N_f})=\frac{1}{\sqrt{N_f!}}\left[\begin{array}{cccc}\mathrm{\Psi }_1(r_1)& \mathrm{\Psi }_1(r_2)& \mathrm{}& \mathrm{\Psi }_1(r_{N_f})\\ \mathrm{\Psi }_2(r_1)& \mathrm{\Psi }_2(r_2)& \mathrm{}& \mathrm{\Psi }_2(r_{N_f})\\ & & \mathrm{}& \\ & & \mathrm{}& \\ & & \mathrm{}& \\ \mathrm{\Psi }_{N_f}(r_1)& \mathrm{\Psi }_{N_f}(r_2)& \mathrm{}& \mathrm{\Psi }_{N_f}(r_{N_f})\end{array}\right],$$
(7)
where $`\mathrm{\Psi }_i(r)`$ is the single particle states determined by Hartree-Fock self-consistent equation

$$[\frac{\mathrm{}^2}{2m_f}^2+\frac{1}{2}m_f\omega _f^2r^2+g_{bf}n_b(r)]\mathrm{\Psi }_i(r)=E_i\mathrm{\Psi }_i(r).$$
(8)
The density of the fermions is given by

$$n_f(r)=|\mathrm{\Psi }(r)|^2.$$
(9)
In the semiclassical (Thomas-Fermi) approximation, the particle are assigned classical position and momenta, but the effects of quantum statistics are taken into account. Under this approximation, the Eqs.(3) and (5) for the boson and fermion wave function are equivalent to

$`{\displaystyle \frac{1}{2}}m_b\omega _b^2r^2+g_{bb}n_b(r)`$ $`=`$ $`\mu _b,`$ (10)
$`{\displaystyle \frac{\mathrm{}^2}{2m_f}}[6\pi ^2n_f(r)]^{\frac{2}{3}}+{\displaystyle \frac{1}{2}}m_f\omega _f^2r^2+g_{bf}n_b(r)`$ $`=`$ $`e_F.`$ (11)
The main conclusion of this equations is discussed in Ref.. We obtain $`n_b(r)=\frac{1}{g_{bb}}(\mu _b\frac{1}{2}m_b\omega _b^2r^2)`$ from the first line of Eqs(7). Substituting $`n_b(r)`$ into the second line of Eqs(7), we yield

$$\frac{\mathrm{}^2}{2m_f}[6\pi ^2n_f(r)]^{\frac{2}{3}}+\frac{1}{2}m_f\omega _f^2r^2+\frac{g_{bf}}{g_{bb}}(\mu _b\frac{1}{2}m_b\omega _b^2r^2)=e_F,$$
(12)
this equation shows that the fermions experience a potential minimum in the center of the trap if $`g_{bf}/g_{bb}<m_f\omega _f^2/m_b\omega _b^2,`$ in this case the entire distribution behaves like a fermionic core within the Bose condensate. The fermion density is a constant throughout the Bose condensate if $`g_{bf}/g_{bb}=m_f\omega _f^2/m_b\omega _b^2.`$ Whereas the fermions are repelled from the center of the trap and localized near the edge of the Bose condensate if $`g_{bf}/g_{bb}>m_f\omega _f^2/m_b\omega _b^2,`$ i.e. a phase separation occurs in this system. We would like to note that the distribution of BEC remains unchanged in the above discussions, since we assume $`N_b>>N_f`$. To drive Eq.(8), we assume that the Thomas-Fermi approximation(TFA) is valid. The coupling constant $`g_{bb}`$ and $`g_{bf}`$ may take any value as long as the TFA is available, and the phase separation depend mainly on ratio $`g_{bf}/g_{bb}`$. In what follows we discuss the separation of the bosonic and fermionic parts from the other aspect for zero temperature by using variation method, the results are indeed different from those under TFA. We note the solution of Eq.(5) requires prior knowledge of the boson density profile $`n_b=|\mathrm{\Phi }(r)|^2`$. To obtain the density profile, we have to solve the Gross-Pitaevskii equation (3). There are a large number of literatures devoted to solve the Gross-Pitaevskii equation, we here use a variation method to solve the problem. For a isotropic trapping potential, we may assume the trial wave function for $`\mathrm{\Phi }(r)`$ in Eq.(3) to be

$$\mathrm{\Phi }(r)=\sqrt{N_b}\omega ^{\frac{3}{4}}(\frac{m_b}{\pi \mathrm{}})^{\frac{3}{4}}e^{m_b\omega r^2/2\mathrm{}},$$
(13)
where $`\omega `$ is the effective frequency and is taken as a variational parameter. Substituting Eq.(9) into Eq.(3), we obtained the ground-state energy

$$E_b[\mathrm{\Phi }]=E_b(\omega )=\frac{3}{4}N_b\mathrm{}\omega +\frac{3}{4}N_b\mathrm{}\frac{\omega _b^2}{\omega }+g_{bb}N_b^2(\frac{\omega m_b}{2\pi \mathrm{}})^{\frac{3}{2}}.$$
(14)
If $`E_b(\omega )`$ is plotted as a function of $`\omega `$ , one sees that a stable local minimum exists only up to a certain maximum number of atoms for $`g_{bb}<0`$. The critical point occurs where

$$\frac{E_b(\omega )}{\omega }|_{(\omega =\omega _c,N_b=N_{bc})}=0,\text{ and }\frac{^2E_b(\omega )}{\omega ^2}|_{\omega =\omega _c,N_b=N_{bc}}>0.$$
(15)
Here, $`\omega _c`$ stands for the variational parameter that minimizes the ground state energy. Using equation (10), for $`g_{bb}<0`$ the critical number of bosons is given by

$$N_b^c<2\mathrm{}\omega _b^2(\omega _c)^{\frac{5}{2}}\frac{1}{|g_{bb}|}(\frac{2\pi \mathrm{}}{m_b})^{\frac{3}{2}}.$$
(16)
where $`\omega _c`$ satisfies

$$\mathrm{}\omega _c^2\mathrm{}\omega _b^2+2g_{bb}N_b(\frac{m_b}{2\pi \mathrm{}})^{\frac{3}{2}}\omega _c^{\frac{5}{2}}=0.$$
(17)
Parameter $`\omega _b166Hz`$ relevant to the experiment gives $`N_b^c1400`$, which is in good agreement with the experiment. The solution $`\omega _c`$ of Eq.(13) against $`g_{bb}`$ is plotted in Fig.1, which shows that as $`|g_{bb}|`$ increases, the variation parameter $`\omega _c`$ decreases, and it has a maximum equal to $`\omega _b`$ at $`g_{bb}=0`$. We will use this solution to study the stability of the mixture at zero temperature below.

We may determine the ground-sate energy functional of the fermions provided $`n_b(r)`$ is known. In terms of the fermion distribution $`n_f(r)`$, the energy functional $`E_f`$ of the fermions is given by

$$E_f=E_f[n_f(r)]=\frac{d^3r}{6\pi ^2}\frac{\mathrm{}^2}{2m_f}[6\pi ^2n_f(r)]^{5/3}+\frac{1}{2}d^3rm_f\omega _f^2r^2n_f(r)+d^3rg_{bf}n_f(r)n_b(r),$$
(18)
since the interaction between the bosons and the fermions is rather week, we may consider a Gaussian function as a trial fermions’ distribution

$$n_f(r)=N_f\mathrm{\Omega }^{3/2}(\frac{m_f}{\pi \mathrm{}})^{3/2}e^{\frac{m_f\mathrm{\Omega }(rr_f)^2}{\mathrm{}}}.$$
(19)
Here $`r_f`$ and $`\mathrm{\Omega }`$ is treated as variation parameters. Substituting this wave functions into Eq.(14), one obtains

$$E_f=E_f(\mathrm{\Omega },r_f,N_f,\omega _c)=P+\frac{1}{2}m_f\omega _f^2r_f^2N_f+hN_bN_f(\pi G)^{\frac{3}{2}}e^{Gr_f^2}$$
(20)
with

$$P=P(\mathrm{\Omega },N_f)=(\frac{3}{5})^{\frac{3}{2}}\frac{1}{\pi }(6\pi ^2)^{\frac{2}{3}}\mathrm{}\mathrm{\Omega }N_f^{\frac{5}{3}}+\frac{3\mathrm{}\omega _f^2}{4\mathrm{\Omega }}N_f$$

$$G=G(\mathrm{\Omega },\omega _c)=\frac{m_fm_b\omega _c\mathrm{\Omega }}{\mathrm{}(m_f\mathrm{\Omega }+m_b\omega _c)}.$$

As known, a physical state corresponds to a stable or metastable point of the energy functional. If a separation of the fermion and boson component occurs, then $`r_{fc}`$ that minimizes the energy $`E_f`$ takes a positive nonzero value. In other words, there are no separations between the two components when energy $`E_f`$ exhibits a minimal value at $`r_f=0`$. When distribution function is restricted to the form of the trial function (15) we may write the conditions of a minimal energy in terms of derivatives of the energy with respect to the adjustable variation parameters of the trial function. We show them as follows

$`{\displaystyle \frac{E_f}{\mathrm{\Omega }}}|_{\mathrm{\Omega }=\mathrm{\Omega }_c}=0,{\displaystyle \frac{E_f}{r_f}}|_{r_f=r_{fc}}=0,`$ (21)
$`{\displaystyle \frac{^2E_f}{\mathrm{\Omega }^2}}{\displaystyle \frac{^2E_f}{r_f^2}}({\displaystyle \frac{^2E_f}{\mathrm{\Omega }r_f}})^2>0,`$ (22)
the stationary conditions (there are no separation between the boson and fermion ) are

$`Y`$ $`=`$ $`Y(N_b,\omega _c,g_{bf})=[\mathrm{}{\displaystyle \frac{\omega _f^2}{\mathrm{\Omega }_c^3}}+g_{bf}N_b\pi ^{\frac{3}{2}}\sqrt{G(\mathrm{\Omega }_c,\omega _c)}{\displaystyle \frac{^2G(\mathrm{\Omega }_c,\omega _c)}{\mathrm{\Omega }^2}}+{\displaystyle \frac{g_{bf}N_b\pi ^{3/2}}{2G(\mathrm{\Omega }_c,\omega _c)}}({\displaystyle \frac{G(\mathrm{\Omega }_c,\omega _c)}{\mathrm{\Omega }}})^2]`$ (23)
$`\times `$ $`[m_f\omega _f^22g_{bf}N_b\pi ^{3/2}G^{5/2}(\mathrm{\Omega }_c,\omega _c)]>0,`$ (24)
where $`\mathrm{\Omega }_c`$ determined by $`E_f/\mathrm{\Omega }|_{\mathrm{\Omega }_c}=0`$ satisfies the following equation

$$\frac{2}{3\pi }(6\pi ^2)^{\frac{2}{3}}(\frac{3}{5})^{\frac{3}{2}}\mathrm{}N_f^{\frac{2}{3}}\frac{1}{2}\mathrm{}(\frac{\omega _f}{\mathrm{\Omega }_c})^2+g_{bf}N_b\pi ^{\frac{3}{2}}\sqrt{G(\mathrm{\Omega }_c,\omega _c)}\frac{G(\mathrm{\Omega }_c,\omega _c)}{\mathrm{\Omega }}=0,$$
(25)
the solution of Eq.(19) as a function of $`g_{bf}`$ is shown in Fig.2, a magnification part of the curve near $`g_{bf}=0`$ (but $`g_{bf}>0`$) is give in the inset. This curve indicates that the fermions prefer to occupy the trap centre for $`g_{bf}<0`$, and the larger the coupling constant $`|g_{bf}|(g_{bf}<0)`$, the sharper the distribution of the fermions. However, as we show below, the fermions and the bosons can not always coexist even if $`g_{bf}<0`$. $`\frac{E_f}{r_f}|_{r_f=r_{fc}}=0`$ has two solutions, one solution is $`r_{fc}=0`$ and the another is

$$r_{fc}=\sqrt{\frac{1}{G}ln(2g_{bf}N_b\pi ^{\frac{3}{2}}G^{\frac{5}{2}}(\mathrm{\Omega }_c,\omega _c))ln(m_f\omega _f^2)}.$$
(26)
For $`g_{bf}<0`$ or

$$0<g_{bf}<\frac{m_f\omega _f^2}{2N_b\pi ^{\frac{3}{2}}G^{\frac{5}{2}}(\mathrm{\Omega }_c,\omega _c)}$$

i.e. the interaction between fermion and the boson is attractive or weekly repulsive, the solution $`r_{fc}=0`$ holds, which indicate that there is not separation between the fermions and the bosons. For $`g_{bf}>\frac{m_f\omega _f^2}{2N_b\pi ^{\frac{3}{2}}G^{\frac{5}{2}}(\mathrm{\Omega }_c,\omega _c)}`$ the fermions experience a effective potential minimum at $`r_{fc}=\sqrt{\frac{1}{G}ln(2g_{bf}N_b\pi ^{\frac{3}{2}}G^{\frac{5}{2}}(\mathrm{\Omega }_c,\omega _c))ln(m_f\omega _f^2)}>0,`$ the Bose condensate is surrounded by a shell of fermions in this case. $`Y(N_b,g_{bf},\omega _c)`$ as functions of the coupling constant $`g_{bf}`$ are shown in figure 3. We see that whether the boson-fermion mixture is stable depends not only on the coupling constant $`g_{bf}`$ and $`g_{bb}`$(through $`\omega _c)`$, but also on $`N_b`$ and $`N_f`$(through $`\mathrm{\Omega }_c`$), i.e., the stability of the mixture system depends on the number of both boson and fermion system. For example, in Fig.3-a we show $`Y`$ given by Eq.(8) as a function of the coupling constant $`g_{bf}`$ for fixed $`N_f=100,N_b=10000`$, while Fig.3-b is for the same parameters as in Fig.3-a except for $`N_b=1000`$, it is obvious that the region of $`g_{bf}`$ in which the system has no phase separation has been broadened with $`N_b`$ decreases (for fixed $`N_f`$). The inset present the dependence of $`Y`$ on $`g_{bf}`$ at a larger scale of $`g_{bf}.`$ It is interesting to compare the above mentioned results with those obtained by treating the fermions in the Thomas-Fermi approximation, this is done in Ref., and we note that the semiclassical description gives a qualitatively correct description and it reliably predicts the phase separation.

Now we tune our attention to discuss the above problem at finite temperature. First of all, we consider the homogenous case, for the boson and fermion system, thermodynamical properties are trivial if there are not interaction between them. But in this case the sympathetic cooling scheme does not take any effect and the degenerate fermions in a trapped potential have not been achieved. The thermodynamical properties may be changed when the interaction between the fermions and bosions is turn on, then a new phenomenon, the phase separation, may occur in this system. for a homogeneous fermion and boson mixture system, the Helmholts free energy can be written as

$$\beta F=\frac{V}{\lambda _2^3}f_{\frac{5}{2}}(z_f)+\frac{1}{2}g_{ff}\rho _fN_f\lambda _f^2+ln(1z_b)\frac{V}{\lambda _b^3}g_{\frac{5}{2}}(z_b)+2g_{bb}\rho _bN_b\lambda _b^2+g_{bf}(\lambda _b^2+\lambda _f^2)N_fN_b/V,$$
(27)
where index $`f`$ refers to the fermionic component, whereas index $`b`$ stands for the bosonic one, $`N_i`$ is the number of particles in component $`i`$, $`\lambda _i`$ denotes the thermal wave length of component $`i`$, $`f_n(z)`$ and $`g_n(z)`$ represent the Fermi and Bose integral, respectively.The equation (21) is based on the pseudopotential form of the atom-atom interaction, and may be assumed accurate when the system is dilute. i.e. $`\rho _ig_{ii}^3<<1`$ and $`g_{ii}/\lambda _i<<1`$, where $`\rho _i`$ is the density of the component $`i`$. This condition is well satisfies for the samples of alkali atoms in experiments to date.

From eq.(21) we obtain the chemical potential for each component straightforwardly,

$`\beta \mu _b`$ $`=`$ $`\beta \mu _b^0+4g_{bb}\rho _b\lambda _b^2+g_{bf}(\lambda _b^2+\lambda _f^2)N_f/V,`$ (28)
$`\beta \mu _f`$ $`=`$ $`\beta \mu _f^0+g_{ff}\rho _f\lambda _f^2+g_{bf}(\lambda _b^2+\lambda _f^2)N_b/V,`$ (29)
where $`\mu _i^0`$ are the chemical potentials of ideal gas. There are three terms in each chemical potential, the second term comes from the interaction within the component and the third term is from the interaction between the fermion and boson component. As known, an homogenous binary mixture is stable only when the symmetric matrix $`\widehat{\mu }`$ given by

$$\widehat{\mu }=\left[\begin{array}{cc}\frac{\mu _b}{\rho _b}\hfill & \frac{\mu _b}{\rho _f}\hfill \\ \frac{\mu _f}{\rho _b}\hfill & \frac{\mu _f}{\rho _f}\hfill \end{array}\right]$$
(30)
is non-negatively definite, in other words, all eigenvalues of matrix $`\widehat{\mu }`$ given in Eq.(23) are non-negative. Mathematically, for homogeneous fermion and boson mixture the stability conditions are

$$\frac{\mu _b}{\rho _b}0,\frac{\mu _f}{\rho _f}0,$$
(31)
and

$$det\left[\begin{array}{cc}\frac{\mu _b}{\rho _b}\hfill & \frac{\mu _b}{\rho _f}\hfill \\ \frac{\mu _f}{\rho _b}\hfill & \frac{\mu _f}{\rho _f}\hfill \end{array}\right]0.$$
(32)
For ideal gas, we have $`\rho _b=\frac{1}{\lambda _b^3}g_{\frac{3}{2}}(z_b),`$ $`\rho _f=\frac{1}{\lambda _f^3}f_{\frac{3}{2}}(z_f),`$ this leads to

$$\beta \frac{\mu _f^0}{\rho _f}=\frac{\lambda _f^3}{f_{\frac{1}{2}}(z_f)},\beta \frac{\mu _b^0}{\rho _b}=\frac{\lambda _b^3}{g_{\frac{1}{2}}(z_b)},$$
(33)
It follows from eqs (24) and (25) that

$`4g_{bb}\lambda _b^2+{\displaystyle \frac{\lambda _b^3}{g_{\frac{1}{2}}(z_b)}}0,`$ (34)
$`g_{ff}\lambda _f^2+{\displaystyle \frac{\lambda _f^3}{f_{\frac{1}{2}}(z_f)}}0,`$ (35)
and (36)
$`Z(T,g_{bf},g_{ff},g_{bb})=Z=(4g_{bb}\lambda _b^2+{\displaystyle \frac{\lambda _b^3}{g_{\frac{1}{2}}(z_b)}})(g_{ff}\lambda _f^2+{\displaystyle \frac{\lambda _f^3}{f_{\frac{1}{2}(z_f)}}})g_{bf}^2(\lambda _b^2+\lambda _f^2)^20.`$ (37)
It is well known that a homogeneous imperfect gas with attractive interaction is not stable. The fermions in this kind of gas could form BCS state, which consists two fermionic particles interacting with each other but not with the other fermions from the Fermi gas, whereas bosons with attractive interaction could collapse into liquid. Hence, we here discuss the system with repulsive interactions. It is obvious that the stability condition (27) and (28) hold always for $`g_{bb}>0`$, $`g_{ff}>0`$. We would like to point out that the stability conditions (27-29) do not involve the densities of the both components. At first sight, this seems to be confusion, in fact, there is no contradiction. One can demonstrate that at low density the Helmholtz free energy of the bogoliubov gas reduce to a quadratic form in $`N_b`$ and $`N_f`$. To have a minimum, this form should be positive definite, i.e., $`\text{det}\frac{^2F}{N_bN_f}0.`$ Therefore, the corresponding stability criterion involves only density-independent constants in the order of approximation used. This criterion is similar to the stability conditions for two-component Bose-Einstein condensate in a trapped untracold gas\[19-25\]. When $`T\mathrm{},\lambda _i0`$, hence $`Z\frac{1}{\rho _b\rho _f}`$. Thus at high temperature, the homogeneous binary gas mixture is always stable and no phase separation occur. In the case considered here, Fermi temperature $`T_F=\frac{h^2}{2mk_B}(\frac{3N_f}{8\pi V})^{\frac{2}{3}}`$ is much lower than BEC temperature $`T_c=\frac{h^2}{2\pi mK_B}(\frac{N_b}{2.612V})^{\frac{2}{3}},`$ i.e.,as temperature decreases, it first passes the BEC transition point $`T_c`$. When $`TT_c`$, $`g_{\frac{1}{2}}(1)\mathrm{}`$, so

$$Z(T,g_{bf},g_{bb},g_{ff})4g_{bb}\lambda _b^2(g_{ff}\lambda _f^2+\frac{\lambda _f^3}{f_{\frac{1}{2}}(z_f)})g_{bf}^2(\lambda _b^2+\lambda _f^2)^2.$$

In particular, when $`T<<T_F`$,i.e., the temperature is much smaller than the Fermi temperature of the fermion system, the stability condition becomes (setting $`m_f=m_b`$)

$$g_{bb}g_{ff}g_{bf}^20,$$
(38)
which does not depend on temperature and coincides with the stability conditions of two-component BEC. Although it is difficult to reach this region of very low temperature, yet it attracts much more attention. Because both superfluidity and shell effects are expected to occur at temperature much smaller than the Fermi temperature. $`Z`$ given by Eq.(29) as a function of the temperature is shown in Fig.4, we see that the system is always stable when $`T0`$ and $`T\mathrm{}`$, and the system is unstable for $`T_{c1}<T<T_{c2}`$, where $`T_{c1}`$ and $`T_{c2}`$ are roots of $`Z(T,g_{bb},g_{bf},g_{ff})=0`$. In particular, $`T_{c1}`$ and $`T_{c2}`$ depend on $`g_{bf}`$ $`g_{ff}`$ and $`g_{bb}`$. As $`g_{bf}`$ decreases (for fixed $`g_{bb}`$ and $`g_{ff}`$), $`T_{c1}`$ tends to $`T_{c2}`$ (in Fig.4 going from dotted line to solid line). The critical temperature $`T_{c1}`$ and $`T_{c2}`$ characterize the onset of the phase separation, which is quite different from the Bose-Einstein condensation and the degenerate fermions. The critical temperature of BEC and of the onset of degenerate fermionic gas depend mainly on the density of the system $`N_i/V(i=b,f)`$. Especially, the BEC and the degenerate fermionic gas may happen even if $`g_{bf}=0`$. For the phase separation, however, nothing will happen if $`g_{bf}=0`$. For a fixed temperature and the coupling constant $`g_{bf}`$, $`Z`$ vs. $`g_{bb}`$ and $`g_{ff}`$ is shown in Fig.5, which represents the dependence of the stability on the interaction strength inside each component.

Until now, we considered only a homogeneous Fermi-Bose gas mixture at finite temperature. In reality, however, experiments with ultracold atoms are performed by trapping and cooling in an external potential that can be generally modeled by an isotropic harmonic oscillator $`V(r)=\frac{m}{2}\omega _t^2r^2`$, where $`\omega _t`$ is the trapping frequency. An exact criterion for the stability of an inhomogeneous Bose-fermi mixture should involve calculating the Helmholts free energy as a function at all eigenstates of the trapping potential. Fortunately, in the system considered here it is a good approximation to take use of the local-density approximation, which treats the system as being locally homogeneous. This requires that the level spacing $`\mathrm{}\omega _t`$ of the trapping potential is much smaller than the Fermi energy. Of course, the local density approximation always breaks down at the edge of the gas cloud where the density vanishes and the effective Fermi energy becomes zero. In this approximation, the stability conditions can still be calculated by means of the equations derived above, with the understanding that now the effective chemical potentials are spatially dependent through

$$\mu _b=\mu _b^0\frac{1}{2}m_b\omega _b^2r^2,\mu _f=\mu _f^0\frac{1}{2}m_f\omega _f^2r^2.$$

Thus a local stability condition is the same as given in Eq.(29) but replacing $`z_i(i=1,2)`$ by

$$\stackrel{~}{z}_1=z_1e^{\frac{\beta }{2}m_b\omega _b^2r^2},\text{ and }\stackrel{~}{z}_2=z_2e^{\frac{\beta }{2}m_f\omega _f^2r^2}.$$

As shown in inset of Fig.4, the region of temperature in which the system is unstable decrease for the case of $`r0`$. As compared with the case without trapped potential, the total energy of the system increases for it in a trap. Alternatively, within the TFA, the chemical potentials decrease in this process. So this effect is equal to be that the particle number of the system has a loss. In this sense, the system is more stable than before.

In summary, we considered a dilute Bose-Fermi gas mixture in an isotropic trap. The atom can interact via s-wave scattering except within the fermions. These interactions strongly affect the stability of the system at zero and finite temperature. In addition, the stability conditions depend on the ratio rate $`N_b/N_f`$, the larger the ratio rate, the smaller the region of stability. For finite temperature, however, the stability conditions depends not only on the interactions,but also on the temperature. The region $`T_{c1}TT_{c2}`$ in which the system is unstable depend on the strength of the interaction between and within the bosons and the fermions. For an anisotropic trap, the stability conditions remain unchanged, whereas somewhat would be changed for zero temperature compared with the case of isotropic trap. To study the effects, we should introduced the another variation parameter in Eqs (9) and (5) to characterize the BEC and the fermions in this trap. Consequently, the stability condition (18) for zero temperature changes and the phase separation could different for different orientation. These need further investigations.
ACKNOWLEDGEMENT:
We thank Dr. Li You for his stimulating and helpful discussions.

Figure captions:
Fig. 1:The parameter $`\omega _c`$ which minimizes the energy functional versus the coupling constant $`g_{bb}`$. The trapped frequency $`\omega _b=166Hz`$ and the number of the bosonic atom $`N_b=1000`$.
Fig. 2:The parameter $`\mathrm{\Omega }_c`$ which minimizes the energy functional as a function of the coupling constant $`g_{bf}`$. The parameters chosen are $`g_{bb}=0.05`$ in units of $`\mathrm{}\omega _fa^3`$ ($`a=\sqrt{\mathrm{}/\omega _bm_b}`$) and all the coupling constants are chosen in this units hence forth, $`N_f=100`$, $`\omega _f=166Hz`$. Scatter and solid line correspond to different number of bosonic atom, as specified in the figure. The inset presents the enlarged part of the curve near $`g_{bf}=0(>0)`$.
Fig.3:Plot of $`Y`$ given by Eq.(18) as a function of coupling constant $`g_{bf}`$. The parameter chosen are a:$`N_f=100,N_b=1000.`$ b:$`N_f=100,N_b=10000.`$ The curve for a larger scale of $`g_{bf}`$ is presented as an inset in the figure.
Fig.4:Plot of $`Z`$ given by Eq.(29) as a function of temperature $`T`$. The parameters chosen are $`N_b=1000,N_f=10000,g_{bb}=0.05,g_{ff}=0.01`$. Dashed-dotted line:$`g_{bf}=0.3`$, dotted line $`g_{bf}=0.02`$, solid line $`g_{bf}=0.01`$. The dotted line in the inset is the same as the dotted line in the figure, while the solid line in the inset is for the gases in a trap with trapped frequency $`166Hz`$.
Fig.5:Plot of $`Z`$ as a function of $`g_{ff}`$ and $`g_{bb}`$. The parameters chosen are temperature $`T=0.1T_F`$, $`g_{bf}=0.2`$.