text:
# Density Fluctuation Mediated Superconductivity

## I Introduction

Soon after the development of the BCS theory of superconductivity it was realized that Cooper instabilities could arise not only from the exchange of phonons, but also from overscreening of the Coulomb interaction in a model Fermi liquid without ion dynamics. Due to the sharpness of the Fermi surface, the generalized spin and charge susceptibilities exhibit Friedel oscillations in space. Kohn and Luttinger argued that oscillations of a similar origin can show up in the effective interaction between quasiparticles. Keeping all diagrams up to second order in a model bare fermion-fermion interaction, they also found an induced attraction which does not rely on the presence of the Friedel oscillations. They demonstrated that it was possible in principle to construct an anisotropic Cooper state that sampled mainly the attractive regions of the effective pairing interaction.

The Kohn-Luttinger model applies where the susceptibilities or response functions are adequately represented by the bare Lindhard function for a homogeneous electron gas for which the oscillations are weak and hence the calculated superconducting transition temperatures $`T_c`$ turn out to be typically well below the experimentally accessible range. As a charge or spin instability is approached, however, we expect the corresponding response function and its real space attractive regions to become enhanced. Provided that it is possible to match the Cooper pair state to the attractive regions, $`T_c`$ may rise to the experimentally accessible range. An extension of the Kohn-Luttinger theory would naturally lead to spin dependent quasiparticle interactions. Current models of magnetically mediated superconductivity focus on such magnetic interactions which are expected to dominate on the border of magnetic long-range order.

It has been shown that this magnetic interaction treated at the mean-field level can produce anomalous normal state properties and superconducting instabilities to anisotropic pairing states. It correctly predicted the symmetry of the Cooper state in the copper oxide superconductors and is consistent with spin-triplet p-wave pairing in superfluid $`{}_{}{}^{3}He`$ \[for a recent review see, e.g., ref. \]. One also gets the correct order of magnitude of the superconducting and superfluid transition temperature $`T_c`$ when the model parameters are inferred from experiments in the normal state of the above systems. There is growing evidence that the magnetic interaction model may be relevant to other materials on the border of magnetism.

In our previous work we focussed on clarifying the general features of the magnetic interaction model. The latter may be relevant to understanding the superconductivity recently discovered, for example, on the border of antiferromagnetism in systems such as cubic $`CeIn_3`$ and its tetragonal counterpart $`CeRhIn_5`$ and on the border of ferromagnetism in $`UGe_2`$, $`URhGe`$ and $`ZrZn_2`$.

In contrast to the conventional phonon-mediated interaction, which is usually taken to be local in space but non-local in time, the magnetic interaction is non-local in both space and time. The non-locality in space leads to anisotropic pairing states whose nature can be acutely sensitive to details of the lattice and electronic structure, and the form of the quasiparticle interaction . For simplicity, in references we considered only a simple cubic or tetragonal crystal structure, a single tight-binding energy band and a magnetic interaction treated at the mean-field level.

One crucial aspect of the magnetic interaction is the vector nature of the spin degree of freedom. At first sight it might appear that the longitudinal and transverse phonons that mediate the usual lattice interaction would be analogous to the longitudinal and transverse spin fluctuations that mediate the magnetic interaction. However, the latter interaction depends on the relative spin orientation of the interacting particles and hence has a different sign and magnitude for the spin-singlet and spin-triplet Cooper states. By contrast, the conventional phonon mediated interaction is oblivious to the spin degree of freedom of the quasiparticles.

One of the consequences is that on the border of ferromagnetism, the magnetic interaction is typically only attractive in the spin-triplet channel. In that case, only the longitudinal spin fluctuations contributing to pairing while all three contribute to the self-interaction that tends to be pair breaking. This disadvantage can be mitigated in systems with strong magnetic anisotropy in that the effect of the transverse spin fluctuations on the self-interaction would be suppressed while the strength of the pairing interaction arising from longitudinal spin fluctuations need not be reduced. By contrast to isotropic magnetic pairing in the spin-triplet channel, in conventional phonon mediated superconductivity and in magnetic pairing in the spin-singlet channel, all modes contribute both to the pairing and to the self-interaction effects.

We have also found that for the model considered in references the robustness of magnetic pairing increases gradually as one goes from a cubic to a more and more anisotropic tetragonal structure under otherwise similar conditions. This is due to the increase with growing anisotropy of the density of states of both the quasiparticles and of the magnetic fluctuations that mediate the quasiparticle interaction. One expects and calculations presented in this paper show that this result carries over to other pairing mechanisms treated at the same level of approximation as in references.

To further our understanding of the conditions favorable to robust pairing it would seem natural to carry out similar types of analyses of superconductivity on the border of other types of instabilities. We consider the possibility of pairing near instabilities signalled by the divergence of a particle density response function. This could include in principle structural instabilities characterized by the softening of phonons in some regions of the Brillouin zone. The induced interaction produced by these soft phonons, in contrast to conventional phonons, is non-local in space. Therefore, one could expect some similarities to the magnetic pairing problem studied in references.

A density response function may also be expected to be strongly enhanced on the border of a charge density wave (CDW) transition, a stripe instability and an $`\alpha \gamma `$ or valence instability. The appropriate density response function may be expected to become large at a wavevector near the Brillouin zone boundary for a CDW, at small but finite wavevectors for stripes and at zero wavevector near the $`\alpha \gamma `$ transition (at which the structure of the unit cell remains the same, but its volume changes).

We note that if the density transition happens to be strongly first order, the appropriate density response function and hence the associated quasiparticle interaction may not be sufficiently enhanced to lead to an observable superconducting phase. This is particularly relevant to the $`\alpha \gamma `$ transition commonly found in heavy fermion systems, which is similar to the liquid-gas transition in that it is first order except at the critical end point. When the latter is at a temperature well above the expected temperature scale for pairing, the enhanced density fluctuations associated with the $`\alpha \gamma `$ transition are unlikely to produce superconductivity. In the temperature region near the critical end point when density fluctuations are strong, superconductivity would be suppressed by thermal fluctuations, while in the low temperature regime the density fluctuations are too weak because of the strong first order character of the density transition. This could explain the absence of superconductivity in $`CeNi`$ where the critical end point is around room temperature, but the existence of superconductivity in $`CeCu_2Si_2`$ and $`CeCu_2Ge_2`$ at high pressures where a corresponding critical end point is believed to exist at low temperatures or may be just suppressed.

## II Model

We consider quasiparticles in a simple tetragonal lattice described by a dispersion relation

$`ϵ_𝐩`$ $`=`$ $`2t(\mathrm{cos}(p_x)+\mathrm{cos}(p_y)+\alpha _t\mathrm{cos}(p_z))`$ (1)
$``$ $`4t^{}(\mathrm{cos}(p_x)\mathrm{cos}(p_y)+\alpha _t\mathrm{cos}(p_x)\mathrm{cos}(p_z)+\alpha _t\mathrm{cos}(p_y)\mathrm{cos}(p_z))`$ (2)
with hopping matrix elements $`t`$ and $`t^{}`$. $`\alpha _t`$ represents the electronic structure anisotropy along the z direction. $`\alpha _t=0`$ corresponds to the 2D square lattice while $`\alpha _t=1`$ corresponds to the 3D cubic lattice. For simplicity, we measure all lengths in units of the respective lattice spacing. In order to reduce the number of independent parameters, we take $`t^{}=0.45t`$ and a band filling factor $`n=1.1`$ as in our earlier work.

The effective interaction between quasiparticles is taken to be the induced density-density interaction and is defined in terms of a coupling constant g and a generalized density susceptibility, which is assumed to have a simple analytical form consistent with the symmetry of the lattice,

$$\chi (𝐪,\omega )=\frac{1}{N_{𝐪_\mathrm{𝟎}}}\underset{𝐪_\mathrm{𝟎}}{}\frac{\chi _0\kappa _0^2}{\kappa ^2+\mathrm{\Delta }(𝐪)i\frac{\omega }{\eta (\widehat{q})}}$$
(3)
where $`\kappa `$ and $`\kappa _0`$ are the correlation wavevectors or inverse correlation lengths in units of the lattice spacing in the basal plane, with and without strong density correlations, respectively. The function $`\mathrm{\Delta }(𝐪)`$, in Eq. (3), is defined as

$$\mathrm{\Delta }(𝐪)=(4+2\alpha _d)2(\mathrm{cos}(q_xq_{0x})+\mathrm{cos}(q_yq_{0y})+\alpha _d\mathrm{cos}(q_zq_{0z}))$$
(4)
where $`\alpha _d`$ parameterises the density anisotropy. $`\alpha _d=0`$ corresponds to quasi-2D density correlations and $`\alpha _d=1`$ corresponds to 3D density correlations. The sum in Eq. (3) is over all the symmetry related vectors $`𝐪_\mathrm{𝟎}`$, with $`N_{𝐪_\mathrm{𝟎}}`$ the number of such vectors. In the following, we only explicitly write one of the vectors. It should be understood that when we say that the incipient wavevector is, for example, $`𝐪_\mathrm{𝟎}=[\pi /4,0]`$, it is implied that the density response function peaks at the four wavevectors $`[\pm \pi /4,0],[0,\pm \pi /4]`$. The parameter $`\eta (\widehat{q})`$ in Eq. (3) is defined as

$`\eta (\widehat{q})`$ $`=`$ $`T_{DF}\widehat{q}^n`$ (5)
$`\widehat{q}^2`$ $`=`$ $`(4+2\alpha _d)2(\mathrm{cos}(q_x)+\mathrm{cos}(q_y)+\alpha _d\mathrm{cos}(q_z))`$ (6)
where $`T_{DF}`$ is a characteristic density fluctuation temperature. In Eq. (5), the exponent $`n=1`$ if the density fluctuations are such that the total density is conserved and $`n=0`$ otherwise. We note that the pole of the density response function Eq. (3) is purely imaginary and therefore the density fluctuations we consider are overdamped. This is believed to apply on the border of CDW, stripe and valence instabilities, but not typically for lattice density fluctuations for which the poles of the density response must have a non-negligible real component. The latter would require the inclusion of an $`\omega ^2`$ term with real coefficient in the denominator of Eq. (3).

In addition to the induced density interaction, we include an on-site Coulomb repulsion $`I`$. In the large $`I`$ limit, the Cooper pair state vanishes when the interacting quasiparticles are on the same site and thus conventional isotropic s-wave pairing is excluded.

We note that in the corresponding problem of magnetic pairing the effective interaction is repulsive when the two interacting quasiparticles are on the same site in the spin-singlet channel. It is, however, attractive in the spin-triplet channel, but this is irrelevant since the required spatial antisymmetry of the pair state means that the two quasiparticles have zero probability of occupying the same site simultaneously.

A complete description of the model, the Eliashberg equations for the superconducting transition temperature and their method of solution can be found in the appendix.

## III Comparison of the Density and Magnetic Pairing Interactions

