text:
# Light front field theory of quark matter at finite temperature11footnote 1Based on invited talks presented at the Workshop on “Effective quark quark interactions” in Bled (Slovenia) July 7-14, 2003, and the 310th WE-Heraeus Seminar “Quarks in Hadrons and Nuclei II”, Rothenfels Castle, Oberwölz (Austria) September 15-20, 2003.

## 1 Introduction

Lattice calculations of quantum chromodynamics (QCD) give firm evidence that nuclear matter undergoes a phase transition to a plasma state at a certain temperature $`T_c`$ of about 170 to 180 MeV. Calculations have been performed at a chemical potential $`\mu =0`$. Recent results are for staggered fermions  and renormalization group improved Wilson fermions . The low density region reflects, e.g., the scenario during the evolution of the early universe. To achieve information from lattice calculations at small $`\mu `$ several methods have recently been developed, i.e., multiparameter reweighting , Taylor expansion at $`\mu 0`$ , imaginary $`\mu `$ . The region of validity is approximately $`\mu T`$ . Effective approaches to QCD indicate an extremely rich phase diagram also for $`\mu >T`$ . Experimentally the QCD phase diagram is accessible to heavy ion collisions. In particular relativistic heavy ion collision at SPS/CERN and RHIC/BNL explore the region where hadronic degrees of freedom are expected to be dissolved. Some results of RHIC are now available that give hints of a non-hadronic state of matter .

On the other hand light-front quantization of QCD can provide a rigorous alternative to lattice QCD . Although the calculational challenge in real QCD of 3+1 dimension seems large (as does lattice QCD) it starts also from the fundamental QCD. The light-front quantization of QCD has the particular advantage that it is completely formulated in physical degrees of freedom. It has emerged as a promising method for solving problems in the strong coupling regime. Light front quantization makes it possible to investigate quantum field theory in a Hamiltonian formulation . This makes it well suited for its application to systems of finite temperature (and density). The relevant field theory has to be quantized on the light front as well, which is presently being developed . I present the light-front field theory at finite temperatures and densities in the next section. For the time being it is applied to the Nambu-Jona-Lasinio model (NJL) model  that is a powerful tool to investigate the non-perturbative region of QCD as it exhibits spontaneous breaking of chiral symmetry and the appearance of Goldstone bosons in a transparent way. Finally, going a step further I shall give the general in-medium light cone time ordered Green functions that allow us to treat quark correlations that lead to hadronization.

## 2 Light front thermal field theory

The four-momentum operator $`P^\mu `$ on the light-front is given by (notation of Ref. )

$$P^\mu =𝑑\omega _+T^{+\mu }(x),$$
(1)
where $`T^{\mu \nu }(x)`$ denotes the energy momentum tensor defined through the Lagrangian of the system and $`S_\mu `$ is the quantization surface. The Hamiltonian is given by $`P^{}`$. To investigate a grand canonical ensemble we need the number operator,

$$N=𝑑\omega _+j^+(x),$$
(2)
where $`j^\nu (x)`$ is the conserved current. These are the necessary ingredients to generalize the covariant partition operator at finite temperature  to the light-front. The grand canonical partition operator on the light-front is given by

$$Z_G=\mathrm{exp}\left\{𝑑\omega _+[\beta _\nu T^{+\nu }(x)+\alpha J^+(x)]\right\},$$
(3)
where $`\alpha =\mu /T`$, with the Lorentz scalars temperature $`T`$ and chemical potential $`\mu `$. The velocity of the medium is given by the time-like vector $`u_\nu u^\nu =1`$ , and $`\beta _\nu =u_\nu /T`$. We choose the medium to be at rest, $`u^\nu =(u^{},u^+,\stackrel{}{u}^{})=(1,1,0,0)`$. The grand partition operator then becomes

$$Z_G=e^{K/T},K\frac{1}{2}(P^{}+P^+)\mu N$$
(4)
with $`P^\pm `$ and $`N`$ defined in (1) and (2). The density operator for a grand canonical ensemble  in equilibrium follows

$$\rho _G=(\mathrm{Tr}e^{K/T})^1e^{K/T}.$$
(5)
The corresponding Fermi distribution functions of particles $`f^+f`$ and antiparticles $`f^{}`$ are given by

$$f^\pm (k^+,\stackrel{}{k}_{})=\left[\mathrm{exp}\left\{\frac{1}{T}\left(\frac{1}{2}k_{\mathrm{on}}^{}+\frac{1}{2}k^+\mu \right)\right\}+1\right]^1$$
(6)
and $`k_{\mathrm{on}}^{}=(\stackrel{}{k}_{}^2+m^2)/k^+`$. This fermionic distribution function (for particles) on the light-front has first been given in . The fermi function for the canonical ensemble can be achieved by simply setting $`\mu =0`$. This then coincides with the distribution function given recently in Ref.  (up to different metric conventions).

The light-front time-ordered Green function for fermions is

$$i𝒢_{\alpha \beta }(xy)=\theta (x^+y^+)\mathrm{\Psi }_\alpha (x)\overline{\mathrm{\Psi }}_\beta (y)\theta (y^+x^+)\overline{\mathrm{\Psi }}_\beta (y)\mathrm{\Psi }_\alpha (x).$$
(7)
We note here that the light-cone time-ordered Green function differs from the Feynman propagator $`S_F`$ in the front form by a contact term $`\gamma ^+/2k^+`$ and therefore coincides with the light-front propagator given previously in Ref. . To evaluate the ensemble average $`\mathrm{}=\mathrm{Tr}(\rho _G\mathrm{})`$ of (7), we utilize the imaginary time formalism . We rotate the light-front time of the Green function to imaginary value. Hence the $`k^{}`$-integral is replaced by a sum of light-front Matsubara frequencies $`\omega _n`$ according to ,

$`{\displaystyle \frac{1}{2}}k^{}i\omega _n{\displaystyle \frac{1}{2}}k^++\mu {\displaystyle \frac{1}{2}}k_n^{}{\displaystyle \frac{1}{2}}z,`$ (8)
where $`\omega _n=\pi \lambda T`$, $`\lambda =2n+1`$ for fermions \[$`\lambda =2n`$ for bosons\]. In the last step we have performed an analytic continuation to the complex plane. For noninteracting Dirac fields the (analytically continued) imaginary time Green function becomes

$`G(z,\overline{k})`$ $`=`$ $`{\displaystyle \frac{\gamma k_{\mathrm{on}}+m}{zk_{\mathrm{on}}^{}+i\epsilon }}{\displaystyle \frac{\theta (k^+)}{k^+}}(1f^+(\overline{k}))+{\displaystyle \frac{\gamma k_{\mathrm{on}}+m}{zk_{\mathrm{on}}^{}i\epsilon }}{\displaystyle \frac{\theta (k^+)}{k^+}}f^+(\overline{k})`$
$`+{\displaystyle \frac{\gamma k_{\mathrm{on}}+m}{zk_{\mathrm{on}}^{}+i\epsilon }}{\displaystyle \frac{\theta (k^+)}{k^+}}f^{}(\overline{k})+{\displaystyle \frac{\gamma k_{\mathrm{on}}+m}{zk_{\mathrm{on}}^{}i\epsilon }}{\displaystyle \frac{\theta (k^+)}{k^+}}(1f^{}(\overline{k})),`$
where $`\overline{k}=(k^+,\stackrel{}{k}_{})`$. For equilibrium the imaginary time formalism and the real time formalism are linked by the spectral function . For $`\mu =0`$ this propagator coincides with that of , but differs from that of .

## 3 Spontaneous symmetry breaking and restoration

### 3.1 NJL model on the light-front

The Nambu-Jona-Lasinio (NJL) originally suggested in  has been reviewed in Ref.  as a model of quantum chromo dynamics (QCD), where also a generalization to finite temperature and finite chemical potential has been discussed. Its generalization to the light-front including a proper description of spontaneous symmetry breaking, which is not trivial, has been done in Ref. , which we use here. The Lagrangian is given by

$$=\overline{\psi }(i\gamma m_0)\psi +G\left((\overline{\psi }\psi )^2+(\overline{\psi }i\gamma _5\tau \psi )^2\right).$$
(10)
In mean field approximation the gap equation is

$$m=m_02G\overline{\psi }\psi =m_0+2iG\lambda \frac{d^4k}{(2\pi )^4}\mathrm{Tr}S_F(k),$$
(11)
where $`\lambda =N_fN_c`$ in Hartree and $`\lambda =N_fN_c+\frac{1}{2}`$ in Hartree-Fock approximation, $`N_c`$ ($`N_f`$) is the number of colors (flavors). For the isolated case $`S_F(k)`$ is the Feynman propagator. Taking only the lowest order in $`1/N_c`$ expansion of the 1-body or 2-body operators the light front gap equation can be achieved by a $`k^{}`$ integration, where in addition $`m_0\stackrel{~}{m}_0`$ and $`G\stackrel{~}{G}`$ have to be renormalized to accommodate the expansion. For details see . The propagator to be used in (11) is given in (2). The gap equation becomes

$$m(T,\mu )=\stackrel{~}{m}_0+2\stackrel{~}{G}\lambda \frac{dk^+d^2k_{}}{2k^+(2\pi )^3}\mathrm{\hspace{0.33em}4}m(T,\mu )(1f^+(k^+,\stackrel{}{k}_{})f^{}(k^+,\stackrel{}{k}_{})).$$
(12)
To regularize (12) we require $`k_{\mathrm{on}}^{}+k^+<2\mathrm{\Omega }`$. As a consequence $`k_1^+<k^+<k_2^+`$ and

$`\stackrel{}{k}_{}^2`$ $`<`$ $`2\mathrm{\Omega }k^+(k^+)^2m^2,`$ (13)
$`k_{1,2}^+`$ $`=`$ $`\mathrm{\Omega }\sqrt{\mathrm{\Omega }^2m^2}.`$ (14)
For the isolated case this regularization is fully equivalent to the Lepage-Brodsky one and the three-momentum cut-off. For the in medium case this $`\mathrm{\Omega }`$ regularization leads to analytically the same expressions as given in  for the instantaneous case . The calculation of the pion mass $`m_\pi `$, the pion decay constant $`f_\pi `$, and the condensate value are also available on the light-front .

### 3.2 Results

The model parameters are adjusted to the isolated system. We use the Hartree approximation, i.e. $`\lambda =N_cN_f=6`$. Parameter values are chosen to reproduce the pion mass $`m_\pi =140`$ MeV, the decay constant $`f_\pi =93`$ MeV, and to give a constituent quark mass of $`m=336`$ MeV, i.e. $`\stackrel{~}{G}=5.51\times 10^6`$ MeV, $`\stackrel{~}{m}_0=5.67`$ MeV, and $`\mathrm{\Omega }=714`$ MeV. The parameters are reasonably close to the cases used in the review Ref. . The difference between the bare mass $`\stackrel{~}{m}_0`$ and the constituent mass is due to the finite condensate, which is $`\overline{u}u^{1/3}=247`$ MeV.

In hot and dense quark matter the surrounding medium leads to a change of the constituent quark mass due to the quasiparticle nature of the quark. The constituent mass as solution of (12) is plotted in Fig. 2 as a function of temperature and chemical potential. The fall-off is related to chiral symmetry restoration, which would be complete for $`m_0=0`$. It is related to the QCD phase transition. For $`T60`$ MeV the phase transition is first order, which is reflected by the steep change of the constituent mass. To keep close contact with the 3M results we have chosen for the $`\mathrm{\Omega }`$ in-medium regulator mass $`\mathrm{\Omega }^2(T,\mu )=\mathrm{\Lambda }_{3\mathrm{M}}^2+m^2(T,\mu )`$ with $`\mathrm{\Lambda }_{3\mathrm{M}}=630`$ MeV fixed for all $`T`$ and $`\mu `$.

We define the phase transition to occur at a temperature at which $`m(T,\mu )`$ is half of the isolated constituent quark mass . The phase diagram is shown in Fig. 2. The line indicates the phase boundary separating the hadronic phase from the quark gluon plasma phase. Results presented in this section are based on an effective interaction in the $`q\overline{q}`$ channel. They have to be supplemented by the medium dependence of $`f_\pi `$ and $`m_\pi `$ that are currently underway.

## 4 Few-particle correlations

