format: text

text:
# Thermodynamic potential of interacting Bose-Einstein gas confined in harmonic potential

## I Introduction

Bose-Einstein condensation (BEC) was first observed in dilute ultracold alkali atoms of rubidium , lithium  and sodium . Furthermore, Mewes et al. and Ensher et al. measured the condensate fraction and the energy of rubidium atoms . It revealed that the transition temperature $`T_0`$ is shifted not only by the interaction effect but also by the finite size effect by a few percent, where $`T_0`$ denotes the transition temperature of the noninteracting Bose gas within the external field in the thermodynamic limit. This means that both the interaction and the finite size effects play an important role in the real Bose gases. In fact, the number of trapped atoms is typically $`N=10^310^6`$, and this may not be sufficiently large to take the thermodynamic limit. For the noninteracting Bose gas, the condensate fraction confined in a spherical harmonic potential can be analytically given as ,

$`{\displaystyle \frac{N_0}{N}}`$ $`=`$ $`1({\displaystyle \frac{T}{T_0}})^3{\displaystyle \frac{3}{2}}{\displaystyle \frac{\zeta (2)}{\zeta (3)^{2/3}}}{\displaystyle \frac{1}{N^{1/3}}}({\displaystyle \frac{T}{T_0}})^2,`$ (1)
where the first two terms in the right hand side correspond to the condensate fraction in the thermodynamic limit and the third term is the finite size correction.

On the other hand, these experimental findings have stimulated much interest in the theory of the interacting Bose gas. Since BEC occurs when atoms are dilute and cold, we can treat the interaction of atoms as the two-body interacting Bose gas. In this case, the s-wave scattering length characterizes the strength of the two-body interaction. Under these conditions, the Gross-Pitaevskii (GP) equation can well describe the behavior of the interacting Bose gas at zero temperature . This is a mean field approach for the order parameter associated with the condensate. Using the GP equation, several authors  studied the ground state and the excitation properties of the condensate.

To study BEC at finite temperature, the Bogoliubov approximation at zero temperature was extended by Griffin . It is called the Hartree-Fock-Bogoliubov (HFB) theory. In particular, Popov approximation to the HFB theory has been employed to explain the experimental results . This approximation neglect the anomalous average of the fluctuating field operator in the HFB theory. With the Popov approximation, thermodynamic quantities, such as the condensate fraction and the internal energy, agree with the experimental results at low temperature. However, recent experimental data of the excitation frequencies of the condensates indicate a discrepancy with these approximate results near the transition temperature $`T_0`$ $`(T/T_0>0.6)`$ .

In this paper, we derive the thermodynamic potential for the interacting Bose gas within a mean field approximation which is different from Bogoliubov approximation. With this thermodynamic potential, we obtain the thermodynamic quantities; the condensate fraction, the entropy and the internal energy. These quantities have the same form as that of the ideal Bose gas, but the energy $`E_𝐢`$ contains the interaction effect and varies as a function of the temperature $`T`$ and the total number $`N`$. Therefore, we must determine not only the chemical potential $`\mu `$ but also the energy $`E_𝐢`$ as a function of $`T`$ and $`N`$ in a self-consistent fashion. In particular, the chemical potential determines the transition temperature on the condition that the chemical potential $`\mu `$ is equal to the ground state energy $`E_0`$ at $`T=T_0`$. Here, we choose the repulsive interaction where the scattering length $`a`$ is positive, and numerically solve the chemical potential and the energy $`E_𝐢`$. By determining the chemical potential near $`TT_0`$, we find that the interaction effect shifts $`T_0`$ by a few percent toward the lower temperature which should be compared with the previous estimations of the semiclassical approximation .

This paper is organized in the following way. In the next section, we derive the thermodynamic potential for the interacting Bose gas at finite temperature. Then, we obtain thermodynamic quantities in a self-consistent way. In section 3, we carry out the numerical evaluation of the chemical potential, the condensate fraction and the internal energy, in particular near $`T_0`$. Section 4 summarizes what we have clarified in this paper.

## II Thermodynamics for Interacting Bose Gas

We derive the thermodynamic potential for the interacting Bose gas within a mean field approach. To begin with, we construct the effective Hamiltonian in this system. The Hamiltonian for the interacting Bose gas confined in an external potential $`V_{\mathrm{ext}}`$ can be written as

$`\widehat{H}`$ $`=`$ $`{\displaystyle 𝑑𝐫\widehat{\mathrm{\Psi }}^{}(𝐫)\left[\frac{\mathrm{}^2}{2m}^2+V_{\mathrm{ext}}(𝐫)\right]\widehat{\mathrm{\Psi }}(𝐫)}`$ (3)
$`+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\widehat{\mathrm{\Psi }}^{}(𝐫)\widehat{\mathrm{\Psi }}^{}(𝐫^{})V(𝐫𝐫^{})\widehat{\mathrm{\Psi }}(𝐫^{})\widehat{\mathrm{\Psi }}(𝐫)},`$
where $`\widehat{\mathrm{\Psi }}(𝐫)`$ is the boson field operator and $`V(𝐫𝐫^{})`$ is the two-body atomic potential. Here, we choose the spherical harmonic potential for the $`V_{\mathrm{ext}}`$ for simplicity. Thus, the $`V_{\mathrm{ext}}`$ is given as

$`V_{\mathrm{ext}}(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{2}}m\omega ^2r^2.`$ (4)
Since we are interested in the dilute and cold Bose gas, we replace $`V(𝐫𝐫^{})`$ by an effective interaction of the $`\delta `$-function type,

$`V(𝐫𝐫^{})`$ $`=`$ $`g\delta ^{(3)}(𝐫𝐫^{}),`$ (5)
where the coupling constant $`g`$ is related to the s-wave scattering length $`a`$ through

$`g`$ $`=`$ $`{\displaystyle \frac{4\pi \mathrm{}^2a}{m}}.`$ (6)
The scattering length is an important length scale in this system. A typical value of the scattering length $`a`$ is of the order of $`nm`$. The boson field operator for the ideal system can be written as

$`\widehat{\mathrm{\Psi }}(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{V}}}{\displaystyle \underset{𝐤}{}}a_𝐤e^{i𝐤𝐫}.`$ (7)
But for the general case, the corresponding field operator is a sum over all normal modes

$`\widehat{\mathrm{\Psi }}(𝐫)`$ $`=`$ $`{\displaystyle \underset{𝐧}{}}a_𝐧\chi _𝐧(𝐫),`$ (8)
when $`\chi _𝐧(𝐫)`$’s define any complete set of normalized single-particle wave function $`a_𝐧`$ is a bosonic annihilation operator for the single-particle state $`𝐧`$. In our case, $`\chi _𝐧(𝐫)`$ is given as

$`\chi _𝐧(𝐫)`$ $``$ $`u_{n_xn_yn_z}(x,y,z)`$ (9)
$`=`$ $`N_{n_x}N_{n_y}N_{n_z}e^{\frac{1}{2}(\frac{r}{a_{ho}})^2}`$ (11)
$`H_{n_x}(x/a_{ho})H_{n_y}(y/a_{ho})H_{n_z}(z/a_{ho}),`$
where $`H_n(\xi )`$ is the Hermite polynomial, and $`N_n`$ is the normalization constant which can be written as

$`N_n`$ $`=`$ $`\sqrt{{\displaystyle \frac{a_{ho}^1}{\pi ^{1/2}2^nn!}}},`$ (12)
where $`a_{ho}`$ denotes the harmonic oscillator length

$`a_{ho}`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mathrm{}}{m\omega }}}.`$ (13)
This is also an important length scale in this system. It is typically of the order of $`\mu m`$. From Eqs. (3) and (8), the Hamiltonian can be written as

$`\widehat{H}`$ $`=`$ $`{\displaystyle \underset{𝐧}{}}\epsilon _𝐧a_𝐧^{}a_𝐧`$ (15)
$`+{\displaystyle \frac{2\pi \mathrm{}^2a}{m}}{\displaystyle \underset{\{𝐧_𝐢\}}{}}𝐧_\mathrm{𝟏}𝐧_\mathrm{𝟐}𝐧_\mathrm{𝟑}𝐧_\mathrm{𝟒}a_{𝐧_\mathrm{𝟏}}^{}a_{𝐧_\mathrm{𝟐}}^{}a_{𝐧_\mathrm{𝟑}}a_{𝐧_\mathrm{𝟒}},`$
with

$`\epsilon _𝐧`$ $`=`$ $`\mathrm{}\omega \left(n_x+n_y+n_z+{\displaystyle \frac{3}{2}}\right),`$ (16)
$`𝐧_\mathrm{𝟏}𝐧_\mathrm{𝟐}𝐧_\mathrm{𝟑}𝐧_\mathrm{𝟒}`$ $`=`$ $`{\displaystyle 𝑑𝐫\chi _{𝐧_\mathrm{𝟏}}^{}(𝐫)\chi _{𝐧_\mathrm{𝟐}}^{}(𝐫)\chi _{𝐧_\mathrm{𝟑}}(𝐫)\chi _{𝐧_\mathrm{𝟒}}(𝐫)}.`$ (17)
Now we consider the first order perturbation theory for the thermodynamic potential . The first order correction is given by the diagonal part of the interacting Hamiltonian in Eq. (15). This corresponds to ignoring the anomalous average. Therefore, the effective Hamiltonian can be written as

$`\widehat{H}_{eff}`$ $`=`$ $`{\displaystyle \underset{𝐢}{}}\epsilon _𝐢n_𝐢+\mathrm{}\omega {\displaystyle \frac{4\pi a}{a_{ho}}}{\displaystyle \underset{𝐢,𝐣}{}}U_{\mathrm{𝐢𝐣}}n_𝐢n_𝐣,`$ (18)
where $`n_𝐢=a_𝐢^{}a_𝐢`$ and

$`U_{\mathrm{𝐢𝐣}}`$ $`=`$ $`F_{i_x,j_x}F_{i_y,j_y}F_{i_z,j_z},`$ (19)
$`F_{i,j}`$ $`=`$ $`{\displaystyle \frac{1}{\pi 2^i2^ji!j!}}{\displaystyle 𝑑\xi e^{2\xi ^2}H_i^2(\xi )H_j^2(\xi )}.`$ (20)
Now we rewrite Eq. (18) in the following way

$`\widehat{H}_{eff}`$ $`=`$ $`\widehat{H}_0+\widehat{H}^{},`$ (21)
with

$`\widehat{H}_0`$ $`=`$ $`{\displaystyle \underset{𝐢}{}}E_𝐢n_𝐢\mathrm{}\omega {\displaystyle \frac{4\pi a}{a_{ho}}}{\displaystyle \underset{𝐢,𝐣}{}}U_{\mathrm{𝐢𝐣}}\nu _𝐢\nu _𝐣,`$ (22)
$`\widehat{H}^{}`$ $`=`$ $`\mathrm{}\omega {\displaystyle \frac{4\pi a}{a_{ho}}}{\displaystyle \underset{𝐢,𝐣}{}}U_{\mathrm{𝐢𝐣}}(n_𝐢\nu _𝐢)(n_𝐣\nu _𝐣),`$ (23)
$`E_𝐢`$ $`=`$ $`\epsilon _𝐢+\mathrm{}\omega {\displaystyle \frac{8\pi a}{a_{ho}}}{\displaystyle \underset{𝐣}{}}U_{\mathrm{𝐢𝐣}}\nu _𝐣.`$ (24)
Here, we minimize $`\widehat{H}^{}`$ by taking $`\nu _𝐢`$ the mean occupation numbers

$`\nu _𝐢=n_𝐢`$ $`=`$ $`{\displaystyle \frac{\mathrm{Tr}n_𝐢e^{(\widehat{H}_0\mu \widehat{N})/k_BT}}{\mathrm{Tr}e^{(\widehat{H}_0\mu \widehat{N})/k_BT}}}`$ (25)
$`=`$ $`{\displaystyle \frac{1}{e^{(E_𝐢\mu )/k_BT}1}}.`$ (26)
Then, the thermodynamic potential $`\mathrm{\Omega }=pV`$ can be given as

$`\mathrm{\Omega }`$ $`=`$ $`k_BT\mathrm{ln}\mathrm{Tr}e^{(\widehat{H}_0\mu \widehat{N})/k_BT},`$ (27)
$``$ $`\mathrm{}\omega {\displaystyle \frac{4\pi a}{a_{ho}}}{\displaystyle \underset{𝐢,𝐣}{}}U_{\mathrm{𝐢𝐣}}n_𝐢n_𝐣`$ (29)
$`+k_BT{\displaystyle \underset{𝐢}{}}\mathrm{ln}(1e^{(E_𝐢\mu )/k_BT}),`$
where the first term in the right hand side is the first order correction. We note that Eq. (26) can be also determined by the condition

$`0`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{n_𝐢}}\right)_{T,\mu ,V,n_𝐣n_𝐢}`$ (30)
$`=`$ $`{\displaystyle \underset{𝐣}{}}U_{\mathrm{𝐢𝐣}}n_𝐣+{\displaystyle \underset{𝐣}{}}U_{\mathrm{𝐢𝐣}}{\displaystyle \frac{1}{e^{(E_𝐣\mu )/k_BT}1}}.`$ (31)
This condition means that the thermodynamic potential $`\mathrm{\Omega }`$ is a minimum with respect to any change of state at constant $`T`$, $`V`$ and $`\mu `$ in a state of thermal equilibrium.