Our assumed form of the density response function is similar to that of the generalized magnetic susceptibility used in our previous papers. However, there is a crucial difference in that the effective magnetic interaction depends on the relative orientation of the spins of the two quasiparticles through the factor $`\sigma _\mathrm{𝟏}\sigma _\mathrm{𝟐}`$. In the spin-singlet state the expectation value of $`\sigma _\mathrm{𝟏}\sigma _\mathrm{𝟐}`$ gives a factor of -3. When the interaction is oscillatory in real space, this sign change leads to an interchange of attractive and repulsive regions. Since one must choose a pair state in which the quasiparticles mainly sample the attractive region of the interaction, the sign inversion implies a change in the symmetry of the Cooper pair state as illustrated in Figs. 1a and 2a for the cases of incipient ordering wavevectors $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$ and $`[\pi ,0]`$ in a square lattice. For the case of small $`𝐪_\mathrm{𝟎}`$, where the oscillations are essentially irrelevant in our model, the density interaction is attractive in real space for both the spin-singlet and spin-triplet states, but the magnetic interaction is attractive solely for the spin-triplet state for which the expectation value of $`\sigma _\mathrm{𝟏}\sigma _\mathrm{𝟐}`$ is +1.

In our model for the generalized magnetic susceptibility we have assumed that the overall magnetization is conserved and hence $`\eta (\widehat{q})`$ vanishes as $`𝐪0`$. This leads to greater incoherent scattering for a nearly ferromagnetic than antiferromagnetic metal, and hence to a reduction of $`T_c`$ on the border of ferromagnetism. If the fluctuations of the density are quasi-local as in some models of valence fluctuations, then $`\eta (\widehat{q})`$ does not vanish at small $`𝐪`$. This corresponds to the case $`n=0`$ in Eq. (5). If $`𝐪_\mathrm{𝟎}`$ is sufficiently far away from the origin in the Brillouin zone, the precise value of $`n`$ is not expected to affect the calculated $`T_c`$. Since we are not going to consider the limit $`𝐪_\mathrm{𝟎}=0`$, for simplicity we take $`n=1`$ as in the case of the magnetic interaction.

## IV Results

### A Quasi 2D: $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$ and $`[\pi /m,0]`$ with $`m=1,2,4`$

The dimensionless parameters at our disposal are $`g^2\chi _0/t`$, $`T_{DF}/t`$, $`\kappa _0`$ and $`\kappa `$. For comparison with results of our earlier work for the case of the magnetic interaction, we take $`T_{DF}=2t/3`$ and $`\kappa _0^2=12`$. In 2D, this $`T_{DF}`$ corresponds to about 1000 K for a bandwidth of 1 eV, while our choice of $`\kappa _0^2`$ is a representative value. We note that $`\kappa _0^2/\kappa ^2`$ represents the density susceptibility enhancement factor, analogous to the Stoner factor in the case of the magnetic interaction.

The results of our numerical calculations of the mean-field critical temperature $`T_c`$ as a function of $`g^2\chi _0/t`$ and of $`\kappa ^2`$ is shown in Fig. 1 for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$ in which the Cooper pair state has $`d_{xy}`$ symmetry. The nodal lines of this state in real space are illustrated in Fig. 1a, which also depicts the static density interaction seen by one of the quasiparticles given that the other is at the origin. For values of the dimensionless coupling parameter $`g^2\chi _0/t`$ corresponding to the Random Phase Approximation (of order 10), $`T_c`$ is found to drop very rapidly as one goes away from the instability, i.e with increasing $`\kappa ^2`$.

The corresponding plots for the cases $`𝐪_\mathrm{𝟎}=[\pi /m,0]`$, where $`m`$ = 1, 2 and 4 are shown in Fig. 2, 3 and 4. In contrast to the case $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$, the next-nearest-neighbor interaction for $`𝐪_\mathrm{𝟎}=[\pi ,0]`$ is repulsive. This requires nodal lines along the diagonal, and hence the $`d_{x^2y^2}`$ instead of $`d_{xy}`$ symmetry. As shown in Fig. 2a, the nearest-neighbor interaction vanishes for the special case $`𝐪_\mathrm{𝟎}=[\pi ,0]`$ and the leading attraction comes from third-nearest-neighbors. This explains why pairing is not as robust in this case compared with the case $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$.

As seen from Figs 3a and 4a, the strength of the nearest-neighbor attraction increases as $`𝐪_\mathrm{𝟎}`$ gets smaller, which correlates with the increased robustness of $`T_c`$.

As $`𝐪_\mathrm{𝟎}`$ decreases the density interaction can also be attractive for other pairing states. In order to avoid the on-site Coulomb interaction, one could use the $`d_{xy}`$ state since the next-nearest-neighbor interaction is attractive for sufficiently small $`𝐪_\mathrm{𝟎}`$. But since the $`d_{x^2y^2}`$ state picks the nearest-neighbor attraction, which is dominant, it is expected to be the favored state. For small $`𝐪_\mathrm{𝟎}`$, the density interaction is also attractive in the spin-triplet channel for a $`p_x`$ or $`p_y`$ Cooper state. This state picks two out of the four nearest-neighbor attractive sites, instead of all four for the $`d_{x^2y^2}`$ state. However, the $`p_x`$ or $`p_y`$ state also picks the attraction on all four next-nearest-neighbor sites where the $`d_{x^2y^2}`$ state vanishes. It is thus not immediately obvious in that case which of the two pairing states has the highest $`T_c`$. Fig. 5 shows the Eliahsberg superconducting transition temperature one obtains for the spin-triplet $`p_x`$ and spin-singlet $`d_{x^2y^2}`$ states as a function of the correlation wavevector $`\kappa ^2`$ for $`g^2\chi _0/t=10`$. The plot shows that the $`d_{x^2y^2}`$ is the favored case and we have found this to be true for the range of values of $`\kappa ^2`$ and $`g^2\chi _0/t`$ studied in this paper.

### B 3D: $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$, $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$

The results of the numerical calculations in 3D are shown in Figs 6 and 7 for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$ in the $`d_{xy}`$ Cooper state and $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$ in the $`d_{x^2y^2}`$ Cooper state, respectively. The pairing for $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$ is more robust than for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$ since the dominant attraction comes from nearest neighbor in the former case rather than next-nearest-neighbor as in the latter case. For $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$, pairing is less robust in 3D than for the corresponding quasi-2D case shown in Fig. 1 for all coupling constants. In the case $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$, pairing is more robust in quasi-2D (Fig. 4) than in the corresponding 3D case for weak to intermediate coupling. At strong coupling, however, pairing is more robust in 3D, but coupling constants $`g^2\chi _0/t`$ in the 20-60 range are less physically realistic.

### C Crossover from 3D to quasi-2D: Tetragonal lattice with $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$, $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$

The calculated $`T_c`$ as a function of the electronic and density response anisotropy parameters $`\alpha _t`$ and $`\alpha _d`$, respectively, are shown in Figs. 8 and 9 for representative values of the parameters $`\kappa ^2`$ and $`g^2\chi _0/t`$. The results reported in sections A and B correspond to the quasi-2D case $`\alpha _t=\alpha _d=0`$ and to the 3D case $`\alpha _t=\alpha _d=1`$.

For $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$, shown in Fig. 8, we find that $`T_c`$ increases gradually and monotonically as the system becomes more and more anisotropic in the density interaction. We also note that the effect of the electronic anisotropy is much less pronounced. In the case of an incipient ordering wavevector $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$, Fig. 9 shows that $`T_c`$ is maximum for an anisotropy parameter $`\alpha _d`$ between 0 and 1, namely for an anisotropic albeit not quasi-2D in the density interaction. Also note that in the $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$ case, $`T_c`$ depends more strongly on the electronic anisotropy parameter $`\alpha _t`$ than for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$.

## V Discussion

### A Role of Real Space Oscillations in the Quasiparticle Interaction

When the wavevector $`𝐪_\mathrm{𝟎}`$ at which the density response function is a maximum lies near the Brillouin zone boundary the quasiparticle interaction has short-range real-space oscillations. As a consequence, the robustness of the pairing depends sensitively on whether one can construct a Cooper pair state from quasiparticle states near the Fermi surface such that given one quasiparticle is located at the origin, the probability of finding the second one in regions where the interaction is repulsive is minimized. For the case $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$, this forces us to consider a Cooper state with nodes along the principal (x and y) axes (see Fig. 1a).

In the density interaction channel, the dominant attraction comes from the next-nearest-neighbor sites and is typically much weaker than the dominant nearest-neighbor attraction for spin-singlet magnetic pairing for the same wavevector $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$. This explains why for this wavevector $`𝐪_\mathrm{𝟎}`$, pairing is not as robust for the density interaction as for the magnetic interaction under otherwise similar conditions.

One might think that if there were a wavevector $`𝐪_\mathrm{𝟎}`$ such that the interaction is attractive at the nearest-neighbor sites, one could achieve pairing in the density channel to the same degree of robustness as in the spin-singlet magnetic channel for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$. A potential candidate wavevector is $`𝐪_\mathrm{𝟎}=[\pi ,0]`$ since by rotating the wavevector one would rotate the oscillation pattern in real space. However, the oscillations one obtains via Eq. (3) are superpositions of oscillations running along the x and y directions coming from the symmetry related components with wavevectors $`[\pi ,0]`$ and $`[0,\pi ]`$. These oscillations perfectly cancel at the odd sites (see Fig. 2a), and in particular at nearest-neighbor sites. The dominant attraction arises from the third-nearest-neighbors, and thus contrary to naive expectations the case with $`𝐪_\mathrm{𝟎}=[\pi ,0]`$ leads to even weaker pairing than with $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$. Note that for the corresponding spin-singlet magnetic pairing for $`𝐪_\mathrm{𝟎}=[\pi ,0]`$, because of the inversion of the sign of the interaction due to the spin factor $`\sigma _\mathrm{𝟏}\sigma _\mathrm{𝟐}`$, the dominant attraction would now come from the next-nearest-neighbor sites. In this case, pairing would be more robust in the magnetic than in the density channel for $`𝐪_\mathrm{𝟎}=[\pi ,0]`$, but still not as favorable as the magnetic spin singlet channel for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$.

The robustness of density pairing for the simple tetragonal lattice is optimized for $`𝐪_\mathrm{𝟎}`$ close to the Brillouin zone center since in that case the interaction at all neighboring sites is maximally attractive (see Fig. 4a). In order to avoid the on-site Coulomb repulsion, the pairing state which is the solution of the gap equation (Appendix) vanishes at the origin and its symmetry is of the form $`p_x`$, $`p_y`$, $`d_{xy}`$ or $`d_{x^2y^2}`$ . Since the $`d_{x^2y^2}`$ state has maximum amplitude at nearest-neighbor sites, it has the highest $`T_c`$.

By contrast to the case of the magnetic interaction where the most robust pairing was shown in our previous work to arise for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$, in the density channel our results indicate that the optimal case is for $`𝐪_\mathrm{𝟎}`$ near the Brillouin zone center. Since the symmetry of the Cooper state is the same in both cases, this would suggest that still stronger pairing should arise when the system is on the border of both a magnetic instability with $`𝐪_\mathrm{𝟎}`$ near $`[\pi ,\pi ]`$ and a density instability with low $`𝐪_\mathrm{𝟎}`$. In that case, the two pairing mechanisms would reinforce each other rather than compete.