We are now interested in the $`qq`$ channel. Since we are going up to the three-particle system, we presently approximate the spin structure. To solve the full three-fermion problem on the light front even for the isolated case is quite a challenge. The spin structure is already rather complex, see e.g. . Therefore the elementary spins are averaged Tr$`\gamma =0`$ (in the medium) and hence, for the time being, we are only dealing with bose type particles however subject to Fermi-Dirac statistics. Our main focus here is to see how such a three-particle system is dynamically influenced by a medium of finite temperature and density; ultimately, how nucleons are formed in the hot and dense environment of a plasma of quarks and gluons as the temperature and the density becomes smaller and how the relevant degrees of freedom in the Fermi function change as the many-particle system undergoes a change to hadronic degrees of freedom. To this end we need to formulate suitable few-body equations that describe clusters of quarks in a medium. In addition, because of the drastic mass change, see Fig. 2, these equations have to be relativistic ones. The equations derived here are based on a systematic quantum statistical framework formulated on the light front using a cluster expansion for the Green functions. The formalism has been given elsewhere . We repeat here the basics to make a connection to the previous sections. The light-front time ordered cluster Green function is defined by

$$i𝒢_{\alpha \beta }(xy)=\theta (x^+y^+)A_\alpha (x)\overline{A}_\beta (y)\theta (y^+x^+)\overline{A}_\beta (y)A_\alpha (x).$$
(15)
where all particles $`A_\alpha (x)A_\alpha (x^+,\underset{¯}{x})=\mathrm{\Psi }_{\alpha _1}(x^+,\underset{¯}{x_1})\mathrm{\Psi }_{\alpha _2}(x^+,\underset{¯}{x_2})\mathrm{\Psi }_{\alpha _3}(x^+,\underset{¯}{x_3})\mathrm{}`$ are all taken at the same light front time $`x^+`$ and $`\underset{¯}{x}=(x_+,\stackrel{}{x}_{})`$. The upper (lower) sign stands for fermion (boson) type clusters. Because of the global light-cone time introduced, the dynamical equation for a cluster is equivalent to a Dyson equation with a complicated mass operator that contains an instantaneous part and a memory (or retardation) part. For the time being we neglect the memory term. This is equivalent to a mean field approximation for clusters leading to Faddeev-type three-body equations.

For a simple zero range interaction the $`t`$ matrix, Fig. 4, separates and is given by the propagator $`t(M_2)`$, i.e.

$$t(M_2)=\left(i\lambda ^1B(M_2)\right)^1.$$
(16)
The expression for $`B(M_2)`$ is represented by the loop diagram of Fig. 4 and, in the rest system of the two-body system $`P^\mu =(M_2,M_2,\stackrel{}{0}_{})`$, given by

$$B(M_2)=\frac{i}{(2\pi )^3}\frac{dxd^2k_{}}{x(1x)}\frac{1f(x,\stackrel{}{k}_{}^2)f(1x,\stackrel{}{k}_{}^2)}{M_2^2M_{20}^2},$$
(17)
where $`M_{20}^2=(\stackrel{}{k}_{}^2+m^2)/x(1x)`$ and $`ff^{}`$ given in (6) with $`x=k^+/P_2^+`$ where $`\overline{P}_2=\overline{k}_1+\overline{k}_2`$. For a fermi system there are two important effects occurring due to the blocking factors of (17). One is the dissociation limit (Mott effect) where $`M_2(T_d,\mu _d)=2m(T_d,\mu _d)`$. Above a certain temperature and density no bound states can be formed. The second effect is related to the appearance of a bose pole in the $`t`$ matrix. This happens for $`M_2(T_c,\mu _c)=2\mu _c`$ and defines the critical temperature below which the system becomes unstable and forms a new vacuum consisting of Cooper pairs or a condensate. This is related to superconductivity or superfluidity. In this case for $`M_2^2M_{02}^2`$, we get

$$f(x,\stackrel{}{k}_{}^2)|_{M_{02}^2M_2^2=4\mu ^2}=f(1x,\stackrel{}{k}_{}^2)|_{M_{02}^2M_2^2=4\mu ^2}=\frac{1}{2},$$
(18)
i.e. both nominator and denominator of (17) are zero.

The three-particle case is driven by the Fadeev-type in medium equation

$$\mathrm{\Gamma }(y,\stackrel{}{q}_{})=\frac{i}{(2\pi )^3}t(M_2)\frac{dxd^2k_{}}{x(1yx)}\frac{1f(x,\stackrel{}{k}_{}^2)f(1xy,(\stackrel{}{k}+\stackrel{}{q})_{}^2)}{M_3^2M_{03}^2}\mathrm{\Gamma }(x,\stackrel{}{k}_{}),$$
(19)
where we have introduced vertex functions $`\mathrm{\Gamma }`$ and $`t(M_2)`$ given before, and an invariant cut-off $`M_{30}^2<\mathrm{\Lambda }^2`$. Here the mass of the virtual three-particle state (in the rest system $`P^\mu =(M_3,M_3,\stackrel{}{0}_{})`$ is

$$M_{03}^2=\frac{\stackrel{}{k}_{}^2+m^2}{x}+\frac{\stackrel{}{q}_{}^2+m^2}{y}+\frac{(\stackrel{}{k}+\stackrel{}{q})_{}^2+m^2}{1xy},$$
(20)
which is the sum of the on-shell minus-components of the three particles. The nucleon scale is introduced by setting $`M_3=938`$ MeV. The isolated quark mass used in these calculations is $`m=386`$ MeV.

Fig. 6 shows a shaded area that reflect the region where the transition from baryons to quarks (or quark diquarks) occur. The area is defined by use of different regularization masses. The chiral phase transition given before is indicated by the solid line. Fig. 6 shows the possible transition of quark matter to a superconducting phase.

## 5 Conclusion and Outlook

We have given a relativistic formulation of field theory at finite temperatures and densities utilizing the light front form. The proper partition operator (and the statistical operator) have been given for the grand canonical ensemble. The special case of a canonical ensemble is given for $`\mu =0`$. The resulting Fermi function depends on transverse and also on the $`k^+`$ momentum components. The $`k^+`$ components emerge in a natural way in a covariant approach. As an application we have revisited the NJL low energy model of QCD. We reproduce the phenomenology of the NJL model, in particular the gap-equation and the chiral phase transition. We have further given consistent relativistic three-quark equations valid in a dense medium of finite temperature. We find that the dissociation transition and the critical temperature for the color superconductivity agree qualitatively with results expected from other sources. However, the latter results are by no means final. We have shown that it is possible to write down meaningful consistent equations to solve the relativistic in-medium problem on the light front. The next steps would be to use the NJL model all through to give a consistent picture for the $`q\overline{q}`$ and the $`qq`$ channel. Further insight into this just emerging possibilities of treating relativistic many-particle systems on the light front might be provided by other theories, like 1+1 QCD, the Yukawa model, and finally real QCD.

#### Acknowledgment:

I gratefully acknowledge the fruitful collaboration with S. Mattiello, T. Frederico, and H.J. Weber, who have substantially contributed to the work I had the pleasure to present in this talk. I also thank S. Brodsky for his interest and discussion on this new approach. In particular, I would like to thank the organizers of the workshop for a fruitful and well organized meeting. Work supported by Deutsche Forschungsgemeinschaft, grant BE 1092/10.


format: text

text:
# Morphology of Rising Hydrodynamic and Magneto-hydrodynamic Bubbles from Numerical Simulations

## 1 INTRODUCTION

Cooling flows have been known for some time to exist at the centers of some clusters of galaxies (see Fabian 1994 for a review). The X-ray emissivity of the intracluster medium (ICM) increases with the square of the plasma density, and as the ICM is very centrally concentrated, radiative energy losses are greatest at the center of a cluster. In rich clusters these losses over the life of the cluster can be significant enough, in the absence of thermal conduction, to cause the plasma to contract gradually within the cluster potential and flow inward with velocities of order a few tens of km s<sup>-1</sup>. It has been understood that — again, in the absence of conduction – the cooling gas must become thermally unstable to the formation of a multiphase medium below about 0.1 $`\mathrm{keV}`$, but the fate of this cold gas – as stars, cold X-ray-absorbing clouds, or even brown dwarfs – has been a mystery. The possibility that some energetic feedback mechanism could shut off the cooling flows has been subject to debate, as well as the role played by thermal conduction and magnetic fields.

Recent high-resolution X-ray observations by the Chandra and XMM-Newton satellites have dramatically changed this picture (Böhringer et al., 2002). First, no evidence from spatially resolved spectroscopy has been observed for the presence of gas colder than about 1–2 $`\mathrm{keV}`$ (e.g., Schmidt et al. 2001; Peterson et al. 2001; Kaastra et al. 2001; Tamura et al. 2001; Matsushita et al. 2002). Second, high-resolution imaging has provided evidence for large-scale motions that can heat the ICM (e.g., McNamara et al. 2000; McNamara et al. 2001; Forman et al. 2002; Blanton et al. 2003). For years it has been known that many of the same clusters that harbor cooling flows also contain large central galaxies with active nuclei. However, the extent to which these active galactic nuclei (AGN) influence the dynamical state of the cooling ICM has only become apparent through the new X-ray observations. These observations show that AGN produce massive outflows of magnetized plasma that displace the cooling gas. The ‘bubbles’ thus produced are known to be magnetized because radio observations show regions of synchrotron emission coinciding with the regions of low X-ray emission, and the polarization of this radiation shows Faraday rotation effects consistent with dynamically important magnetic fields (e.g., Allen et al. 2001; Nulsen et al. 2002). The emerging picture is that bubbles represent the late stages of propagation of the magnetized, relativistic jets produced by AGN into the ICM, after they have slowed and reached approximate pressure equilibrium with the ICM (Reynolds et al., 2002).

The influence of these magnetized bubbles on the cooling ICM is still poorly understood. How efficiently does the bubble plasma mix with the ICM? Does it heat the cooling flow sufficiently to avoid the formation of multiphase gas and the possible formation of stars? If so, how – by radiative heating, magnetic reconnection, or some other process? What fraction of the total energy budget of a cluster is contributed by magnetic fields and cosmic rays? Could these active regions be sites for the acceleration of the ultra-high-energy cosmic rays seen in some air-shower experiments?

Owing to the geometrical and physical complexity of the AGN-cooling flow environment, numerical simulations are the best theoretical tools available for addressing these questions. While bubbles in liquids are well-studied (see, e.g., Magnaudet and Eames 2000; Harper 1972), bubbles in cluster cooling flows should differ from bubbles in liquids in several important ways. In particular, molecular forces present in liquids (e.g., surface tension, viscosity) will not play a significant role, and magnetic fields may be important. These differences will affect the dynamics and stability properties of bubbles.

Several numerical studies have been performed recently of rising bubbles in the context of cluster cooling flows. Churazov et al. (2001) used 2-D hydrodynamic simulations in a cooling-flow model atmosphere to gain insight into the bubbles observed in M87. Brüggen and Kaiser (2002, 2001) studied 2-D spherical and elliptical magnetohydrodynamic (MHD) bubbles rising in a hydrostatic background medium with an isothermal $`\beta `$-model density profile (Cavaliere and Fusco-Femiano, 1976),

$$\rho (r)[1+(r/r_c)^2]^{3\beta /2},$$
(1)
where $`r_c200`$ $`\mathrm{keV}`$ is the core radius and $`\beta `$ was taken to be 0.5. Brüggen et al. (2002) performed 3-D hydrodynamic bubble simulations with continuous energy injection. More recently, Basson and Alexander (2003) have performed 3-D hydrodynamic simulations of both the active and inactive phases of AGN jets propagating into $`\beta `$-model density profiles, including radiative cooling. In addition, several papers have addressed the detectability in radio and X-rays of cluster bubbles, including Churazov et al. (2001), Brüggen and Kaiser (2001), Ensslin and Heinz (2002), and Soker et al. (2002). Finally, bubbles in stratified atmospheres have been studied in the context of contact binary stars by Brandenburg and Hazlehurst (2001).

In this paper, we study the general question of the behavior of a buoyant gas bubble rising in a denser gas, investigating the relative importance of the effects of geometry, stratification, density contrast, and thermal and magnetic pressure on the stability of the bubble. Most of the previous studies (both 2-D and 3-D) have used finite-difference methods to evolve the gas on relatively low-resolution grids, up to about 400 zones on a side. Brüggen (2003) and Brüggen and Kaiser (2002) used the Flash code (Fryxell et al., 2000), which has an adaptive mesh (MacNeice et al., 2000) to achieve an effective resolution of 4,000$`\times `$2,000 zones for purely hydrodynamic calculations. We also employ Flash, but in addition to hydrodynamic bubbles we study the effects of magnetic fields with different spatial configurations, exploiting a new MHD module we have developed for Flash. In addition, because we do not treat the bubbles as self-gravitating, we are able to use for the hydrodynamic case a version of the piecewise parabolic method (PPM) that has been modified to handle nearly hydrostatic flows (Zingale et al., 2002). Our initial model is deliberately simplified as a first step in a research program to better understand the dynamics of bubbles in clusters; thus we postpone detailed observational comparisons to a later paper.

The paper is organized as follows: in § 2 we describe the construction of our initial models and the numerical methods used. In § 3 we discuss the results of our purely hydrodynamic calculations. In § 4 we discuss results of simulations of hot bubbles rising in background magnetic fields, and in § 5 we consider magnetized bubbles initially supported partly by force-free magnetic fields. In § 6 we summarize our conclusions.

## 2 INITIAL MODELS

We consider as initial models an under-dense circular bubble in an isothermal, stratified atmosphere with an imposed constant gravitational acceleration in the $`y`$-direction. An ideal equation of state is used, with a ratio of specific heats $`\gamma =5/3`$. The sound speed was chosen, through the temperature and mean molecular weight, so that our standard simulation domain would span approximately three pressure scale heights. For most of our simulations, we simulate 2-D boxes in planar geometry. We vary the stratification by varying the gravitational acceleration. We also vary the size of the bubble.

We perform simulations in which the density contrast between the bubble and the background is caused by a higher temperature in the bubble in horizontal pressure equilibrium with its surroundings (the ‘hot bubble’ case) and also in which the pressure support for the under-dense bubble comes from an azimuthal magnetic field (the ‘magnetized bubble’ case). In the ‘hot bubble’ case, we perform simulations with and without a background magnetic field. The bubble is in horizontal pressure equilibrium at every point inside the bubble; since the bubbles considered here are significant fractions of a pressure scale height, the bubble’s pressure structure is itself somewhat stratified.

The physics in the simulations presented is all scale-free. For numerical reasons, we choose scales for the simulations so that numbers simulated are all near unity. Thus, our simulations were performed with a domain size of $`40\times 50`$ cm, and the bubble radius of $`3`$ cm. The base density and pressure are taken to be of order unity ($`\rho _0=1`$, $`p_0.64`$) and gravitational acceleration is $`0.04\mathrm{cm}\mathrm{s}^2`$. For comparison to the astrophysical systems of interest, however, we rescale the results when plotting figures. We rescale by scaling the length scales by a factor $`L`$, the gravitational acceleration by $`G`$, and the density by $`R`$. In that case, the hydrodynamic quantities scale as

$`x`$ $`=`$ $`Lx^{}`$ (2)
$`t`$ $`=`$ $`\sqrt{{\displaystyle \frac{L}{G}}}t^{}`$ (3)
$`v`$ $`=`$ $`\sqrt{LG}v^{}`$ (4)
$`g`$ $`=`$ $`Gg^{}`$ (5)
$`\rho `$ $`=`$ $`R\rho ^{}`$ (6)
$`p`$ $`=`$ $`\left(RGL\right)p^{},`$ (7)
where the primed quantities represent the simulated values, and unprimed represent the scaled values. Here, $`x`$ is a length, $`t`$ time, $`v`$ velocity, $`g`$ the gravitational acceleration, $`\rho `$ density, and $`p`$ pressure. The values presented here are all scaled by factors to be roughly consistent with the model of Abel 2052 used in Soker et al. (2002): $`L=(3\mathrm{k}\mathrm{p}\mathrm{c})/(1\mathrm{cm})`$, $`G=(7\times 10^9\mathrm{cm}\mathrm{s}^2)/(.04\mathrm{cm}\mathrm{s}^2)`$, and $`R=(4\times 10^{26}\mathrm{g}\mathrm{cm}^3)/(1\mathrm{g}\mathrm{cm}^3)`$. With these scalings, the physical background temperature is approximately $`1.6\times 10^7`$ K, and the sound speed is approximately $`0.41\mathrm{kpc}\mathrm{Myr}^1`$.

We perform our simulations with the Flash code (Fryxell et al., 2000). The Flash code is an adaptive-mesh reactive hydrodynamics code that has both hydrodynamic and magneto-hydrodynamic solvers. The Flash code’s main hydrodynamic solver has undergone a rigorous validation and verification process (Calder et al., 2002). The method we use for performing our magnetic calculations is described in Powell et al. (1999). The two solvers use similar methods, although the MHD solver uses lower-order reconstructions than the PPM solver, meaning for the same number of grid points, it is ‘lower resolution’. We use the techniques described in Zingale et al. (2002) to maintain hydrostatic equilibrium; in particular, we use the initialization procedure described therein and the boundary conditions. The modification to PPM to provide more accurate stable solutions was not found to significantly effect the very dynamic solutions described here, and was not used with either solver. Unless otherwise noted, the simulations described here were performed with an effective grid size of $`512\times 640`$ — with an unrefined grid of $`128\times 160`$ and two further levels of mesh refinement applied where there are large second derivatives of density or pressure, as described in Fryxell et al. (2000).

We use reflecting boundaries on the left and right, and the hydrostatic boundaries described in Zingale et al. (2002) with ‘outflow’ velocities. In simulations with magnetic fields, the magnetic boundary conditions are reflecting on the left and right, and zero-gradient at the top and bottom.

Because the gas is miscible, there is no surface tension. We can estimate the magnitudes of other unmodelled effects such as viscosity, conduction, and cooling effects (Spitzer, 1962). If the gas is a completely ionized hydrogen gas, scaling by the typical values relevant to this problem allows us to write

$$\mathrm{Re}=UL/\nu =400T_7^{5/2}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{32}\right)\left(\frac{n_i}{10^2\mathrm{cm}^3}\right)\left(\frac{U}{0.1\mathrm{kpc}\mathrm{Myr}^1}\right)\left(\frac{L}{18\mathrm{kpc}}\right),$$
(8)
where the Reynolds number, Re, expresses the relative importance of inertial to viscous forces, $`T_7`$ is the temperature in units of $`10^7K`$, $`\mathrm{ln}\mathrm{\Lambda }`$ is the Coulomb logarithm, measuring the range of length scales over which collisions are important, $`n_i`$ is the number density of ions, and $`\nu `$ is the kinematic viscosity. The length scale used, $`L=18\mathrm{kpc}`$, is a typical size of the bubbles we study, and the velocity scale, $`U=0.1\mathrm{kpc}\mathrm{Myr}^1`$ is, as we will see, a typical bubble velocity.