Now we derive the thermodynamic quantities from Eqs. (29) and (31). First, we obtain the total number $`N`$ from Eq. (29) with the condition of Eq. (31)

$`N`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{\mu }}\right)_{T,V}`$ (32)
$`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{\mu }}\right)_{T,V,\mu ,n_𝐢}{\displaystyle \underset{𝐢}{}}\left({\displaystyle \frac{\mathrm{\Omega }}{n_𝐢}}\right)\left({\displaystyle \frac{n_𝐢}{\mu }}\right)`$ (33)
$`=`$ $`{\displaystyle \underset{𝐢}{}}{\displaystyle \frac{1}{e^{(E_𝐢\mu )/k_BT}1}}.`$ (34)
We can determine the chemical potential $`\mu `$ and the energy $`E_𝐢`$ from Eqs. (24) and (34) which should be solved self-consistently. Then, we obtain the thermodynamic quantities through Eq. (29) with Eq. (31). For example, the entropy $`S`$ can be given by

$`S`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{T}}\right)_{V,\mu }`$ (35)
$`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{T}}\right)_{V,\mu ,n_𝐢}{\displaystyle \underset{𝐢}{}}\left({\displaystyle \frac{\mathrm{\Omega }}{n_𝐢}}\right)\left({\displaystyle \frac{n_𝐢}{T}}\right)`$ (36)
$`=`$ $`k_B{\displaystyle \underset{𝐢}{}}\left[(1+n_𝐢)\mathrm{ln}(1+n_𝐢)n_𝐢\mathrm{ln}n_𝐢\right],`$ (37)
where we also use the condition Eq. (31). Finally, the internal energy $`E`$ can be given as

$`E`$ $`=`$ $`{\displaystyle \underset{𝐢}{}}{\displaystyle \frac{E_𝐢}{e^{(E_𝐢\mu )/k_BT}1}}.`$ (38)
We obtain the thermodynamic quantities from Eq. (29) by solving $`E_𝐢`$ and $`\mu `$ from Eqs. (24) and (29) self-consistently. This is a mean field approach within the first order perturbation theory. The chemical potential cannot exceed the energy $`E_𝐢`$ at any temperature. At the transition temperature, the chemical potential is equal to the ground state energy. Therefore, the transition temperature has to be determined by the chemical potential. We note that the second term of Eq. (24) shifts the transition temperature. If we choose $`U_{𝐢,𝐣}`$ as a constant i.e. $`U`$, then it only causes the constant shift for all $`E_𝐢`$, but the transition temperature does not change. Therefore, the second term of Eq. (24) plays an important role to understand the interacting Bose gas. We also note that, when we put $`g=0`$, Eq. (34) reproduces the total number of the noninteracting Bose gas confined in the spherical harmonic potential, because there are no constraints in between the chemical potential $`\mu `$ and the energy $`E_𝐢`$. Therefore, when we consider the interaction effect, we can also investigate the finite size effect at the same time.

## III Numerical Results

We present the numerical results of thermodynamic quantities which appear in Eqs. (34) and (38). First, we numerically solve the chemical potential $`\mu `$ and the energy $`E_𝐢`$ by using Eqs. (24) and (34) with the iteration method. Then, we obtain $`\mu `$ and $`E_𝐢`$ as a function of $`T`$ and $`N`$. If the temperature of this system is lowered at fixed $`N`$, the chemical potential $`\mu `$ given by Eq. (34) increases. It finally reaches the ground state energy $`E_0`$, since the chemical potential cannot exceed the ground state energy. Therefore, we can define the transition temperature $`T_c`$ as $`\mu =E_0`$ at $`T=T_c`$. We note that the transition temperature $`T_c`$ for the noninteracting Bose gas confined in the spherical harmonic potential can be given as

$`T_c`$ $`=`$ $`T_0\left(1{\displaystyle \frac{\zeta (2)}{2\zeta (3)^{2/3}}}{\displaystyle \frac{1}{N^{1/3}}}\right),`$ (39)
where the second term is the finite size effect, and $`T_0`$ denote the transition temperature in the thermodynamic limit which can be written as

$`T_0`$ $`=`$ $`{\displaystyle \frac{\mathrm{}\omega }{k_B}}\left({\displaystyle \frac{N}{\zeta (3)}}\right)^{1/3}.`$ (40)
Because the interaction effect also shifts the transition temperature $`T_c`$, we have to redefine the transition temperature.

For $`T>T_0`$, we solve Eqs. (24) and (34) with fixed $`N`$ and given $`T`$. Here, the total number $`N`$ and the strength of the interaction with the repulsive force are chosen to be $`N=1000`$ and $`a/a_{ho}=7\times 10^3`$, respectively. In Fig. 1, we show the calculated values of the chemical potential $`\mu /k_BT_0`$ as a function of $`T/T_0`$. As can be seen, when $`T/T_00.91`$, $`\mu /k_BT_0`$ approaches to a finite value. This means that the chemical potential $`\mu `$ reaches the ground state energy $`E_0`$. Therefore, the transition temperature becomes $`T/T_00.91`$ (see, Fig. 3). Comparing this result with the transition temperature $`T_c/T_00.93`$ given by Eq. (39), the interaction effect further shifts the transition temperature about $`2`$ percent toward the lower temperature.

For $`T<T_0`$, the chemical potential $`\mu `$ is equal to the ground state energy $`E_0`$. In this case, since the ground state occupation numbers $`N_0`$ diverges, we have to rewrite Eq. (24) as follows

$`E_𝐢`$ $`=`$ $`\epsilon _𝐢+\mathrm{}\omega {\displaystyle \frac{8\pi a}{a_{ho}}}U_{𝐢,0}N`$ (42)
$`+\mathrm{}\omega {\displaystyle \frac{8\pi a}{a_{ho}}}{\displaystyle \underset{𝐣0}{}}n_𝐣(U_{𝐢,𝐣}U_{𝐢,0}),`$
where we use the relation $`N=N_0+_{𝐣0}n_𝐣`$. Under this condition, we solve Eq. (42) with $`\mu /k_BT_0=0.178`$. We note that $`\mu /k_BT_0`$ is $`0.159`$ for the noninteracting model at the transition temperature, since $`\mu `$ is $`3/2\mathrm{}\omega `$ in our approach.

In Fig. 2, we show the condensation fraction as a function of $`T/T_0`$. The dotted line refers to the condensation fraction of the noninteracting Bose gas in the thermodynamic limit, and the solid line represents the condensation fraction including the finite size effect (see, Eq. (1)). Finally, the condensation fraction with the interaction and the finite size effects is shown as solid circles. These correspond to each other at $`T>T_0`$. This means that we can ignore the atom-atom interaction because the number density is low at $`T>T_0`$. On the other hand, at near the transition temperature, atoms must occupy the ground state energy $`E_0`$. Since the number density becomes high at this temperature, the atom-atom interaction becomes strong, and shifts the transition temperature. Although the finite size effect can also shift the transition temperature, it disappears at large $`N`$ limit. Therefore, the interaction effect is more important than the finite size effect.

In Fig. 4, we show the results for the internal energy $`E/Nk_BT_0`$ as a function of $`T/T_0`$. It reveals that the internal energy curve changes slope near the transition temperature. We note that the specific heat can be calculated by differentiating the internal energy with respect to the temperature $`T`$. Therefore, this result means that the specific heat is discontinuous at the transition temperature. As can be seen from the result of Ensher et al. , the derivative of the internal energy is discontinuous at the transition temperature. This behavior of our calculated internal energy is in a good agreement with the experimental result.

Finally, we compare our results with the semiclassical results. In the semiclassical approximation, the dimensionless parameter $`E_{int}/E_{kin}=Na/a_{ho}`$ is important to discuss the two-body atomic interaction energy $`E_{int}`$ compared to the kinetic energy $`E_{kin}`$. For the Thomas-Fermi (TF) limit $`(Na/a_{ho}1)`$, the chemical potential $`\mu _0`$ can be written in terms of this parameter

$`\mu _0`$ $`=`$ $`{\displaystyle \frac{\mathrm{}\omega }{2}}\left(15N{\displaystyle \frac{a}{a_{ho}}}\right)^{2/5}.`$ (43)
In our case, $`\mu _0/k_BT_0`$ is 0.342. This is about twice as large as the $`\mu /k_BT_0=0.178`$. Since we take the parameter $`Na/a_{ho}1`$, it does not achieve the TF limit. Therefore, the kinetic term plays an important role in this value. Further, the interaction term $`U_{\mathrm{𝐢𝐣}}`$ in Eq. (24) suppresses the total number $`N`$ of $`Na/a_{ho}`$. The shift of the transition temperature $`\delta T_0/T_0`$ in the semiclassical approximation for the large $`N`$ limit can be given as

$`{\displaystyle \frac{\delta T_0}{T_0}}`$ $`=`$ $`1.3{\displaystyle \frac{a}{a_{ho}}}N^{1/6}.`$ (44)
In our case, the shift of the transition temperature $`\delta T_0/T_0`$ is about 3 percent. Because of the interaction term $`U_{\mathrm{𝐢𝐣}}`$, the transition temperature is also suppressed by its term in our model.

## IV Conclusions

We have presented the thermodynamic potential of the interacting Bose gas confined in the spherical harmonic potential within a mean field approach. In this approach, we minimize the interacting Hamiltonian $`\widehat{H}^{}`$ under the condition that the mean occupation numbers $`\nu _𝐢`$ should be the Bose distribution function. Then, we can derive the thermodynamic potential in a self-consistent fashion. With this thermodynamic potential, we obtain the thermodynamic quantities which have the same form as that of the ideal Bose gas. Further, we carry out the numerical calculations of the chemical potential, the condensate fraction and the internal energy. In particular, because the energy $`E_𝐢`$ also varies as a function of $`T`$ and $`N`$ in our approach, we have to solve $`E_𝐢`$ and $`\mu `$ with the conditions that $`\mu `$ cannot exceed $`E_𝐢`$ at any temperature. Further, since we take the parameter $`Na/a_{ho}1`$ which corresponds neither to the TF limit nor to the noninteracting limit for the semiclassical point of view, we cannot expect $`\mu `$ from the semiclassical approximation. Then, we find that the interaction effect shifts the transition temperature $`T_0`$ given by $`\mu =E_0`$ at $`T=T_c`$ about 2 percent toward the lower temperature compared with the noninteracting model including the finite size effect.

###### Acknowledgements.

I thank T. Fujita for careful reading.

## Figure Captions


text:
# Asymmetric Nuclear Matter from Extended Brueckner-Hartree-Fock Approach

## 1 Introduction

Within the general interest for the equation of state ( EOS ) of nuclear matter in nuclear physics as well as in nuclear astrophysics, increasing attention is currently paid to the isospin degree of freedom.

The EOS of isospin asymmetric nuclear matter plays a central role for our understanding of astrophysical phenomena like supernova explosions, neutron stars structure, X-ray bursts, neutron stars merging and possibly $`\gamma `$-ray bursts. The study of asymmetric nuclear matter represents also the first step for a microscopic theory of the structure of nuclei far from the valley of beta stability. This “terra incognita” is going to be explored in the near future thanks to a new generation of experimental facilities with high intensity radioactive ion beams. Moreover, dynamical simulations of collisions between neutron-rich nuclei show that the main reaction mechanisms including fragmentation are quite sensitive to the density dependence of the nuclear symmetry energy . Such calculations mainly make use of phenomenological Skyrme-like forces where the symmetry energy at high density can also be in strong disagreement with the one extracted from the microscopic predictions.

On a microscopic basis the EOS of asymmetric nuclear matter has been studied within the variational approach  as well as relativistic  and non-relativistic  Brueckner-Bethe-Goldstone (BBG) theory. Within the Brueckner-Hartree-Fock (BHF) approximation to the BBG theory a systematic study of isospin effects on the EOS of asymmetric nuclear matter has been carried out in Ref. , where a separable version  of the Paris potential  was adopted to describe the two-body nuclear force.

Beside the bulk properties (EOS), the authors of Ref.  focused also on the single-particle (s.p.) properties of neutrons and protons in isospin-asymmetric nuclear medium. The neutron and proton s.p. potentials were calculated  to the lowest order in the Brueckner reaction matrix (BHF approximation), using the so-called continuous choice .