This observation may be very relevant to the superconductivity in f-electron compounds such as $`CeCu_2Si_2`$ and $`CeCu_2Ge_2`$. In these systems the superconductivity extends over a region in pressure containing both an antiferromagnetic and a valence instability. What is special about these materials is that the critical end point of the latter instability lies at unusually low temperatures or is incipient, which means that density fluctuations are expected to be important in the temperature regime where superconductivity is observed. Since the two instabilities do not occur at the same pressure one would expect that near the magnetic instability the pairing would be dominated by the magnetic channel and as the pressure is increased that it would cross over to a regime dominated by the density channel.

When the two instabilities are sufficiently widely separated, one might expect to see two distinct superconducting domes, one centred near the magnetic instability and the other near the density instability. A double domed superconducting temperature-pressure phase diagram has in fact been observed in $`CeCu_2Si_2`$ and $`CeCu_2Ge_2`$ systems and in $`CeNi_2Ge_2`$. Some of these experimental findings have been interpreted in terms of the effects of magnetic and valence fluctuations.

The overall scale of $`T_c`$ is set by the characteristic temperature of magnetic and density fluctuations which tends to be below 100 K in the above f-systems. One way to increase the value of $`T_c`$ is to increase these characteristic temperatures. This could be achieved by looking for analogous d-metal systems with broader electron bands. The antiferromagnetic and stripe fluctuations in the cuprates may be an example where magnetic and density fluctuations with high characteristic temperature scales reinforce to produce high temperature superconductivity.

### B Role of Crystalline Anisotropy

The numerical results show that the robustness of density mediated superconductivity increases gradually and monotonically as one goes from a cubic to a more anisotropic tetragonal structure for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$ and that $`T_c`$ is optimum for an anisotropic albeit not quasi-2D system for $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$. One can partly understand this result by looking at the evolution of the density interaction in real space with increasing anisotropy as illustrated qualitatively in Figs. 10 and 11 for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$ and $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$, respectively. We see that the attraction in the basal plane gets enhanced as one goes from the cubic to a more anisotropic tetragonal lattice. This enhancement is the consequence of the increase of the phase space of soft density fluctuations as one goes from a cubic to a quasi-2D structure. Note that for $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$, the model pairing potential is not continuous at $`\alpha _d=0`$ since the number of peaks of the density response goes from four (at $`[\pm \pi /4,0],[0,\pm \pi /4]`$) in stricly 2D to six (at $`[\pm \pi /4,0,0],[0,\pm \pi /4,0],[0,0,\pm \pi /4]`$) for $`\alpha _d>0`$. Other than that, our model potential varies smoothly with the tetragonal distortion, parameterised by $`\alpha _d`$ in Figs. 10 and 11, and it is clear that this effect grows gradually with increasing separation between the basal planes. In our Eliashberg calculations, mass renormalization effects, which tend to suppress $`T_c`$, also increase as one goes to a more and more anisotropic crystal structure. Our results thus depend on the interplay between the strengths of the pairing interaction and mass renormalization, and the fact that the maximum $`T_c`$ in the case $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$ occurs for anisotropic but not quasi-2D systems reflects the delicate balance between these opposing effects. The above given phase space argument is similar to that used to explain the increased robustness of magnetic pairing with increasing lattice anisotropy and, hence, as anticipated in Ref., carries over to other pairing mechanisms treated at the one-loop mean-field level. Another potential benefit of going to a more anisotropic crystal structure is the narrowing of the electronic band and the associated increase in the electronic density of states. Our results show that in the case of $`𝐪_\mathrm{𝟎}=[\pi ,\pi ,\pi ]`$ and the model parameters considered, this does not play the dominant role. However, for an incipient ordering wavevector $`𝐪_\mathrm{𝟎}=[\pi /4,0,0]`$, the increase in the electronic density of states with increased lattice anisotropy plays a more important role. This effect could also be sensitive to details of the electronic and crystal structure not considered here.

The calculations presented in this paper and in our previous work show that, in the majority of cases considered, the lattice anisotropy increases the robustness of magnetic and density pairing in the mean-field approximation. Superconducting phase fluctuations which are not included in this approximation may be expected to suppress $`T_c`$ in the 2D limit. Therefore, in practice, one would think that the most favorable case for magnetic or density pairing is that of strong, but not extreme, anisotropy where the effect of the superconducing phase fluctuations are typically weak. This is to be contrasted with the effect of order parameter fluctuations on magnetic and density transitions that can be large in metals even in 3D and more so in 2D. In the case of the density transition, even a small lattice anisotropy and the resulting increase in the order parameter fluctuations can lower the critical end point significantly. By weakening the first order transition at low temperatures, this would enhance the density fluctuations that mediate the pairing on the border of the density instability and lead to a superconducting phase.

The importance of crystalline anisotropy in enhancing the superconducting $`T_c`$ on the border of antiferromagnetism has been dramatically demonstrated in going from the simple cubic system $`CeIn_3`$ to the related tetragonal compounds $`CeMIn_5`$ where $`M=Co`$, $`Rh`$ and $`Ir`$, as correctly anticipated by our earlier model calculations of magnetic pairing.

In addition to an antiferromagnetic instability at relatively low pressure, $`CeIn_3`$ is also thought to have a strongly first order $`\alpha \gamma `$ transition at high pressures. Superconductivity is only observed in a narrow range of pressure and temperature around the antiferromagnetic quantum critical point. Because of the wide separation in pressure between the magnetic and density transitions and the strongly first order nature of the latter, one would expect the observed superconductivity to be magnetically mediated. In the tetragonal compounds $`CeMIn_5`$, however, superconductivity is observed over a wide range of pressures. Were an $`\alpha \gamma `$ transition present in these compounds, the critical end point would be expected to be at much lower temperatures than in $`CeIn_3`$ due to the role of anisotropy as discussed above. This would result in stronger density fluctuations in the neighborhood of the $`\alpha \gamma `$ instability. Could this be another example where antiferromagnetic and low $`𝐪_\mathrm{𝟎}`$ density fluctuations both contribute to the attractive pairing interaction in the $`d_{x^2y^2}`$ Cooper state? This would explain the unusually wide extent of the superconducting domes observed in these materials.

It would not be surprising that such a density transition has not been reported because we expect its signature to be weak. Moreover, it is likely to be observable as a well-defined transition only over a very narrow range in pressure in the temperature-pressure phase diagram and would require very careful examination pressure scans at fixed temperatures in order to detect it.

## VI Outlook

One can expect that the total effective interaction between particles in a strongly correlated electron system to be very complex. The interaction will clearly depend on the charge, but also more generally on the spin and current carried by the particles. The border of a density or spin or current instability is characterized by strongly enhanced order-parameter fluctuations and it is therefore plausible that the dominant interaction channel is mediated by the density, spin or current fluctuations, respectively.

In this paper we have shown how the framework developed for systems on the border of magnetism can be translated to describe systems on the border of density instabilities. A striking feature of the model we have considered is that the most robust pairing is obtained in the spin-singlet $`d_{x^2y^2}`$ Cooper state on the border of both the density and spin instabilities. However, crucially the wavevector $`𝐪_\mathrm{𝟎}`$ at which the response function is most enhanced is different in the two cases. Density fluctuations give rise to the highest superconducting $`T_c`$ for $`𝐪_\mathrm{𝟎}`$ near the center of the Brillouin zone while magnetic pairing is strongest for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$. While it is possible to construct a Cooper pair state that samples mainly the most attractive regions of the density and magnetic interaction for $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$, the attraction is weaker in the density channel because the minimum separation of the two interacting particles is larger in the $`d_{xy}`$ state for the density interaction than in the $`d_{x^2y^2}`$ state for the magnetic interaction (Fig. 1a). For low $`𝐪_\mathrm{𝟎}`$, however, the density interaction is mostly attractive provided that the particles are not on the same site (Fig. 4a) and thus the most favored state is $`d_{x^2y^2}`$ in which the two interacting particles can take advantage of the strong nearest-neighbor attraction.

It would seem that the same argument could apply to the low $`𝐪_\mathrm{𝟎}`$ magnetic interaction. However, in contrast to the density interaction which has the same sign in the spin-singlet and spin-triplet channels, the magnetic interaction depends on the relative spin orientation of the two interacting particles and thus has a different sign for the two cases. Magnetic pairing in the spin-singlet state is only possible if the real space interaction has sufficiently short-wavelength oscillations. Therefore, when $`𝐪_\mathrm{𝟎}`$ is near the center of the Brillouin zone magnetic pairing in the spin-singlet state is not possible, but is allowed in principle in the spin-triplet state for which the magnetic interaction has the opposite sign. However, magnetic pairing in this state has the disadvantage that only the exchange of spin fluctuations polarized along the direction of the interacting spins, i.e., the longitudinal fluctuations contribute to the particle interaction. For a spin rotationally invariant system, both the longitudinal and transverse spin fluctuations contribute to pairing only for a spin-singlet Cooper state. For the model considered in Refs this effect results in much weaker pairing on the border of ferromagnetism ($`𝐪_\mathrm{𝟎}=0`$) than antiferromagnetism with $`𝐪_\mathrm{𝟎}=[\pi ,\pi ]`$.

Another disadvantage of being on the border of ferromagnetism is that for otherwise similar conditions the suppression of $`T_c`$ due to the self-interaction arising from the exchange of magnetic fluctuations is stronger than in the corresponding case on the border of antiferromagnetism. This disadvantage can be mitigated in systems with strong magnetic anisotropy in that the effect of the transverse magnetic fluctuations on the self-interaction would be suppressed while the strength of the pairing interaction arising from the longitudinal magnetic fluctuations need not be reduced. This may apply in systems with strong spin-orbit interactions or in the spin-polarized state close to the border of ferromagnetism.

These arguments  have stimulated a new search for evidence of superconductivity on the border of itinerant electron ferromagnetism in cases where spin anisotropy is expected to be pronounced, such as $`UGe_2`$. This search has proved fruitful because it led to the first observation of the coexistence of superconductivity and itinerant electron ferromagnetism in $`UGe_2`$ and shortly thereafter in $`ZrZn_2`$ and $`URhGe`$.

In the previous section and in Figs. 10 and 11 we gave simple arguments to explain how the pairing effect of the interactions are strengthened by a tetragonal distortion in our model. However, the same effects also contribute to an enhanced self-interaction which acts to suppress $`T_c`$. The relative importance of the pair-forming and pair-breaking effects of the effective interaction cannot be inferred solely from the above physical picture for the density channel and the analogous arguments given in Ref. for the magnetic channel. The numerical calculations show that for most cases considered here and in Ref. the pair-forming effects dominate.

A most striking manifestation of the interplay between the pair-forming and pair-breaking tendency of both the density and magnetic interactions is the breakdown of the McMillan-style expression for $`T_c`$. This was noted in Refs. and has been interpreted in Ref. in terms of the important role played by the incoherent part of the Green function which is ignored in the simplest treatments, but is included in the present and earlier work where the full momentum and frequency dependence of the self-energy is taken into account.