Considerable effort has gone into investigating the relatively low numerical dissipation in PPM and other shock-capturing schemes as well as the issue of convergence of solutions of these methods (Sytine et al., 2000; Garnier et al., 1999; Porter and Woodward, 1994). While there is not complete agreement about these results and their interpretation, it is apparent that the intrinsic viscosity of methods such as PPM, even at the highest attainable resolutions, prevents simulation of high Reynolds number flows. Because the Reynolds number of the flow is certainly larger than can be obtained in these simulations (the 18 kpc bubble is resolved with only 77 points), the omission of an explicit viscosity is appropriate.

Similarly, one can compute a Peclet number, Pe, which expresses the relative importance of heat transport by convection to that of thermal diffusion:

$$\mathrm{Pe}=UL/D_{\mathrm{th}}=9T_7^{5/2}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{32}\right)\left(\frac{n_i}{10^2\mathrm{cm}^3}\right)\left(\frac{U}{0.1\mathrm{kpc}\mathrm{Myr}^1}\right)\left(\frac{L}{18\mathrm{kpc}}\right).$$
(9)
Here $`D_{\mathrm{th}}`$ is the thermal diffusivity, which is approximately

$$D_{\mathrm{th}}=0.209T_7^{5/2}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{32}\right)^1\left(\frac{n_i}{10^2\mathrm{cm}^3}\right)^1\mathrm{kpc}^2/\mathrm{Myr}.$$
(10)
Here we see that thermal diffusion may be significant; it will certainly be significant at the interface between the bubble and the ambient medium. We do not include thermal conductivity in the simulations, and we discuss heat diffusion out of the bubble in §6.

Another neglected effect is radiative cooling of any hot bubbles. If the radiative cooling is dominated by free-free emission, then we will have

$$\frac{dT}{dt}=2300T_7^{1/2}\left(\frac{\overline{g}_{\mathrm{ff}}}{1.103}\right)\left(\frac{n_i}{10^2\mathrm{cm}^3}\right)K\mathrm{Myr}^1,$$
(11)
where $`dT/dt`$ is the change in temperature per unit time. In this case, gas at 1 $`\mathrm{keV}`$ $`10^7K`$ will take time scales of order $`5000\mathrm{M}\mathrm{y}\mathrm{r}`$ to radiate away significant amounts of its internal energy, meaning that such radiative cooling can be neglected for the current simulations. However, using more realistic cooling rates such as in Tucker and Rosner (1983), we have

$$\frac{dT}{dt}=48000T_7^{1/2}\left(\frac{n_i}{10^2\mathrm{cm}^3}\right)K\mathrm{Myr}^1,$$
(12)
so that the cooling time scale at $`10^7K`$ is 210 $`\mathrm{Myr}`$, meaning that cooling would be starting to play a role during these simulations were it included.

## 3 HYDRODYNAMIC BUBBLES

We first examine our reference case – a simulation with $`g=7\times 10^9\mathrm{cm}\mathrm{s}^2`$ and a density at the bottom of the simulation domain of $`\rho =4\times 10^{26}\mathrm{g}\mathrm{cm}^3`$. With these parameters, our box (150 $`\mathrm{kpc}`$ tall) contains $`3.1`$ pressure scale heights. The density contrast between the bubble and the maximum density in the box is 100:1, and there are no magnetic fields. The bubble radius is $`r=9\mathrm{kpc}`$, making the whole bubble a little more than one-third of a pressure scale height in size. This simulation is in 2D planar coordinates. The evolution of density is shown in Fig. 1.

As we have scaled our results to Abell 2052, our bubble size (in units of pressure scale heights) is consistent with that of bubbles observed in that cluster (approximately 0.35). In other clusters, bubbles are seen with sizes of approximately 0.1–0.9 pressure scale heights: examples are 0.08 in Abell 478(Sun et al., 2003); 0.25–0.37 in Abell 2597(McNamara et al., 2001); 0.64–0.9 in Hydra A(McNamara et al., 2000).

In the absence of magnetic effects, nothing prevents the complete disruption of the bubble. The bubble’s motion itself generates vortical motion in the surrounding fluid of the same size as the bubble and with rise velocity of the bubble, as shown in Fig. 2. These motions will then be able to rip the bubble apart on order of a bubble rise time, $`t_R=\sqrt{4R/|g|}`$, with $`R`$ being the bubble’s radius, and $`g`$ being the local gravitational acceleration. These effects have been studied previously (*e.g.*, Harper 1972 and Layzer (1955)).

The bubble will also be susceptible to growth of the Rayleigh-Taylor instability at the top. Since we are dealing with large density contrasts, the Rayleigh-Taylor growth time (from linear theory) will be $`t_{RT}=1/\sqrt{kg}=\sqrt{R/(2\pi gn)}`$ for the mode that is $`1/n`$ the size of the bubble; however, modes of size on order the bubble radius, or small modes near the sides of the bubble will be suppressed because of the geometry of the bubble. Comparing time scales, the Rayleigh-Taylor instability will disrupt the bubble on small scales near the center of the bubble quickly, but the disruption of the entire bubble will be caused by induced vortical motions in the fluid. The shear along the side of the bubble will also be unstable to the Kelvin-Helmholtz instability; this will generate secondary instabilities as the evolution progresses.

The large scale motions of the fluid can be seen in Fig. 2 at time $`t70\mathrm{Myr}`$. Two large rolls extend out, and their convergence under the bubble results in an up-welling and compression of material, forming a slight high-density tail in the wake of the bubble.