Motivated by the renewed interest in this subject, in the present paper, we report an extension of the calculations of Ref.  along the following lines. First, in the calculations we make use of a different realistic nucleon-nucleon ( $`NN`$ ) potential, i.e., the full Argonne $`V_{14}`$ potential , which enables us to take into account a larger number of partial waves with respect to the calculation  with the separable Paris potential. These additional partial waves ( $`3L6`$ ) give a non-negligible contribution both to the EOS and the nucleon mean field, expecially in the high-density region, which is relevant for applications in astrophysics as well as heavy ion physics.

Second, the Bethe-Goldstone equation is now solved for the complex $`G`$-matrix. This enables us to calculate the complex nuclear mean field and some closely related quantities such as the optical potential and the mean free path.

Third, according to the Landau definition of quasi-particle energy (for an extended discussion see Ref. ) in the calculations of the mass operator (nucleon self-energy) and single-particle properties, we go beyond the BHF approximation by including some higher-order correlation contributions. In particular, in the present work, we include the so-called rearrangement term $`M_2`$ which is a second order diagram in the $`G`$-matrix and accounts for particle-hole excitations in nuclear matter ground state. Next we consider also the renormalization contributions of the third and forth order in the $`G`$-matrix, which account for the partial depletion of the neutron and proton Fermi seas due to the nuclear correlations . It has been shown, in the case of pure neutron matter  and also symmetric nuclear matter  that the new terms give a large contribution to s.p. properties like the mean field and the nucleon effective mass. We will refer to the present approach to compute nuclear s.p. properties as extended Brueckner–Hartree–Fock ( EBHF ) approximation .

As is well known, the BHF approximation largely violates the Hugenholtz-Van Hove ( HVH ) theorem , which basically measures the consistency of a given order of approximation in a perturbative approach. In symmetric nuclear matter, the inclusion of the rearrangement contribution greatly improves the fulfillment of the HVH theorem . In the present paper, we study this problem in the case of asymmetric nuclear matter within the EBHF approximation.

## 2 EBHF and Nucleon Self-energy for Asymmetric Nuclear Matter

In this section the formalism of the Brueckner-Bethe-Goldstone ( BBG ) theory is described for the case of asymmetry nuclear matter . The proton and the neutron Fermi momenta are related to their corresponding densities $`\rho _\mathrm{p}`$ and $`\rho _\mathrm{n}`$ through the relations

$`k_F^\mathrm{p}=[{\displaystyle \frac{3\pi ^2}{2}}(1\beta )\rho ]^{1/3},`$
$`k_F^\mathrm{n}=[{\displaystyle \frac{3\pi ^2}{2}}(1+\beta )\rho ]^{1/3},`$ (1)
where $`\rho =\rho _\mathrm{p}+\rho _\mathrm{n}`$ is the total density, and $`\beta =(\rho _\mathrm{n}\rho _\mathrm{p})/\rho `$ the asymmetry parameter determining the neutron excess (from now on we assume $`\rho _\mathrm{n}\rho _\mathrm{p}`$).

The starting point in BBG theory, is the Brueckner reaction matrix $`G`$, which in the case of asymmetric nuclear matter depends also on the isospin components of the two colliding nucleons. The $`G`$-matrix satisfies the Bethe-Goldstone equation,

$`G(\rho ,\beta ;\omega )=v_{NN}+v_{NN}{\displaystyle \underset{k_1k_2}{}}{\displaystyle \frac{|k_1k_2Q(k_1,k_2)k_1k_2|}{\omega ϵ(k_1)ϵ(k_2)+i\eta }}G(\rho ,\beta ;\omega ),`$ (2)
where $`v_{NN}`$ is the two-body nuclear interaction and $`\omega `$ the starting energy. Here $`k(\stackrel{}{k},\sigma ,\tau )`$ denotes s.p. momentum, $`z`$-components of spin and isospin, respectively.

The $`G`$-matrix can be considered as an in-medium effective interaction between two nucleons. The surrounding nucleons renormalize the bare $`NN`$ interaction via the Pauli blocking and the nuclear mean field. The Pauli operator, defined as

$`Q(k_1,k_2)=[1n(k_1)][1n(k_2)],`$ (3)
prevents two nucleons in intermediate states from scattering into states inside their respective Fermi seas. By $`n(k)`$ we denote the Fermi distribution function, which at zero temperature is given by the step function $`\theta (kk_F^\tau )`$ ( uncorrelated ground state ). The s. p. energy 222In the present work, we assume the neutron and proton rest masses equal to their average value $`m`$

$`ϵ(k)={\displaystyle \frac{\mathrm{}^2k^2}{2m}}+U(k),`$ (4)
appearing in the energy denominator of Eq. (2), involves the auxiliary potential $`U(k)`$, which controls the convergence rate of the hole-line expansion. Within the BHF approximation the neutron and proton s.p. auxiliary potentials are calculated from the real part of the on-shell antisymmetrized $`G`$-matrix, via the relation

$`U(k)={\displaystyle \underset{k^{}}{}}n(k^{})\mathrm{Re}kk^{}|G(ϵ(k)+ϵ(k^{}))|kk^{}_A.`$ (5)
Here we adopt the continuous choice  for the auxiliary s.p. potential. In this context, the auxiliary potential has the physical meaning of the mean field that each nucleon feels during its propagation between two successive scatterings.

In the BHF approximation, Eqs. (2),(4), and (5) are solved self-consistently for given total density $`\rho `$ and asymmetry $`\beta `$. Then the energy per particle is evaluated at the lowest order ( two hole-line diagrams ) of the BBG hole-line expansion ( see Ref.  for the case of asymmetric matter ).

### 2.1 Mass operator and quasi-particle energy

One of the main purposes of the present paper is to calculate s.p. properties of neutrons and protons in asymmetric matter going beyond the BHF approximation. To this end we introduce the mass operator

$`M^\tau (k,\omega )=V^\tau (k,\omega )+iW^\tau (k,\omega ),`$ (6)
which is a complex quantity and can be identified with the potential energy felt by a neutron ($`\tau =n`$) or a proton ($`\tau =p`$) with momentum $`\stackrel{}{k}`$ and energy $`\omega `$ in asymmetric nuclear matter ( hereafter, we will write out explicitly the isospin index $`\tau `$ ). In the same spirit of the BBG theory, the mass operator $`M^\tau (k,\omega )`$ can be expanded in a perturbation series according to the number of hole lines  and various terms of this expansion can be represented by means of Goldstone diagrams a few of which are shown in Fig. 1.

In analogy with the case of symmetric nuclear matter, the neutron and proton quasi-particle energies $`E^\tau (k)`$ are the solutions of the energy-momentum relation

$`E^\tau (k)={\displaystyle \frac{\mathrm{}^2k^2}{2m}}+V^\tau (k,E^\tau (k)),`$ (7)
i.e., $`E^\tau (k)`$ is obtained from the on-shell values of the real part of the mass operator.

To the lowest order in the hole-line expansion the mass operator is given by ( diagrams of Fig. 2 )

$`M_1^\tau (k,\omega )={\displaystyle \underset{\tau ^{}}{}}{\displaystyle \underset{\stackrel{}{k}^{}\sigma ^{}}{}}n^\tau ^{}(k^{})kk^{}|G^{\tau \tau ^{}}(\omega +ϵ^\tau ^{}(k^{}))|kk^{}_A{\displaystyle \underset{\tau ^{}}{}}M_1^{\tau \tau ^{}}(k,\omega ).`$ (8)
In this approximation the quasi-particle energy $`E_1^\tau (k)`$ coincides with the BHF s.p. energy given by Eqs. (4) and (5), i.e., $`E_1^\tau (k)=ϵ^\tau (k)`$.

### 2.2 The rearrangement contribution to the s.p. energy

The next contribution to the perturbative expansion of the mass operator is given by the so-called rearrangement term $`M_2^\tau (k,\omega )`$ . The associated Goldstone diagrams are shown in Fig. 3. $`M_2^\tau `$ is a second-order diagram in the $`G`$-matrix and accounts for particle-hole excitations in nuclear matter. Its expression, extended to asymmetric nuclear matter, reads

$`M_2^\tau (k,\omega )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\tau ^{}}{}}{\displaystyle \underset{\stackrel{}{k}^{}\sigma ^{}}{}}(1n^\tau ^{}(k^{})){\displaystyle \underset{k_1k_2}{}}n^\tau (k_1)n^\tau ^{}(k_2){\displaystyle \frac{|kk^{}|G^{\tau \tau ^{}}(ϵ^\tau (k_1)+ϵ^\tau ^{}(k_2))|k_1k_2_A|^2}{\omega +ϵ^\tau ^{}(k^{})ϵ^\tau (k_1)ϵ^\tau ^{}(k_2)i\eta }}`$ (9)
$``$ $`{\displaystyle \underset{\tau ^{}}{}}M_2^{\tau \tau ^{}}(k,\omega ),`$ (10)
where $`ϵ^\tau (k)`$ is the s.p. spectrum in BHF approximation, given by Eqs. (4) and (5). In this approximation for the mass operator \[ i.e., $`M^\tau (k,\omega )M_1^\tau (k,\omega )+M_2^\tau (k,\omega )`$ \], the quasi-particle energy (7) is given by the approximate relation