In this and our earlier work we deliberately avoided modelling a specific system since our main goal is to gain insights into the nature of the pairing problem on the border of a density and spin instabilities. We have focussed on understanding trends and certain general factors affecting the robustness of the pairing mechanism. Even the simplest models considered display surprising sensitivity to factors such as the nature of the instability, the wavevector $`𝐪_\mathrm{𝟎}`$ at which it occurs, the total spin of the Cooper pair, details of the electronic and lattice structure as well as the form of the relevant response function. Therefore, one should exercise caution in making quantitative comparisons between the results of our calculations and experiment.

In particular, our model may not apply to situations where there is a large, local in space, contribution to the dynamical response function. This would not contribute to the pairing interaction for anisotropic Cooper states, but could greatly enhance the self-interaction effect that is pair-breaking. This could for instance greatly increase the sensitivity of $`T_c`$ to lattice anisotropy as observed in $`CeMIn_5`$ systems and to the correlation length ($`1/\kappa `$) characterising the relevant response function as indirectly seen in the strong pressure dependence of $`T_c`$ in, for example, $`CeIn_3`$. Such a local contribution to the magnetic response function has been observed in heavy fermion systems.

The results of the calculations would be very sensitive to the particular choice of the wavevector dependence of the response function. In cases where it falls off in $`𝐪`$ faster than in our model, the response is appreciably enhanced in a smaller portion of the Brillouin zone and one would then expect the effect of the density or magnetic interaction on the thermal, transport and superconducting properties to be reduced. This could explain the surprisingly weak effects on these properties of the CDW fluctuations in systems such as $`NbSe_2`$.

At first sight, our results seem to imply that anisotropic forms of superconductivity should be a generic property of systems on the border of density and magnetic instabilities. It may seem surprising therefore that there are still relatively few observations of this phenomenon. In addition to the sensitivity of $`T_c`$ to details of the system as discussed above, in many cases the multiplicity of bands and lattice structure may be unfavorable for pairing to such an extent that quenched disorder may completely suppress superconductivity. An illustration of this latter point is the dramatic collapse of the spin-triplet superconducting $`T_c`$ in $`Sr_2RuO_4`$ in the presence of Al impurity concentrations as low as 0.1% . Another factor that may explain the absence of superconductivity is the common occurrence of first order rather than continuous magnetic as well as density instabilities. Our results show that in many cases one has to be close to the instability. A first order transition may make this region of the phase diagram inaccessible.

The results of the mean-field calculations presented here and in our earlier papers show that robust pairing can occur in both density and magnetic channels under suitable conditions. Therefore, it would seem that one should not favor one mechanism over another in the search for new examples of high temperature superconductivity. This conclusion may turn out to be incorrect when corrections to the one-loop mean-field calculations are taken into account. In contrast to the conventional electron-phonon pairing theory where corrections to the Eliashberg value of $`T_c`$ are small, it has been argued for many years that this may not be the case for other types of pairing mechanisms.

It has been shown that the mean-field approximation of the kind we have considered here qualitatively breaks down in a half-filled Hubbard model in 2D which is a Mott insulator with long-range antiferromagnetic order at absolute zero. This breakdown has been interpreted in terms of the effect of thermal magnetic fluctuations in the renormalized classical regime. Thermal density fluctuations near Peirels CDW transition in 2D also lead to qualitative changes to the electronic spectrum that are not captured in the present model.

Recent non-perturbative calculations have shown that dynamical fluctuations even at the Gaussian level are sufficient to cause a breakdown of the present mean-field model for sufficiently small $`\kappa ^2`$. In this and our earlier work on the magnetic interaction model, we had to solve the Eliashberg equations for the superconducting transition temperature $`T_c`$ for very many choices of model parameters. Even with the best numerical algorithms, this is only practically feasible, especially in 3D, if the overall scale of the $`T_c`$’s obtained is rather high (say of the order of 0.1t). Our choice of the charactersitic density-fluctuation temperature $`T_{DF}=2t/3`$ (or spin-flucutation temperature $`T_{SF}=2t/3`$ in our earlier work) was therefore dictated by such considerations. We now know that for these values of $`T_{DF}`$ (or $`T_{SF}`$) and the range of other model parameters considered here and in our earlier papers that vertex corrections are important. Our results are nevertheless useful if interpreted in the following way. We expect, and have checked in a number of cases, that the trends in $`T_c/T_{DF}`$ (or $`T_c/T_{SF}`$) are weakly dependent on the abolute scale of the characteristic temperatures $`T_{DF}`$ and $`T_{SF}`$. Therefore the conclusions drawn from our calculations are expected to remain qualitatively correct for smaller values of $`T_{DF}`$ or $`T_{SF}`$ and hence $`T_c`$, values for which the mean-field theory of superconductivity is likely to be more accurate.

It has generally been believed that the most important factor in determining the applicability of Migdal’s theorem is the form and parameters entering the relevant fluctuation spectrum. Therefore, a surprising finding was that in the range of model parameters studied in Ref., the vertex corrections to the Eliashberg self-energy led to qualitatively different electron spectral functions for a coupling to magnetic and density fluctuations with identical fluctuation spectra. In those calculations, it was found that the corrections to the Eliashberg theory enhanced the magnetic interaction, but suppressed the density interaction. This effect can readily be seen at the two-loop level.

The contribution of the transverse magnetic fluctuations to the first order vertex correction is opposite in sign to that of longitudinal ones and dominates. On the other hand, in the density channel one has essentially only ’longitudinal’ fluctuations, which as in the magnetic case suppress the interaction at this level of approximation. It is also known that the two fluctuation exchange processes lead to the enhancement of the spin-singlet magnetic pairing interaction. While detailed calculations of $`T_c`$ beyond the single-fluctuation exchange approximation have yet to be carried out, the above findings suggest that spin-singlet magnetic pairing may turn out to be more robust than density pairing under otherwise equivalent conditions.

## VII Acknowledgments

We would like to thank A.V. Chubukov, P. Coleman, F.M. Grosche, S.R. Julian, P.B. Littlewood, A.J. Millis, A.P. Mackenzie, D. Pines, D.J. Scalapino, M. Sigrist and H. Yuan, for discussions on this and related topics. We acknowledge the support of the EPSRC, the Newton Trust and the Royal Society.

## VIII Appendix

We consider quasiparticles on a cubic or tetragonal lattice. We assume that the dominant interaction is in the density channel and postulate the following low-energy effective action for the quasiparticles:

$`S_{eff}`$ $`=`$ $`{\displaystyle \underset{𝐩,\alpha }{}}{\displaystyle _0^\beta }𝑑\tau \psi _{𝐩,\alpha }^{}(\tau )\left(_\tau +ϵ_𝐩\mu \right)\psi _{𝐩,\alpha }(\tau )+{\displaystyle \frac{I}{N}}{\displaystyle \underset{𝐪}{}}{\displaystyle _0^\beta }𝑑\tau \rho _{}(𝐪,\tau )\rho _{}(𝐪,\tau )`$ (8)
$`{\displaystyle \frac{g^2}{2N}}{\displaystyle \underset{𝐪}{}}{\displaystyle _0^\beta }𝑑\tau {\displaystyle _0^\beta }𝑑\tau ^{}\chi (𝐪,\tau \tau ^{})\rho (𝐪,\tau )\rho (𝐪,\tau ^{})`$
where $`N`$ is the number of allowed wavevectors in the Brillouin zone and the carrier density $`\rho _\sigma (𝐪,\tau )`$ is given by

$$\rho _\sigma (𝐪,\tau )\underset{𝐩}{}\psi _{𝐩+𝐪,\sigma }^{}(\tau )\psi _{𝐩,\sigma }(\tau )$$
(9)
and $`\rho (𝐪,\tau )=_\sigma \rho _\sigma (𝐪,\tau )`$. The quasiparticle dispersion relation $`ϵ_𝐩`$ is defined in Eq. (2), $`\mu `$ denotes the chemical potential, $`\beta `$ the inverse temperature, $`g^2`$ the coupling constant and $`\psi _{𝐩,\sigma }^{}`$ and $`\psi _{𝐩,\sigma }`$ are Grassmann variables. We measure temperatures, frequencies and energies in the same units. Our effective density interaction consists of an induced part, the last term in Eq. (8), and a local on-site Coulomb repulsion, the second term in Eq. (8).

The retarded generalized susceptibility $`\chi (𝐪,\omega )`$ that defines the effective interaction, Eq. (8), is defined in Eq. (3).

The density-fluctuation propagator on the imaginary axis, $`\chi (𝐪,i\nu _n)`$ is related to the imaginary part of the response function $`Im\chi (𝐪,\omega )`$, Eq. (3), via the spectral representation

$$\chi (𝐪,i\nu _n)=_{\mathrm{}}^+\mathrm{}\frac{d\omega }{\pi }\frac{Im\chi (𝐪,\omega )}{i\nu _n\omega }$$
(10)
To get $`\chi (𝐪,i\nu _n)`$ to decay as $`1/\nu _n^2`$ as $`\nu _n\mathrm{}`$, as it should, we introduce a cutoff $`\omega _0`$ and take $`Im\chi (𝐪,\omega )=0`$ for $`\omega \omega _0`$. A natural choice for the cutoff is $`\omega _0=\eta (\widehat{q})\kappa _0^2`$.

The Eliashberg equations for the critical temperature $`T_c`$ in the Matsubara representation reduce, for the effective action Eq. (8), to

$$\mathrm{\Sigma }(𝐩,i\omega _n)=g^2\frac{T}{N}\underset{\mathrm{\Omega }_n}{}\underset{𝐤}{}\chi (𝐩𝐤,i\omega _ni\mathrm{\Omega }_n)G(𝐤,i\mathrm{\Omega }_n)$$
(11)
$$G(𝐩,i\omega _n)=\frac{1}{i\omega _n(ϵ_𝐩\mu )\mathrm{\Sigma }(𝐩,i\omega _n)}$$
(12)
$`\mathrm{\Lambda }(T)\mathrm{\Phi }(𝐩,i\omega _n)`$ $`=`$ $`{\displaystyle \frac{T}{N}}{\displaystyle \underset{\mathrm{\Omega }_n}{}}{\displaystyle \underset{𝐤}{}}\left(g^2\chi (𝐩𝐤,i\omega _ni\mathrm{\Omega }_n)I\right)|G(𝐤,i\mathrm{\Omega }_n)|^2\mathrm{\Phi }(𝐤,i\mathrm{\Omega }_n)`$ (13)
$`\mathrm{\Lambda }(T)`$ $`=`$ $`1T=T_c`$ (14)
where $`\mathrm{\Sigma }(𝐩,i\omega _n)`$ is the quasiparticle self-energy, $`G(𝐩,i\omega _n)`$ the one-particle Green’s function and $`\mathrm{\Phi }(𝐩,i\omega _n)`$ the anomalous self-energy. The Hartree terms coming from the on-site Coulomb repulsion and induced density interaction have been absorbed in the definition of the chemical potential, which is adjusted to give an electron density of $`n=1.1`$. $`N`$ is the total number of allowed wavevectors in the Brillouin Zone. Eq. (14) has been written for spin-singlet Cooper pairs. In the spin-triplet channel, the sign and coefficient of the term $`g^2\chi (𝐩𝐤,i\omega _ni\mathrm{\Omega }_n)`$ remains unchanged. The on-site Coulomb interaction $`I`$ in Eq. (14) can be ignored since in the for spin-triplet Cooper pairs, the amplitude for the two particles to be on the same site simultaneously vanishes.