The position of the bubble top as a function of time is shown in Fig. 3. This position was calculated by horizontally averaging density throughout the domain and finding the first position from the top of the domain of a density jump; this local drop in density represents the presence of an under-dense fluid. Because our underdensities were so large, the position measured this way was insensitive to the precise value of the threshold we chose to mark the density jump. This measure allows us to track the top of a ‘bubble’ even when no coherent bubble exists. The bubble slows down over time primarily due to the entrainment of material from the environment as the bubble is disrupted.

As resolution varies, the small-scale structure varies, but the overall behavior remains the same, as is shown in Fig. 4 and Fig. 5. In these simulations, there is no modeled dissipation mechanism such as viscosity or thermal or material diffusion. In this case, only the numerical dissipation on small scales limits the amount of small scale structure. The small scale structure affects entrainment of ambient (denser) material, so that as resolution increases and numerical dissipation decreases, late-time heights change, but the early behavior is nearly identical.

A typical rise velocity for the still largely intact bubble can be seen from looking at the early-time slope in Fig. 3; it is approximately $`0.1\mathrm{kpc}\mathrm{Myr}^1`$. That this is near the sonic velocity ($`0.41\mathrm{kpc}\mathrm{Myr}^1`$) is not a coincidence. The density contrasts we consider here are large enough that the buoyant acceleration felt by the bubble is approximately $`g`$; since in an isothermal atmosphere the scale height is $`H=c_s^2/(\gamma g)`$, where $`c_s`$ is the sound speed, the velocity that the bubble picks up as it rises a significant fraction of a scale height is on the order $`c_s`$. The disruption of the bubble and entrainment of denser material before a scale height is reached prevents the bubble from accelerating past the sound speed.

### 3.1 Effects of Density Contrast

We ran simulations with density contrasts between the bubble and the maximum density in the box of 100:1 and 10:1. While there are some differences between the resulting structures (compare Figs.1 and 6), caused partly by the different growth of the Rayleigh-Taylor instability at the top of the bubble, the underlying dynamics are broadly similar. One way of understanding this is that in both cases, the bubble densities are very much less than the surrounding environment, so that it is the momentum of the fluid in the environment, not the bubble, which dominates the dynamics. Quantitatively, the Atwood number, $`A=(\rho _1\rho _2)/(\rho _1+\rho _2)`$ is $`0.82`$ in the 10:1 case, and $`0.98`$ in the 100:1 case – so while the bubble density changes by a factor of 10, the Atwood number, which controls both the buoyant force and the evolution of the Rayleigh-Taylor instability, changes only by 20%.

In particular, looking at the position of the top of the bubble versus time (in Figs. 3 and 7) we see a very similar evolution after the initial acceleration phase. Perhaps surprisingly, the bubble with the larger density (the 10:1 case) actually rose slightly further during the course of the simulation. In the late-time evolution, it is the entrainment of higher-density fluid which determines the buoyancy, and thus the acceleration, of the bubble. The bubble in the 10:1 case entrained less ambient fluid, thus its buoyancy ended up being greater than in the nominally more buoyant 100:1 case.

### 3.2 Effects of Geometry

In the purely hydrodynamic case it is possible to perform the same simulations described above in 2-D axisymmetric cylindrical coordinates, making the bubble a sphere instead of a cylinder. In this more physically relevant case, the bubble deforms into tori (Fig. 8), rather than the wing-like shape found in planar coordinates, as the shedded vorticity from the bubble rising becomes a vortex ring rather than two parallel lines of vorticity.

Because in this geometry there is less volume near the center of the bubble than there is towards the edges, the vortices that form can advect essentially all of the material from within the bubble into the vortex structure, so that the bubble becomes a rising vortex ring which then undergoes secondary instabilities. Also because of the geometry, the effect of the initial Rayleigh-Taylor instability in disrupting the bubble is diminished, although it does not vanish; the material that is disrupted by the Rayleigh-Taylor instability forms a smaller, secondary vortex ring that leads the larger ring.

This simulation was performed with several density differences which all exhibited similar behavior. Fig. 8 shows the evolution of one of these cylindrical cases with a density contrast of 100:1 and a 9 $`\mathrm{kpc}`$ radius bubble. Fig. 9 shows the final frame of this simulation along with the corresponding frame of a simulation performed with a 10:1 density contrast. As with the planar case, although there are differences in the small-scale structure, the medium- and large-scale behavior is essentially unchanged.

### 3.3 Effects of Bubble Size

Simulations were also performed with a larger and a smaller bubble in the planar case. In the case of the larger $`r=15\mathrm{kpc}`$ bubble of Fig. 10, the Raleigh-Taylor instability plays a much larger role in the evolution of the bubble than in the previous case, as there are more modes that can grow on the larger top of the bubble. Conversely, the $`r=10\mathrm{kpc}`$ bubble (Fig. 11) tends to remain flat across the top, and is primarily disrupted by the induced vortical motions.

### 3.4 Effects of Stratification

The isothermal atmosphere surrounding the bubble has an exponential density profile provided by a constant gravitational field in the negative $`y`$-direction. Changing the stratification, through the gravitational acceleration, has two major effects. Immediately, it changes the time scale on which the dynamics occurs, as shown in the scaling relations in §2. The stratification also tends to elongate vortical motions in the atmosphere, due to lower densities at higher positions.

In Fig. 12 we show the evolution of a bubble moving through a less stratified domain, with the gravity set to be one quarter of its value in our reference case. In Fig. 13, we compare the final frame of the evolution with the frame from the corresponding time (taking into account the change of the buoyant time scale) of our reference case.

Some of the differences between the two simulations can be understood in terms of the height of the vortical motions induced by the rising bubble. The large-scale motions, if they are to conserve momentum, must be elongated in atmospheres where the density is stratified, as there is less mass at higher positions in the atmosphere. Fig. 14 shows total velocity contours on top of a zoom-in of density for the same simulations as in Fig. 13, but at earlier times. The vortical motions extend noticeably higher in the more stratified case, meaning that in this case the bubble will be disrupted over longer length scales.

## 4 MAGNETIC FIELDS: HOT BUBBLE IN A BACKGROUND FIELD

From purely hydrodynamic bubbles, we move on to examining hot bubbles rising in a fluid with an initially constant background field running through the entire domain. In the rest of this paper, we describe solutions performed with the MHD solver in the Flash code. The rest of the code, including the boundary conditions and initialization, remained unchanged; the MHD solver was simply compiled into the simulation code in place of the PPM solver. This MHD solver is based on the code described in Powell et al. (1999), with two methods of ‘$`𝐁`$-cleaning’ – a diffusive method due to Marder (1987), and a projection method due to Brackbill and Barnes (1980). We used the diffusive mechanism for the results described here due to the more modest computational cost of this method and the simple magnetic dynamics in these simulations. Our implementation of this MHD solver works only in Cartesian coordinates, so we consider only the 2-D planar cases.

Since we are changing the numerical solver as well as adding magnetic fields to the initial conditions, it is important to ensure that we get the same results in the purely hydrodynamic case with both solvers. Calculating our reference case without magnetic fields with both solvers, we find noticeable, but understandable, differences in the final results. Shown in Fig. 15 is our reference calculation at the final time calculated with our MHD solver and the same result with the PPM solver. The MHD result looks something like the lower-resolution versions of the PPM solution (Fig. 4), which can be understood to be a consequence of the the lower spatial order of accuracy of the MHD solver. Fig. 15 also shows the results of using the purely hydrodynamic solver with reconstruction functions more similar to those in the MHD solver, which improves the agreement between the simulations. Understanding the differences between the results of the different solvers gives us some confidence that the results obtained with the magnetic field solver can be meaningfully compared with those calculated using the purely hydrodynamic solver.

In the results that follow, a density ratio of 10:1 was used. This avoids numerical issues which can arise if there is too large a jump in Alfvén velocity across a small region (eg., between the bubble and the ambient material.)

### 4.1 Horizontal Field

#### 4.1.1 z-field

The simplest case to consider in 2-D is an initial magnetic field in the $`\widehat{𝐳}`$-direction, out of the plane of the simulation domain. In this case, the magnetic field lines can slip around the cylinder as it rises, enhancing the initial production of large-scale vorticity around the bubble. The magnetic field’s contribution to the dynamics in this case is simply as an additional, constant, pressure term $`B^2/(8\pi )`$, reducing the effective relative pressure stratification $`p/p`$ by a factor of $`\beta _0/(\beta _0+1)`$ while leaving the density structure unchanged. Here $`\beta _0`$ is the plasma $`\beta `$ parameter (gas pressure divided by magnetic pressure) at the base of the atmosphere. This effective pressure scale height will be one characteristic length scale for eddies, so that decreasing $`\beta `$ will, up to a limit, decrease the characteristic vortical motion size without the elongation that increasing the density stratification would produce.

The evolution of a bubble in a very weak $`\widehat{𝐳}`$-field is shown in Fig. 16. In Fig. 17 is shown the final frame for simulations with $`\beta _0`$ of 462, 0.19, 0.046, and 0.012. As $`\beta _0`$ decreases, the vortical structures are seen to ‘tighten’.

At late times, the high-density wake that the bubble leaves begins to contain low-density structures; one explanation for this is that the large-scale motions induced by the flow sweep up magnetic field lines, resulting in a higher magnetic field density in the wake of the bubble; the resulting increase in magnetic pressure means a lower mass density under these constant-pressure conditions, leaving a different tail in this case than in the purely hydrodynamic case.

#### 4.1.2 x-field

As opposed to the effects of a field in the $`\widehat{𝐳}`$-direction, even an extremely weak field in the $`\widehat{𝐱}`$ direction in this geometry can strongly contain the rising bubble in our planar coordinates. We see from Chandrasekhar (1981) that the magnetic field suppresses the Kelvin-Helmholtz instability completely unless the relative velocity between the bubble and the ambient medium exceeds the root-mean square Alfvén speed in the two media; or, as written in Vikhlinin et al. (2001),

$$4\pi \frac{\rho _0\rho _b}{\rho _0+\rho _b}v^2>B_0^2+B_b^2.$$
(13)
At the top of the bubble, near the stagnation point, the relative velocity is nearly zero and of course the instability is suppressed. Even near the sides, however, because the Alfvén speed in the bubble is so high (due to the low density), the shear velocity would have to be many times the ambient Alfvén speed to for Kelvin-Helmholtz instabilities to occur.

Initially, since $`B_0=B_b`$, this reduces to $`\beta _0>(4/\gamma )(\rho _0/\rho _b+1)^2`$ where $`\gamma `$ is the ratio of specific heats for the gas, and $``$ is the flow Mach number. Since we have a density ratio of 10:1, a Mach number of 0.25, and a 5/3-law ideal gas, this suggests that for $`\beta _0>420`$, the bubble should be unstable to shear. However, we show the evolution of a bubble with $`\beta _0=462`$ in Fig. 18, and there is clearly no shear instability present. In fact, as the bubble rises, not only does shear occur but the magnetic field lines are wrapped around the rising bubble, increasing the field density. Further, the shear at the bubble interface is overestimated by the net bubble velocity upwards. Applying Eq. 13 to the data of Fig. 18, we find that the shear everywhere is several orders of magnitude below the threshold for instability. This is tabulated in Table 1.

In this case of a horizontal magnetic field, the wake left behind by the bubble is an unstructured high-density entrained tail; this will be true in the spherical bubble case, as it is caused by the stagnation line caused by the induced vortical motions of the rising bubble.

The drastically different effects of the two horizontal orientations of the uniform magnetic field is partly an artifact of the geometry of the simulation. For a spherical bubble, rather than a cylindrical bubble, a uniform horizontal field would act on the bubble asymmetrically – as a tension in one plane, and as a pressure component in the other. In particular, the strong containing effect of even a weak $`\widehat{𝐱}`$-field is a geometrical artifact of these simulations.

### 4.2 Vertical Field

A vertical field has a more direct dynamical effect on the rising bubble. Adding a weak vertical field does not prevent the bubble from breaking up, but somewhat suppresses the cascade of vortical motions at late times by suppressing horizontal motions as compared to vertical motions (Fig. 19). If we increase the strength of the magnetic field, however, the suppression of horizontal motions can be so strong that the bubble is deformed into a chevron and travels upwards intact (Fig. 20.) A sufficiently strong vertical magnetic field can completely contain the bubble, although it is difficult to imagine a cluster having a magnetic field of this magnitude over the entire size of the observed bubbles.