$`E_2^\tau (k)`$ $`=`$ $`E_1^\tau (k)+Z_2^\tau (k)V_2^\tau (k,E_1^\tau (k))`$ (11)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2k^2}{2m}}+V_1^\tau (k,E_1^\tau (k))+Z_2^\tau (k)V_2^\tau (k,E_1^\tau (k)),`$ (12)
where

$`Z_2^\tau (k)=\left\{1{\displaystyle \frac{}{\omega }}\left[V_1^\tau (k,\omega )+V_2^\tau (k,\omega )\right]\right\}_{\omega =E_1^\tau (k)}^1`$ (13)
is an approximation of the quasi-particle strength for asymmetric nuclear matter

$`Z^\tau (k)=\left\{1{\displaystyle \frac{}{\omega }}\left[V^\tau (k,\omega )\right]\right\}_{\omega =E^\tau (k)}^1.`$ (14)
### 2.3 The renormalization contributions to the s.p. energy

Due to many-body correlations the two Fermi seas are partially depleted, and the correlated momentum distributions $`\stackrel{~}{n}^\tau (k)`$ differ from the uncorrelated ones $`n^\tau (k)=\theta (kk_F^\tau )`$. To account for this physical effect, one considers the contribution $`M_3^\tau (k,\omega )`$ ( last diagram of Fig. 1 ) given by

$`M_3^\tau (k,\omega )={\displaystyle \underset{\tau ^{}}{}}{\displaystyle \underset{\stackrel{}{h^{}}\sigma ^{}}{}}\kappa _2^\tau ^{}(h^{})kh^{}|G^{\tau \tau ^{}}(\omega +ϵ^\tau ^{}(h^{}))|kh^{}_A,`$ (15)
where $`h^{}`$ refers to “hole” state with momentum smaller than $`k_\mathrm{F}^\tau `$, and

$`\kappa _2^\tau ^{}(h^{})=\left[{\displaystyle \frac{}{\omega }}M_1^\tau ^{}(h^{},\omega )\right]_{\omega =ϵ^\tau ^{}(h^{})}`$
is at the lowest order the depletion of neutron (proton) Fermi sea , i.e., $`\kappa _2^\tau ^{}(h^{})`$ is the probability that a neutron (proton) hole-state $`(|\stackrel{}{h}^{}|k_F^\tau ^{})`$ is empty. Let us consider now the sum

$`\stackrel{~}{M}_1^\tau (k,\omega )`$ $``$ $`M_1^\tau (k,\omega )+M_3^\tau (k,\omega )`$ (16)
$`=`$ $`{\displaystyle \underset{\tau ^{}}{}}{\displaystyle \underset{\stackrel{}{h}^{}\sigma ^{}}{}}\left[1\kappa _2^\tau ^{}(h^{})\right]kh^{}|G^{\tau \tau ^{}}(\omega +ϵ^\tau ^{}(h^{}))|kh^{}_A`$ (17)
$`=`$ $`{\displaystyle \underset{\tau ^{}}{}}{\displaystyle \underset{\stackrel{}{h}^{}\sigma ^{}}{}}\stackrel{~}{n}_2^\tau ^{}(h^{})kh^{}|G^{\tau \tau ^{}}(\omega +ϵ^\tau ^{}(h^{}))|kh^{}_A,`$ (18)
$`\stackrel{~}{n}_2^\tau ^{}(h^{})=\left[1\kappa _2^\tau ^{}(h^{})\right]`$ being the second-order approximation for the correlated momentum distribution. $`\stackrel{~}{M}_1^\tau (k,\omega )`$ is the so-called renormalized BHF approximation for the off-shell mass operator \[ compare to Eq. (8) \].

An accurate approximation consists in using the average value of the depletion, which is

$`\kappa ^\tau ^{}=\kappa _2^\tau ^{}\left(h^{}=0.75k_\mathrm{F}^\tau ^{}\right).`$ (19)
Then Eqs. (15) and (18) yield

$`M_3^\tau (k,\omega ){\displaystyle \underset{\tau ^{}}{}}\kappa ^\tau ^{}M_1^{\tau \tau ^{}}(k,\omega ),`$ (20)
$`\stackrel{~}{M}_1^\tau (k,\omega ){\displaystyle \underset{\tau ^{}}{}}\left[1\kappa ^\tau ^{}\right]M_1^{\tau \tau ^{}}(k,\omega ).`$ (21)
From the similar considerations, a renormalization correction should also be brought from the four hole-line terms to the second-order contribution to $`M_2^\tau `$ in order to take into account the fact that the hole-state $`k_1`$ in Eq. (10) is partially empty ( see also Ref.  for symmetric nuclear matter ). Along the same line of the previous correction one gets the renormalized $`M_2`$, which is approximately given by

$`\stackrel{~}{M}_2^\tau (k,\omega )={\displaystyle \underset{\tau ^{}}{}}\left[1\kappa ^\tau ^{}\right]M_2^{\tau \tau ^{}}(k,\omega ).`$ (22)
The renormalized contributions can also be traced to the functional dependence of the $`G`$-matrix on the quasi-particle occupation numbers within the Landau theory of Fermi liquids. It can be shown, in fact, that taking the functional derivative of the binding energy (at two hole-line level) includes also the terms of the third and fourth order in the self-energy, the effect of which has just been discussed. Taking into account all the corrections discussed above, from Eq. (7) one can get the following expression for the quasi-particle energy ,

$`E^\tau (k)`$ $``$ $`{\displaystyle \frac{\mathrm{}^2k^2}{2m}}+V_1^\tau (k,E_1^\tau (k))`$ (23)
$`+`$ $`Z^\tau (k){\displaystyle \underset{\tau ^{}}{}}\left[\kappa ^\tau ^{}V_1^{\tau \tau ^{}}(k,E_1^\tau (k))+(1\kappa ^\tau ^{})V_2^{\tau \tau ^{}}(k,E_1^\tau (k))\right],`$ (24)
where

$`Z_3^\tau (k)=\left\{1{\displaystyle \underset{\tau ^{}}{}}(1\kappa ^\tau ^{}){\displaystyle \frac{}{\omega }}\left[V_1^{\tau \tau ^{}}(k,\omega )+V_2^{\tau \tau ^{}}(k,\omega )\right]\right\}_{\omega =E_1^\tau (k)}^1.`$ (25)
In the following we refer to this approximation for the quasi-particle energy as the extended Brueckner-Hartree-Fock (EBHF) approximation .

### 2.4 Partial wave expansion and angular averaging

After the usual angular averaging on the Pauli operator and the energy denominator , the Bethe-Goldstone equation can be expanded in partial waves,

$`G_{\alpha LL^{}}^{\tau \tau ^{}}(q,q^{},P,\omega )=v_{\alpha LL^{}}(q,q^{})`$ (26)
$`+`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle \underset{L^{\prime \prime }}{}}{\displaystyle q^{\prime \prime 2}dq^{\prime \prime }v_{\alpha LL^{\prime \prime }}(q,q^{\prime \prime })\frac{Q^{\tau \tau ^{}}(q^{\prime \prime },P)}{\omega e_{12}^{\tau \tau ^{}}(q^{\prime \prime },P)+i\eta }G_{\alpha L^{\prime \prime }L^{}}^{\tau \tau ^{}}(q^{\prime \prime },q^{},P,\omega )},`$ (27)
where $`\stackrel{}{q}=(\stackrel{}{k}_1\stackrel{}{k}_2)/2`$ and $`\stackrel{}{P}=\stackrel{}{k}_1+\stackrel{}{k}_2`$ are the relative momentum and total momentum, respectively. $`e_{12}^{\tau \tau ^{}}(q^{\prime \prime },P)=ϵ^\tau (k_1)+ϵ^\tau ^{}(k_2)`$ is the angle average of the energy denominator. The angular-averaged Pauli operator is

(i) for $`\tau =\tau ^{}`$ (neutron-neutron or proton-proton),

$`Q^{\tau \tau }(q,P)`$ $`=`$ $`\{\begin{array}{ccc}\mathrm{min}(1,\xi _\tau )\hfill & & \mathrm{if}\xi _\tau 0\\ 0\hfill & & \mathrm{otherwise},\end{array}`$ (30)
(ii) for $`\tau \tau ^{}`$ (neutron-proton or proton-neutron),

$`Q^{\tau \tau ^{}}(q,P)`$ $`=`$ $`\{\begin{array}{ccc}{\displaystyle \frac{1}{2}}[\mathrm{min}(1,\xi _\mathrm{p})+\mathrm{min}(1,\xi _\mathrm{n})]\hfill & & \mathrm{if}\xi _\mathrm{n}\xi _\mathrm{p},\xi _\mathrm{n}1\\ 0\hfill & & \mathrm{otherwise},\end{array}`$ (33)
where

$`\xi _\tau ={\displaystyle \frac{P^2/4+q^2(k_\mathrm{F}^\tau )^2}{Pq}}.`$
The mass operators $`M_1`$ and $`M_2`$ become

$`M_1^{\tau \tau ^{}}(k,\omega )`$ $`=`$ $`{\displaystyle \frac{1+\delta _{\tau ,\tau ^{}}}{2\pi }}{\displaystyle \underset{\alpha L}{}}(2J+1){\displaystyle _0^{k_\mathrm{F}^\tau ^{}}}k^2dk^{}\mathrm{sin}\theta \mathrm{d}\theta G_{\alpha LL}^{\tau \tau ^{}}(q,q,P,\omega +ϵ^\tau ^{}(k^{})),`$ (34)
$`M_2^{\tau \tau ^{}}(k,\omega )`$ $`=`$ $`{\displaystyle \frac{2(1+\delta _{\tau ,\tau ^{}})}{\pi ^2k}}{\displaystyle \underset{\alpha LL^{}}{}}(2J+1){\displaystyle qdqPdP\left[1n^\tau ^{}\left(\sqrt{P^2/2+2q^2k^2}\right)\right]}`$ (35)
$`\times `$ $`{\displaystyle q^2dq^{}R^{\tau \tau ^{}}(q^{},P)\frac{|G_{\alpha LL^{}}^{\tau \tau ^{}}(q,q^{},P,e_{12}^{\tau \tau ^{}}(q^{},P))|^2}{\omega +ϵ^\tau ^{}\left(\sqrt{P^2/2+2q^2k^2}\right)e_{12}^{\tau \tau ^{}}(q^{},P)i\eta }}.`$ (36)
The integrations of $`q`$ and $`P`$ in the expression of $`M_2`$ are limited to

$`q_{\mathrm{min}}`$ $`=`$ $`\{\begin{array}{ccc}\mathrm{max}[0,k{\displaystyle \frac{1}{2}}(k_\mathrm{F}^\tau +k_\mathrm{F}^\tau ^{}),{\displaystyle \frac{1}{2}}(k_\mathrm{F}^\tau ^{}k)]\hfill & \mathrm{if}\hfill & kk_\mathrm{F}^\tau \\ & & \\ \mathrm{max}[{\displaystyle \frac{1}{2}}\sqrt{2k^22(k_\mathrm{F}^\tau )^2+(k_\mathrm{F}^\tau k_\mathrm{F}^\tau ^{})^2},{\displaystyle \frac{1}{2}}(k+k_\mathrm{F}^\tau ^{})]\hfill & \mathrm{if}\hfill & k>k_\mathrm{F}^\tau \end{array}`$ (39)
$`q_{\mathrm{max}}`$ $`=`$ $`k+{\displaystyle \frac{1}{2}}(k_\mathrm{F}^\tau +k_\mathrm{F}^\tau ^{})`$ (40)
and

$`P_{\mathrm{min}}`$ $`=`$ $`\{\begin{array}{ccc}\mathrm{max}[2(kq),\sqrt{2k^2+2(k_\mathrm{F}^\tau ^{})^24q^2}]\hfill & \mathrm{if}\hfill & q\frac{1}{2}(k+k_\mathrm{F}^\tau ^{})\\ & & \\ 2|qk|\hfill & \mathrm{if}\hfill & k>\frac{1}{2}(k+k_\mathrm{F}^\tau ^{})\end{array}`$ (43)
$`P_{\mathrm{max}}`$ $`=`$ $`\mathrm{min}[2(k+q),(k_\mathrm{F}^\tau +k_\mathrm{F}^\tau ^{})].`$ (44)
The angular averaging of the anti-Pauli operator $`R^{\tau \tau ^{}}(k_1,k_2)n^\tau (k_1)n^\tau ^{}(k_2)`$ can be written as

(i) for $`\tau =\tau ^{}`$, i.e., neutron-neutron or proton-proton

$`R^{\tau \tau }(q,P)`$ $`=`$ $`\{\begin{array}{ccc}\mathrm{min}(1,\eta _\tau )\hfill & & \mathrm{if}\eta _\tau 0\\ 0\hfill & & \mathrm{otherwise},\end{array}`$ (47)
(ii) for $`\tau \tau ^{}`$, i.e., neutron-proton or proton-neutron

$`R^{\tau \tau ^{}}(q,P)`$ $`=`$ $`\{\begin{array}{ccc}{\displaystyle \frac{1}{2}}[\mathrm{min}(1,\eta _\mathrm{p})+\mathrm{min}(1,\eta _\mathrm{n})]\hfill & & \mathrm{if}\eta _\mathrm{p}\eta _\mathrm{n},\eta _\mathrm{p}1\\ 0\hfill & & \mathrm{otherwise},\end{array}`$ (50)
where $`\eta ^\tau =\xi ^\tau `$.

## 3 Results

We have performed a set of nuclear matter calculations for the asymmetric case within the EBHF approximation. Three different densities have been selected: $`\rho _0/2`$, $`\rho _0`$ and $`2\rho _0`$, being $`\rho _0=0.17fm^3`$ the saturation density of symmetric nuclear matter. For each density the whole range of asymmetry parameter ( $`0\beta 1`$ ) has been spanned. The self-consistent solution of the Bethe-Goldstone equation yielding simultaneously $`G`$-matrix and auxiliary potential $`U^\tau (k)`$ needed five iterations to reach a satisfactory convergence. The bare potential adopted as input in the calculation was the Argonne $`V_{14}`$  with 24 channels up to $`L=6`$.

### 3.1 Symmetry Energy

In Fig. 4 (left panel) we report the results (symbols) for the energy per nucleon $`B(\rho ,\beta )`$, calculated self-consistently within the BHF approximation . $`B(\rho ,\beta )`$ is plotted as a function of $`\beta ^2`$, for three values of density. The numerical results lie on a linear fit performed with only the first three values of the asymmetry parameter. This proves that the empirical parabolic law

$`B(\rho ,\beta )=B(\rho ,0)+E_{sym}(\rho )\beta ^2,`$ (51)
taken from the nuclear mass table can be extended up to the highest asymmetry of nuclear matter, in good agreement with our previous BHF calculation with the separable Paris potential .

Equation (51) can be considered as the $`\beta ^2`$ expansion of the binding energy truncated at the lowest order. Only even powers of the asymmetry parameter $`\beta `$ may occur in the expansion for charge-independent $`NN`$ interactions, such as the Argonne $`V_{14}`$ used in the present work. A $`\beta ^4`$ contribution might arise at the three hole-line order of the BBG expansion. Unfortunately, no such calculation for $`B(\rho ,\beta )`$ has been done yet. However, it has been shown recently that the three hole-line contribution to the binding energy of symmetric nuclear matter  is rather small within the continuous choice. Therefore, we do not expect a large deviation from the parabolic law after including the three hole-line contribution. A deviation from the parabolic law could be expected at densities higher with respect to those considered in the present paper .

The symmetry energy is defined as

$`E_{sym}(\rho )={\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{^2B(\rho ,\beta )}{\beta ^2}}\right]_{\beta =0}.`$ (52)
Due to the simple $`\beta ^2`$-law the symmetry energy can be equivalently calculated as the difference between the binding energy of pure neutron matter and symmetric nuclear matter: $`E_{sym}(\rho )=B(\rho ,1)B(\rho ,0)`$, but one would refrain from using that recipe at very high density. The results of our BHF calculations for $`E_{sym}(\rho )`$ are depicted by the continuous curve in the right panel of Fig. 4. In the same figure, we also show the results from the variational approach using the same Argonne $`V_{14}`$ potential .