The momentum convolutions in Eqs. (11,14) are carried out with a Fast Fourier Transform algorithm on a $`128\times 128`$ lattice for calculations in two dimensions and $`48\times 48\times 48`$ lattice for three dimensional calculations. The frequency sums in both the self-energy and linearized gap equations are treated with the renormalization group technique of Pao and Bickers. We have kept between 8 and 16 Matsubara frequencies at each stage of the renormalization procedure, starting with an initial temperature $`T_0=0.4t`$ in two dimensions and $`T_0=0.6t`$ in three dimensions and cutoff $`\mathrm{\Omega }_c30t`$. The renormalization group acceleration technique restricts one to a discrete set of temperatures $`T_0>T_1>T_2\mathrm{}`$. The critical temperature at which $`\mathrm{\Lambda }(T)=1`$ in Eq. (14) is determined by linear interpolation.


format: text

text:
# Hot Gas in the Circumstellar Bubble S 308

## 1 Introduction

Fast winds from massive stars can interact with their ambient medium and blow bubbles. The physical structure of a bubble in a homogeneous interstellar medium (ISM) has been modeled by Weaver et al. (1977). In their model of a pressure-driven bubble, the shocked fast wind reaches temperatures of 107–108 K and forms a contact discontinuity with the cool (104 K), dense shell of swept-up ambient medium. Heat conduction and mass evaporation across this interface lowers the temperature and raises the density of the hot interior gas, producing 105–106 K gas that can be observed at ultraviolet and X-ray wavelengths.

Of all known bubbles blown by single massive stars, only two show detectable diffuse X-ray emission: NGC 6888 and S 308 (Bochkarev, 1988; Wrigge et al., 1994; Wrigge, 1999). $`ROSAT`$ PSPC and $`ASCA`$ SIS observations of NGC 6888 showed a limb-brightened X-ray morphology and spectra that were best fitted by a dominant component at 1.5$`\times `$106 K and a weaker component at 8$`\times `$106 K (Wrigge et al., 1994, 1998, 2003). $`ROSAT`$ PSPC observations of S 308 detected X-ray emission near the shell rim, but detailed analysis was hampered by the low surface brightness and the obstruction of the PSPC window support structure projected near the shell rim (Wrigge, 1999).

Both NGC 6888 and S 308 have Wolf-Rayet (WR) central stars and show enhanced N and He abundances in the nebulae, suggesting that they consist of stellar material ejected by the central stars during the red supergiant (RSG) phase (Esteban et al., 1992). If the RSG mass loss rate was constant, the density of the circumstellar material would fall off with a $`r^2`$ dependence, where $`r`$ is the radius. The formation and evolution of a bubble blown in such a circumstellar medium has been calculated analytically by García-Segura & Mac Low (1995), and numerically simulated by García-Segura et al. (1996) specifically for WR central stars with RSG progenitors.

Wrigge et al. (1998, 2003) modeled the diffuse X-ray emission from NGC 6888 by incorporating Weaver et al.’s (1977) heat conduction in the analytical solution of García-Segura & Mac Low (1995), and found that the model could not explain simultaneously the observed limb-brightened morphology, surface brightness, and the presence of a high temperature component. By reducing the efficiency of heat conduction at the interface, it was possible to reproduce the limb-brightened morphology but not the surface brightness.

The recently launched XMM-Newton Observatory has unprecedented sensitivity in the soft X-ray band, and thus provides an excellent opportunity to observe the hot gas in bubble interiors. We have obtained XMM-Newton observations of S 308. This paper reports our analysis of this new X-ray observation and discusses the implications of our results on the physical structure of S 308’s hot interior.

## 2 XMM-Newton Observations of S 308

S 308 was observed with the XMM-Newton Observatory on 2001 October 23-24, in Revolution 343, using the EPIC/MOS1, EPIC/MOS2, and EPIC/pn CCD cameras (Observation ID 79570201). The two EPIC/MOS cameras were operated in the full-frame mode for a total exposure time of 47.0 ks, and the EPIC/pn camera in the extended full-frame mode for 43.0 ks. For all observations, the medium optical blocking filter was used. The angular resolution of the EPIC cameras at energies below 1 keV is $``$6″, and the energy resolutions of the EPIC/MOS and EPIC/pn CCDs are 60 eV and 100 eV at 1 keV, respectively (Dahlem, 1999; Strüder et al., 2001).

We received the XMM pipeline products, and processed them further using the XMM-Newton Science Analysis Software (SAS version 5.3.3) and the calibration files from the Calibration Access Layer available on 2002 October 22. The event files were screened to eliminate bad events, such as those due to charged particles. For the EPIC/MOS observations, only events with CCD patterns 0–12 (similar to ASCA grades 0–4) were selected; for the EPIC/pn observation, only events with CCD pattern 0 (single pixel events) were selected. As the background directly affects the detectability of diffuse X-ray emission, periods of high background have to be excluded. We assessed the background by binning the data over 100 s time intervals for each instrument in the 10–12 keV energy band; the background was considered high if the count rates were $``$0.3 cnts s-1 for the EPIC/MOS or $``$1.5 cnts s-1 for the EPIC/pn. Unfortunately, 15 ks after the observation of S 308 started, a high-background flare occurred and persisted throughout the rest of the observation. The useful exposure times were thus reduced to 22.7 ks and 16.9 ks for the EPIC/MOS and EPIC/pn observations, respectively.

## 3 Spatial Distribution of Diffuse X-ray Emission from S 308

We have observed only the bright northwest quadrant of S 308 because the field-of-view of the EPIC cameras is smaller than the angular extent of the bubble. To produce an X-ray image of S 308, we first extract images in the 0.25–1.15 keV energy band from the EPIC/pn and EPIC/MOS observations individually, then mosaic these three images together to increase the S/N ratio and to reduce the null exposure in the gaps between CCDs. This energy band includes all photons detected from S 308, as shown later in §4. The images are sampled with a pixel size of 2″, which is adequate for the resolution of EPIC cameras. The mosaicked raw EPIC image is presented in Fig. 1a. This image is adaptively smoothed with Gaussian profiles of FWHM ranging from 2″ to 50″, and divided by a normalized exposure map to remove the telescope and instrument sensitivity variations across the field. The resultant image is presented in grey scale with contours in Fig. 1b.

Diffuse X-ray emission from S 308 is clearly detected, and a limb-brightened morphology is seen. To compare the distribution of diffuse X-ray emission to the cool nebular shell, we present an \[O III\] $`\lambda `$5007 line image of S 308 from Gruendl et al. (2000) in Fig. 1c, and plot the X-ray contours extracted from the smoothed EPIC image over the \[O III\] image in Fig. 1d. The diffuse X-ray emission is completely interior to the optical shell. An apparent gap exists between the outer edge of the diffuse X-ray emission and the outer rim of the optical shell. Close inspection shows bright optical filaments delineating the brightest X-ray-emitting regions. It is not clear whether the gap between the outer edge of the X-ray and optical emission regions is caused by a lack of X-ray emission or additional absorption from the nebular shell.

The optical shell of S 308 shows a protrusion in the northwest corner, suggesting a blowout. This is a very mild case of blowout because its X-ray edge recedes the farthest from the optical shell edge, unlike the superbubble N44’s blowout where the diffuse X-ray emission extends beyond the optical shell (Chu et al., 1993).

## 4 Physical Properties of the Hot Gas in S 308

The physical properties of the hot gas in S 308 can be determined by modeling the observed X-ray spectrum. We have extracted spectra from the event files of the three EPIC cameras separately, using a source aperture outlined in Fig. 1a, and a background aperture exterior to S 308 but with a similar area, $``$210 arcmin2. Small regions encompassing X-ray point sources are excised from the source and background apertures. For illustration, regions around the point sources within the source aperture are marked in Fig. 1a. The background-subtracted, vignetting-corrected spectra binned by a constant channel width of 15 eV are displayed in Fig. 2. The EPIC/pn and EPIC/MOS spectra show similar shapes, but the EPIC/pn camera is more sensitive and has detected a larger number of counts, 7,900 $`\pm `$ 150 cnts. We have thus concentrated on the EPIC/pn spectrum for further analysis and discussion.

The EPIC/pn spectrum of S 308 is very soft, peaking at the He-like triplet of N VI at 0.42–0.43 keV and declining sharply toward higher energies. Above 0.7 keV, the spectrum shows only possibly weak emission at $``$0.92 keV, corresponding to the He-like triplet of Ne IX; no emission is detected beyond 1.0 keV. The presence of line features suggests thermal plasma emission; therefore, we have used the MEKAL optically-thin plasma emission model (Kaastra & Mewe, 1993; Liedahl et al., 1995) to simulate the X-ray spectrum of S 308. The observed spectrum is simulated by passing the MEKAL model spectrum through an interstellar absorption with cross-sections from Balucińska-Church & McCammon (1992) and convolving it with the EPIC/pn response matrices. The simulated spectra are compared to the observed spectrum in the 0.3–2.0 keV range. The best-fit model, judged by the $`\chi ^2`$ statistics, gives the plasma temperature, volume emission measure, and foreground absorption.

The input parameters for the spectral fitting include the temperature and abundances of the emitting plasma, and the column density and abundances of the foreground absorbing medium. It is impractical to allow all these parameters to vary freely in the spectral fits because the spectral resolution is limited and the number of photons detected is modest. We have therefore used independently determined information to constrain the model parameters and fix them whenever possible. First, we adopt the solar abundances for the foreground absorbing material, since S 308 is at a distance of $`1.5\pm 0.3`$ kpc.111Several distances to HD 50896 have been reported, ranging from 0.6+0.4-0.2 kpc measured by Hipparcos (Perryman et al., 1997) to $``$1.8 kpc determined from interstellar absorption lines (Howarth & Schmutz, 1995). We have critically examined the data from the latter work and conclude that the interstellar absorption convincingly demonstrates that HD 50896 cannot be closer than 1 kpc, ruling out the near distance by Hipparcos. The systemic velocity of S 308, derived from the average velocity near the nebular center, and the average velocities of the neighboring H II regions S 303 and S 304 are all near $`V_{\mathrm{LSR}}`$ = $`20\pm 2`$ km s-1 (Chu et al., 1982). The kinematic distance estimated from this radial velocity and the Galactic rotation is $`1.5\pm 0.2`$ kpc. The photometric distance to HD 50896 is $`1.0\pm 0.2`$ kpc (van der Hucht, 2001). Combining the latter three distances, we decide to adopt a distance of $`1.5\pm 0.3`$ kpc to HD 50896. Next we consider the abundances of the X-ray-emitting plasma. While the hot gas contains shocked fast stellar wind, its mass may be dominated by nebular material evaporated across the conduction front or ablated by the fast stellar wind (Weaver et al., 1977; Pittard et al., 2001a, b). Thus, we start the spectral fits using the nebular abundances of S 308. For He, O, N, and Ne, we adopt the abundances determined from optical spectrophotometry by Esteban et al. (1992). For Mg and Fe, we assume that their abundances relative to the solar values are tied to that of O. The C abundance of the nebula is unknown, although the N/C abundance ratio of the fast wind has been determined to be 14 within a factor of 3 uncertainty (Hillier, 1988). We adopt a 0.1 solar abundance for C, which is somewhat arbitrary but compatible with the stellar wind abundances. The final abundances of C, N, O, Ne, Mg, and Fe we have adopted for the MEKAL model are 0.1, 1.6, 0.13, 0.22, 0.13, and 0.13 times the solar values (Anders & Grevesse, 1989), respectively.