## 5 MAGNETIC FIELDS: MAGNETIZED BUBBLE IN A BACKGROUND FIELD

In the previous sections, under-dense bubbles were supported by thermal pressure. This isotropic thermal pressure can support the bubble against collapse, but not against being torn apart by the fluid motions induced by its own rising. The ambient magnetic fields could only prevent this if the fields were strong enough (or oriented fortuitously) to suppress those large-scale fluid motions. We now consider bubbles which are supported instead by magnetic fields, which provide an anisotropic ‘pressure’ support that also acts as a tension; this tension can be enough to maintain the bubble even in the presence of such motions.

One can construct families of solutions to the MHD equations with the property that $`\stackrel{}{B}`$ is force-free (i.e. $`\stackrel{}{J}\times \stackrel{}{B}=0`$) everywhere in the domain except at the boundary between the under-dense bubble and the ambient medium, where the Lorentz force balances the gas pressure jump. We choose a field used in Cargill and Chen (1996) for flux tubes:

$`B_\theta `$ $`=`$ $`\{\begin{array}{cc}B_b\frac{r}{a_0}\hfill & ra_0\text{,}\hfill \\ 0\hfill & r>a_0\text{.}\hfill \end{array}`$ (14)
$`B_z`$ $`=`$ $`\{\begin{array}{cc}B_b\sqrt{42\left(\frac{r}{a_0}\right)^2}\hfill & ra_0\text{,}\hfill \\ B_{z0}\hfill & r>a_0\text{.}\hfill \end{array}`$ (15)
$`B_b`$ $`=`$ $`\sqrt{{\displaystyle \frac{8\pi }{3}}p_0\left(1{\displaystyle \frac{\rho _b}{\rho _0}}\right)+B_{z0}^2}.`$ (16)
Here $`p_0`$ is the ambient gas pressure, $`\rho _0`$ is the ambient gas density, $`\rho _b`$ is the gas density inside the bubble, $`B_{z0}`$ is the ambient magnetic field in the $`\widehat{𝐳}`$ direction, and $`a_0`$ is the radius of the bubble. Thus we have a uniform magnetic field in the $`\widehat{𝐳}`$ direction in the ambient medium, as discussed in our previous section, but the bubble is also supported by an azimuthal field and an enhanced $`\widehat{𝐳}`$field.

The hydrodynamic initial conditions are computed as in previous sections, except that in this case the underdense region is at the same temperature as the background material, and a constant Lorentz force at the bubble interface supports the bubble against the larger external ambient pressure. The evolution is shown in Fig. 21. In this case, the bubble maintains its form as it rises.

Because the ambient pressure changes by approximately 30% over the height of the bubble, while the bulk of the pressure support of the bubble is fixed by the magnetic field jump, a transient occurs over the course of a few sound-crossing times of the bubble (that is, a few times 44 $`\mathrm{Myr}`$) as the bubble equilibrates; the short-time evolution of gas pressure is is shown in Fig. 22. The magnetic field is enough to maintain the bubble’s shape until gas pressure and density inside the bubble is redistributed to maintain equilibrium; we see in Fig. 22 that the bubble expands overall, and gas slightly settles towards the bottom.

The bubble adiabatically expands as it rises, and eventually overshoots the height of neutral buoyancy, decelerates, and falls back again (Fig. 23). We also see again at late times the development of an under-dense wake as in the uniform $`\widehat{𝐳}`$-field case.

Also notable is that the bubble carries a ring of higher density material (Fig. 24). Such a high-density ring has also been noted in Chandra observations of bubbles. This collar of material appears to be “grabbed” by the azimuthal field at time t=0 when the entire bubble is surrounded by higher density material. As shown in Fig. 25, the high density collar surrounding the low-density bubble does not change in density significantly as the bubble moves upwards, suggesting that the bubble simply advects ambient material.

## 6 CONCLUSIONS

This work focused on examining the primary physical conditions that have to occur for the continued existence of the dark spots in galaxy clusters found by the Chandra X-ray observatory. We confirm that bubbles without supporting magnetic fields are torn apart by instabilities and vortical motions before they can move an entire bubble height. This is surprising, as ‘ghost’ bubbles — not radio bright, and thus presumably no longer powered — are observed (Mazzotta et al., 2002; Fabian et al., 2000; McNamara et al., 2000) which are a significant distance in bubble lengths away from the radio source which presumably formed them. The existence of these bubbles can be explained if they are supported by an internal magnetic field.

Furthering this argument is the importance of thermal diffusion absent any magnetic field effects. One can numerically integrate the diffusion equation using the thermal diffusivity quoted in §2 given by Spitzer (1962). On doing so, one finds that an $`r=9\mathrm{kpc}`$ bubble at a temperature of 100 $`\mathrm{keV}`$ in an ambient gas of 1 $`\mathrm{keV}`$ is largely diffused away after 1 $`\mathrm{Myr}`$; a 10 $`\mathrm{keV}`$ bubble survives for perhaps 150 $`\mathrm{Myr}`$. Indeed, if the Spitzer diffusivity is appropriate, there are constraints much tighter than this; the Chandra observations show a sharp edge to the bubble, and this edge would be blurred by thermal diffusion on time scales orders of magnitude shorter than those required to completely dissipate the bubble.

Thermal arguments alone, however, are not enough to require a magnetic field to support the bubble; a weak tangled magnetic field may reduce the conductivity without being strong enough to support the bubble, for instance, or one could require more sophisticated diffusivities in the presence of such large gradients than the Spitzer model. In addition, the synchrotron emission in the bubble may suggest that a 10-100 $`\mathrm{keV}`$ thermal gas is an insufficient model for the gas inside the bubble.

This research also suggests that the ring of brighter material surrounding these bubbles may be caused by magnetic diffusion of the field that maintains them. There is also numerical evidence suggesting a wake left behind the bubble as it moves. Searching the radio emissions for such a wake would be a good indicator as to whether or not these bubbles are moving.

Since magnetic fields may be necessary to keep the bubbles intact as they travel, future work should focus on MHD simulations, and in particular performing simulations with more physically meaningful geometry than 2D planar symmetry. As was seen in the hydrodynamic case, the difference between cylindrical and planar geometry significantly altered the morphology, if not the timescale, of the disruption of the bubble; the difference in geometries would be even more pronounced in the presence of magnetic fields, since the orientation of the field also plays a role. The hot bubble in magnetic field simulations generalize easily to 3D; the magnetized bubble simulation will be more complicated, as one can not write down an analytical magnetic field in 3D which would support a spherical bubble analogous to the cylindrical flux-tube presented here. Such a field, numerically generated, however, should also be able to support the bubble without disruption. Including the divergent geometry appropriate to the center of a cluster will also be a necessary step.

Support for this work was provided by DOE grant number B341495 to the ASCI/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago. K. Robinson was supported by the NSF REU program at the University of Chicago. LJD was supported by the Department of Energy Computational Science Graduate Fellowship Program of the Office of Scientific Computing and Office of Defense Programs in the Department of Energy under contract DE-FG02-97ER25308. MZ is supported by the Scientific Discovery through Advanced Computing (SciDAC) program of the DOE, grant number DE-FC02-01ER41176 to the Supernova Science Center/UCSC. All simulations were performed with Flash 2.2. The authors thank T. Emonet for very useful discussions, and the anonymous referee whose comments improved this work. Flash is freely available at http://flash.uchicago.edu/.


text:
# Energy dependence of current noise in superconducting/normal metal junctions

## Abstract

Interference of electronic waves undergoing Andreev reflection in diffusive conductors determines the energy profile of the conductance on the scale of the Thouless energy. A similar dependence exists in the current noise, but its behavior is known only in few limiting cases. We consider a metallic diffusive wire connected to a superconducting reservoir through an interface characterized by an arbitrary distribution of channel transparencies. Within the quasiclassical theory for current fluctuations we provide a general expression for the energy dependence of the current noise.

Interference of electronic waves in metallic disordered conductors is responsible for weak localization corrections to the conductance Weakloc . If these are neglected, the probability of transferring an electron through the diffusive medium is given by the sum of the modulus squared of the quantum probability amplitudes for crossing the sample along all possible paths. This probability is denoted as semiclassical, since quantum mechanics is necessary only for establishing the probability for following each path independently of the phases of the quantum amplitudes. In superconducting/normal metal hybrid structures, interference contributions are not corrections, they may actually dominate the above defined semiclassical result for temperatures and voltages smaller than the superconducting gap. This is seen experimentally as an energy dependence of the conductance on the scale of the Thouless energy. Indeed, the energy dependence comes from the small wavevector mismatch, linear in the energy of the excitations, between the electron and the Andreev reflected hole. This is responsible for the phase difference in the amplitudes for two different paths leading to interference. The effect is well known and explicit predictions and measurements exist for a number of systems VZK93 ; Hekking93 ; Beenakker97 .

Interference strongly affects the current noise too BlanterButtiker . The largest effects are expected in the tunneling limit, when the transparency of the barrier is small and its resistance is much larger than the resistance of the diffusive normal region. Then, the conductance has a strong non linear dependence at low bias (reflectionless tunneling) VZK93 ; Hekking93 . This is actually the case, but the zero-temperature noise (or shot noise) does not give any additional information on the system since it is simply proportional to the current, as shown numerically in a specific example in Ref. SH  and quite generally in Ref. PBH . The double tunnel barrier system has been considered in Ref. samuelsson . In the case of a diffusive metal wire in contact with a superconductor through an interface of conductance $`G_B`$ much larger than the wire conductance $`G_D`$, Belzig and Nazarov BN01  found that the differential shot noise, $`dS/dV`$, shows a reentrant behavior, as a function of the voltage bias, similar, but not identical, to the conductance one. (The extension of the Boltzman-Langevin approach to the coherent regime in Ref. houzet  neglects this difference.) In order to compare quantitatively with actual experiments Jehl ; KSP ; Lefloch  and to gain more insight in the interference phenomenon, it is necessary to obtain the energy dependence of noise in more general situations. The numerical method used in Ref. BN01  is, in principle, suitable to treat more general cases, notably the case when $`G_DG_B`$, but only if all channel transparencies, $`\{\mathrm{\Gamma }_n\}`$, that characterize the interface are small. When arbitrary transparencies are present one has to solve numerically an additional self-consistent equation private . This appears particularly heavy numerically if one is interested to treat a distribution of transparencies. We are not aware of results in this direction. In this paper we present an analytical solution for the diffusion-type differential equation for the noise within the theory of current fluctuations FCS  in the quasiclassical dirty limit BN01 . It allows to treat the general case of arbitrary values for $`\{\mathrm{\Gamma }_n\}`$ and $`G_B/G_D`$. We express the noise in terms of a function that satisfies a linear differential equation to be solved numerically once the channel distribution is given. We can thus isolate the energy-independent (semiclassical) contribution and the interference contribution to both the conductance and noise. We discuss their relation in the following.

We begin by stating the framework for the theory of current fluctuations FCS ; BN01 . It relies on the evaluation of the quasiclassical Green’s functions $`\stackrel{ˇ}{g}(x,\epsilon ,\chi )`$ in the Nambu$`(\widehat{})`$-Keldysh$`(\overline{})`$ space, at a given energy $`\epsilon `$ and counting field $`\chi `$. The counting field appears as a gauge transformation of the Green’s function in the normal reservoir

$$\stackrel{ˇ}{g}_N(\chi )=e^{\frac{i}{2}\chi \stackrel{ˇ}{\tau }_K}\stackrel{ˇ}{g}_N^0e^{\frac{i}{2}\chi \stackrel{ˇ}{\tau }_K},$$
(1)
where $`\stackrel{ˇ}{g}_N^0=\widehat{\tau }_3\overline{\sigma }_3+(f_{T0}+f_{L0}\widehat{\tau }_3)(\overline{\sigma }_1+i\overline{\sigma }_2)`$ is the standard Green’s function for a metallic reservoir, $`\overline{\sigma }_i,\widehat{\tau }_{j(i,j=1,2,3)}`$ are Pauli matrices, $`\stackrel{ˇ}{\tau }_K=\widehat{\tau }_3\overline{\sigma }_1`$, $`f_{L0}=1f_+f_{}`$, $`f_{T0}=f_{}f_+`$, $`f_\pm (\epsilon )=f(\epsilon \pm eV)`$, $`f`$ is the Fermi function at temperature $`T`$, and $`V`$ is the voltage bias. Given $`\stackrel{ˇ}{g}(x)`$, one can calculate the spectral matrix current in the wire