The systematic disagreement displayed by the two many-body approaches has been believed to be a shortcoming of the Brueckner approach in view of the fact that the BHF result lies above the variational one. However, in Ref.  (and similar works), the variational expectation value $`E_{var}`$ of the Hamiltonian is calculated in a diagrammatic cluster expansion (FHNC-SOC), which is of course truncated to some order. To estimate the convergence of this diagrammatic cluster expansion, we plot, in the same figure, the results of a lowest-order constrained variational calculation , which includes only two-body cluster contributions to $`E_{var}`$. Moreover, the variational trial wave function used in Ref.  does not contain the correlations which arise from $`\stackrel{}{L}^2`$, $`\stackrel{}{L}^2(\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2)`$, and $`(\stackrel{}{L}\stackrel{}{S})^2`$ terms of the nucleon-nucleon potential. Finally, spin-orbit correlations are not treated accurately, as discussed in the same paper . All these $`NN`$ correlations are included in a self-consistent way in the BHF approach. All the above-mentioned approximations could give large uncertainties in the calculated expectation value of the energy in the high density region. The same discrepancy has also been observed in our previous calculations for asymmetric nuclear matter  and also in neutron matter calculations .

From the previous discussion we guess that the nice agreement between our calculation and a lowest order constrained variational calculation  is fortuitous. On the other hand, the agreement up to $`\rho 0.24`$ fm-3 with the variational calculation including three-body force , also plotted in Fig. 4, is hardly understandable.

More recent versions of the $`NN`$ potential do not provide any appreciable difference of the symmetry energy from the present calculation except the CD-Bonn potential as discussed in Ref. . All Brueckner calculations predict the symmetry energy to increase with the nucleon density and no saturation is observed up to $`\rho =0.5fm^3`$ at variance with the preceeding variational results . In the relativistic mean-field theory this behavior is easily understood in terms of the $`\rho `$-meson exchange, which leads to a repulsive symmetry potential at all densities . In order to try to explain what happens in the nonrelativistic case, we report in Fig. 5 the different contributions to the symmetry energy, plotted as a function of density. The kinetic contribution monotonically increases as $`\rho ^{\frac{2}{3}}`$ according to the free Fermi-gas model. In the figure (right two panels), the isoscalar and isovector contributions of potential part are plotted separately. The density dependence of the symmetry energy is dominated at high density by the kinetic contribution, where the opposite behavior vs density of the potential contributions from $`T=0`$ and $`T=1`$ channels results in a very flat density dependence of the symmetry potential. As already found in the previous paper , the most important contribution to the $`T=0`$ component is due to the deuteron $`{}_{}{}^{3}S_{1}^{}`$-$`{}_{}{}^{3}D_{1}^{}`$ coupled channels of the interaction, which exhibits a maximum at $`\rho 0.3fm^3`$. This peak can be traced back to the behavior of the two components: the attractive $`{}_{}{}^{3}S_{1}^{}`$ channel dominates at low energy whereas the repulsive $`{}_{}{}^{3}D_{1}^{}`$ dominates at high energy. Two terms compensate each other at the energy $`E4E_F200MeV`$, where $`E_F`$ is the Fermi energy corresponding to $`\rho 0.3fm^3`$.

### 3.2 Single-particle energy

For asymmetric nuclear matter the neutron mass operator $`M^n`$ is different from the proton mass operator $`M^p`$. Moreover, as shown in Figs. 2 and 3 \[see also Eqs.(7) and (8)\], both of them can be split into two components: $`M^p=M^{pp}+M^{pn}`$ for protons and $`M^n=M^{nn}+M^{np}`$

for neutrons. In Fig. 6 the on-shell values of the real part (upper panels) of $`M_1^\tau `$ are reported as a function of the s.p. momentum, for different values of the asymmetry parameter $`\beta `$ at fixed density $`\rho =0.17fm^3`$. The proton mean field $`V_1^p(k)\text{Re}M_1^p(k)`$ becomes more attractive, while the neutron mean field $`V_1^n(k)\text{Re}M_1^n(k)`$ becomes more repulsive going from symmetric ($`\beta =0`$) to neutron ($`\beta =1`$) matter. The $`\beta `$ dependence of $`V_1^n`$ and $`V_1^p`$ is almost linear and nearly symmetric with respect to their common value at $`\beta =0`$. This result supports from a microscopic point of view the validity of the so-called Lane potential .

It is worth noticing that a crossing point occurs for both $`V_1^p`$ and $`V_1^n`$, where the isospin effect on neutron and proton mean fields versus $`\beta `$ is inverted. This behavior of the neutron and proton mean field can be understood in terms of phase-space arguments, as pointed out Ref.. To this end, we write the s.p. potentials $`V_1^n`$ and $`V_1^p`$ in terms of their components $`V_1^{\tau \tau ^{}}`$ \[defined according to Eq.(7)\]:

$`V_1^p(k){\displaystyle \frac{1}{2}}(1\beta )\rho G^{pp}+{\displaystyle \frac{1}{2}}(1+\beta )\rho G^{pn}`$ (53)
$`V_1^n(k){\displaystyle \frac{1}{2}}(1\beta )\rho G^{np}+{\displaystyle \frac{1}{2}}(1+\beta )\rho G^{nn}`$ (54)
where $`G^{pp}`$ is the average value of the real part of the matrix $`G^{pp}`$ in the proton Fermi sphere ($`|\stackrel{}{h}^{}|k_F^p`$), and $`G^{pn}`$ the average value of the real part of the matrix $`G^{pn}`$ in the neutron Fermi sphere ($`|\stackrel{}{h}^{}|k_F^n`$). $`G^{nn}`$ and $`G^{np}`$ have similar definitions. This approximation is suggested by the almost linear dependence of $`V_1^n`$ and $`V_1^p`$ on $`\beta `$ and, in fact, is numerically fulfilled with a good accuracy (see also Fig. 8). The crossing point in momentum space is determined by the occurrence of $`G^{pp}=G^{pn}`$ for $`V_1^p`$ and $`G^{np}=G^{nn}`$ for $`V_1^n`$ at a certain value of the momentum which does not depend upon $`\beta `$. A signature of the inversion of the isospin effect at the crossing point could be found in those collective observables measured in heavy ion collisions which are sensitive to the momentum dependence of the mean field.

The imaginary part $`W_1^\tau `$ of the mass operator $`M_1^\tau `$ is due to the virtual collisions of a single nucleon with a neutron or a proton of the background, promoting it to a particle state. $`W_1^\tau `$ is vanishing below the Fermi momentum $`k_F^\tau `$ due to the Pauli blocking. It is worth noticing that reducing the proton Fermi momentum implies a less Pauli blocking for protons. This means that high asymmetric nuclear matter is less transparent to the proton propagation.

The second-order terms of the on-shell mass operator $`M_2^\tau `$ are plotted in Fig. 7. The real part $`V_2^\tau `$ (upper panels) gives the contribution to the mean field due to the coupling of the single-particle motion with the ground state particle-hole excitations. As is well known, $`V_2^\tau `$ is repulsive and reduces to a large extent the pure BHF mean field $`V_1^\tau `$ which is too attractive compared with the phenomenological optical potential .

The imaginary part $`W_2^\tau `$ plays a role complementary to $`W_1^\tau `$: it describes the virtual collisions of a single nucleon of the background with an excited neutron or proton, making it to decay into a hole state. $`W_2^\tau `$ is vanishing above the Fermi momentum $`k_F^\tau `$.

In order to focus on only the isospin dependence, we plot in Fig. 8 the mass operator as a function of $`\beta `$ at $`k=0fm^1`$ except for the imaginary part of $`M_1`$ for which a value of $`k`$ above the Fermi momentum has to be taken.

The first-order contribution has the linear behavior for the real part as well as for the imaginary part as expected from phase-space arguments. The slope of $`|W_1^{pn}|`$ is more pronounced than that of $`|W_1^{np}|`$ since the neutron particle-hole excitations coupled to a proton in a particle state are more favored than the proton particle-hole excitations coupled to a neutron (see also Fig. 6).

The isospin dependence of the second-order contribution $`M_2^\tau `$ is affected by the coupling between the nucleon hole stats and particle-hole excitations \[see the bubble in Fig.3 and Eq.(10)\], which yields a nonlinear variation of the mixed components $`M_2^{pn}`$ and $`M_2^{np}`$ vs $`\beta `$. The nonlinearity is much more sizeable for $`V_2^{pn}`$ and $`W_2^{np}`$, which can be easily explained as a phase-space effect as well, i.e., of the interplay between the neutron and proton phase-spaces as increasing neutron excess.

### 3.3 Fermi Energy and Hugenholtz-Van Hove theorem

The EBHF approximation basically relies on the Landau definition of quasi-particle energy as showed in sec.I, whose relation to the Brueckner theory has been well established . Study of the HVH theorem within the EBHF approximation could provide an additional support to a proper definition of the quasi-particle energy and, at the same time, a more realistic evaluation of the Fermi energy.

Strictly speaking, the HVH theorem concerns only symmetric nuclear matter at saturation point ($`P=\mathrm{\hspace{0.17em}0}`$), and it states that the energy per nucleon must be exactly equal to the Fermi energy. In the case of asymmetric nuclear matter (two-component system) at zero temperature, the HVH theorem can be generalized via the thermodynamic relation

$`{\displaystyle \frac{E(\rho ,\beta )}{A}}+{\displaystyle \frac{P(\rho ,\beta )}{\rho }}=Y^pE_F^\mathrm{p}(\rho ,\beta )+Y^nE_F^\mathrm{n}(\rho ,\beta ),`$ (55)
$`P(\rho ,\beta )`$ being the pressure, $`Y^p=\rho _p/\rho `$ and $`Y^n=\rho _n/\rho `$ the proton and neutron fractions, respectively. The Fermi energy is calculated from the quasi-particle energy spectrum at Fermi surface according to Eq. (6).
TAB. I As a function of the asymmetry parameter (first comumn) are reported (in MeV) the physical quantities involved in the Hugenholtz-Van Hove theorem: pressure over $`\rho `$ (second column), energy per nucleon (third column) and “weighted” chemical potentials of asymmetric nuclear matter in different approximations, as discussed in the text. $`Y^p=Z/A`$ and $`Y^n=N/A`$ are the proton and neutron fractions, respectively. The total density is $`\rho =0.17fm^3`$.

| $`\beta `$ | $`P/\rho `$ | $`E/A`$ | $`P/\rho +E/A`$ | $`Y^\mathrm{p}E_F^p+Y^\mathrm{n}E_F^n`$ | $`Y^\mathrm{p}E_F^p+Y^\mathrm{n}E_F^n`$ | $`Y^\mathrm{p}E_F^p+Y^\mathrm{n}E_F^n`$ |
| --- | --- | --- | --- | --- | --- | --- |
|  |  |  |  | BHF | BHF+$`M_2`$ | EBHF |
| 0.0 | $`5.02`$ | $`15.92`$ | $`20.94`$ | $`34.27`$ | $`28.50`$ | $`19.28`$ |
| 0.2 | $`4.40`$ | $`14.73`$ | $`19.13`$ | $`32.29`$ | $`26.43`$ | $`17.35`$ |
| 0.4 | $`3.27`$ | $`11.36`$ | $`14.63`$ | $`26.28`$ | $`20.74`$ | $`12.42`$ |
| 0.6 | $`0.08`$ | $`5.75`$ | $`5.67`$ | $`16.56`$ | $`11.44`$ | $`4.41`$ |
| 0.8 | $`3.76`$ | $`2.24`$ | 6.00 | $`2.40`$ | $`2.07`$ | 7.42 |
| 1.0 | $`8.91`$ | $`12.83`$ | $`21.74`$ | 16.28 | 19.67 | 22.39 |