Using the above nebular abundances for single-temperature MEKAL models, it is possible to reproduce the bulk spectral shape with a plasma temperature of $`kT0.1`$ keV (or $`T1.1\times 10^6`$ K). To improve the spectral fits and to search for the possible existence of hotter gas expected from a shocked fast wind, we have also fit the observed spectrum with two temperature components. Only a small improvement is achieved, as the reduced $`\chi ^2`$ of the single-temperature and two-temperature best fits are $`\chi ^2/\mathrm{DoF}`$ = 1.12 and 1.02, respectively. The best-fit single-temperature model and the low-temperature component of the best-fit two-temperature model are virtually indistinguishable from each other; thus only the best-fit two-temperature MEKAL model is overplotted on the EPIC/pn spectrum in Fig. 3. The absorption column density is $`N_\mathrm{H}=`$ 1.1$`\times `$1021 cm-2; the temperatures, normalization factors, observed fluxes, and unabsorbed fluxes of these two temperature components are summarized in Table 1. The low temperature component has an accurate temperature and dominates the X-ray emission, while the high temperature component is not well constrained in temperature because it is weak and contributes only 6% of the observed flux, or 1.5% of the unabsorbed flux.

We have also tried spectral fits using MEKAL models with two temperature components and freely varying abundances, and find that the EPIC/pn spectrum can be reasonably well described by temperature components at $`kT_1=0.1`$ keV and $`kT_2=0.7`$ keV with N, O, and C abundances at 5.6, 0.37, and 0.3 times solar values, respectively. These abundances are 3–4 times as high, but the N/C ratio ($``$18 times the solar value) and the temperatures of the two thermal components are similar to those of our best-fit model with fixed nebular abundances. Furthermore, both fixed abundance and freely varying abundance model fits show that the low-temperature component dominates the observed X-ray flux. For a $`1\times 10^6`$ K plasma with collisional ionization equilibrium, N and C exist mainly in He-like ions, and thus the N/C abundance ratio can be constrained by the relative strengths of N VI lines at 0.43 keV and C V lines at 0.3 keV, or the spectral shape below 0.5 keV. The absolute abundances, on the other hand, are not well constrained because the spectral resolution does not allow unambiguous separation of line emission and bremsstrahlung emission. Therefore, we conclude that the diffuse X-ray emission of S 308 originates mostly from N-enriched gas at a temperature of $`1.1\times 10^6`$ K. This temperature is noticeably lower than the plasma temperatures observed in planetary nebulae or superbubbles (Chu et al., 2001; Guerrero, Gruendl, & Chu, 2002; Kastner et al., 2000, 2001; Dunne et al., 2003; Townsley et al., 2003).

Finally, we have compared our EPIC/pn spectrum of S 308 to Wrigge’s (1999) best-fit model for the ROSAT PSPC spectrum. Wrigge’s model consists of two components at temperatures of $`(1.5\pm 0.1)\times 10^6`$ K and $`(2.8\pm 0.4)\times 10^7`$ K, with the high-temperature component contributing about 1/3 of the total X-ray flux. Our EPIC/pn spectrum clearly rules out the presence of this high-temperature component, which may be caused by a combination of a low S/N in the PSPC data and contaminations from the numerous point sources projected within the boundary of S 308. These point sources are resolved by the XMM-Newton EPIC cameras but not the ROSAT PSPC detector.

## 5 Discussion

### 5.1 Dynamical Parameters of the S 308 Bubble

S 308 is a bubble blown by the WN4 star HD 50896 (Johnson & Hogg, 1965) in a circumstellar medium that was ejected by its progenitor during the RSG phase (Esteban et al., 1992). At a distance of $`1.5\pm 0.3`$ kpc, its 20 radius corresponds to $`9\pm 2`$ pc. The nebular shell of S 308 is photoionized, but additional shock heating is present as the nebular \[O III\]$`\lambda `$5007/H$`\beta `$ ratios are as high as $``$20 (Esteban et al., 1992). The presence of shocks is also evidenced in the nebular morphology at the shell rim where the \[O III\] emission leads the H$`\alpha `$ emission by 16–20′′ (Gruendl et al., 2000). The S 308 shell resides inside an H I cavity evacuated by the stellar wind of HD 50896’s main-sequence progenitor, with the optical nebulae S 303 and S 304 (Fig. 1a of Chu et al., 1982) as part of the fossil swept-up interstellar shell (Arnal & Cappa, 1996). The overall nebular structure and environment of S 308 agree well with the expectations of García-Segura et al.’s (1996) models for WR bubbles.

Assuming that the RSG progenitor had a constant mass loss rate $`\dot{M}_{\mathrm{RSG}}`$ and a constant wind velocity $`V_{\mathrm{RSG}}`$, the circumstellar medium of HD 50896 would have a radial density profile $`r^2`$. Assuming further that the WR star has a constant mass loss rate $`\dot{M}`$ and a constant terminal wind velocity $`V_{\mathrm{}}`$, the bubble expansion velocity will be constant with

$$V_{\mathrm{exp}}=(\dot{M}V_{\mathrm{}}V_{\mathrm{RSG}}/\dot{M}_{\mathrm{RSG}})^{1/2}$$
(1)
(García-Segura & Mac Low, 1995). For a shell expansion velocity of $`63\pm 3`$ km s-1 (Chu et al., 1982) and a radius of $`9\pm 2`$ pc, the dynamic age of S 308 is thus $`(1.4\pm 0.3)\times 10^5`$ yr. The terminal velocity of HD 50896’s fast stellar wind has been well measured to be $`V_{\mathrm{}}`$ = 1,720 km s-1 (Prinja et al., 1990). The mass loss rate of HD 50896’s fast wind has been determined by Nugis et al. (1998) from IR and radio continuum measurements assuming a smooth wind or a clumped wind for a distance of 1.8 kpc. With the distance adjusted from 1.8 kpc to $`1.5\pm 0.3`$ kpc, the smooth-wind and clumping-corrected mass loss rates of HD 50896 are $`(3.6\pm 1.1)\times 10^5`$ M yr-1 and $`(1.4\pm 0.5)\times 10^5`$ M yr-1, respectively. If we assume that the high electron temperature measured in the nebular shell, $`(1.41.8)\times 10^4`$ K (Esteban et al., 1992), is caused by shock heating, then the shock velocity has to be $``$30 km s-1. The RSG wind velocity can be approximated by the difference between the shell expansion velocity and the shock velocity, thus $`V_{\mathrm{RSG}}`$ 30 km s-1; from Eq. (1) we then find the RSG mass loss rate to be $`\dot{M}_{\mathrm{RSG}}(4.7\pm 1.4)\times 10^4`$ M yr-1 for the case of smooth fast wind, or $`(1.9\pm 0.6)\times 10^4`$ M yr-1 for the case of clumped fast wind.

The nebular shell of S 308 has a sharp rim, indicating that it is still surrounded by RSG wind material; however, the RSG wind material must not extend much further in radius, as a breakout has already occurred at the northwestern part of the shell. The unperturbed RSG wind material just exterior to the nebular shell was ejected at least $`(2.9\pm 0.6)\times 10^5`$ years ago. As the fast WR wind started $`(1.4\pm 0.3)\times 10^5`$ yr ago, the RSG wind lasted $`(1.5\pm 0.3)\times 10^5`$ yr, and the total RSG mass loss would be $`70_{30}^{+45}`$ M and $`29_{13}^{+15}`$ M for the cases of smooth and clumped fast winds, respectively. Evidently, the former value is unrealistically large for a normal massive star; thus the total RSG mass loss confirms Nugis et al.’s (1998) suggestion that the clumping-corrected mass loss rates of WR stars are more accurate. Therefore, we will use only the clumping-corrected mass loss rate for the rest of the discussion. The physical parameters of the RSG progenitor and the nebular shell derived above are summarized in Table 2. Note that the velocity, mass loss rate, and total amount of mass loss of the RSG wind derived from bubble dynamics are all well within the canonical ranges of these parameters, and therefore provide a very satisfactory consistency check for the bubble dynamics.

### 5.2 The Hot Interior of S 308

The XMM-Newton observations of S 308 show an X-ray spectrum dominated by emission from a N-enriched medium at a temperature of $`1.1\times 10^6`$ K, with less than 6% of the observed flux contributed by hotter gas. This dominant temperature is almost two orders of magnitude lower than the post-shock temperature of the 1,720 km s-1 fast stellar wind, $`8.3\times 10^7`$ K, for a H/He ratio of 0.2 in the WR wind (Nugis et al., 1998). The cooling time scale for such hot gas is much longer than the dynamic age of the bubble; therefore, the observed low plasma temperature indicates that cool nebular material has been mixed into the shocked WR wind to lower the temperature. This mixing is supported by the limb-brightened X-ray morphology of S 308.

Two known mechanisms can mix nebular material with the hot gas: mass evaporation due to thermal conduction and mass loading due to dynamic ablation. The former mechanism has been incorporated into bubble models of Weaver et al. (1977) and Pittard et al. (2001a), and the latter mechanism has been used by Pittard et al. (2001b). To determine whether evaporation or ablation is responsible for the mixing in S 308, we need to know whether the WR wind is shocked and where the shock front is. In the bubble model of Weaver et al. (1977), the fast stellar wind encounters an inner, stagnation shock and the shocked stellar wind provides the thermal pressure to drive the expansion of the nebular shell. Thermal conduction at the interface between the hot interior gas and the cool nebular shell evaporates nebular material into the hot interior, but the bulk of the shocked fast wind remains hot and occupies most of the volume in the bubble interior. We may estimate the X-ray emission expected from the “uncontaminated” shocked stellar wind of HD 50896 and compare it to the observed diffuse X-ray emission from S 308. The WR wind of HD 50896 has injected $`2\pm 1`$ M into the bubble interior during the lifetime of the bubble. If the shocked fast wind is uniformly distributed within the bubble interior (the worst-case scenario for emissivity and detectability), the electron density would be $`0.015\pm 0.003`$ cm-3, for a H/He ratio of 0.2. The total volume emission measure of this shocked wind at $`8.3\times 10^7`$ K can be scaled to determine its expected normalization factor $`A`$ (see Table 1 for the definition of the normalization factor) for the aperture with which our EPIC spectrum was extracted; we find an expected normalization factor to be $`(1.2\pm 0.1)\times 10^4`$ cm-5. Compared with the two temperature components from the best spectral fit listed in Table 1, the normalization factor expected for the X-ray emission from the uncontaminated shocked wind is fortuitously close to that of the $`T_2`$ component. If the hot shocked wind occupies only a fraction of the total volume, the expected normalization factor will be higher and would have been detected by our XMM-Newton observation if its useful exposure time was not shortened by a high background. Deeper X-ray observations are needed to place more stringent limits on the existence or absence of the shocked stellar wind at $`8.3\times 10^7`$ K.