$$\stackrel{ˇ}{J}(x)=LG_D\stackrel{ˇ}{g}(x)_x\stackrel{ˇ}{g}(x),$$
(2)
with $`x`$ the coordinate along the wire of length $`L`$. The current and the zero-frequency noise are given respectively by BN01

$$I=J(\chi =0)\text{and}S=2ie\frac{J(\chi )}{\chi }|_{\chi =0},$$
(3)
where $`J(\chi )=1/(8e)𝑑\epsilon \text{Tr}\left[\stackrel{ˇ}{\tau }_K\stackrel{ˇ}{J}(x)\right]`$ and $`J(\chi )`$ does not depend on the position $`x`$ where $`\stackrel{ˇ}{J}`$ is evaluated.

The Green’s function $`\stackrel{ˇ}{g}(x)`$ in the wire is determined by the diffusion-like Usadel equation

$$\frac{\mathrm{}D}{LG_D}_x\stackrel{ˇ}{J}(x)=i\epsilon [\stackrel{ˇ}{\tau }_3,\stackrel{ˇ}{g}(x)],$$
(4)
(where $`\stackrel{ˇ}{\tau }_3=\widehat{\tau }_3\overline{1}`$ and $`D`$ is the diffusion constant in the wire) and by boundary conditions at the two extremities. We assume a good contact on the normal side at $`x=L`$. This implies $`\stackrel{ˇ}{g}(L)=\stackrel{ˇ}{g}_N`$ as defined in (1). On the superconducting side, at $`x=0`$, the boundary condition expresses the conservation of the spectral matrix current (2) through the interface zaitsev ; CircuitTheory :

$$\stackrel{ˇ}{J}(0)=G_Q\underset{n}{}\frac{2\mathrm{\Gamma }_n[\stackrel{ˇ}{g}(0),\stackrel{ˇ}{g}_S]}{4+\mathrm{\Gamma }_n(\{\stackrel{ˇ}{g}(0),\stackrel{ˇ}{g}_S\}2)},$$
(5)
where the eigenvalues $`\mathrm{\Gamma }_n`$ of the transmission matrix through the interface appear explicitly and $`G_B=G_Q_n\mathrm{\Gamma }_n`$, with $`G_Q=e^2/(\pi \mathrm{})`$ the quantum of conductance. The Thouless energy, $`E_T=\mathrm{}D/L^2`$, is the relevant energy scale in Eq. (4). We restrict to the case where it is much smaller than the superconducting gap. With the same restrictions on the voltage bias and temperature, the Green’s function in the superconducting reservoir is $`\stackrel{ˇ}{g}_S=\widehat{\tau }_2\overline{1}`$. Extension to higher energies is straightforward, but it will not be considered to keep the presentation simple. The quasiclassical Green’s functions obey the normalization condition $`\stackrel{ˇ}{g}^2=1`$ and the symmetry property

$$\stackrel{ˇ}{g}(x)^{}=\stackrel{ˇ}{\tau }_L\stackrel{ˇ}{g}(x)\stackrel{ˇ}{\tau }_L,\mathrm{with}\stackrel{ˇ}{\tau }_L=\widehat{\tau }_3\overline{\sigma }_2.$$
(6)
Finding $`J(\chi )`$ for all values of $`\chi `$ is equivalent to calculate the full counting statistics of charge transfer FCS  and it may be a formidable task. One of the main difficulties comes from the normalization condition $`\stackrel{ˇ}{g}^2=1`$. As a matter of fact, for $`\chi =0`$ the Green’s functions have a triangular structure in Keldysh space:

$$\stackrel{ˇ}{g}_0=\stackrel{ˇ}{g}(\chi =0)=\left(\begin{array}{cc}\widehat{g}^R& \widehat{g}^K\\ 0& \widehat{g}^A\end{array}\right).$$
(7)
In the absence of supercurrent it exists a simple parametrization fulfilling the normalization condition and the symmetry property (6): $`\widehat{g}^R=\widehat{\tau }_3\mathrm{cosh}\theta +i\widehat{\tau }_2\mathrm{sinh}\theta `$, $`\widehat{g}^K=\widehat{g}^R\widehat{f}\widehat{f}\widehat{g}^A`$, and $`\widehat{f}=f_L+\widehat{\tau }_3f_T`$, with $`\theta =\theta _1+i\theta _2`$ and $`f_L`$, $`f_T`$, $`\theta _1`$, and $`\theta _2`$ real. When $`\chi 0`$ the triangular structure is lost and we are not aware of a simple parametrization for $`\stackrel{ˇ}{g}`$ in this case.

A less ambitious program is to calculate only the noise. In this case it is clear from Eq. (3) that the full dependence of $`\stackrel{ˇ}{g}(\chi )`$ is not needed: $`\stackrel{ˇ}{g}(\chi )/\chi |_{\chi =0}`$ is enough. As suggested in Ref. BN01  one can thus develop $`\stackrel{ˇ}{g}(\chi )`$ in powers of $`\chi `$:

$$\stackrel{ˇ}{g}(x,\chi )=\stackrel{ˇ}{g}_0(x)i(\chi /2)\stackrel{ˇ}{g}_1(x)+𝒪(\chi ^2)$$
(8)
and solve the problem order by order in $`\chi `$. For the leading order Green’s function, $`\stackrel{ˇ}{g}_0`$, the above parametrization is appropriate. Eq. (4) for $`\widehat{g}^R`$ gives the non-linear equation

$$\mathrm{}D\theta ^{\prime \prime }(x)+2i\epsilon \mathrm{sinh}\theta (x)=0,$$
(9)
with boundary conditions $`\theta (L)=0`$ and

$$L\theta ^{}(0)=\frac{i}{r}\frac{\mathrm{cosh}\theta (0)}{1+\frac{\mathrm{\Gamma }}{2}\left(i\mathrm{sinh}\theta (0)1\right)}.$$
(10)
We defined $`\psi (\mathrm{\Gamma })_n\mathrm{\Gamma }_n\psi (\mathrm{\Gamma }_n)/_n\mathrm{\Gamma }_n`$ and $`r=G_D/G_B`$. Eq. (4) for $`\widehat{g}^K`$ gives $`f_L(x)=f_{L0}`$ and the differential equation for $`f_T`$: $`(\mathrm{cosh}^2\theta _1(x)f_T^{}(x))^{}=0`$ with boundary conditions $`f_T(L)=f_{T0}`$ and

$$f_T^{}(0)=\frac{f_T(0)\theta _1^{}(0)}{\mathrm{cosh}\theta _1(0)\mathrm{sinh}\theta _1(0)}.$$
(11)
Then the current is $`I=1/(2e)𝑑\epsilon 𝒢(\epsilon )f_{T0}(\epsilon )`$ with VZK93

$$𝒢(\epsilon )=G_D\left[𝒟^1(\epsilon )+\frac{\mathrm{tanh}\theta _1(0)}{L\theta _1^{}(0)}\right]^1$$
(12)
and $`𝒟^1(\epsilon )=1/L_0^L𝑑s/\mathrm{cosh}^2\theta _1(s).`$ At low temperatures, $`k_BTeV`$, $`G(V)=𝒢(eV)`$ is the differential conductance. These equations have been recently exploited in Ref. tanaka  to discuss the conductance.

At the next order in $`\chi `$, we get a linear matrix differential equation for $`\stackrel{ˇ}{g}_1(x)`$:

$$\mathrm{}D_x\stackrel{ˇ}{J}_1(x)=i\epsilon LG_D[\stackrel{ˇ}{\tau _3},\stackrel{ˇ}{g}_1(x)]$$
(13)
with $`\stackrel{ˇ}{J}_1(x)=LG_D[\stackrel{ˇ}{g}_0(x)_x\stackrel{ˇ}{g}_1(x)+\stackrel{ˇ}{g}_1(x)_x\stackrel{ˇ}{g}_0(x)]`$. The boundary conditions read $`\stackrel{ˇ}{g}_1(L)=[\stackrel{ˇ}{\tau }_K,\stackrel{ˇ}{g}_N^0]`$ on the normal side and $`\stackrel{ˇ}{J}_1(0)=2G_B\stackrel{ˇ}{A}\stackrel{ˇ}{B}\stackrel{ˇ}{A}`$ where $`\stackrel{ˇ}{A}=[4+\mathrm{\Gamma }(\{\stackrel{ˇ}{g}(0),\stackrel{ˇ}{g}_S\}2)]^1`$ and $`\stackrel{ˇ}{B}=[4+2\mathrm{\Gamma }(\stackrel{ˇ}{g}_S\stackrel{ˇ}{g}_0(0)\stackrel{ˇ}{g}_1(0)\stackrel{ˇ}{g}_S\stackrel{ˇ}{g}_0(0)\stackrel{ˇ}{g}_1(0)[\stackrel{ˇ}{g}_1(0),\stackrel{ˇ}{g}_S])]`$ on the superconducting side. Finally, the normalization of $`\stackrel{ˇ}{g}`$ gives the condition $`\{\stackrel{ˇ}{g}_0(x),\stackrel{ˇ}{g}_1(x)\}=0`$ that can be fulfilled with the change of variable $`\stackrel{ˇ}{g}_1(x)=[\stackrel{ˇ}{g}_0(x),\stackrel{ˇ}{\varphi }(x)]`$. The matrix $`\stackrel{ˇ}{\varphi }(x)`$ is now constrained by the symmetry properties (6) only. We find that Eq. (13) can be conveniently solved with the following parametrization for $`\stackrel{ˇ}{\varphi }`$:

$$\stackrel{ˇ}{\varphi }=\left(\begin{array}{cc}af_{T0}\widehat{\tau }_1c\widehat{f}\widehat{\tau }_3& b\widehat{\tau }_3+d\\ c\widehat{\tau }_3& a^{}f_{T0}\widehat{\tau }_1+c\widehat{f}\widehat{\tau }_3\end{array}\right)$$
(14)
with $`a=a_1+ia_2`$ and $`a_1`$, $`a_2`$, $`b`$, $`c`$, $`d`$ real functions of $`x`$. Substituting $`\stackrel{ˇ}{g}_1`$ in terms of $`\stackrel{ˇ}{\varphi }`$ into (13) after straightforward but lengthly calculations we obtain the set of equations and boundary conditions for the four parameters $`a`$, $`b`$, $`c`$, and $`d`$. The parameter $`a`$ plays the role of $`\theta `$ in the lowest order equation (9):

$$\mathrm{}Da^{\prime \prime }(x)+2i\epsilon a(x)\mathrm{cosh}\theta (x)=2E_T\frac{\mathrm{sinh}\theta _1(x)}{\mathrm{cosh}^3\theta _1(x)}\frac{G(\epsilon )^2}{G_D^2}.$$
(15)
The boundary conditions are $`a(L)=0`$ and $`La^{}(0)=\alpha a(0)/r+\beta /r`$ with

$$\alpha =\frac{i\mathrm{sinh}\theta \mathrm{\Gamma }(i\mathrm{sinh}\theta 1)/2}{\left[1+\mathrm{\Gamma }(i\mathrm{sinh}\theta 1)/2\right]^2},\beta =\frac{ic^2}{8}\frac{2\mathrm{\Gamma }^2\mathrm{cosh}\theta ^{}+8(\mathrm{\Gamma }1)\mathrm{cosh}\theta 2i\mathrm{\Gamma }(\mathrm{\Gamma }2)\mathrm{sinh}\theta \mathrm{cosh}\theta ^{}}{\left|1+\mathrm{\Gamma }(i\mathrm{sinh}\theta 1)/2\right|^2\left(1+\mathrm{\Gamma }(i\mathrm{sinh}\theta 1)/2\right)},$$
(16)
both evaluated at $`x=0`$. The equation for $`b`$ resembles that for $`f_T`$: it is solved analytically in terms of an integral that enters the final expression for the noise. The parameter $`c`$ turns out to be simply proportional to $`f_T(x)`$: $`c(x)=f_T(x)/f_{T0}`$ with $`c(0)=1𝒢(\epsilon )/[G_D𝒟(\epsilon )]`$. Finally we find $`d(x)=2f_{L0}f_{T0}[1+c^3(x)2\mathrm{tan}\theta _2(x)]`$. Its knowledge is not necessary to obtain the noise.