In Tab. I it is numerically shown to what extent the HVH theorem is fulfilled by the EBHF approximation. The pressure has been calculated using the relation $`P(\rho ,\beta )=\rho ^2[E_A(\rho ,\beta )/\rho ]`$, where $`E_A(\rho ,\beta )E(\rho ,\beta )/A`$ being the energy per nucleon. In the forth column the left-hand side of Eq. (55) is calculated for several asymmetries (density fixed at $`\rho =\mathrm{\hspace{0.17em}0.17}fm^3)`$. One would notice that, despite the fact that the total density is fixed at the empirical saturation value, our calculated saturation point lies at higher density, because, as is well known, Brueckner theory with two-body force misses the empirical saturation point. The last three columns provide different approximations for the right-hand side of Eq. (55). The pure BHF approximation by itself is far from fulfilling the HVH theorem. Including the unrenormalized ground-state correlations (indicated by BHF$`+M_2`$ in the table), where the Fermi energy is calculated according to Eq. (9), provides some improvement but it is not enough to fulfill the HVH theorem. One needs to include both the rearrangement and the renormalized contributions (EBHF) if a satisfactory agreement within less than $`10\%`$ is to be attained (last column of Tab. I). This result is in keeping with the uncertainty in the calculation of the pressure because the binding energy curve is rather flat as a function of density.

## 4 Applications

### 4.1 Effective mass

The effective mass incorporates the non-local part of the mean field which makes the local part less attractive for a nucleon travelling with momentum $`k>0`$. It is defined as

$`{\displaystyle \frac{m_\tau ^{}(k)}{m}}={\displaystyle \frac{k}{m}}\left({\displaystyle \frac{dE^\tau (k)}{dk}}\right)^1.`$ (56)
The momentum dependence of $`m^{}`$ is characterized by the wide bump inside the Fermi sphere due to the high probability amplitude for particle-hole excitations near the Fermi surface . The effect of correlations is a flattening of the slope of the mean field around the Fermi energy, which implies an enhancement of the effective mass at $`k_F`$ with respect to BHF value . This result is shown in Fig. 9, where an increase from 0.8 to 0.92 is observed for symmetric nuclear matter at the empirical saturation density. Also shown in the figure is the isospin dependence of the neutron (upper curve) and proton (lower curve) effective masses. In both the BHF and EBHF calculations, $`m_n^{}`$ increases and $`m_p^{}`$ decreases as increasing $`\beta `$. Compared to the BHF approximation, the corrections of EBHF shift $`m_n^{}`$ and $`m_p^{}`$ to higher values, a feature which can be traced to the depletions of the proton and neutron Fermi surfaces due to the ground-state correlations. The value of $`m_p^{}`$ calculated from EBHF approaches its BHF value as increasing $`\beta `$ since the correlations become smaller.

### 4.2 Mean free path

Information on the in-medium cross section or, equivalently, on the mean free path of a nucleon travelling inside a nuclear medium can be obtained from the transparency of a nucleus measured in $`(e,e^{}p)`$ reactions  and, in general, from nucleon-induced reactions at low energy . The underlying assumption is that the behavior of a nucleon located at the position $`\stackrel{}{r}`$ in a nucleus is the same as a nucleon in nuclear matter at density $`\rho (\stackrel{}{r})`$. Such an assumption is the well known local-density approximation (LDA) . The mean free path is intimately related to the imaginary part of the optical potential or, equivalently, to the imaginary part of the mean field. The latter comes from the collisions of a single nucleon with the background of neutrons and protons: a nucleon with momentum $`kk_F`$ can collide with a neutron or proton of its Fermi sea and promote it to a particle state, or a nucleon with momentum $`kk_F`$ interacting with an excited neutron or proton can make it decay into a hole state. The first process is related to the imaginary part of $`M_1`$, the second one to the imaginary part of $`M_2`$, both of which have been plotted in the lower panels of Figs. 6 and 7. But, in the case of asymmetric matter the collisions between like and unlike nucleons yield contributions to the mean field which are very different.

The mean free path $`\lambda _\tau `$ is given by

$`\lambda _\tau (E)={\displaystyle \frac{\mathrm{}^2k(E)}{2\stackrel{~}{m}_\tau }}{\displaystyle \frac{1}{|\mathrm{Im}M_\tau (k(E),E)|}},`$
where $`\stackrel{~}{m}_\tau `$ is the so-called $`k`$-mass, and $`E`$ is the single-particle energy . In Fig. 10 the proton (upper panels) and neutron (lower panels) mean free paths calculated within the EBHF approximation are shown for three values of the total density. In each panel the values of $`\lambda _\tau `$ for several asymmetries are plotted as a function of single-particle energy. The most relevant effect of the isospin asymmetry is the increasing deviation from the symmetric values (solid lines), upward for $`\lambda _n`$ and downward $`\lambda _p`$, as increasing asymmetry. The nonvanishing values of neutron and proton mean free paths below their respective Fermi energies are effects of ground-state correlations which prevents a full occupancy of the Fermi spheres. Comparing with the BHF calculation it turns out that the correlation effects tend to rise the asymptotic value of the mean free path from about $`3`$fm up to about $`4`$fm at the saturation density of symmetric nuclear matter .

In Fig. 11 it is shown how the isospin dependence of the inverse $`\lambda `$ develops as increasing nuclear matter densities at a fixed value of the energy. Except for very small asymmetries the shift of $`\lambda _p`$ and $`\lambda _n`$ is not symmetric with respect to their common value at $`\beta =0`$. At any density the slope of the neutron inverse $`\lambda `$ is less than the proton one. This effect can be traced to the reduction, as increasing neutron excess, of the proton particle-hole excitations contributing to the neutron optical potential (see Figs. 6 and 7). Moreover, the EBHF $`\lambda _n`$ seems to reach the asymptotic value of pure neutron matter much faster than in the uncorrelated case.

The most striking effect of the isospin-asymmetry is the sizeable reduction of the proton mean free path at high asymmetry. Accordingly, the nuclear surface would become more transparent to neutrons than protons in nucleon-induced reactions on nuclei near the neutron drip-line. This effect would be more pronounced at higher density as shown in Fig. 11.

### 4.3 Proton fraction in $`\beta `$-equilibrium matter

The core of a neutron star is expected to be formed by an uncharged mixture of neutrons, protons, electrons and muons in equilibrium with respect to the weak interactions ($`\beta `$-stable matter). The concentrations of different particles are then obtained under the requirements

$`\mu _n\mu _p=\mu _e,\mu _\mu =\mu _e.`$ (57)
$`\rho _p=\rho _e+\rho _\mu .`$ (58)
The difference between the neutron and proton chemical potentials can be expressed as

$`\mu _n\mu _p={\displaystyle \frac{B}{Y^p}}|_\rho =2{\displaystyle \frac{B}{\beta }}|_\rho .`$ (59)
In the parabolic approximation, Eq. (51), for the energy per particle of asymmetric nuclear matter, one has

$`\mu _n\mu _p=4E_{sym}(\rho )(12Y^p).`$ (60)
Therefore, the composition of $`\beta `$-stable matter, and in particular, the proton fraction $`Y^p`$ present at a given density, is strongly dependent on the nuclear symmetry energy. The proton fraction plays also a crucial role in the thermal evolution of neutron stars. In fact, if the proton fraction in the core of a neutron star, is above a critical value $`Y_{Urca}^p`$, the so-called direct Urca processes can occur . If they occur, the direct Urca processes enhance the neutrino emission and neutron star cooling rate by a large factor compared to the standard cooling scenario. The critical proton fraction has been estimated  to be in the range 11 – 15%. In a recent paper , based on microscopic EOS of dense matter, it has been found that the onset of the direct Urca processes occurs at densities $`\rho >0.54`$–0.65 $`fm^3`$, depending on the nuclear interaction used to get the EOS ( see Ref.  for more details ).

In Tab. II, we report our present calculations of the proton fraction $`Y^p(\rho )`$ for $`\beta `$-stable matter in comparison with the one obtained with the separable Paris ( see Ref.  ) and variational calculation of Ref.  using the Argonne $`V_{14}`$ potential plus the Urbana model ( UVII ) three-body force. In the calculations reported in Tab. II, muons have been not included.
Tab. II Proton fraction in $`\beta `$-stable nuclear matter (no muons) versus the total baryonic density from different forces. The values reported are $`10^2Y`$. The results in the second column are taken from Ref. . Those in the third column have been given by A. Fabrocini (private communication).

| $`\rho _B`$ | Paris | AV14+UVII | present |
| --- | --- | --- | --- |
| 0.038 | 2.75 |  |  |
| 0.076 | 2.80 | 1.85 | 2.40 |
| 0.11 | 3.09 | 2.48 | 2.74 |
| 0.14 | 3.48 | 2.96 | 3.03 |
| 0.17 | 3.70 | 3.37 | 3.32 |
| 0.20 | 4.10 | 3.74 | 3.50 |
| 0.30 | 4.90 | 3.67 | 4.07 |
| 0.40 | 5.79 | 3.56 | 4.57 |
| 0.50 |  | 3.63 | 5.01 |

Our purpose, in the present paper, is not an accurate determination of the proton fraction in dense stellar matter. Here, we aim to study how the inclusion of contributions beyond the BHF to the chemical potentials could alter the proton fraction in $`\beta `$-stable matter. In fact, to solve the $`\beta `$-equilibrium conditions (57) and (58), the shift between neutron and proton chemical potentials $`\widehat{\mu }\mu _n\mu _p`$ has to be evaluated. In Tab. III the neutron and proton chemical potentials and their difference $`\widehat{\mu }`$, are reported for the different approximations used in the present work. From the results reported in Tab. III we see that the chemical potential, approximated by the Fermi energy, in the EBHF are noticeably affected by the rearrangement and renormalization contributions. However, their difference and consequently the proton fraction is almost unchanged with respect to the BHF approximation. The EBHF approximation provides neutron and proton Fermi energies, which are in better agreement with the empirical values extracted from the mass table of atomic nuclei  than the BHF approximation .
Tab. III Proton and neutron chemical potentials (Fermi energies) calculated in different approximations and compared with the symmetry energy. The results corresponding to four values of asymmetry parameter $`\beta `$ for each of the three densities are reported.

| $`\rho `$ | $`\beta `$ | $`\mu ^p`$ | $`\mu ^n`$ | $`\widehat{\mu }`$ | $`\mu ^p`$ | $`\mu ^n`$ | $`\widehat{\mu }`$ | $`\mu ^p`$ | $`\mu ^n`$ | $`\widehat{\mu }`$ | $`4\beta E_{sym}`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`fm^3`$ |  | BHF | | | BHF+$`M_2`$ | | | EBHF | | |  |
|  | 0.2 | $`33.99`$ | $`18.05`$ | 15.94 | $`30.44`$ | $`14.18`$ | 16.26 | $`23.00`$ | $`6.92`$ | 16.08 | 16.22 |
| 0.085 | 0.4 | $`42.98`$ | $`10.77`$ | 32.20 | $`39.75`$ | $`7.01`$ | 32.74 | $`32.60`$ | $`0.16`$ | 32.44 | 32.45 |
|  | 0.6 | $`51.71`$ | $`3.31`$ | 48.40 | $`49.26`$ | $`0.15`$ | 49.41 | $`43.24`$ | $`5.67`$ | 48.91 | 48.67 |
|  | 0.8 | $`61.32`$ | $`3.68`$ | 65.00 | $`59.63`$ | $`6.58`$ | 66.21 | $`54.85`$ | $`10.49`$ | 65.34 | 64.90 |
|  | 0.2 | $`45.94`$ | $`23.19`$ | $`22.75`$ | $`40.45`$ | $`17.08`$ | 23.37 | $`31.01`$ | $`8.24`$ | 22.77 | 23.00 |
| 0.170 | 0.4 | $`58.08`$ | $`12.65`$ | $`45.43`$ | $`53.46`$ | $`6.72`$ | 46.74 | $`44.30`$ | $`1.25`$ | 45.55 | 46.00 |
|  | 0.6 | $`71.73`$ | $`2.77`$ | $`68.96`$ | $`68.10`$ | $`2.72`$ | 70.82 | $`59.60`$ | $`9.39`$ | 68.99 | 69.00 |
|  | 0.8 | $`86.22`$ | $`6.91`$ | $`93.13`$ | $`83.99`$ | $`11.63`$ | 95.62 | $`76.42`$ | $`16.73`$ | 93.15 | 92.00 |
|  | 0.2 | $`47.53`$ | $`16.11`$ | 31.42 | $`38.89`$ | $`5.75`$ | 33.14 | $`24.04`$ | $`7.12`$ | 31.16 | 32.30 |
| 0.340 | 0.4 | $`64.84`$ | $`1.79`$ | 63.05 | $`57.48`$ | $`8.88`$ | 66.36 | $`42.31`$ | $`20.22`$ | 62.53 | 64.59 |
|  | 0.6 | $`82.75`$ | $`12.21`$ | 94.96 | $`77.17`$ | $`22.85`$ | 900.02 | $`62.44`$ | $`32.35`$ | 94.79 | 96.89 |
|  | 0.8 | $`103.15`$ | $`25.48`$ | 128.63 | $`99.87`$ | $`35.52`$ | 135.39 | $`86.32`$ | $`42.95`$ | 129.27 | 129.18 |

## 5 Conclusions

In this paper we have reported the study of asymmetric nuclear matter within the Brueckner-Bethe-Goldstone approach. The isospin effect on the equation of state has been investigated by performing a set of calculations at the two hole-line level of the BBG expansion for the energy per particle $`B(\rho ,\beta )`$. The Bethe-Goldstone equation has been solved with the Argonne $`V_{14}`$ interaction. The continuous choice has been adopted for the auxiliary potential since it makes the convergence of the hole-line expansion faster than the gap choice . Ranging the asymmetry parameter from $`\beta =0`$ (symmetric nuclear matter) to $`\beta =1`$ (pure neutron matter) it was possible to check that $`B(\rho ,\beta )`$ exhibits a linear dependence on $`\beta ^2`$ for baryonic densities as large as at least two times the saturation density. This result confirms the empirical law introduced in the mass formula of atomic nuclei and also extends its validity up to the highest asymmetries. As a consequence, the entire isospin effect is incorporated in the symmetry energy. The calculation of the symmetry energy in the BHF approximation shows a monotonic increase as a function of baryonic density. Its value calculated at the saturation density is about 28.7MeV, in agreement with the empirical one. The comparison with the variational prediction is made rather difficult due to the contradictory results still existing in this approach. An accurate determination of the symmetry energy is required for dynamical simulations of collisions between neutron-rich nuclei, where the collective observables including collective flows, balance energy and other quantities are expected to be sensitive to the isospin degree of freedom . Ground state correlations were included in the mass operator up to the four hole-line order contributions. Their effect on the single particle properties has been investigated. The first-order contribution to the mass operator displays a linear dependence on the asymmetry parameter confirming a long-standing analysis by Lane . A new effect of the isospin degree of freedom appears when the ground-state correlations induced by the second-order contribution are introduced in the mass operator. That is a nonlinear effect due to the particle-hole excitations of, say, protons induced by the propagation of a neutron in the nuclear medium. This new feature affects the isospin dependence of single-particle properties such as mean field, effective mass and mean free path. Along with the symmetry energy the heavy ion collisions with asymmetric nuclei could also probe the isospin dependence of mean free path and effective mass, which play also an important role in the collision dynamics.

The EBHF approximation for asymmetric matter results in a satisfactory fulfillment of the Hugenholtz-Van Hove theorem in all asymmetry range $`0\beta 1`$. This property makes us more confident of the hole-line expansion of the mass operator for calculating the single-particle properties including the Fermi energy. We found that the neutron and proton chemical potentials are largely affected by contributions beyond the BHF approximation. This could have far reaching consequences for the physics of the neutron star crust. In fact, the proton chemical potential in asymmetric nuclear matter is a very important ingredient in locating the inner boundary of the neutron star crust. However the difference $`\widehat{\mu }\mu _n\mu _p`$, and consequently the proton fraction in $`\beta `$-stable matter, is almost unchanged in the EBHF approximation with respect to the BHF approximation.

## Acknowledgment

The authors are indebted to Professor A. Fabrocini for valuable discussions.


text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The Navier-Stokes equations are notoriously difficult to solve. However, from the viewpoint of differential topology, the Navier-Stokes equations may be viewed as a statement of cohomology: the difference between two non-exact 1-forms is exact. Abstractly, the idea is similar to the cohomology statement of the first law of thermodynamics.

$$QW=dU$$
(1)
For the Navier-Stokes case, define the two inexact 1-forms in terms of the dissipative forces

$$W_D=𝐟_Dd𝐫=\rho \{\nu ^2𝐕\}d𝐫$$
(2)
and in terms of the advective forces of virtual work

$$W_V=𝐟_Vd𝐫=\rho \{𝐕/t+grad(𝐕𝐕/2)𝐕\times curl𝐕\}d𝐫$$
(3)
Then the abstract statement of cohomology, formulated as $`W_VW_D=dP`$ , when divided by the common function, $`\rho ,`$ is precisely equivalent to an exterior differential system whose coefficients are the partial differential equations defined as the Navier-Stokes equations,

$$\{𝐕/t+grad(𝐕𝐕/2)𝐕\times curl𝐕\}\{\nu ^2𝐕\}=gradP/\rho $$
(4)
The cohomological constraint on the velocity field, $`𝐕`$, is such that the kinematically defined vector, $`𝐟`$,

$$𝐟=𝐟_V𝐟_D$$
(5)
is a vector field that satisfies the Frobenius integrability theorem . That is,

$$𝐟curl𝐟=0\text{ even though }𝐯curl𝐯0.$$
(6)
The meaning of the Frobenius criteria is that the vector $`𝐟`$ has a representation in terms of only two independent functions of $`\{x,y,z,t\}`$. The Navier-Stokes equations makes this statement obvious. One of these functions has a gradient, $`gradP,`$ in the direction of the tangent vector to $`𝐟,`$ and the other function, $`\rho ,`$ is a renormalization, or better, a reparametrization factor for the dynamical system represented by $`𝐟`$.

These observations suggest that there must exist certain constraint relationships on the functional forms that make up the components of any solution vector field, $`𝐕`$, (which usually does not satisfy the Frobenius condition in general) such that the covariant kinematic vector, $`𝐟`$, is decomposable in terms of at most two functions. If such a constraint equation can be found in terms of the component functions that represent $`𝐕`$, then its solutions may be easier to deduce than the direct solutions of the Navier-Stokes equations. For example, the constraint relation may involve only 1 partial differential equation rather than 3. In fact such a single constraint relation can be found by imposing a type of symmetry condition on the system, a symmetry condition that expresses the existence of a two dimensional (functional) representation for the vector field,$`𝐟.`$ In this article attention will be focused on the two spatial variables, $`r`$ and $`z`$, such that the solution examples will have a certain degree of cylindrical symmetry. As these solutions involve dissipative terms with a kinematic viscosity coefficient, $`\nu ,`$ they are not necessarily equilibrium solutions of an isolated thermodynamic system.

Closed form solutions are few in number , but it appears that many of the known steady-state solutions to the Navier-Stokes equations fall into the following class of systems: Consider a variety$`\{x,y,z,t\}`$ with $`r^2=x^2+y^2`$. Consider three arbitrary functions, $`\mathrm{\Theta }(r,z)`$ and $`\mathrm{\Phi }(r,z,t)`$, and $`\mathrm{\Lambda }(r,z)`$ which are defined in terms of two independent variables spatial variables, $`(r,z),`$ and time. Define the flow field, $`𝐕`$ in cylindrical coordinates as,

$$𝐕=\mathrm{\Lambda }(r,z)𝐮_z+\mathrm{\Theta }(r,z)𝐮_r+\mathrm{\Phi }(r,z,t)𝐮_\varphi /r,$$
(7)
where $`𝐮_\varphi `$ is a unit vector in the azimuthal direction. Note that this vector field does not necessarily satisfy the Frobenius theorem. Note that for simplicity, the only time dependence permitted is in the azimuthal direction.

Substitution of this format for $`𝐕`$ into the equation for $`𝐟`$ will yield a vector equation of the form

$$𝐟=\alpha (r,z)𝐮_z+\beta (r,z)𝐮_r+\gamma (r,z,t)𝐮_\varphi /r.$$
(8)
The Pfaffian form $`W=𝐟d𝐫`$ will become an expression in two variables if the azimuthal factor $`\gamma (r,z,t)`$ is constrained to the value zero. In other words, a single constraint on the functions, $`\mathrm{\Theta }(r,z)`$ and $`\mathrm{\Phi }(r,z,t)`$, and $`\mathrm{\Lambda }(r,z)`$, defined by the equation $`\gamma (r,z,t)=0,`$can be used to reduce the Pfaffian form to the expression

$$\alpha (r,z)dz+\beta (r,z)dr=dP/\rho (r,z)$$
(9)
The left hand side represents a Pfaffian form in two variables, and therefore always admits an integrating factor. It is this idea that is used to find new solutions to the Navier-Stokes equations. First a solution to constraint equation is determined. Then the Cartan 1-form of total work is computed. The 1-form is either exact, or can be made exact by an appropriate integrating factor. If the 1-form is exact then the Pressure is obtained by integration. It the 1-form is not exact a suitable integrating factor is found, and that integrating factor represents a variable fluid density, $`\rho .`$ For a given choice of integrating factor, the Pressure is again obtained by integration.

It is also useful to consider a rotating frame of reference defined by the equation

$$𝛀=\omega 𝐮_z.$$
(10)
It is the choice of rotational axis that defines the cylindrical symmetry. For such rotating systems the same technique will insure that the flow field, $`𝐕`$, is a solution of the Navier-Stokes equations in a rotating frame of reference,

$`𝐕/t+grad(𝐕𝐕/2)𝐕\times curl𝐕`$ (11)
$`=`$ $`gradP/\rho +\nu ^2𝐕2𝛀\times 𝐕𝛀\times (𝛀\times 𝐫)`$ (12)
By direct substitution, into the Navier-Stokes equation above, of the presumed format for the velocity field $`𝐕`$ yields an expression for $`\gamma (r,z)`$ in terms of the three functions $`\mathrm{\Theta }(r,z)`$ and $`\mathrm{\Phi }(r,z)`$, and $`\mathrm{\Lambda }(r,z):`$

$`\gamma (r,z)`$ $`=`$ $`\{\mathrm{\Phi }/t+\mathrm{\Lambda }(r,z)\mathrm{\Phi }/z+\mathrm{\Theta }(r,z)(\mathrm{\Phi }/r2\omega r)`$ (14)
$`\nu \{^2\mathrm{\Phi }/z^2+^2\mathrm{\Phi }/r^2(\mathrm{\Phi }/r)/r\}.`$
Similar evaluations of the standard formulas of vector calculus in terms of the assumed functional forms for the velocity field lead to the useful expressions:

$$div𝐕=\mathrm{\Theta }/r+\mathrm{\Theta }/r+\mathrm{\Lambda }/z$$
(15)
$$curl𝐕=\{\mathrm{\Phi }/r𝐮_z\mathrm{\Phi }/z𝐮_r\}/r+\{\mathrm{\Theta }/z\mathrm{\Lambda }/r\}𝐮_\varphi $$

$`curlcurl𝐕`$ $`=`$ $`\{^2\mathrm{\Lambda }/r^2+^2\mathrm{\Theta }/zr\}𝐮_z+`$
$`\{^2\mathrm{\Theta }/z^2+^2\mathrm{\Lambda }/zr\}𝐮_r+`$
$`\{^2\mathrm{\Phi }/z^2^2\mathrm{\Phi }/r^2+(\mathrm{\Phi }/r/r)\}𝐮_\varphi /r`$
$`𝐕\times curl𝐕`$ $`=`$ $`\mathrm{\Theta }\{\mathrm{\Theta }/z\mathrm{\Lambda }/r\}𝐮_z+`$
$`\mathrm{\Lambda }\{\mathrm{\Lambda }/r\mathrm{\Theta }/z\}𝐮_r+`$
$`\{(1/r^2)grad(\mathrm{\Phi }^2/2)\{\mathrm{\Lambda }\mathrm{\Phi }/z+\mathrm{\Theta }\mathrm{\Phi }/r\}𝐮_\varphi /r`$
$`grad(𝐕𝐕)/2`$ $`=`$ $`\{\mathrm{\Theta }\mathrm{\Theta }/z+\mathrm{\Lambda }\mathrm{\Lambda }/z+\mathrm{\Phi }(\mathrm{\Phi }/z)/r^2\}𝐮_z+`$
$`\{\mathrm{\Theta }\mathrm{\Theta }/r+\mathrm{\Lambda }\mathrm{\Lambda }/r+\mathrm{\Phi }(\mathrm{\Phi }/r)/r^2\mathrm{\Phi }^2/r^3\}𝐮_r`$
$`grad(div𝐕)`$ $`=`$ $`\{^2\mathrm{\Theta }/zr+(\mathrm{\Theta }/z)/r+^2\mathrm{\Lambda }/z^2\}𝐮_z+`$
$`\{\mathrm{\Theta }^2/r^2+^2\mathrm{\Lambda }/rz+(\mathrm{\Theta }/r)/r\mathrm{\Theta }/r^2\}𝐮_r`$
It is remarkable that many solutions to the Navier-Stokes equations then can be found by using the following algorithm: Choose a functional form $`\mathrm{\Phi }(r,z)`$ of interest and then deduce functions $`\mathrm{\Lambda }(r,z)`$ and $`\mathrm{\Theta }(r,z)`$ to satisfy the azimuthal constraint,

$$\gamma (r,z,t)=0.$$
(20)
The flow field $`𝐕`$ so obtained is therefore a candidate solution to the compressible, viscous, three dimensional Navier-Stokes equations for a system with a density distribution, $`\rho `$ and a pressure, $`P`$. The components of flow field so determined then permit the evaluation of the coefficients of the Pfaffian form

$$W=\alpha (r,z)dz+\beta (r,z)dr$$
(21)
If the expression is not a perfect differential, then use the standard methods of ordinary differential equations to find an integrating factor, $`\rho (r,z).`$ The integrating factor represents the density distribution of the resulting Navier-Stokes solution. The Pressure follows by integration.

This method is demonstrated in the next section for the known viscous vortex examples reported in Lugt. In addition, several new closed form exact solutions are generated by the technique. Among these closed form solutions are exact solutions to the Navier Stokes equations (in a rotating frame of reference) that exhibit the bifurcation classifications for N = 3 as given by Langford . In particular, exact, non-truncated solutions are given that represent the trans-critical Hopf bifurcation, the saddle-node Hopf bifurcation, and the hysteresis Hopf bifurcation. It has been long suspect that many phenomena in hydrodynamics exhibit Hopf bifurcation; now these exact solutions to the Navier-Stokes equation formally justify this position, and are especially interesting for the understanding of slightly perturbed Poiseuille flow and the onset of turbulence in a pipe.

## 2 EXAMPLES

In the following examples, the vector field specified has been used to compute the various terms in the Navier-Stokes equations. The algebra has been simplified by use of a symbolic computation program written in the Maple syntax. For each example, the two vector components that make up the work one form have been evaluated and are displayed with the solution. For the divergence free cases, the pressure function also has been computed. First, known solutions are exhibited, and are shown to be derived from the above technique. Then a few new solutions are exhibited.

### 2.1 Old solutions

#### 2.1.1 Example 1. The Rankine Vortex

$$\mathrm{\Phi }(r,z)=a+(b+\omega )r^2,\mathrm{\Theta }=0,\mathrm{\Lambda }=1$$
(22)
$$𝐟_V=\{(a+br^2)^2/r^3\}𝐮_r+\{0\}𝐮_\varphi /r+\{0\}𝐮_z$$
(23)
$$𝐟_D=\{0\}𝐮_r+\{0\}𝐮_\varphi /r+\{0\}𝐮_z$$
(24)
This flow is a solution independent of the kinematic viscosity coefficient (the velocity field is harmonic, as $`𝐟_D`$ $`=0`$) and therefore could be construed as an equilibrium solution. This solution, for a and b equal to piecewise constants, will generate the Rankine vortex.

As the flow is isochoric ($`divV=0`$), the steady pressure can be determined by quadrature, and is given by the expression,

$$P=1/2(b^2r^4+4abr^2ln(r)a^2)/r^2$$
(25)
#### 2.1.2 Example 2. Diffusion Balancing Advection.

$$\mathrm{\Phi }(r,z)=a+br^{2+m/\nu },\mathrm{\Theta }(r,z)=m/r,\mathrm{\Lambda }=1,\omega =0$$
(26)
$`𝐟_V`$ $`=`$ $`\{m^2(a+br^{(2\nu +m)/\nu })^2/r^3\}𝐮_r+`$
$`\{br^{(2\nu +m)/\nu }m(2\nu +m)/\nu r^2\}𝐮_\varphi /r+\{0\}𝐮_z`$
$$𝐟_D=\{0\}𝐮_r+\{br^{(2\nu +m)/\nu }m(2\nu +m)/\nu r^2\}𝐮_\varphi /r+\{0\}𝐮_z$$
(28)
In this case the Laplacian of the vector field is not zero, but the dissipative parts exactly cancel the advective parts in the coefficient of the azimuthal field, thereby satisfying the constraint condition. As the functions depend only on r, the integrability (gradient) condition is satisfied, and these solutions obey the Navier-Stokes equations for a system of constant density. The Pressure function may be computed as

$`P`$ $`=`$ $`(\nu b(4a(m+\nu )+bmr^{(2\nu +m)/\nu )})r^{(2\nu +m)/\nu )}`$
$`(m(\nu +m)(a^2+m^2))/(2m(\nu +m)r^2)`$
The solutions are cataloged in Lugt. As these solutions explicitly involve the kinematic viscosity, $`\nu ,`$ they cannot be equilibrium solutions to isolated systems. Instead they represent steady state solutions, far from equilibrium. A special case exists for $`m/\nu =2`$.

#### 2.1.3 Example 3. Burger’s Solution, but with Helicity and Zero Divergence.

$$\mathrm{\Phi }(r,z)=k(1e^{ar^2/2\nu }),\mathrm{\Theta }(r,z)=ar,\mathrm{\Lambda }=U+2az,\omega =0$$
(30)
$`𝐟_V`$ $`=`$ $`\{(ke^{(ar^2/2\nu )}+r^2ak)(ke^{(ar^2/2\nu )}r^2ak)/r^3\}𝐮_r+`$
$`\{kra^2/\nu e^{(1/2ar^2/\nu )}\}𝐮_\varphi /r+`$
$`\{2(U+2az)a\}𝐮_z`$
$$𝐟_D=\{0\}𝐮_r+\{kra^2/\nu e^{(1/2ar^2/\nu )}\}𝐮_\varphi /r+\{0\}𝐮_z$$
(31)
This solution corresponds to a modification of Burger’s solution and exhibits a 3-dimensional flow (in 2-variables) in which the diffusion is balanced by convection to give azimuthal cancellation. The Burgers solutions has been modified to exhibit zero divergence. This flow in a non-rotating frame of reference exhibits a helicity.

$$Helicity=(U+2az)(ka/\nu )e^{(1/2ar^2/\nu )}$$
(32)
### 2.2 New Solutions

#### 2.2.1 Example 4. A Beltrami Type Solution

$$\mathrm{\Phi }(r,z)=r^2\mathrm{cos}(z/a),\mathrm{\Theta }(r,z)=r\mathrm{sin}(z/a),\mathrm{\Lambda }(r,z)=2a\mathrm{cos}(z/a),\omega =0$$
(33)
$$𝐟_V=\{r\}𝐮_r+\{0\}𝐮_\varphi /r+\{4a\mathrm{cos}(z/a)sin(z/a)\}𝐮_z$$
(34)
$$𝐟_D=\nu /a^2[\{rsin(z/a)\}𝐮_r+\{rcos(z/a)\}𝐮_\varphi /r+\{2acos(z/a)\}𝐮_z]$$
(35)
This solution is a Beltrami-like solution, has zero divergence, and can be made time harmonic by multiplying the velocity field by any function of t. The flow exhibits Eckman pumping and has a superficial resemblance to a hurricane. The time independent steady flow is a strictly Beltrami $`(curl𝐯=a𝐯)`$with the vorticity proportional to the velocity field. In all cases the helicity is given by the expression,

$$Helicity:=(r^2+4a^2\mathrm{cos}(z/a)^2)/a.$$
(36)
The kinetic energy is a/2 times the helicity, which is a times the enstrophy. The Pressure generated from the Navier Stokes equation is given by the expression

$$P=1/2(r^2+(r^2(\nu /a^2)4\nu )sin(z/a)+4a^2\mathrm{sin}(z/a)^2)$$
(37)
#### 2.2.2 Example 5. A Saddle Node Hopf Solution

$$\mathrm{\Phi }(r,z)=\omega r^2,\mathrm{\Theta }(r,z)=r(a+bz),\mathrm{\Lambda }(r,z)=Udr^2+Bz^2$$
(38)
The components of the advective force and dissipative force are given by the expressions,

$`𝐟_V`$ $`=`$ $`\{r(a+bz)^2+(Udr^2+Bz^2)rb\}𝐮_r+`$
$`\{0\}𝐮_\varphi /r+`$
$`\{2r^2(a+bz)d+2(Udr^2+Bz^2)Bz\}𝐮_z`$
and

$$𝐟_D=\{0\}𝐮_r+\{0\}𝐮_\varphi /r+\nu \{4d+2B]\}𝐮_z$$
(40)
The divergence of the velocity field is given by the expression:

$$div𝐕:=2\{a+(b+B)z\}$$
(41)
The helicity of the flow depends upon the rotation, $`\omega ,`$

$$Helicity:\omega (+r^2b+2U2bz^2)$$
(42)
but remarkably changes for finite values of $`r`$ and $`z`$, depending on mean flow speed, $`U.`$

Note that when $`b=0,B=0,a=0,`$ the solution is equivalent to the standard incompressible Poiseuille solution for flow down a pipe. The vector velocity field is not harmonic, but vector Laplacian of the velocity field is a constant.

Without these constraints, it is remarkable that the ordinary differential equations that represent the components of the velocity field are in one to one correspondence with the saddle node - Hopf bifurcation of Langford. That is, the ODE,s representing the Langford format for the SN-Hopf are given by the expressions:

$`dz/dt`$ $`=`$ $`\mathrm{\Lambda }(r,z)=Udr^2+Bz^2`$ (43)
$`dr/dt`$ $`=`$ $`\mathrm{\Theta }(r,z)=r(a+bz)`$
$`d\theta /dt`$ $`=`$ $`\omega `$
.This first order system which exhibits tertiary bifurcation is associated with an exact solution of the Navier Stokes partial differential system in a rotating frame of reference. In principle, the method also relaxes the constraint on incompressibility, and allows a density distribution, or integrating factor, to be computed for an exact solution to the compressible Navier-Stokes equations which can be put into correspondence with saddle node-Hopf bifurcation process.