We have used the best-fit spectral model to determine the unabsorbed X-ray flux from S 308, $`7.2\times 10^{12}`$ ergs cm-2 s-1 within the 210 arcmin2 aperture in the 0.25–1.5 keV energy band. This X-ray flux can be scaled to the entire area covered by the bubble to estimate the total X-ray luminosity. As the observed aperture contains the brightest patch of diffuse X-ray emission within S 308, the X-ray luminosity of S 308 should be $`(1.2\pm 0.5)\times 10^{34}`$ ergs s-1. The rms electron density of the X-ray-emitting gas in S 308 can be determined from the normalization factor $`A`$ given in Table 1, $`N_\mathrm{e}=(0.20\pm 0.03)ϵ^{1/2}`$ cm-3, where $`ϵ`$ is the volume filling factor of the hot gas. For $`ϵ`$ = 0.5, the rms $`N_\mathrm{e}=0.28\pm 0.04`$ cm-3 and the total mass of the hot gas is $`11\pm 5`$ M; while for $`ϵ`$ = 0.1, the rms $`N_\mathrm{e}=0.63\pm 0.09`$ cm-3 and the total mass of the hot gas is $`5\pm 2`$ M. As the limb-brightened morphology of the diffuse X-ray emission indicates a thick shell structure, the filling factor must be closer to 0.5 than 0.1. Therefore, the hot gas mass is much larger than the total WR wind mass that has been injected into the bubble interior, $`2\pm 1`$ M, suggesting that nebular material contributes to the hot gas mass and justifying the fixed nebular abundances we used for the spectral fits in §4.

Our present XMM-Newton observations do not detect enough photons to warrant a spatially resolved spectral analysis of S 308. Future deeper observations are needed to determine the radial temperature structure of the X-ray emitting gas. A complete mapping of S 308 is also needed to determine whether the varying morphology of the optical shell rim coincides with spatial and physical variations in the adjacent hot, X-ray emitting gas in the bubble interior.

### 5.3 Interface between the Hot Interior and the Cool Shell

UV observations in conjunction with X-ray observations of S 308 allow us to probe the interface between the hot gas in the bubble interior and the cool gas in the swept-up shell. High-resolution UV spectra of HD 50896 have detected a nebular N V absorption line component associated with the approaching side of the S 308 shell, and its large column density suggests an overabundance in nitrogen (Boroson et al., 1997). As the ionization potential of N IV is 77.5 eV, N V must be produced by collisional ionization in gas at temperatures of a few $`\times 10^5`$ K, which is expected at the interface between the hot interior gas and the cool nebular shell of a bubble. The detection of N V absorption from the S 308 shell has provided the strongest evidence of the existence of an interface layer.

Our XMM-Newton X-ray image of S 308 reveals the relative location of the hot interior gas and the cool nebular shell (see Fig. 1d). The projected outer edge of the diffuse X-ray emission is offset from the projected outer edge of the \[O III\] emission by 90′′ to over 200′′, corresponding to 0.5 – 1.7 pc. As the \[O III\] rim depicts the outer shock associated with the swept-up shell advancing into the unperturbed RSG wind, this “gap” contains the swept-up shell and the interface layer where the temperature drops radially from 106 K to 104 K. The large angular size of this gap makes it the most promising site where an interface layer may be spatially resolved and studied in detail. Future long-slit high-dispersion spectra of the H$`\alpha `$ \+ \[N II\] lines are needed to determine kinematically the exact width of the swept-up shell. We have been awarded Far Ultraviolet Spectroscopic Explorer (FUSE) observations of the O VI emission for three positions in S 308. While these FUSE observations of O VI emission may reveal the existence of $`3\times 10^5`$ K gas at the interface, the relatively large aperture ($`30^{\prime \prime }\times 30^{\prime \prime }`$) provides only modest spatial resolution. Future high-resolution, long-slit Hubble Space Telescope STIS UV spectra of the N V and C IV emission lines at the gap are needed to probe the temperature structure in the interface layer. Only through the analysis of these detailed multi-wavelength observations can we finally learn empirically how fast stellar winds interact with the ambient medium and test the various theoretical models of bubbles.

## 6 Summary and Conclusions

We have obtained XMM-Newton EPIC observations of the northwest quadrant of the circumstellar bubble S 308 blown by the WN4 star HD 50896. Diffuse X-ray emission is unambiguously detected, showing a limb-brightened morphology with a clear gap extending from the outer edge of the diffuse X-ray emission to the outer rim of the nebular shell. The X-ray spectrum of the diffuse emission is very soft, and is well fitted by an optically thin plasma model for a plasma at temperatures of $`1.1\times 10^6`$ K. For such a low temperature, the N/C abundance ratio is well constrained by the spectral shape below 0.5 keV; spectral fits with freely varying abundances and fixed nebular abundances show consistently that the N/C ratio of the hot gas in S 308 is 15–20 times the solar value. A hotter component may exist, but its temperature is not well constrained, $`8_6^{+17}\times 10^6`$ K, and contributes less than 6% of the observed X-ray flux.

We have applied García-Segura & Mac Low’s (1995) analytical model for a WR bubble blown in the progenitor’s RSG wind to S 308, using observed bubble dynamics and stellar wind parameters, and derived the masses lost by HD 50896 via the RSG wind and the WR wind. We have also used the X-ray data to determine the mass of the X-ray-emitting gas in the bubble interior. The results are summarized in Table 2. It can be seen that the mass in the hot gas is indeed much larger than the total mass injected by the WR wind, indicating that the hot gas is dominated by the nebular material mixed into the bubble interior. It is not clear whether the mixing is through conductive evaporation or dynamic ablation. Deep X-ray observations are needed to detect hotter gas and to resolve the temperature structure in the bubble interior. Spatially resolved UV spectroscopic observations of the interface between the X-ray-emitting gas and the cool nebular gas are also needed to probe the physical mechanisms that mix the nebular material into the hot interior.

This research was supported by the NASA grant NAG 5-10037.


text:
# MODEST: modeling stellar evolution and (hydro)dynamics

## 1. Introduction

Large-scale computer simulations in astrophysics require teamwork. Gone are the days that an individual graduate student could write from scratch all the software needed to model a complex astrophysical system. Whether a simulation requires a state-of-the-art stellar evolution code, hydrodynamics code or stellar dynamics code, it has become standard practice to use a ‘legacy code’. One may write some modifications and additions, but it is hard to compete with the tens of person-years that have gone into some of the standard codes.

As a result, codes have steadily grown in complexity within separate disciplines in astrophysics, but until recently little effort has been invested to make these legacy codes compatible across different areas of astrophysics. In the area of dense stellar systems, the topic of this meeting, the result has been that we have access to very detailed codes that can model the evolution of individual stars, other codes that can handle collisions between stars, and yet other codes that can follow the motion of the whole ensemble of stars. However, none of these codes come even close to being able to talk with each other.

To let different codes talk together, similar requirements hold as for humans trying to communicate between cultures and across language barriers: each person needs to adapt to some extent to the customs of the other culture, and there needs to be a dictionary to translate between the languages.

The first requirement reads like a medical oath: first of all, do no harm, i.e. do not cause the combined codes to halt. It is better that a large-scale simulation is allowed to come to completion, with messages in a log file indicating where serious errors may have been made, than to have the simulation halt at each bend of the road, when one of the modules is not optimally happy.

The second requirement boils down to the definition of specific interfaces, and a willingness in the community to adopt a standard definition. Producing such a definition requires a fair amount of care, flexibility and vision, to avoid the danger of saddling a field with a standard that can hinder future progress.

In the next few pages, I summary recent efforts to let codes cross the boundaries of stellar evolution, stellar dynamics and stellar hydrodynamics, in order to enable next-generation large-scale simulations of dense stellar systems.

## 2. MODEST-1: New York

It was the goal of the MODEST-1 meeting to begin addressing this problem of letting codes talk to each other. The workshop was held at the Hayden Planetarium in the Museum for Natural History in New York City, in June 2002. The MODEST acronym was coined during this meeting, and it can stand not only for MOdeling DEnse STellar systems, but also for MODifying Existing STellar codes. The latter description stresses the desirability to start with what is already available, and to find ways to put it all together, rather than to try to write a kitchen-sink type over-arching super code from scratch.

The format of the meeting was unconventional. The organizers, Piet Hut and Mike Shara, polled whether there was enough interest to warrant a workshop, and then sent out an email with the following note: “this will be a work-workshop, not a lectures-workshop. Believe it or not, we have not a single scheduled lecture! We are curious to see how this format will work out. Hopefully we will actually get a lot of work done, together as well as in smaller groups.”

We did. The 34 participants met for five days, splitting up in topical groups and reporting their results back to the group on a daily basis. Major progress was made in defining and starting to tackle the two main requirements listed above, that of culture and translation. Within a few weeks of the meeting, 10 of the participants had written a review paper, in lieu of proceedings, which was posted on astro-ph, and subsequently published (Hut et al. 2003).

Also immediately following the workshop, a web site was established, at http://www.manybody.org/modest.html , as a place to accumulate general information, demos, toy models, links to more detailed models and simulation software, etc. At the same time an email list was started, with occasional announcements of major events and progress reports. The web site contains instructions for subscribing to this email list.

As a result of the great enthusiasm displayed during the first workshop, it was decided to hold similar meetings frequently enough so as to enable ongoing projects to be critiqued at each meeting. Such peer-review feedback is essential for the health of any really large-scale project, and therefore a frequency of two meetings per year was chosen as an optimal compromise. While individual projects could of course use more frequent reviews, getting the main players together every half year was already seen as quite a challenge. It was a challenge that was met, it turns out, judging from the attendance of each of the following two workshops, between 30 and 40 people, comparable to that of MODEST-1.

## 3. MODEST-2: Amsterdam

MODEST-2 was organized by Simon Portegies Zwart en Piet Hut, at the University ofu Amsterdam, Holland, in December 2002. This workshop followed a similar format: there were no scheduled talks. Instead, individual participants sometimes showed a couple view graphs, to illustrate a particular point. Also, during the first day a number of specific topics of interest were formulated, and speakers were invited to give an impromptu brief introduction to each topic. In this way, eight short talks were delivered by participants outlining how their own work was fitting into the MODEST framework, what they wanted to get out of participation in MODEST, what the most relevant questions were for their area, and what they had accomplished since MODEST-1.