The expression for the noise is obtained by evaluating Eq. (3):

$$S=𝑑\epsilon 𝒢(\epsilon )\left\{1f_{L0}^2(\epsilon )[1(\epsilon )]f_{T0}^2(\epsilon )\right\},$$
(17)
where

$`(\epsilon )={\displaystyle \frac{2}{3}}(1+c(0)^3)+{\displaystyle \frac{2𝒢(\epsilon )}{G_D}}{\displaystyle _0^1}{\displaystyle \frac{\mathrm{sinh}\theta _1a_1}{\mathrm{cosh}^3\theta _1}}𝑑s`$ (18)
$`c(0)\left({\displaystyle \frac{G_Da_1^{}(0)c(0)}{𝒢(\epsilon )\mathrm{tanh}\theta _1(0)}}+{\displaystyle \frac{2a_1(0)}{\mathrm{sinh}2\theta _1(0)}}\right).`$
In equilibrium, Eq. (17) yields the fluctuation dissipation relation $`S=4k_BTG`$, as expected. Out of equilibrium, in the shot-noise regime $`k_BTeV`$, Eq. (18) defines the experimentally accessible differential Fano factor $`F(V)(dS(V)/dV)/(2eG(V))=(eV)`$.

Eq. (18) is the central result of the paper. Once the problem for the conductance has been solved and $`\theta (x)`$ and $`f_T(x)`$ are known VZK93 ; tanaka , we provide a simple and efficient way to calculate the noise. It suffices to solve the linear differential equation (15) with the given boundary conditions, and substitute the result into Eq. (18). This program has to be followed numerically in most situations, but its implementation is straightforward and allows to obtain quantitative predictions for a wide range of realizable experiments. The previous approach to the same problem has been to discretize Eq. (4) from the outset and solve it numerically CircuitTheory ; BN01 . That technique is particularly efficient to treat the simplified boundary conditions $`\stackrel{ˇ}{J}(0)=2G_B[\stackrel{ˇ}{g}(0),\stackrel{ˇ}{g}_S]`$ KL , valid only when all $`\mathrm{\Gamma }_n1`$. Its application to the general case is numerically much more difficult.

We can now discuss the crossover from the coherent ($`\epsilon E_T`$) to the incoherent semiclassical ($`\epsilon E_T`$) regime quite generally.

In the fully coherent regime we recover known expressions obtained with random-matrix theory dJB94 ; Beenakker97  or quasiclassical Green’s functions Nazarov94 . In the opposite limit, propagation is incoherent: due to the large value of $`\epsilon `$ the phases accumulated along two different paths are always uncorrelated. In order to verify this point we first repeated the calculation when both reservoirs are normal. It suffices to substitute $`\stackrel{ˇ}{g}_S`$ with $`\stackrel{ˇ}{g}_N`$ into (5) with $`eV=0`$. (The procedure greatly simplifies since $`\stackrel{ˇ}{g}`$ becomes diagonal in Nambu indexes.) As expected, we find no energy dependence for $`G`$ and $`F`$. For the conductance the usual Ohm’s law holds: $`G^1=G_D^1+G_B^1`$ while for the Fano factor we find

$$F=\frac{1}{3}\left[1+\left(23\frac{\underset{n}{}\mathrm{\Gamma }_n^2}{_n\mathrm{\Gamma }_n}\right)\frac{G_D^3}{(G_D+G_B)^3}\right].$$
(19)
Eq. (19) coincides with the semiclassical result obtained in Ref. deJongBen . In the superconducting case, Eq. (12) and Eq. (18) can be evaluated analytically for $`\epsilon E_T`$ defining $`G_{\mathrm{inc}}`$ and $`F_{\mathrm{inc}}`$, respectively. It turns out that the conductance and the Fano factor are given by the same expressions for the normal case multiplied by a factor $`2`$ and where the substitutions $`G_DG_D/2`$ and $`\mathrm{\Gamma }_n\mathrm{\Gamma }_n^A=\mathrm{\Gamma }_n^2/(2\mathrm{\Gamma }_n)^2`$ have been operated (the last one implies $`G_BG_B\mathrm{\Gamma }/(2\mathrm{\Gamma })^2`$). This result can be understood Bezuglyi ; SB ; BelzigSamuelsson  in terms of semiclassical probability for charge transfer, that is the probability for an electron to cross the distance $`L`$ to the barrier, to be Andreev reflected, and then to come back as a hole on the same distance. Therefore, the wire in contact with the superconductor can be described by two normal wires, one for the electrons and the other for the holes, coupled through a fictitious barrier that converts electrons into holes with probability $`\mathrm{\Gamma }_n^A`$ and to which a voltage bias of $`2V`$ is applied Bezuglyi .

Quantum corrections to both $`G(\epsilon )`$ and $`F(\epsilon )`$ can be evaluated as an expansion for large $`\epsilon `$. We give the result for all $`\mathrm{\Gamma }_n=\mathrm{\Gamma }`$: $`G=G_{\mathrm{inc}}+G_1(\epsilon )`$ , with $`G_1=2\sqrt{E_T/\epsilon }(44\mathrm{\Gamma }\mathrm{\Gamma }^2)/(r(\mathrm{\Gamma }2)^2+2\mathrm{\Gamma })^2`$, and

$$F(\epsilon )=F_{\mathrm{inc}}r\gamma (\mathrm{\Gamma },r)[G(\epsilon )G_{\mathrm{inc}}]$$
(20)
with $`\gamma =[r(2\mathrm{\Gamma })^2(64+(2\mathrm{\Gamma })\mathrm{\Gamma }(64+r(2\mathrm{\Gamma })(1212\mathrm{\Gamma }\mathrm{\Gamma }^2)))]/[(r(2\mathrm{\Gamma })^2+2\mathrm{\Gamma })^2(44\mathrm{\Gamma }\mathrm{\Gamma }^2)]`$. $`\gamma `$ is positive for most values of $`\mathrm{\Gamma }`$. Eq. (20) relates the interference contributions to $`G`$ with those to $`F`$ for large $`\epsilon `$. When the barrier dominates ($`r1`$) the relation (20) holds for any value of $`\epsilon `$ up to first order in $`1/r`$ with $`\gamma (\mathrm{\Gamma })=\gamma (\mathrm{\Gamma },r\mathrm{})`$ and $`G_1(\epsilon )=2\varphi (2\sqrt{\epsilon /E_T})(44\mathrm{\Gamma }\mathrm{\Gamma }^2)/[r^2(\mathrm{\Gamma }2)^4]`$ and $`\varphi (x)=(\mathrm{sinh}x+\mathrm{sin}x)/[x(\mathrm{cosh}x+\mathrm{cos}x)]`$. This suggests that the simple relation (20) stands, though approximately, beyond the range of validity for which it has been proved. A possible interpretation is that interference modifies the effective transparency of the whole system. The linear relation between the quantum corrections to $`G`$ and $`F`$ then corresponds to the single channel quantum result: $`G\mathrm{\Gamma }`$ and $`F1\mathrm{\Gamma }`$ BlanterButtiker .

For intermediate values of the parameters $`r`$, $`\epsilon `$, and $`\mathrm{\Gamma }_n`$, quantum interference contribution can be studied numerically. We report the results in Fig. 1 for $`G_D=G_B`$ and different values of $`\mathrm{\Gamma }=\mathrm{\Gamma }_n`$ for all $`n`$.

The transparency $`\mathrm{\Gamma }`$ drives a smooth crossover from a reflectionless tunneling to a reentrant behavior which can be seen both in the conductance and the Fano factor. This proves that the energy dependence of both may be strongly affected by the precise set $`\{\mathrm{\Gamma }_n\}`$ of transparencies of the interface. As a relevant example we plot the result for the universal distribution for disordered interfaces $`_n\delta (\mathrm{\Gamma }\mathrm{\Gamma }_n)1/(\mathrm{\Gamma }^{\frac{3}{2}}\sqrt{1\mathrm{\Gamma }})`$ SchepBauer . Fig. 1 also confirms that Eq. (20) is actually qualitatively satisfied in the whole range of parameters.

In the limit $`G_DG_B`$ ($`r0`$) the actual values of $`\{\mathrm{\Gamma }_n\}`$ drop from the boundary conditions. Indeed Eqs. (10,11,16) force $`\theta (0)=i\pi /2`$, $`f_T(0)=0`$ (and thus $`c(0)=0`$), and $`a(0)=0`$. The conductance becomes $`𝒢(\epsilon )=𝒟(\epsilon )G_D`$ and the Fano factor is given by Eq. (18) where only the $`2/3`$ and the integral terms survive. The first term is the semiclassical incoherent value NagaevButtiker , which also coincides with the fully coherent one Beenakker97 . The integral term singles out the interference contribution to the Fano factor. The energy dependence coincides with that obtained previously in a different way BN01  which, in turn, agrees qualitatively with the experimental result KSP . However, a broader voltage range for the reentrance in $`G(V)`$ and $`S(V)`$ is predicted impr . A possible explanation is that the barrier resistance is not negligible. Assuming a disordered interface, we find that the reasonably small value of $`r=0.3`$ allows to fit the conductance. It also improves the fit for the noise, but the agreement is not perfect. This may be due either to a different distribution of transparencies at the barrier, or, as suggested in Ref. BN01 , to heating or interaction effects.

In conclusion, we provide a framework to calculate both conductance and noise in a normal metallic wire connected to a superconducting lead through an arbitrary interface. We predict a strong energy dependence of the coherent contribution to the Fano factor. We suggest to exploit this dependence to experimentally characterize the transparency of interfaces.

We thank V. P. Mineev, M. Sanquer, and F. W. J. Hekking for useful discussion. F.P. acknowledges financial support from CNRS through contract ATIP-JC 2002.


text:
# Cosmological dynamics of scalar field with non-minimal kinetic term

## I Introduction

In this paper we investigate dynamics of scalar field theory with the kinetic term dependent on the field value. This hypothesis provides an additional resources in an explanation of observational data. The small value of cosmological term, consistent with recent experimental data Knop03  can be explained using a non-canonical form of the kinetic term in the scalar field (like in the quintessential model Damour99 ; Chiba99 ). Substantial review of different ways of origin of the cosmological term can be found in Sahni00 . A non-trivial kinetic term could be responsible for a new coupling between adiabatic and entropy perturbations Marco02 .

Such a model may be considered as the restricted scalar-tensor theory (STT) of gravity Bergmann . There is much current interest in the STT, because it allows for explanations of variety of cosmological phenomena Boisseau00 ; Wang ; Esp03  like, e.g., inflation at a low energy scale Morris01  and dark matter caused by phantom fields Singh03 . This theory, being a generalization of general relativity, possesses a much richer structure then the latter. The model has a sector of a scalar field and it seems natural to relate this component to the inflation. Historically, the first internally self-consistent inflationary model Star79  appears to be equivalent to the special limiting case of STT Equiv .

The exact form of STT action has not been derived yet from a more fundamental theory. Thus a large number of phenomenological potentials and kinetic terms are considered in current models. Some of those can be excluded being inconsistent with astrophysical data. Still a substantial number of them remains topical. In this connection two questions can be raised and need to be answered: (i) What shape of potential is realized in our universe? Finding an answer would require that the theoretical prediction from such model potential would fit most observational data. Suppose one would find answer, then the next question is: (ii) Why is this particular shape of potential (or kinetic term) realized in nature? What are the underlying reasons?

Some theoretical hints on the form of the potential have been given by supergravity, which predicts an infinite power series expansion in the scalar field potential Nilles . Its minima, if they exist, correspond to the stationary states of the field. In the low energy regime it is reasonable to retain only a few terms (lowest powers in the Taylor expansion) of the scalar field Lyth96 . However, the potential caused, e.g., by supergravity could correspond to a function with infinite set of the potential minima. There is still no physical law which limits the potential to possess only a finite number of minima. In the vicinity of each of the minima the potential has an individual form and hence the universe associated with such minimum may be individually different from others. Our own universe is associated with a particular potential minimum, not necessarily located at $`\phi =0`$. The supposition that the potential possesses infinite number of randomly distributed minima is self-consistent Ru42a . A similar behavior may hold also for the kinetic term of STT.