This example exhibits isochoric ($`divV=0`$) flow for $`B+b=0,a=0.`$ The steady isochoric pressure is then determined by quadrature, and is given by the expression,

$$P=b(dr^4/2(r^22z^2)Ubz^4)/2\nu (4d2b)z$$
(44)
where the constant U can be interpreted as the mean flow down the pipe. Part of the pressure is due to geometry, and part is due to the kinematic viscosity. Note that the pressure is independent from the viscosity coefficient when the velocity field is harmonic; e.g. when $`(2db)=0.`$ As the vector Laplacian of the velocity field determines the dissipation in the system, intuition would say that the harmonic solution is some form of a limit set for the otherwise viscous flow.

#### 2.2.3  Example 6. A Transcritical Hopf Bifurcation

$$\mathrm{\Phi }(r,z)=\omega r^2,\mathrm{\Theta }(r,z)=r(Aa+cz),\mathrm{\Lambda }(r,z)=br^2+Az+Bz^2$$
(45)
$`𝐟_V`$ $`=`$ $`\{r(Aa+cz)^2+(br^2+Az+Bz^2)rc\}𝐮_r+`$
$`\{0\}𝐮_\varphi /r+`$
$`\{2r^2(Aa+cz)b+(br^2+Az+Bz^2)(A+2Bz)\}𝐮_z`$
$$𝐟_D=\{0\}𝐮_r+\{0\}𝐮_\varphi /r+\nu \{4b+2B\}𝐮_z$$
(47)
This example exhibits isochoric ($`divV=0`$) flow for $`a=3A/2`$ and $`B=c.`$ The steady isochoric pressure is then determined by quadrature, and is given by the expression,

$$P=1/4cbr^41/8A^2r^2+1/2A^2z^2Az^3c+1/2c^2z^4\nu (4b2c)z.$$
(48)
Again it is apparent that the pressure splits into a viscous and a non-viscous component, and when the flow is harmonic ($`2bc=0),`$ the pressure is independent from viscosity, and there is no dissipation in the flow.

The transcritical Hopf bifurcation is represented by the Langford system

$`dz/dt`$ $`=`$ $`\mathrm{\Lambda }(r,z)=br^2+Az+Bz^2`$ (49)
$`dr/dt`$ $`=`$ $`\mathrm{\Theta }(r,z)=r(Aa+cz)`$
$`d\theta /dt`$ $`=`$ $`\omega `$
#### 2.2.4  Example 7. A Hysteritic Hopf Bifurcation

$$\mathrm{\Phi }(r,z)=\omega r^2,\mathrm{\Theta }(r,z)=r(a+bz),\mathrm{\Lambda }(r,z)=Udr^2+Az+Bz^3$$
(50)
$`𝐟_V`$ $`=`$ $`\{r(a+bz)^2+(Udr^2+Az+Az^3)rb\}𝐮_r+`$
$`\{0\}𝐮_\varphi /r+`$
$`\{2r^2(a+bz)d+(Udr^2+Az+Az^3)(A+3Az^2)\}𝐮_z`$
$$𝐟_D=\{0\}𝐮_r+\{0\}𝐮_\varphi /r+\nu \{4d+6Az\}𝐮_z$$
(52)
This system has the remarkable property that the vector Laplacian changes sign at a position $`z=2d/3A`$ down stream. There is no global way of making this solution isochoric, for the divergence is equal to

$$div𝐕=(A+2a)+2bz+3Az^2.$$
(53)
The hysteretic Hopf bifurcation exhibits what has been called intermittency. The Langford system is

$`dz/dt`$ $`=`$ $`\mathrm{\Lambda }(r,z)=Udr^2+Az+Bz^3`$ (54)
$`dr/dt`$ $`=`$ $`\mathrm{\Theta }(r,z)=r(a+bz)`$
$`d\theta /dt`$ $`=`$ $`\omega `$
## 3 Acknowledments

This work was supported in part by the Energy Lab at the University of Houston in 1989, and was discussed at the Permb Conference in 1990

## 4 References

 FLANDERS, H. (1963) ”Differential Forms”. Academic Press, New York.

 LANGFORD, W. F. (1983) in ”Non-linear Dynamics and Turbulence” Edited by G.I. Barrenblatt, et. al. Pitman, London.

 LUGT, H. J. (1983) ”Vortex flow in Nature and Technology” Wiley, New York, p.33.