During the previous workshop in New York, a broad discussion had been started concerning the question when to use full-blown stellar evolution codes and when to use approximate recipes. This discussion was continued in Amsterdam. The conclusions reached were similar, but more refined in detail: to use recipes for single stars; a mix of recipes and codes for interacting binaries; and ‘life’ stellar evolution codes for merger products. Given the very many ways that stars can form from single or even repeated mergers, there is just no way that we can expect to construct a grid of model tracks to anticipate the specific needs for evolving such merger products. Not only will the chemical compositions be different from standard values, but incomplete mixing of the progenitor stars will add a full functional amount of degrees of freedom.

Another outcome of this workshop was the organization of eight working groups, listed below (the ninth working group was added during MODEST-3). The main MODEST web site now contains pointers to these working groups. While Steve McMillan is the webmaster for the main web site, the following people are the contact persons for the web sites of the individual working groups:

| Working Group |  | Contact Person |
| --- | --- | --- |
| Star Formation |  | Ralf Klessen |
| Stellar Evolution |  | Onno Pols |
| Stellar Dynamics |  | Rainer Spurzem |
| Stellar Collisions |  | Marc Freitag |
| Simulating Observations of Simulations |  | Simon Portegies Zwart |
| Data Structures |  | Peter Teuben |
| Validation |  | Douglas Heggie |
| Literature |  | Melvyn Davies |
| Observations |  | Giampaolo Piotto |

The highlights of the workshop appeared on astro-ph within a few weeks after concluding the meeting, and were published as a review (Sills et al. 2003).

## 4. MODEST-3: Melbourne

The third workshop was held at Monash University in Melbourne, Australia, and was organized by Rosemary Mardling and Piet Hut. The meeting was held in early July 2003, in the three days before the General Assembly assembled in Sydney. In contrast to the previous two workshops, we experimented with a more traditional style of scheduled presentations (see the MODEST-3 web site http://www.manybody.org/modest/Workshops/modest-3.html for the detailed program with speakers, titles, and abstracts).

These talks filled the first two days of the meeting. In the morning of the third day, each of the working groups reported about their progress, and a ninth working group was added to focus on observations. The afternoon featured an open discussion. One of the outcomes was a series of proposals for extending the web site, by including a page for job opportunities and another page for project proposals. Examples of the latter can be notes from people who have observational data sets, for which they invite theoreticians to join them in simulations of the systems observed; or simulators who have large data sets for which they invite students to join them in analyzing the results.

After all the informal discussions during the first two workshops, this meeting with more formal talks played a complementary role, as a one-time occasion to inform people from all three fields (evolution, dynamics, hydrodynamics) about the main activities in the other fields. Subsequent workshops will return to the old format of more free-wheeling discussions, starting with MODEST-4, in January 2004 at Geneva Observatory, located between Geneva and Lausanne in Switzerland.

## 5. Future Meetings and Activities

The MODEST initiative has grown in one year from an initial informal meeting into a broad framework to facilitate collaborations across boundaries between various disciplines within computational astrophysics. So far, these collaborations have taken on the form of sharing codes, defining and building code interfaces, developing demos, writing review papers, applying for joint grants, and starting joint research projects. Future workshops may coincide with the organization of MODEST summer schools and other teaching activities, as well as outreach projects aimed at the general public. MODEST is an open forum: we invite anyone interested in simulations of dense stellar systems to join us.

## References

Hut, P., Shara, M.M., Aarseth, S.J., Klessen, R.S., Lombardi, J.C., Makino, J., McMillan, S., Pols, O.R., Teuben, P.J., & Webbink, R.F. 2003, New Astronomy 8, 337

Sills, A., Deiters, S., Eggleton, P., Freitag, M., Giersz, M., Heggie, D., Hurley, J., Hut, P., Ivanova, N., Klessen, R.S., Kroupa, P., Lombardi, J.C., McMillan, S., Portegies Zwart, S., & Zinnecker, H. 2003, New Astronomy 8, 605


text:
# Lattice QCD at finite isospin density and/or temperatureTalk presented by D. K. Sinclair, Lattice2003, Tsukuba.

## 1 Introduction

QCD at finite temperature and/or densities is relevant to the physics of the early universe, neutron stars and relativistic heavy-ion collisions — RHIC and the CERN heavy-ion program.

Because of the difficulties involved with simulating QCD at finite chemical potential $`\mu `$ for quark number with its complex fermion determinant, we are simulating $`N_f=2`$ lattice QCD with a finite chemical potential $`\mu _I`$ for isospin ($`I_3`$), which has a real positive fermion determinant . We are currently studying the $`\mu _I`$ dependence of the finite temperature transition . At small $`\mu `$, the Bielefeld-Swansea collaboration has observed that the phase of the fermion determinant is sufficiently well behaved that the $`\mu `$ and $`\mu _I`$ dependence of the transition are identical . Our predictions are in good agreement with those of de Forcrand and Philipsen .

At large values of $`\mu _I`$ ($`\mu _I>m_\pi `$), the low temperature phase is characterized by a pion condensate. The finite temperature transition is now a true phase transition, which appears to be second order for lower $`\mu _I`$ values and first order at higher $`\mu _I`$ values.

We also simulate lattice QCD with an irrelevant chiral 4-fermion term which allows us to simulate at zero quark mass, giving us direct access to the critical exponents at the finite temperature transition . We are using $`16^3\times 8`$ and $`24^3\times 8`$ lattices for these simulations.

In section 2 we present preliminary results of our finite $`\mu _I`$ and temperature simulations. Section 3 gives some preliminary graphs from our $`N_t=8`$ finite-temperature simulations using our modified action. Finally we present our conclusions and indicate further avenues of research in section 4.

## 2 QCD at finite $`\mu _I`$ and temperature

The quark part of our lattice action is

$$S_f=\underset{sites}{}\left[\overline{\chi }[D/(\frac{1}{2}\tau _3\mu _I)+m]\chi +i\lambda ϵ\overline{\chi }\tau _2\chi \right].$$
(1)
For simulations at $`\mu _I<m_\pi `$, we set the symmetry breaking parameter $`\lambda =0`$.

We perform simulations on an $`8^3\times 4`$ lattice with $`m=0.05`$, $`\lambda =0`$ and $`0\mu _I0.55`$, for a set of $`\beta `$s covering the crossover region for each $`\mu _I`$. We measure the chiral condensate, the plaquette, the Wilson Line and the isospin density and their corresponding susceptibilities for each set of parameters. The position of the crossover was obtained as the peak of the susceptibilities. This was determined using Ferrenberg-Swendsen reweighting techniques . Figure 1 shows the chiral susceptibilities obtained from such reweightings for 3 of the 7 $`\mu _I`$ values we use. The multiple values for each $`\mu _I`$ are the results of using distributions from several $`\beta `$s close to the peak.

The crossover $`\beta `$ values, $`\beta _c`$, obtained from the peaks in the susceptibilities for each of the 4 observables, are plotted in figure 2 as functions of $`\mu _I^2`$, since the leading term is expected to be quadratic in $`\mu _I`$. We note that the 4 predictions for each $`\mu _I`$ appear consistent within errors, indicating that this is a reasonable definition of the position of the crossover. The straight line

$$\beta _c=5.3220.143\mu _I^2$$
(2)
in this figure is only meant as a rough guide. Using the relation $`\beta _c(\mu )=\beta _c(\mu _I=2\mu )`$, which should hold for small $`\mu ,\mu _I`$ our results are consistent with those of de Forcrand and Philipsen .

We are repeating these simulations for $`m=0.1`$ and $`m=0.2`$ where the critical $`\mu _I`$s, $`\mu _c=m_\pi `$ are larger. In no case does there appear to be a critical end point, beyond which the transition becomes first order, for $`\mu _I<m_\pi `$.

We are extending our $`m=0.05`$ simulations to $`\mu _I>m_\pi `$. While the thermal phase transition where the pion condensate evaporates appears to be second order for $`\mu _I=0.6`$, simulations on $`16^3\times 4`$ lattices show indications of the metastability expected for a first order transition, at $`\mu _I=0.8`$.

## 3 $`\chi `$QCD” at finite temperature

We simulate lattice QCD with 2 flavours of massless staggered quarks and an irrelevant chiral 4-fermion interaction which allows us to run at zero quark mass, at finite temperature on $`16^3\times 8`$ and $`24^3\times 8`$ lattices. Using $`N_t=8`$ should free us from the lattice artifacts which were present for $`N_t=4,6`$ .

Having massless quarks should enable us to extract the critical exponent $`\beta _m`$ which describes the vanishing of the chiral condensate as the chiral transition is approached from below, and the critical $`\beta `$ ($`\beta _c`$). Running at $`\beta _c`$ with small quark masses will give us the critical exponent $`\delta `$.

Preliminary results for the chiral condensate and Wilson Line are given in figure 3. These clearly indicate that the transition occurs in the range $`5.530<\beta _c<5.545`$ and appears sufficiently smooth to be second order.

## 4 Conclusions

We have determined the $`\mu _I`$ dependence of the finite temperature transition for 2-flavour lattice QCD. For small $`\mu `$, the fluctuations of the phase of the fermion determinant are small enough that the $`\mu `$ and $`\mu _I`$ dependence of this transition are the same . The $`\mu _I`$ dependence we measure predicts a $`\mu `$ dependence consistent with that obtained from simulations with imaginary $`\mu `$ . We will check if it is also consistent with series expansions around $`\mu =0`$ . Our simulations at 3 different quark masses indicate that the transition is a crossover for $`\mu _I<m_\pi `$. At higher $`\mu _I`$, where the the pion condensate evaporates at the transition, we see evidence for a change to first order behaviour, but it is unclear if this is related to the critical endpoint expected at finite $`\mu `$.

Others have estimated the position of the critical endpoint at finite $`\mu `$ for the more physical $`2+1`$-flavour QCD , the closely related 3-flavour QCD  and 4-flavour QCD . We are extending our simulations to 3-flavours, where one can tune the quark mass to make the critical endpoint as close to $`\mu =0`$ as desired. If the endpoint is at small enough $`\mu `$, we expect a critical endpoint at the corresponding $`\mu _I`$, giving another estimate for the position of this endpoint.

We are using the “$`\chi `$QCD” action to determine the critical exponents for the chiral transition of finite temperature QCD, on $`N_t=8`$ lattices.

We thank F. Karsch, P. de Forcrand, O. Philipsen and D. Toublan for helpful discussions. These simulations were performed on the IBM SP at NERSC and the Jazz cluster at Argonne.


text:
# Acknowledgements

## Acknowledgements

This collaboration has been encouraged by the stimulating discussion in the Summer Institutes 2003. We would like to thank to A. Sugamoto, T. Kugo and N. Okamura and T. Kurimoto who encouraged us very much. We are also indebted to SI2003 by many helpful discussions. M. Bando and M. Tanimoto are supported in part by the Grant-in Aid for Scientific Research No.12047225 and 12047220.