The observed smallness of the value of the $`\mathrm{\Lambda }`$ term is explained usually on the basis of a more fundamental theory like supergravity or the anthropic principle. Our point of view is that we have to merge these approaches. The more fundamental theory supplies us with an infinite set of minima of the potential. These minima having an individual shape are responsible for the formation of those universes used in the anthropic picture.

Actually, any way of introducing a scalar field leads to a theory similar to STT provided that quantum corrections are taken into account. Indeed, STT arises necessarily in ordinary quantum field theory with a scalar field. For example, it is a standard result that the kinetic term of the effective action acquires a multiplier which is a function of the field Itzykson .

Here we consider an action with a non-trivial kinetic term of the following form

$`S={\displaystyle d^4x\sqrt{g}\left[\frac{R}{16\pi G}+\frac{1}{2}K(\phi )_\mu \phi ^\mu \phi V(\phi )\right]}.`$ (1)
In this paper we concentrate on the influence of singular points of the kinetic term on the scalar field dynamics. The singular points are not exceptional cases. E.g., the well known Brans - Dicke model Brans  does contain a singularity at zero value of the field $`\phi `$. If a multiplier of the Ricci scalar equals to zero at some point in the Jordan frame, the kinetic term will be singular in the Einstein frame (see e.g. Morris01 ). Quintessential models with negative power law Peebles  are another well known example for the potential of such a sort.

By a suitable change of variable the function $`K`$ can be transformed to become $`K(\phi )=\pm 1`$. Such behavior is physical during some inflationary period Ru  but not at some recent epoch where the field fluctuates around the singular point. The similar problem is discussed in refs. Star81 , Bronnikov02  and Esp  in the framework of more general case of STT.

The general form of the kinetic term considered in this paper is

$`K(\phi )=M^n/(\phi \phi _s)^n.`$ (2)
This form is correct at least in the vicinity of the singularity if $`n>0`$, or zero if $`n<0`$. Here $`M`$ is some parameter (its value will be discussed below).

The equation of motion has the form

$`K(\phi )[\ddot{\phi }+3H\dot{\phi }]+{\displaystyle \frac{1}{2}}K(\phi )^{}\dot{\phi }^2+V(\phi )^{}=0`$ (3)
in the Friedmann-Robertson-Walker universe where $`H`$ denotes the Hubble parameter. Keeping in mind expression (2) we have

$`\ddot{\phi }+3H\dot{\phi }{\displaystyle \frac{n}{2(\phi \phi _s)}}\dot{\phi }^2+V(\phi _s)^{}(\phi \phi _s)^n/M^n=0.`$ (4)
We observe that $`\phi =\phi _s`$ is a stationary solution for any smooth potential $`V`$ and $`n>0`$. The cosmological energy density of the vacuum is connected usually with one of the potential minima. In the case considered here it is not like this - it could be located at the singular point of the kinetic term $`K(\phi )`$.

In the following we limit ourselves by considering the first two terms in Taylor series of the potential and shift the field variable so that the potential acquires the simple form

$`V(\phi )=V_0+m^2\phi ^2/2.`$ (5)
To facilitate a more detailed analysis, a new auxiliary variable $`\chi `$ is useful. We suggest the substitution of variables $`\phi \chi `$ in the following manner

$`d\chi =\pm \sqrt{K(\phi )}d\phi ,K(\phi )>0.`$ (6)
An action in terms of the auxiliary field $`\chi `$ has the form

$`S={\displaystyle d^4x\sqrt{g}\{\frac{R}{16\pi G}+sgn(\phi \phi _s)\frac{1}{2}_\mu \chi ^\mu \chi U(\chi )\}},`$ (7)
where the potential $`U(\chi )V(\phi (\chi ))`$ ’partly discontinues’ function. Its form depends on the form of initial potential $`V(\phi )`$, the form of kinetic term and a position of the singularities. Now let us consider particular cases of $`K(\phi )`$.

## II The case $`n=1`$

The kinetic term (2) is $`K(\phi )=M/(\phi \phi _s)`$ and Eq.(6) reads

$`\phi =\phi _s+sgn(\phi \phi _s)\chi ^2/4M`$ (8)
with the potential

$`U(\chi )V(\phi (\chi ))=V_0+{\displaystyle \frac{m^2}{2}}(\phi _s+sgn(\phi \phi _s){\displaystyle \frac{\chi ^2}{4M}})^2;\phi _s>0;|\chi |<\mathrm{}.`$ (9)
The situation is illustrated by Figure 1. If the physical field $`\phi <\phi _s`$, than the auxiliary field $`\chi `$ behaves like the phantom field Phantom ; if the the field $`\phi `$ fluctuates around $`\phi _s`$ the substitution (6) is meaningless.

The inflationary scenario on the basis of potential (9) is well known and we discuss only its final stage. If $`\phi _s=0`$ and $`V_0=0`$ exactly, we obtain a potential $`U\chi ^4`$. If $`\phi _s>0`$, the auxiliary field $`\chi `$ oscillates, according to classical equation, with damping around its final position at $`\chi =0`$ \- see Figure 2.

For the physical field $`\phi `$, the motion looks like oscillations around the point $`\phi _s`$. Finally, the field is captured in the vicinity of the singular point, supplying an energy density of the vacuum equal to

$`\mathrm{\Lambda }=V_0+{\displaystyle \frac{1}{2}}m^2\phi _s^2,`$ (10)
rather than $`\mathrm{\Lambda }=V_0`$ as it could be expected from expression (5).

The problem of smallness of the vacuum energy density, $`\mathrm{\Lambda }=10^{123}M_P^4`$, is still topical. In other words, what physical reason adjusts the location of the singular point $`\phi _s`$ such that it lies in the extremely tiny interval $`\sqrt{2V_0/m^2}÷\sqrt{2V_0/m^2+2\mathrm{\Lambda }/m^2}`$? In the view of the above discussion, we could expect that some minima $`\phi _m`$ of the potential could incorporate critical point(s) $`\phi _s`$ at specific distance from this minimum. Now the problem is reformulated as: “what part of an infinite amount of the minima contains singular points located so closely to the minima?” It seems plausible that this part is very small, but not zero, due to an infinite number of minima. Only this part is important - it represents those vacua where galaxies could be formed Weinberg .

## III The case $`n=2`$

In this case, the inflationary period is more interesting than in the previous case, so we will discuss it more thoroughly. The kinetic term has the form $`K(\phi )K_s(\phi )=M^2/(\phi \phi _s)^2`$. The new auxiliary field is connected with the physical one by

$$\chi =M\mathrm{ln}\left|\frac{\phi \phi _s}{\phi _s}\right|,$$

so that the potential has the form

$`U(\chi )={\displaystyle \frac{1}{2}}m^2\phi _s^2(1+sgn(\phi _s)sgn(\phi \phi _s)e^{\chi /M})^2+V_0.`$ (11)
If $`\phi _s>0`$ and $`\phi >\phi _s`$, the potential mimics the quintessential model with nonzero vacuum energy density $`\mathrm{\Lambda }=\frac{1}{2}m^2\phi _s^2+V_0`$.

The case $`\phi _s<0`$ and $`\phi >\phi _s`$ is much more interesting. This potential is highly asymmetric so that a behavior of inflaton is rather different at $`\chi <0`$ and at $`\chi >0`$

Let us suppose that inflation starts with $`\chi _{in}>0`$, i.e. point $`1`$ in Figure 3. The picture is similar to the improved quintessence potential Albrecht , without problems of the radiation-dominated stage during Big Bang nucleosynthesis. A modern epoch is characterized by a large negative value of the field $`\chi `$. It slowly varies along the flat part of the potential with exponentially slow variation of the vacuum energy density around the value (10).

It is worth to discuss another possibility, when the inflation starts at $`\chi =\chi _{in}<0`$, see Figure 3, point 2. The form of the potential $`\chi \mathrm{}`$ is rather plain being suitable for slow rolling, so it is reasonable to expect some new results. Indeed, our estimations indicate that the inflation could take place at the parameter values $`MM_P`$ (Planck scale), $`m\phi _s10^3M_P`$ (GUT scale) and $`|\chi |M_p\mathrm{ln}5`$. Evidently, this inflationary model does not require fine tuning of its parameters. The inflation ends at the minimum of the potential, $`\chi =0`$ and the final energy density of the vacuum becomes $`\mathrm{\Lambda }=V_0`$.

## IV The case $`n=1`$

A nontrivial situation takes place if the kinetic function has zero value at some point, i.e. $`K(\phi )=(\phi \phi _s)/M`$. Eq.6 leads to the classical equation

$`(\phi _s\phi )(\ddot{\phi }+3H\dot{\phi }){\displaystyle \frac{1}{2}}\dot{\phi }^2+MV(\phi )^{}=0`$ (12)
This equation is quite uncommon. Indeed, if the zero point $`\phi _s`$ of the kinetic term does not coincide exactly with the position of the potential minimum at $`\phi =0`$, the point $`\phi =\phi _s`$ is not a stationary solution of this equation. On the other hand, the point $`\phi _s`$ is some kind of attractor. Namely, if the field value is larger than $`\phi _s`$, then $`K(\phi )>0`$ and we have standard rolling of the field down to the singular point $`\phi _s`$. If the field value is smaller than $`\phi _s`$, then $`K(\phi )<0`$. The field behaves like the phantom field Phantom , climbing up to the potential and thus tending toward the point $`\phi _s`$. Classically, the situation looks very strange - the singular point attracts the solution, but forbids it to stay there forever. Evidently, the field fluctuates stochastically around the singular point. Additional discussion on this problem can be found in the paper Frolov02 .

Let the initial field value $`\phi =\phi _{in}>\phi _s`$. An appropriate variable substitution (6) looks like

$$\phi =\phi _s+sgn(\phi \phi _s)\gamma |\chi |^{2/3},\gamma (3\sqrt{M}/2)^{2/3}\phi >\phi _s.$$

The potential of the auxiliary field $`\chi `$ becomes

$`U(\chi )={\displaystyle \frac{1}{2}}m^2(\phi _s+sgn(\phi \phi _s)\gamma |\chi |^{2/3})^2+V_0.`$ (13)
$`U(\chi )`$ is finite at $`\chi =0`$ but its derivative is singular,

$$U_{\chi +0}^{}=\frac{2}{3}m^2\phi _s^2\gamma |\chi |^{1/3}.$$

The potential (13) behaves like $`\chi ^{4/3}`$ at large field values. It leads to standard inflation with moderate fine tuning of the parameters. If $`\phi _s>0`$, the field $`\phi `$ will oscillate around critical point and energy density (10).

## V Discussion

It is shown that the region where the kinetic term changes its sign gives new possibilities for the scalar field dynamics. It takes place even for the simplest form of potential. Depending on the position of the singular point of the kinetic term specific forms of the potential of the auxiliary field can be obtained. One of the main result is that the scalar field could be located in the singular points of the kinetic term rather than in minima of the potential. Another interesting result is that if the kinetic term becomes zero at the singular point, the behavior of scalar field acquires stochastic character.

The problem of the cosmological $`\mathrm{\Lambda }`$ \- term acquires another sense. One has to explain the extremal proximity of the singular point of the kinetic term and zero point of the potential. Such situation seems absolutely accidental. Moreover, the probability of such event to occur is very small. The problem can be solved in the framework of the random potential Ru42a  and a kinetic term of the scalar field discussed in the Introduction. In this case we have infinite number of singular points of the kinetic term. Fluctuations of the scalar field which were created at high energies moves classically to stationary points. Those of them who reach stationary points with appropriate energy density could form a universe similar to our Universe. This energy density ($`10^{123}M_P^4`$) is the result of small value of a concrete potential minimum or small value of the difference $`\phi _s\phi _0`$, $`\phi _0`$ is a zero of the potential. The fraction of such universes is extremely small, but nevertheless is finite because of infinite number of stationary states.

The abundance of the variants discussed in the paper is based on the simplest forms of the potential and the kinetic terms. Nevertheless, we reproduced the quintessence feature, proposed a new possibility for inflation without fine tuning and discussed a fluctuating state of the scalar field as a final state which be could realized in the Universe.

One of us (S.G.R) is grateful to K. A. Bronnikov and V. M. Zhuravlev for discussion. The work was partially performed in the framework of Russian State contract $`\mathrm{40.022.1.1.1106}`$, RFBR grant $`020217490`$ and of grant of Russian Universities UR.02.01.026.