text:
# Quantum Monte Carlo study of confined fermions in one-dimensional optical lattices.

## I Introduction

The realization of Bose-Einstein condensation (BEC) of trapped atomic gases anderson ; bradley ; davis  has generated in the last years a huge amount of experimental and theoretical research bose ; dalfovo . BEC was achieved by confining atoms in magnetic traps and lowering the temperature of the system via the evaporative cooling technique. In this technique, the hottest atoms are selectively removed from the system and the remaining ones rethermalize via two-body collisions. A common feature of these experiments is that the trapped gases are dilute; mean-field theory then provides a useful framework to study the role of the interaction between particles.

Recently, a new feature has been added to the experiments: the magnetically trapped condensate is transferred into an optical lattice generated by interfering laser beams greiner1 . With this new experimental setup it is possible to access the strongly correlated regime. The superfluid-Mott-insulator transition was studied greiner1  and other interesting physical phenomena such as the collapse and revival of the condensate have been found greiner2 . The presence of the optical lattice and the fact that the particles interact only via contact interaction lead in a natural way to the Hubbard model as a paradigm for these systems. A theoretical work proposing such experiments jaksch  and recent quantum Monte Carlo (QMC) simulations batrouni ; kaskurnikov  have investigated these systems beyond the mean-field approximation. It has been found that the incompressible Mott-insulating phase always coexists with compressible phases so that a local order parameter batrouni  has to be defined to characterize the system. Also, the phase diagram of trapped bosons was found to be more complex than the one in the homogeneous system batrouni .

The experimental realization of ultracold fermionic gases is more difficult than the bosonic case, and has been achieved only recently demarco ; truscott ; schreck ; granade ; hadzibabic ; roati ; ohara . Unlike bosons, single species fermions cannot be directly evaporatively cooled to very low temperatures because the s-wave scattering that could allow the gas to rethermalize during the evaporation is prohibited for identical fermions. This problem has been overcome by simultaneously trapping and evaporatively cooling two-component Fermi gases demarco ; granade , and introducing mixed gases of bosons and fermions in which bosons enable fermions to rethermalize through their elastic interactions truscott ; schreck ; hadzibabic ; roati . More recently, rapid forced evaporation employing a Feshbach resonance of two different spin states of the same fermionic atom has been used ohara . Now that it is possible to go well below the degeneracy temperature in the experiments and superfluidity appears within reach ohara , it is expected that the metal-Mott-insulator transition (MMIT) could also be realized for fermions on an optical lattice.

Motivated by this expectation we study the ground state of the one-dimensional (1D) fermionic Hubbard model with a harmonic trap and with repulsive contact interaction, using QMC simulations. Like in the bosonic case batrouni , Mott-insulating domains appear over a continuous range of fillings and always coexist with compressible phases, such that global quantities are not appropriate to characterize the system. Instead, we define a local-order parameter (a local compressibility) in order to characterize the local phases present in the system. We also analyze the generic features that are valid for any kind of confining potential, and not only for the harmonic one.

It is well known that the Hubbard model in the periodic case displays a MMIT phase transition at half filling and at a finite value of the on-site repulsive interaction $`U`$. (In the case of perfect nesting the transition occurs at $`U=0`$.) The two routes to this MMIT are the filling-controlled MMIT and the bandwidth-controlled MMIT imada . In the confined case, it is possible to drive the transitions by changing the total filling of the trap, the on-site repulsive interaction, or varying the curvature of the confining potential. In the following section, we study a number of local quantities as a function of the filling and the strength of the interaction. Although they signal the appearance of a Mott-insulating phase, such quantities provide no rigorous criterion to characterize the Mott-insulating region. We therefore define a local compressibility that acts as a genuine local-order parameter. In Sec. III, we study the momentum distribution function for the confined system. In Sec. IV, we discuss a mean-field approach for the 1D trapped system and compare the results with the ones obtained using QMC. In Sect. V, we analyze the phase diagram and determine its generic form, which allows to compare systems with different sizes, number of particles and curvatures of the harmonic confining potential. In this section we also analyze the extension, for arbitrary confining potentials, of the results obtained for the harmonic case. Finally, the conclusions are given in Sec. VI. Some of the results presented here were summarized in a previous publication rigol1 .

## II Local order parameter

Due to the fact that the confining potential leads to an inhomogeneous density profile, we study in the present section how local quantities behave in the trapped system when the parameters at hand are changed. The Hamiltonian of the fermionic Hubbard model with a confining parabolic potential has the form

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(c_{i\sigma }^{}c_{i+1\sigma }+\text{H.c.})+U{\displaystyle \underset{i}{}}n_in_i`$ (1)
$`+V{\displaystyle \underset{i\sigma }{}}x_i^2n_{i\sigma },`$
where $`c_{i\sigma }^{}`$, $`c_{i\sigma }`$ are creation and annihilation operators, respectively, for a fermion with spin $`\sigma `$ at site $`i`$, and $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$, such that $`t`$ is the hopping amplitude, $`U`$ is the on-site interaction that in the present work will be considered repulsive ($`U>0`$), $`V`$ is the curvature of the confining harmonic potential, and $`x_i`$ measures the position of the site $`i`$ ($`x_i=ia`$ with $`a`$ the lattice constant). The number of lattice sites is $`N`$ and is selected so that all the fermions are confined in the trap. We denote the total number of fermions in the trap as $`N_f`$ and consider equal number of fermions with spins up and down ($`N_f=N_f=N_f/2`$). In our simulations, we used the zero-temperature projector method sugiyama ; sorella  adapted from the QMC determinantal algorithm by Blankenbecler, Scalapino, and Sugar scalapino ; blankenbecler ; sugar . The discrete Hubbard-Stratonovich transformation by Hirsch was used hirsch . For details of the algorithm we refer to a number of reviews alejandro ; fakher .

The results for the evolution of the local density ($`n_i=n_i+n_i`$), as a function of the total number of the confined particles, are shown in Fig. 1. For the lowest filling, so that $`n<1`$ at every site, the density shows a profile with the shape of an inverted parabola, similar to that obtained in the non-interacting case vignolo1 , and hence, such a situation should correspond to a metallic phase. Increasing the number of fermions up to $`N_f=60`$, a plateau with $`n=1`$ appears in the middle of the trap, surrounded by a region with $`n<1`$ (metallic). Since in the homogeneous case, a Mott insulator appears at $`n=1`$, it is natural to identify the plateau with such a phase. The Mott-insulating domain in the center of the trap increases its size when more particles are added, but at a certain filling ($`N_f=70`$ here) this becomes energetically unfavorable and a new metallic phase with $`n>1`$ starts to develop in the center of the system. Upon adding more fermions, this new metallic phase widens spatially and the Mott-insulating domains surrounding it are pushed to the borders. Depending on the on-site repulsion strength, they can disappear and a complete metallic phase can appear in the system. Finally, a “band insulator” (i.e., $`n=2`$) forms in the middle of the trap for the highest fillings (after $`N_f=144`$ here). Due to the full occupancy of the sites, it will widen spatially and push the other phases present in the system to the edges of the trap when more fermions are added.

Although the existence of flat regions in the density profile is an indication that there is an insulator there, a more quantitative characterization is needed. As shown in Ref. batrouni , the variance of the local density ($`\mathrm{\Delta }_i=n_i^2n_i^2`$) may be used on a first approach (from here on, we refer to the variance as the variance of the local density). In Fig. 2 we show four characteristic profiles present in Fig. 1 and their respective variances when the number of fermions in the system are $`N_f=`$50 (a), 68 (b), 94 (c), and 150 (d). For the case in which only a metallic phase is present \[Fig. 2(a)\], it is possible to see that the variance decreases when the density approaches $`n=1`$ and has a minimum for densities close to that value. For the Mott-insulating domain \[$`n=1`$ in Fig. 2(b)\], the variance has a constant value that is smaller than that of the metal surrounding it. As soon as the Mott-insulating phase is destroyed in the middle of the trap and a new metallic region with $`n>1`$ develops there \[Fig. 2(c)\], the variance increases in this region. The variance in the metallic region with $`n>1`$ will start to decrease again when the density approaches $`n=2`$ and will have values even smaller than those in the Mott-insulating phase. Finally, when the insulator with $`n=2`$ is formed in the center of the trap the variance vanishes there \[Fig. 2(d)\].

Alternatively, Mott-insulating regions can be obtained by increasing the ratio between the on-site repulsive interaction and the hopping parameter as was done in experiments for bosons confined in optical lattices greiner1 . Fig. 3(a) shows the evolution of the density profiles in a trapped system with $`N=100`$, $`N_f=70`$, and $`Va^2=0.0025t`$ when this ratio is increased from $`U/t=2`$ to $`U/t=8`$ \[for more details of these density profiles see Fig. 7(e)-7(h)\]. It can be seen that for small values of $`U/t`$ ($`U/t=2`$) there is only a metallic phase present in the trap. As the value of $`U/t`$ ($`U/t=4`$) is increased, a Mott-insulating phase tries to develop at $`n=1`$ while a metallic phase with $`n>1`$ is present in the center of the system. As the on-site repulsion is increased even further ($`U/t=6,8`$), a Mott-insulating domain appears in the middle of the trap suppressing the metallic phase that was present there. In Fig. 3(b) we show the variance of the density for the profiles in Fig. 3(a) (from top to bottom, the values presented are for $`U/t=2,4,6,8`$). As expected, the variance decreases in both the metallic and Mott-insulating phases when the on-site repulsion is increased. When the Mott-insulating plateau is formed in the density profile, a plateau with constant variance appears in the variance profile with a value that will vanish only in the limit $`U/t\mathrm{}`$. As shown in Fig. 3(b), whenever a Mott-insulating domain is formed in the trap, the value of the variance in it is exactly the same as the one for the Mott-insulating phase in the homogeneous system for the same value of $`U/t`$ (horizontal dashed lines). This would support the validity of the commonly used local density (Thomas-Fermi) approximation butts97 . However, the insets in Fig. 3(b), show that this is not necessarily the case, since for $`U/t=4`$, the value of the variance in the Mott-insulating phase of the homogeneous system is still not reached in the trap, although the density reaches the value $`n=1`$. Therefore, in contrast to the homogeneous case, a Mott-insulating region is not determined by the filling only. In the cases of $`U/t=6`$ \[inset in Fig. 3(b) for a closer look\] and $`U/t=8`$, the value of the variance in the homogeneous system is reached and then we can say that Mott-insulating phases are formed there.

Although the variance gives a first indication for the formation of a local Mott insulator, an ambiguity is still present, since there are metallic regions with densities very close to $`n=0`$ and $`n=2`$, where the variance can have even smaller values than in the Mott-insulating phases. Therefore, an unambiguous quantity is still needed to characterize the Mott-insulating regions. We propose a local compressibility as a local-order parameter to characterize the Mott-insulating regions, that is defined as

$$\kappa _i^{\mathrm{}}=\underset{j\mathrm{}(U)}{}\chi _{i,i+j},$$
(2)
where

$$\chi _{i,j}=n_in_jn_in_j$$
(3)
is the density-density correlation function and $`\mathrm{}(U)b\xi (U)`$, with $`\xi (U)`$ the correlation length of $`\chi _{i,j}`$ in the unconfined system at half-filling for the given value of $`U`$. As a consequence of the charge gap opened in the Mott-insulating phase at half filling in the homogeneous system, the density-density correlations decay exponentially \[$`\chi _{\left(x\right)}\mathrm{exp}^{x/\xi \left(U\right)}`$\] enabling $`\xi (U)`$ to be determined. The factor $`b`$ is chosen within a range where $`\kappa ^{\mathrm{}}`$ becomes qualitatively insensitive to its precise value. Since there is some degree of arbitrariness in the selection of $`b`$, we will examine its role in the calculation of $`\kappa ^{\mathrm{}}`$ here.

We show in Fig. 4 (in a semi-log plot) how the local compressibility behaves as a function of $`\mathrm{}`$ for a homogeneous and a trapped systems in the Mott-insulating ($`n=1`$) and metallic ($`n0.66`$) phases. In the Mott-insulating phase, it is possible to see that in the homogeneous and confined systems the local compressibility decays exponentially to zero in exactly the same way due to the charge gap present there. In the metallic case there is no charge gap so that density-density correlations decay as a power law and the local compressibility, shown in Fig. 4, decays very slowly. In Fig. 4, a departure in the behavior of the local compressibility for the trapped system from the homogeneous case for large values of $`\mathrm{}`$ can be seen. This occurs because correlations with points very close to the band insulating ($`n=0`$) and the Mott-insulating ($`n=1`$) regions surrounding the metallic phase \[see Fig. 3(a)\] are included in $`\kappa `$. However, this effect does not affect the fact that the value of the local compressibility is clearly different from zero as long as the size of the metallic phase is larger than $`2\mathrm{}`$. We obtain that considering $`b10`$ (no matter what the exact value of $`b`$ is), the local compressibility is zero with exponential accuracy (i.e. $`10^6`$ or less) for the Mott-insulating phase and has a finite value in the metallic phase. For the case in Fig. 4 where $`U=6t`$, the correlation length is $`\xi 0.8`$ so that $`\mathrm{}8`$. Physically, the local compressibility defined here gives a measure of the change in the local density due to a constant shift of the potential over a finite range but over distances larger than the correlation length in the unconfined system.

In Fig. 3(c), we show the profiles of the local compressibility for the same parameters as Figs. 3(a) and (b) (we did not include the profile of the local compressibility for $`U=2t`$ because for that value of $`U`$ we obtain that $`\mathrm{}`$ is bigger than the system size). In Fig. 3(c), it can be seen that the local compressibility only vanishes in the Mott-insulating domains. For $`U=4t`$, it can be seen that in the region with $`n1`$ the local compressibility, although small, does not vanish. This is compatible with the fact that the variance is not equal to the value in the homogeneous system there, so that although there is a shoulder in the density profile, this region is not a Mott insulator. Therefore, the local compressibility defined here serves as a genuine local order parameter to characterize the insulating regions that always coexist with metallic phases. At this point we would like to remark that, as shown in Ref. rigol1 , the local compressibility shows critical behavior on entering the Mott-insulating region. This is the reason to speak about phases, since on passing from one region to the other, critical behavior sets in.

Finally, we discuss in this section the spin-spin correlation function, since in periodic chains, quasi-long-range antiferromagnetic correlations appear in the Mott-insulating phase. In Fig. 5 we show the local $`S_i^zS_j^z`$ correlation function for some points of the profiles presented in Fig. 2. We measured the spin-spin correlations in the trap at the points in the figures where it can be clearly seen that $`S_i^zS_j^z`$ has the maximum value. In Fig. 5(a) it can be seen that in the metallic phase with $`n<1`$, the spin-spin correlations decay rapidly and do not show any clear modulation, which is due to the fact that the density is changing in this region. In the local Mott-insulating phase \[center of the trap in Fig. 5(b)\], short-range antiferromagnetic correlations are present, and they disappear completely only on entering in the metallic regions. For the shoulders with $`n1`$ \[Fig. 5(c)\], the antiferromagnetic correlations are still present but due to the small size of these regions they decay very rapidly. Finally for the metallic regions with $`n>1`$, the spin correlations behave like in the metallic phases with $`n<1`$ as shown in Fig. 5(d).

## III Momentum distribution function

In most experiments with quantum gases carried out so far, the momentum distribution function, which is determined in time-of-flight measurements, played a central role. A prominent example is given by the study of the superfluid-Mott-insulator transition greiner1  in the bosonic case. Also a QMC study relating this quantity to the density profiles and proposing how to determine the point at which the superfluid-Mott-insulator transition occurs was presented in Ref. kaskurnikov . As shown below, we find that this quantity is not appropriate to characterize the phases of the system in the fermionic case, and does not show any clear signature of the MMIT.

In Fig. 6 we show the normalized momentum distribution function ($`n_k`$) for the same density profiles presented in Fig. 2. For the trapped systems, we always normalize the momentum distribution to be unity at $`k=0`$. In the case presented in Fig. 6, transitions between phases occur due to changes in the total filling of the trap. We first notice that $`n_k`$ for the pure metallic phase in the harmonic trap \[Fig. 6(a)\] does not display any sharp feature corresponding to a Fermi surface, in clear contrast to the homogeneous case. The lack of a sharp feature for the Fermi surface is independent of the presence of the interaction and is also independent of the size of the system. In the non-interacting case, this can be easily understood: the spatial density and the momentum distribution will have the same functional form because the Hamiltonian is quadratic in both coordinate and momentum. When the interaction is present, it could be expected that the formation of local Mott-insulating domains generates a qualitatively and quantitatively different behavior of the momentum distribution, like in the homogeneous case where in the Mott-insulating phase the Fermi surface disappears and $`n_k`$ is smoother. In Fig. 6(b), it can be seen that there is no qualitative change of the momentum distribution when the Mott-insulating phase is present in the middle of the trap. Quantitatively $`n_k`$ in this case is very similar to the pure metallic cases Fig. 6(a) (for densities smaller than one) and Fig. 6(c) (where in the middle of the trap the density is higher than one). Only when the insulator with $`n=2`$ appears in the middle of the system we find a quantitative change of $`n_k`$ for $`U=4t`$, as shown in Fig. 6(d).

We also studied the momentum distribution function when the MMIT is driven by the change of the on-site repulsion. We did not observe any clear signature of the formation of the Mott-insulating phase in $`n_k`$. In Fig. 7 we show the normalized momentum distribution function \[Fig. 7(a)-7(d)\] for density profiles \[Fig. 7(e)-7(h)\] in which the on-site repulsion was increased from $`U/t=2`$ to 8. It can be seen that the same behavior present in Fig. 6 and the quantitative changes in $`n_k`$ appear only when the on-site repulsion goes to the strong-coupling regime, but this is long after the Mott-insulating phase has appeared in the system.

At this point one might think that in order to study the MMIT using the momentum distribution function, it is necessary to avoid the inhomogeneous trapping potential and use instead a kind of magnetic box with infinitely high potential on the boundaries. However, in that case one of the most important achievements of the inhomogeneous system is lost, i.e., the possibility of creating Mott-insulating phases for a continuous range of fillings. In the perfect magnetic box, the Mott-insulating phase would only be possible at half filling, which would be extremely difficult (if possible at all) to adjust experimentally. The other possibility is to create traps which are almost homogenous in the middle and which have an appreciable trapping potential only close to the boundaries. This can be studied theoretically by considering traps with higher powers of the trapping potentials. As shown below, already non-interacting systems make clear that a sharp Fermi edge is missing in confined systems.

In Fig. 8(a) we show the density profile of a system with 1000 sites, $`N_f=840`$, and a trapping potential of the form $`V_{10}x_i^{10}`$ with $`V_{10}a^{10}=7\times 10^{27}t`$. It can be seen that the density is almost flat all over the trap with a density of the order of one particle per site. Only a small part of the system at the borders has the variation of the density required for the particles to be trapped. In Fig. 8(b) (continuous line), we show the corresponding normalized momentum distribution. It can be seen that a kind of Fermi surface develops in the system but for smaller values of $`k`$, $`n_k`$ is always smooth and its value starts decreasing at $`k=0`$. In order to see how $`n_k`$ changes when an incompressible region appears in the system, we introduced an additional alternating potential, so that in this case the new Hamiltonian has the form

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(c_{i\sigma }^{}c_{i+1\sigma }+\text{H.c.})+V_{10}{\displaystyle \underset{i\sigma }{}}x_i^{10}n_{i\sigma }`$ (4)
$`+V_a{\displaystyle \underset{i\sigma }{}}\left(1\right)^in_{i\sigma },`$
where $`V_a`$ is the strength of the alternating potential. For the parameters presented in Fig. 8(a), we obtain that a small value of $`V_a`$ ($`V_a=0.1t`$) generates a band insulator in the trap, which extends over the region with $`n1`$ (when $`V_a=0`$). However, the formation of this band insulator is barely reflected in $`n_k`$, as can be seen in Fig. 8(b) (dashed line). Only when the value of $`V_a`$ is increased and the system departs from the phase transition \[$`V_a=0.5t`$, dotted line in Fig. 8(c)\], does a quantitatively appreciable change in $`n_k`$ appear.

In general, it is expected that on increasing the system size, the situation in the homogeneous system is approached. However, in the present case it is not merely a question of boundary conditions, but the whole system is inhomogeneous. Therefore, a proper scaling has to be defined in order to relate systems of different sizes. In the case of particles trapped in optical lattices when there is a confining potential with a power $`\alpha `$ and strength ($`V_\alpha `$), a characteristic length of the system ($`\zeta `$) is given by $`\zeta =(V_\alpha /t)^{1/\alpha }`$, so that a characteristic density ($`\stackrel{~}{\rho }`$) can be defined as $`\stackrel{~}{\rho }=N_fa/\zeta `$. We find that this characteristic density is the one meaningful in the thermodynamic limit. In Fig. 9 we show three systems in which the total number of particles and the curvature of the confining potential $`V_{10}`$ were changed, keeping $`\stackrel{~}{\rho }`$ constant.

We measured the positions in the trap in units of the characteristic length $`\zeta `$. These systems have occupied regions ($`n>0`$) with very different sizes, of the order of 500 lattice sites for the dashed line, of the order of 1000 sites for the continuous line, and of the order of 2000 lattice sites for the dotted line. Fig. 9(a) shows that the density profiles scale perfectly when the curvature of the confining potential and the number of fermions are changed in the system, so that in order to show the changes in the density profile, we introduced an inset that expands a region around a density 0.85. In the case of the momentum distribution function \[Fig. 9(b)\], it is possible to see that this quantity also scales very well so that almost no changes occur in the momentum distribution function when the occupied system size is increased, implying that in the thermodynamic limit the behavior of $`n_k`$ is different from the one in the homogeneous system. An expanded view of the region with $`n_k`$ around 0.87 is introduced to better see the scale of the differences in this region. In Sec. V we show that the characteristic density defined here is also the meaningful quantity to define the phase diagram of the system.

## IV Mean-field approximation

The mean-field (MF) approach has been very useful in the study of dilute bosonic gases confined in harmonic potentials dalfovo . This theory for the order parameter associated with the condensate of weakly interacting bosons for inhomogeneous systems takes the form of the Gross-Pitaevskii theory. When the bosonic condensate is loaded in an optical lattice, it is possible to go beyond the weakly interacting regime and reach the strongly correlated limit for which a superfluid-Mott-insulator transition occurs. Also in this limit a MF study was done jaksch , and the results were compared with exact diagonalization results for very small systems, reporting a qualitatively good agreement between both methods jaksch .

In this section, we compare a MF approximation with QMC results for the system under consideration. It will be shown that the MF approach not only violates the Mermin-Wagner theorem, leading to long-range antiferromagnetic order in one dimension, as expected, but also introduces spurious structures in the density profiles. In order to obtain the MF Hamiltonian, we rewrite the Hubbard Hamiltonian in Eq. (1) in the following form (up to a constant shift in the chemical potential):

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(c_{i\sigma }^{}c_{i+1\sigma }+\text{H.c.})+V{\displaystyle \underset{i,\sigma }{}}x_i^2\widehat{n}_{i\sigma }`$ (5)
$`+{\displaystyle \frac{\left(1\mathrm{\Lambda }\right)}{2}}[\widehat{n}_i\widehat{n}_i\widehat{n}_i^2+\widehat{\delta }_{n_i}^2]`$
$`{\displaystyle \frac{\mathrm{\Lambda }}{2}}[\widehat{\mu }_i\widehat{\mu }_i\widehat{\mu }_i^2+\widehat{\delta }_{\mu _i}^2],`$
where the density and the magnetization in each site are denoted by $`\widehat{n}_i=\widehat{n}_i+\widehat{n}_i`$ and $`\widehat{\mu }_i=\widehat{n}_i\widehat{n}_i`$, respectively. The fluctuations of the density and magnetization are given by $`\widehat{\delta }_{n_i}=\widehat{n}_i\widehat{n}_i`$ and $`\widehat{\delta }_{\mu _i}=\widehat{\mu }_i\widehat{\mu }_i`$, respectively, and $`\mathrm{\Lambda }`$ is an arbitrary parameter that was introduced in order to allow for the most general variation in parameter space. The following relations for fermionic operators were used:

$`\widehat{n}_i\widehat{n}_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}\widehat{n}_i^2{\displaystyle \frac{1}{2}}\widehat{n}_i`$
$`\widehat{n}_i\widehat{n}_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}\widehat{\mu }_i^2+{\displaystyle \frac{1}{2}}\widehat{n}_i.`$ (6)
Neglecting the terms containing the square of the density and magnetization fluctuations in Eq. (5), the following MF Hamiltonian

$`H_{MF}`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(c_{i\sigma }^{}c_{i+1\sigma }+\text{H.c.})+V{\displaystyle \underset{i,\sigma }{}}x_i^2\widehat{n}_{i\sigma }`$ (7)
$`+{\displaystyle \frac{1}{2}}[\left(1\mathrm{\Lambda }\right)\widehat{n}_i\mathrm{\Lambda }\widehat{\mu }_i]\widehat{n}_i`$
$`+{\displaystyle \frac{1}{2}}[\left(1\mathrm{\Lambda }\right)\widehat{n}_i+\mathrm{\Lambda }\widehat{\mu }_i]\widehat{n}_i`$
$`{\displaystyle \frac{\left(1\mathrm{\Lambda }\right)}{2}}\widehat{n}_i^2+{\displaystyle \frac{\mathrm{\Lambda }}{2}}\widehat{\mu }_i^2`$
is obtained. Due to the inhomogeneity of the trapped system, an unrestricted Hartree-Fock scheme is used to determine the local densities $`\widehat{n}_i`$ and local magnetic moments $`\widehat{\mu }_i`$. Given a MF ground state $`|\mathrm{\Psi }_{MF}`$, the minimum energy $`E=\mathrm{\Psi }_{MF}|H|\mathrm{\Psi }_{MF}`$ (not the MF one $`E_{MF}=\mathrm{\Psi }_{MF}|H_{MF}|\mathrm{\Psi }_{MF}`$) is reached when $`\mathrm{\Lambda }=0.5`$, so that all the results that follow were obtained for this value of $`\mathrm{\Lambda }`$.

We first recall some of the discrepancies between MF approximations and the known valid facts for the 1D homogenous system ($`V=0`$) as follows.

(i) At half filling it is known that the Hubbard model exhibits a Mott-insulating phase with a charge gap while the spin sector remains gapless, so that the density-density correlations decay exponentially and the spin-spin correlations decay as a power law. Within the MF approximation, there is a band insulator in the system at half filling so that both the density-density and spin-spin correlations decay exponentially. The MF value of the charge gap is an overestimation of the real one, as can be seen in Fig. 10 where we present the MF results for the density-density correlations as function of the distance, together with the QMC results.

The slopes of the curves are proportional to the charge gap and the MF slope is approximately twice the QMC one (the results were obtained for a system with 102 sites and $`U=6t`$). The correlation length of $`\chi `$ is the inverse of these slopes. Finally, at half filling the MF theory leads to an antiferromagnetic state, while in the Hubbard model the magnetization is always zero although quasi-long-range antiferromagnetic correlations appear in the system.

(ii) At any noncommensurate filling, the Hubbard model in 1D describes a metal (Luttinger liquid), so that no gap appears in the charge and spin excitations, and there is a $`2k_F`$ modulation in the spin-spin correlation function that leads to the well-known $`2k_F`$ singularity. Within the MF approximation, there is always a band insulator at any noncommensurate filling, with a gap that decreases when the system departs from half filling. The appearance of this gap for any density is due to the perfect nesting present in one dimension. The insulating nature of these solutions can be seen in the global compressibility of the system that is always zero, in the behavior of the density-density and the spin-spin correlations which decay exponentially, and in the momentum distribution function where there is no Fermi surface, as shown in Fig. 11 for $`N=102`$, $`N_f=66`$, and $`U=6t`$.

Within this MF approach, there is a $`2k_F`$ modulation of the magnetization only in the $`z`$ component \[the $`SU(2)`$ symmetry was broken\], which leads to a divergence of the Fourier transform of $`S^zS^z`$ at $`k=2k_F`$. In Fig. 12 we compare the MF result (a) for $`S^zS^z_k`$ with the QMC one (b), for $`N=102`$, $`N_f=66`$, and $`U=6t`$, where it can be seen that the $`2k_F`$ peak is one order of magnitude bigger in the MF case compared to the QMC case. This is due to the existence of the magnetization in the MF solution (the values at the peaks will only diverge in the thermodynamic limit).

We consider next the confined inhomogeneous case. In Fig. 13(a) we show the density profiles obtained for a trap with $`N=100`$, $`N_f=70`$, and $`Va^2=0.0025t`$ (like the one presented in Fig. 3) for different values of the on-site repulsion $`U/t=2,3,4,6`$ \[for more details see Figs. 15(e)-15(h)\]. It can be seen that the MF solutions have also density profiles in which “metallic” ($`n1`$) and insulating ($`n=1`$) phases coexist. However, the MF insulating plateaus with $`n=1`$ \[$`U/t=4`$ in Fig. 13(a)\] appear for smaller values of $`U`$ than the ones required in the Hubbard model for the formation of the Mott-insulating plateaus \[$`U/t=6`$ in Fig. 3(a)\]. In the “metallic” phases it is possible to see that the charge density shows rapid spatial variations that are very large for $`n>1,U/t=4`$ \[see also Fig. 15(f)\]. These density variations in the “metallic” phases are reflected in the variance profiles \[Fig. 13(b)\]. The values of the variance in the plateaus are the same as the ones obtained in the homogeneous MF case for the same value of $`U`$ when $`n=1`$ \[Fig. 13 (b)\]. This feature was discussed in Sec. II for the QMC solutions.

Since the MF solution leads to an insulating state for incommensurate fillings in the homogeneous case, we analyze in the following more carefully the regions with $`n1`$ in the presence of the trap. We find that the local compressibility always decays exponentially as a function of $`\mathrm{}`$ over the entire system, although the exponents are different depending on the density of the point analyzed. Some additional modulation appears in $`\kappa ^{\mathrm{}}`$ when $`n1`$, as shown in Fig. 14 where the local compressibility is displayed as a function of $`\mathrm{}`$ for $`n=1`$ and $`n0.66`$ in the trapped and homogeneous cases. Therefore, although the MF approximation leads to a density profile similar to the one obtained in QMC as long as $`n<1`$, it gives a qualitative wrong description of the character of those regions.

The MF results for the normalized momentum distribution function are presented in Figs. 15(a)-15(d) with their corresponding density profiles \[Figs. 15(e)-15(h)\], for the same parameter values as in Fig. 13. It can be seen that the normalized momentum profiles are similar to the ones obtained with QMC when the density profiles are similar, so that the MF results are very similar to the QMC ones for this averaged quantity although the physical situation is very different. (Within MF there is always an insulator in the system and within QMC there are metallic and insulating phases coexisting.) In Figs. 15(e)-15(h), we also show the $`S^z`$ component of the spin on each site of the trap \[$`S^z=\left(\widehat{n}_i\widehat{n}_i\right)/2`$\]. When the density is around $`n=1`$, it can be seen that antiferromagnetic order appears, and that a different modulation exists for $`n1`$. We should mention that we have also found some MF solutions for the trapped system where there was formation of spin domain walls in the insulating regions with $`n=1`$ and there the variance was not constant (a plateau in the density does not correspond to a plateau in the variance), so that care should be taken with the mean field solutions.

## V Phase Diagram

In the present section, we study the phase diagram for fermions confined in harmonic traps with an underlying lattice.

As it can be inferred from Fig. 1, the phases present in the system are very sensitive to the values of the parameters in the model, i.e., to the curvature of the parabolic potential ($`V`$), to the number of fermions present in the trap ($`N_f`$) and to the strength of the on-site repulsion ($`U`$). In order to be able to relate systems with different values of these parameters and different sizes, we use the characteristic density ($`\stackrel{~}{\rho }`$) defined in Sec. III for the definition of the phase diagram. This quantity has been shown to be a meaningful quantity in the thermodynamic limit since density profiles (as function of $`x/\zeta `$) and normalized momentum distributions do not change when $`\stackrel{~}{\rho }`$ is kept constant and the filling in the system (or the occupied system size) is increased. For a harmonic potential, the characteristic length is given by $`\zeta =\left(V/t\right)^{1/2}`$ and the characteristic density is then $`\stackrel{~}{\rho }=N_fa\left(V/t\right)^{1/2}`$. In Fig. 16 we show the phase diagrams of two systems with different strengths of the confining potential ($`Va^2=0.006t`$ and $`Va^2=0.002t`$) and different sizes ($`N=100`$ and $`N=150`$, respectively). This shows that the characteristic density is a meaningful quantity to characterize a generic phase diagram. It allows to compare different systems, and, hence, to relate the results of numerical simulations with larger experimental sizes (although currently the linear dimension of the experiments is $`65a`$).

The phases present in the phase diagram in Fig. 16 can be described as follows: (A) a pure metal without insulating regions, (B) a Mott insulator in the middle of the trap always surrounded by a metal, (C) a metallic region with $`n>1`$ in the center of the Mott-insulating phase, that is surrounded by a metal with $`n<1`$, (D) an insulator with $`n=2`$ in the middle of the trap surrounded only by a metallic phase, and finally (E) a “band insulator” in the middle of the trap surrounded by metal with $`n>1`$, and these two phases inside a Mott insulator that as always is surrounded by a metal with $`n<1`$. As it can be seen, this phase diagram is more complex than the one for the 1D homogeneous Hubbard model.

There are some features that we find interesting in this phase diagram: (i) For all the values of $`U`$ that we have simulated, we see that the lower boundary of phase B with the phase A has a constant value of the characteristic density. (ii) We have also seen that the upper boundary of phase B with phases A and C and the boundary of phases D and E with A and C, respectively, are linear within our errors. (iii) Finally, there is another characteristic of phase B, and of the metallic phase A that is below phase B, that we find intriguing and could be related to point (i). As can be seen in Fig. 3(a) for $`U/t=6`$ and 8, once the Mott-insulating phase is formed in the middle of the trap, a further increase of the on-site repulsion leaves the density profile almost unchanged. (The two mentioned profiles are one on top of the other in that figure.) However, if we look at the variance of these profiles it decreases with increasing $`U`$, which means that the double occupancy $`n_{}n_{}=\left[\mathrm{\Delta }n\left(1n\right)\right]/2`$ is decreasing in the system when $`U`$ is increased, without a redistribution of the density. Less expected for us is that this also occurs for the metallic phase A that is below phase B. In Fig. 17(a) we show four density profiles and in Fig. 17(b) their variances, for a trap with $`Va^2=0.006t`$ and $`N_f=30`$ when the on-site repulsion has values $`U/t=2`$, 4, 6, and 8. If we look at the phase diagram, the point with $`U/t=2`$ is not below phase B whereas the other three values are. In Fig. 17(a) we can see that as the value of $`U/t`$ is increased from 2 to 4 there is a visible change in the density profile, but after a further increase the profile remains almost the same for the other values of $`U`$. The inset shows a magnification of the top of the profiles for a more detailed view. Fig. 17(b) shows that although the density stays almost constant for $`U/t=4,6`$, and 8, the variance decreases. In the other phases (C, D, and E), an increase of the on-site repulsion always changes the local densities pushing particles to the edges of the trap.

Finally, in order to see, whether the features of the phase diagram discussed previously are particular of a harmonic confining potential, we have also performed an analogous study for a quartic confining potential. In Fig. 18 we show how the density profiles evolve in a trap with a quartic confining potential ($`V_4a^4=3\times 10^6t`$) when the total filling is increased from $`N_f=20`$ to 140. It can be seen that the shape of the metallic regions is flatter than in the harmonic case. However, all the features discussed previously for the parabolic confining potential are also present in the quartic case. Local phases appear in the system, the Mott-insulating plateaus with $`n=1`$ have values of the variance equal to the ones in the homogeneous case for the same value of $`U`$, the local compressibility vanishes in the Mott-insulating regions and the spin and momentum distribution exhibit a behavior similar to that of the parabolic case.

In addition, the form of the phase diagram for the quartic confining potential is similar to that of the harmonic trap. In Fig. 19 we show the phase diagram for a system with $`N=100`$ and $`V_4a^4=3\times 10^6t`$, where the characteristic density is given by $`\stackrel{~}{\rho }=N_fa\left(V_4/t\right)^{1/4}`$. We have also calculated some points of the phase diagram for another system with $`V_4a^4=10^6t`$ and $`N=100`$ in order to check that the scaling relation also works properly in this case. The phases are labeled in the same way as in Fig. 16, up to the largest value of $`U`$ that we have studied we do not obtain the phase E which for the quartic confining case seems to be moved further towards the strong coupling regime. Figure 19 shows that the phase diagram in this case is very similar to that of the parabolic case, so that all the features discussed previously apply to the quartic confining potential. The form of the phase diagram in Fig. 16 seems to be of general applicability for any confining potential when the proper characteristic density is considered.

## VI Conclusions

We have studied the ground-state properties of the 1D fermionic Hubbard model in harmonic traps with an underlying lattice, using QMC simulations and a MF approach.

Due to the inhomogeneous density distribution, in general metallic and insulating regions coexist in the trap. Therefore, local quantities have to be used to characterize the system. We considered first the variance of the local density, which, when a plateau develops in the density profile with $`n=1`$, equals the value in the homogeneous system in the Mott-insulating phase. However, it is not an unambiguous quantity, since it can be larger in the Mott-insulating phases than in metallic regions with densities close to $`n=0`$ or $`n=2`$. Therefore, we have defined a local-order parameter (local compressibility) that vanishes in the insulating regions and is finite in the metallic ones. The local compressibility gives a measure of the local change in density due to a constant shift of the potential over finite distances larger than the correlation length of the density-density correlation function for the Mott-insulating phase in the unconfined system. As expected, antiferromagnetic correlations are present in the Mott-insulating phases and they remain in the shoulders of the density profile for $`n1`$.

As opposed to homogeneous or periodic systems, where the appearance of a gap can in general be seen by the disappearance of the Fermi edge in the momentum distribution function, it is easily seen that due to the presence of a confining potential, such a feature is much less evident. Although increasing the power of the confining potential sharpens the features connected with the Fermi edge, it remains qualitatively different from the homogeneous case. Increasing the system size with a proper definition of the thermodynamic limit, does not change at all the behavior of the density and momentum distribution functions in the trapped case, such that, in fermionic systems it does not seem to be the appropriate quantity to look for evidence of gaps opening in the system.

The MF study has shown that within this approach the trapped system is an insulator for all the densities in the trap, a qualitatively wrong picture in contradiction with the QMC results. In addition, the metallic regions showed large oscillations in the density that are not present in the real solutions. However, it is possible to obtain some qualitative information from MF, such as the coexistence of $`n1`$ regions with $`n=1`$ plateaus, the shape of the momentum distribution function and the existence of antiferromagnetic correlations in the insulating regions with $`n=1`$.

Finally, we have determined a generic form for the phase diagram that allows to compare systems with different values of all the parameters involved in the model. It can be used to predict the phases that will be present in future experimental results. The phase diagram also reveals interesting features such as reentrant behavior in some phases when some parameters are changed and phase boundaries with linear forms. A similar form of the phase diagram has been also found for a quartic confining potential, such that the results obtained here are generally applicable to cases beyond a perfect harmonic potential.

###### Acknowledgements.

We gratefully acknowledge financial support from the LFSP Nanomaterialien. We are grateful to G. G. Batrouni and R. T. Scalettar for interesting discussions at the early stages of the project. We are grateful to F. F. Assaad for valuable discussions about the MF approach we used here and also for providing us the core of the code we employed in our MF calculations. We thank R. Noack and G. Modugno for useful comments on the manuscript, and M. Feldbacher and T. Pfau for useful conversations. We thank HLR-Stuttgart (Project DynMet) for allocation of computer time, the calculations were carried out on the HITACHI SR8000.


format: text

text:
# Experimental observation of superluminal group velocities in bulk two-dimensional photonic bandgap crystals

## I Introduction

The superluminal propagation of wave packets with faster-than-$`c`$, infinite and negative group velocities has been observed in a wide range of physical systems, including both passive and active optical media Chiao1997 ; Nimtz1997 ; Boyd2002 ; Rochester ; Wang2001 . Garrett and McCumber Garrett1970  first predicted that a Gaussian wave packet can propagate with negligible distortion through a linear optical medium with a superluminal group velocity, if its bandwidth is restricted to a narrow spectral region of anomalous dispersion within an absorption or gain line. This prediction was later experimentally verified for the propagation of picosecond laser pulses through a GaP:N sample Chu1982 . Such superluminal behavior does not violate relativistic causality. Any analytic function, such as a Gaussian wave packet, contains sufficient information in its early tail to reconstruct the entire wave packet with a superluminal pulse advancement, but little distortion. Relativistic causality forbids only the *front* velocity, i.e., the velocity of any *discontinuity*, from exceeding the vacuum speed of light, $`c`$. It does not forbid the *group* velocity of a wavepacket from being superluminal.

The superluminal behavior of electromagnetic pulses can occur not only in passive, *absorptive* optical media, but also in passive, *dissipationless* optical media, such as periodic dielectric photonic bandgap (PBG) structures Yablonovitch . Optical pulses whose bandwidths are restricted to the PBG are exponentially damped not due to absorption, but interference between multiple Bragg reflections from planes within the transparent periodic dielectric structure. The superluminal propagation of wavepackets through 1D periodic dieletric structures has been observed in the optical spectral region Steinberg . Similarly, superluminal propagation has been observed for classical microwave pulses in a 1D periodic dielectric stack placed inside a microwave waveguide Nimtz1994 . Other work has demonstrated superluminal behavior in the frequency range of tens of MHz in a 1D photonic structure made of a sequence of coaxial cables with alternating impedance values Hache2002 . No reported experiment has yet demonstrated superluminal propagation of electromagnetic wavepackets in 2D periodic structures, although we have recently demonstrated for the first time clear experimental evidence for the existence of exponential wave decay in such structures at microwave frequencies Hickmann2002 . Here we report on experimental results which imply that superluminal group velocities do occur in such structures.

Microwave measurements have certain advantages which we exploit here. First, one can perform direct, precise phase measurements using a microwave vector network analyzer (VNA). Second, the wavelength is sufficiently large so that high-precision, well-fabricated 2D photonic crystals can be relatively easily constructed. Third, one can more easily tailor the dispersive properties of the PBG by changing structural parameters, such as the structure air-filling fraction (AFF).

## II Experiment

The fabrication of our photonic crystals has been described in detail elsewhere Hickmann2002 . In summary, we stack acrylic rods in a hexagonal lattice and glue them in place with an acrylic solvent cement. We constructed and tested crystals with as many as eighteen layers of rods, and AFFs of 0.60, using rods with a 1/2-inch outer diameter. The crystals are “bulk” samples rather than “slabs” in the sense that their extruded (non-periodic) dimension is many times the lattice spacing in length.

In Fig. 1 we show a schematic of the experimental apparatus used for the microwave transmission measurements. We used an HP 8720A VNA, which is capable of measuring both the transmission and phase delay of microwaves detected by a polarization-sensitive receiver antenna (horn). The VNA was swept from 8 to 14 GHz in 15 MHz increments and was connected to identical transmitter and receiver horns. The crystal and receiver horn were placed 1.6 m away in a box measuring 60 cm on a side, with microwave-absorbing walls. Microwaves entered the box through a 14 cm $`\times `$ 17 cm rectangular aperture. These dimensions were chosen to minimize diffraction effects through the aperture without allowing leakage around the crystal. With the crystal removed and the rectangular aperture closed, the microwave signal at the receiver was suppressed by more than 45 dB, indicating good shielding by the box. As previously reported, our bulk crystal with AFF equal to 0.60 displays a strong bandgap at $``$ 11 GHz Hickmann2002 .

To establish the reliability of the phase delay measurements, we placed sheets of ordinary acrylic of varying widths in front of the receiver horn. The phase delay relative to that of free space propagation was linear in the transmitted microwave frequency, with a slope that depended on the thickness of the acrylic. Using a simple model of phase delay given by

$$\mathrm{\Delta }\varphi (\omega )=\frac{\omega }{c}d(n(\omega )1),$$
(1)
where $`d`$ is the sheet thickness, we were able to extract a value of $`n=1.58`$ for the index of refraction of acrylic, in good agreement with the accepted value of $`1.61`$ for this frequency range Gray1972 .

The group velocity measurements were carried out by measuring, for both polarizations, the spectral dependence of the phase delay for crystals with different numbers of layers. The group index was calculated according to the relation Garrett1970

$$n_g(\omega _0)=n(\omega _0)+\omega _0\frac{n}{\omega }|_{\omega _0},$$
(2)
where $`n(\omega )`$ is the normal (phase) index of refraction and the group index $`n_g`$ is defined as the ratio of the speed of light to the group velocity. Thus, superluminal group velocities exist in any spectral region where the group index is less than unity.

In order to calculate the group index, it was necessary to unwrap the phase data obtained from the VNA, since it is sensitive to phase modulo $`2\pi `$. We were able to determine when such phase slips occurred by observing that the unwrapped phase delay at fixed frequencies was a monotonic function of the number of layers except for sharp changes (phase slips) at certain numbers. We adjusted the data by adding $`2\pi \times m`$ to each phase delay spectrum, where $`m`$ is the cumulative number of sharp drops observed over smaller layer numbers. We also confirmed the validity of this phase adjustment by considering $`n/\omega `$, which should go to zero far from the bandgap. We found that it did so if and only if we applied this phase delay adjustment as described.

## III Results and Discussion

In Fig. 2, we show the measured transmission bandgap of a crystal with eighteen layers and AFF equal to 0.60, for both incident microwave polarizations. We consider the plane of incidence to be the plane of periodicity (perpendicular to the rods) and define the TE (TM) polarization to have its electric field perpendicular (parallel) to this plane. We normalized transmission measurements by removing the crystal and measuring the total microwave signal at the receiving horn. We find that this crystal exhibits anomalous dispersion within the bandgap, passing briefly through zero dispersion at the bandgap edges, for both TE and TM polarizations (see Fig. 3). We note that the region of anomalous dispersion becomes increasingly well-defined as the number of layers is increased.

From the phase delay data in Fig. 3 we calculate the group index using Eq. 2. The results are displayed in Fig. 4. The region of superluminal $`n_g`$ extends over a wider spectral band for the TM polarization than the TE polarization, and deepens more rapidly as the number of layers is increased. This result is not surprising given that the TM gap is known to form more rapidly than the TE gap as the number of crystal layers is increased Hickmann2002 .

We did similar measurements (not shown here) using other crystals with AFF equal to 0.32 and obtained qualitatively similar bandgaps, although they are narrower and shallower for both polarizations, as expected Broeng1998 , and have slightly different center frequencies. Both crystal sets display similar superluminal spectral regions, although we found that those with higher AFF typically have greater bandwidth and yield superluminal group velocities which are more superluminal.

## IV Conclusion

We have experimentally demonstrated that superluminal including infinite group velocities can exist for analytic signals whose spectral bandwidth lies within the bandgap of a two-dimensional hexagonal photonic crystal. Superluminal phenomena are exhibited by these crystals for both tested AFFs, although the effects are more significant (i.e., group velocities become more superluminal and exist over a wider spectral range) for structures of higher AFF.

This work was supported by ARO grant number DAAD19-02-1-0276. We thank the UC Berkeley Astronomy Department, in particular Dr. R. Plambeck, for lending us the VNA. JMH thanks the support from Instituto do Milênio de Informação Quântica, CAPES, CNPq, FAPEAL, PRONEX-NEON, ANP-CTPETRO.


text:
# BUBBLES IN PLANETARY NEBULAE AND CLUSTERS OF GALAXIES: JET PROPERTIES

## 1 INTRODUCTION

Chandra X-ray observations of clusters of galaxies reveal the presence of X-ray-deficient bubbles in the inner regions of many (cooling flow) clusters, e.g., Hydra A (McNamara et al. 2000), Perseus (Fabian et al. 2000, 2003), A 2597 (McNamara et al. 2001), RBS797 (Schindler et al. 2001), Abell 4059 (Heinz et al. 2002), and Abell 2052, (Blanton et al. 2003). A nice pair of bubbles is seen in the poor cluster HCG 62 (Vrtilek et al. 2002). X-ray-deficient pairs of bubbles, although less prominent, exist also in elliptical galaxies; e.g., in M84 (Finoguenov & Jones 2001) and NGC 4125 and NGC 4552 (White & Davis 2003). The low X-ray emissivity implies low density inside the bubbles, while the absence of evidence of strong shocks suggests that the bubbles are expanding and moving at subsonic or mildly transonic velocities (Fabian et al. 2000, 2003; McNamara et al. 2000; Blanton et al. 2001). In paper-I (Soker 2003) I pointed out an interesting and not trivial similarity in the morphology and some non-dimensional quantities between pairs of X-ray-deficient bubbles in clusters of galaxies and pairs of optical-deficient bubbles in planetary nebulae (PNs). Examples of PNs with nice pairs of bubbles are the Owl nebula (NGC 3587; PN G148.4+57.0: Guerrero et al. 2003), Cn 3-1 (VV 171; PN G038.2+12.0; Sahai 2000), and Hu 2-1 (PN G051.4+09.6), which possesses two prominent pairs of bubbles, one pair closer to the center and the other farther out, with inclination between the two symmetry axes (Miranda et al. 2001). Many other PNs also possess pairs of low emissivity bubbles, although less prominent, e.g., NGC 2242 (PN G170.3+15.8: Manchado et al. 1996), and a pair of bubbles within bipolar lobes in M 2-46 (PN G024.8-02.7: Manchado et al. 1996). More details of the morphological similarities and the similarity in some non-dimensional variables between pair of bubbles in these two classes of objects is given in paper I (see Table 1 there).

The similarity in morphology and in some relevant non-dimensional variables led me to postulate in paper I a similar formation mechanism. It is commonly accepted that bubbles in clusters are blown by AGN jets (e.g., Brighenti & Mathews 2002; Brüggen 2003; Brüggen et al. 2002; Fabian et al. 2002; Nulsen et al. 2002; Quilis, Bower, & Balogh 2001; Reynolds, Heinz, & Begelman 2001; Soker, Blanton, & Sarazin 2002). This postulate, then, was used to strengthen models for PN shaping by jets (or collimated fast winds: CFW). In particular, the presence of dense material in the equatorial plane observed in the two classes of bubbles constrains the jets and CFW activity in PNs to occur while the AGB star still blows its dense wind, or very shortly after. I argued there that only a stellar companion can account for such fast jets (or CFW). Very recent Chandra observations hint at more similarities between the jets in PNs and AGN jets. Kastner et al. (2003) point out that the blobby X-ray appearance of the jets in the symbiotic bipolar nebula Mz 3, probably to become a PN, is similar somewhat to the blobby X-ray nature of some AGN jets, e.g., Centaurus A (Hardcastle et al. 2003). The similarity is not only of fat bubbles in cluster and bipolar PNs. The general X-ray image of the cooling flow cluster Cygnus A (Smith et al. 2002; see their fig. 2) resembles that of the optical image of many elliptical PNs, e.g., NGC 6826 (Balick 1987; Balick et al. 1998). Cygnus A hot spots, i.e., the heads of the two jets, resemble the structure of pairs of ansae$``$two opposite dense blobs along the major axis of some elliptical PNs$``$suggesting that the ansae are formed by jets (see for example the ansae of the PN NGC 7009; Balick 1987; Balick et al. 1998; Sabbadin et al. 2003).

These similarities motivate me to consider a unified scheme for the formation of bubbles in clusters and PNs, based on the inflation of pairs of bubbles by oppositely ejected jets. Although commonly accepted by people studying clusters, the jet-inflated bubble model is still in dispute among people studying PNs. In particular, it is not widely accepted that jets which shape bipolar PNs$``$those with two large lobes and an equatorial waist between them$``$are blown by a binary companion. As argued in several papers (Paper I and references therein), more and more observations, including the similarity in morphology discussed above, point to jets blown by a binary companion in bipolar PNs (but not in all PNs).

In the present paper I study in more detail the conditions under which, more or less, spherical bubbles are formed, as opposed to the case where a propagating jet inflates a long narrow region. A ‘fat’ low-density inflated region, more or less spherical but not necessarily exactly so, is termed in the present paper a fat bubble (e.g., as in the Owl nebula), while an elongated low density inflated region is termed lobe, e.g., as in OH 231.8+4.2 (Bujarrabal et al. 2002), and He2-115 (Sahai & Trauger 1998).

## 2 JET PROPAGATION IN PLANETARY NEBULAE

The flow structure is of a jet flowing into the dense circumstellar material, which is the expanding AGB wind. The speed of the material inside the jet is $`v_j`$, while the jets’ head propagates at speed $`v_h<v_j`$. The slow wind speed is $`v_s<v_h`$. The aim here is to examine the conditions for the formation of fat bubbles. The condition is that the expansion speed of the bubble surface relative to the slow wind $`v_b`$ be faster than the propagation speed of the jet’s head relative to the slow wind

$$v_{hs}v_hv_s<v_b\mathrm{for}\mathrm{spherical}\mathrm{bubble}\mathrm{formation}.$$
(1)
If, on the other hand, $`v_{hs}>v_b`$, the expanding bubble formed by the jet lags behind the jet’s head, and an elongated lobe, rather than a fat bubble, is formed (Soker 2002; Lee & Sahai 2003). The condition in equation (1) implies that the jet is slow-propagating, with its velocity relative to the slow wind speed given by momentum conservation, which here reads (see eq. 1 in Soker 2002)

$$v_{hs}\left(\frac{\dot{M}_jv_j}{\beta \dot{M}_sv_s}\right)^{1/2}v_s,\mathrm{for}\rho _s\rho _j.$$
(2)
where $`\beta `$ is defined such that the jet expands into a solid angle of $`\mathrm{\Omega }_j=4\pi \beta `$, and $`\rho _s(r)`$ and $`\rho _j(r)`$ are the density of the slow wind and jet at a distance $`r`$ from the center, respectively. $`\dot{M}_s`$ and $`\dot{M}_j`$ are the mass loss rates into the slow wind and one jet, respectively; both are constant, and defined positively.

The jet deposits kinetic energy in a power of $`0.5\dot{M}_jv_j^2`$ into the bubble. The radius of the bubble as a function of time $`t`$, can be approximated by using the results of Castor, McCray, & Weaver (1975),

$$R_b=0.76\left(\frac{0.5\dot{M}_jv_j^2}{\rho _s}\right)^{1/5}t^{3/5},$$
(3)
where $`\rho _s=\dot{M}_s/(4\pi v_sr^2)`$. The expansion velocity of the bubble’s surface relative to the slow wind is

$`v_b=0.6R_b/t=18\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{2/5}\left({\displaystyle \frac{\dot{M}_j}{0.01\dot{M}_s}}\right)^{1/5}\left({\displaystyle \frac{v_s}{10\mathrm{km}\mathrm{s}^1}}\right)^{1/5}`$
$`\times \left({\displaystyle \frac{z}{10^{17}\mathrm{cm}}}\right)^{2/5}\left({\displaystyle \frac{t}{1000\mathrm{yr}}}\right)^{2/5}\mathrm{km}\mathrm{s}^1,`$ (4)
where the expression for $`\rho _s`$ was used, and $`z`$ is the coordinate along the propagation direction of the jet, with the origin at the origin of the jet. With the parameters used in the last equation, equation (2) gives $`v_h10[1+2(\beta /0.1)^{1/2}]\mathrm{km}\mathrm{s}^130\mathrm{km}\mathrm{s}^1`$, where $`\beta =0.1`$ is for a jet with a half opening angle (i.e., measured from its symmetry axis) of $`\alpha =37^{}`$. The distance of the jet’s head from the source at time $`t`$ is given by

$$z_h10^{17}\left[\frac{1}{3}+\frac{2}{3}\left(\frac{\beta }{0.1}\right)^{1/2}\right]\left(\frac{t}{1000\mathrm{yr}}\right)\mathrm{cm}.$$
(5)
This is the reason for the scaling of $`z`$ in equation (4). Note that what matters is the ratio $`z/t`$, which changes slowly, if at all, with time. Substituting $`v_{hs}`$ from equation (2), and $`v_b`$ from equation (4) in condition (1) for inflating a fat bubble, gives the condition of spherical bubble formation in the form

$`1{\displaystyle \frac{v_b}{v_{hs}}}=0.9\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{1/10}\left({\displaystyle \frac{\dot{M}_j}{0.01\dot{M}_s}}\right)^{3/10}\left({\displaystyle \frac{\beta }{0.1}}\right)^{1/2}\left({\displaystyle \frac{v_s}{10\mathrm{km}\mathrm{s}^1}}\right)^{3/10}`$
$`\times \left({\displaystyle \frac{z}{10^{17}\mathrm{cm}}}\right)^{2/5}\left({\displaystyle \frac{t}{1000\mathrm{yr}}}\right)^{2/5}\mathrm{km}\mathrm{s}^1.`$ (6)
From equation (5) the ratio $`z/t`$ does not change much in different relevant cases. Also, the slow wind speed is $`v_s10\mathrm{km}\mathrm{s}^1`$ for all relevant AGB stars. The condition for fat-bubble inflation becomes mainly a condition on the jet’s properties

$`\beta 0.12\left({\displaystyle \frac{v_j}{400\mathrm{km}\mathrm{s}^1}}\right)^{1/5}\left({\displaystyle \frac{\dot{M}_j}{0.01\dot{M}_s}}\right)^{3/5}.`$ (7)
Namely, more or less spherical fat bubbles will be inflated by non-well collimated jets. For these jets, a more appropriate term is CFW (collimated fast wind).

For small angles $`\beta \alpha ^2/4`$, and the last condition can be written for the half opening angle of the jet

$$\alpha 40^{}\left(\frac{v_j}{400\mathrm{km}\mathrm{s}^1}\right)^{1/10}\left(\frac{\dot{M}_j}{0.01\dot{M}_s}\right)^{3/10}.$$
(8)
The scaling of the jet (or CFW) properties used here is appropriate for a jet blown by a main sequence companion, because the jet’s speed is of the order of the escape velocity. The mass loss rate into the jets$``$two oppositely blown jets$``$is quite high, e.g., a companion accreting a fraction of $`0.2`$ of the AGB wind, and blows a fraction of $`0.1`$ of the accreted mass into the two jets. In this case, the last condition implies that the jets, or more appropriately a CFW, will be non-well collimated. This is reasonable for a close main sequence companion, since the accretion disk will not be extended. The main sequence companion can resides at larger distances, but not too large $`a20\mathrm{AU}`$, otherwise the jet will be weak and the bubble will not be inflated much.

At the other extreme the jets may be well collimated, and their mass deposition rate very low $`\dot{M}_j0.01\dot{M}_s`$. For the bubbles to be fat, the jet’s speed must be high. For example, for $`\dot{M}_j=10^4\dot{M}_s`$ and $`v_j=5000\mathrm{km}\mathrm{s}^1`$, condition (8) reads $`\alpha 10^{}`$. These parameters are appropriate for a widely separated, $`a50\mathrm{AU}`$, white dwarf companion.

For a large pair of bubbles to be observed in PNs there are two other conditions, in addition to the condition derived here. First, the jet should be strong enough to inflate a large bubble that expands fast enough, i.e., faster than the expansion speed of the slow wind,

$$v_{hs}v_s.$$
(9)
The condition on the jet’s properties implied by this inequality can be derived from equation (4). Second, the jet’s material should not cool fast. This condition is that the bubble starts to expand close to the center, $`z10^{17}\mathrm{cm}`$. The location where the bubble starts to expand is given by equation (6) in Soker (2002). Basically, keeping the other parameters as used in this section, it constrains the jet’s speed to

$$v_j150\mathrm{km}\mathrm{s}^1.$$
(10)
Finally, I note that the conditions on the jet’s properties given in equations (8) (or 6), (9) and (10), were derived for the flow occurring before ionization and the spherical fast wind from the central post-AGB star start. These are likely to erode the bubble; in particular, after ionization the thermal pressure of the slow wind material can’t be neglected any longer. Therefore, the real constraints on jets to inflate observable fat bubbles in PNs are somewhat stronger than those given by these equations.

## 3 JET PROPAGATION IN CLUSTERS

There are several significant differences between bubble inflation in PNs and clusters. (1) In clusters the thermal pressure of the ambient gas is non-negligible. (2) In clusters the ambient medium does not flow outward. (3) The inflating jets in clusters may be relativistic, and the magnetic pressure inside the bubble can be large. (4) In clusters the bubbles can be observed as they form, unlike in PNs, where they are observed long after the jets have ceased (old bubbles may be observed in clusters$``$termed ghost-bubbles$``$as in the Perseus cluster; Fabian et al. 2000). However, as argued in paper I, these don’t prevent a similar bubble-formation mechanism in PNs and clusters.

Following Soker et al. (2002) I neglect the ambient pressure and use the results of Castor et al. (1975), as in the previous section. The expansion velocity of the bubble is given by equation (20) of Soker et al. (2002), which is written here as

$$v_b(\tau _I)720\left(\frac{\tau }{10^7\mathrm{yr}}\right)^{2/5}\left(\frac{\dot{E}_j}{10^{45}\mathrm{erg}\mathrm{s}^1}\right)^{1/5}\left(\frac{\rho _c}{10^{25}\mathrm{g}\mathrm{cm}^3}\right)^{1/5}\mathrm{km}\mathrm{s}^1.$$
(11)
For the ambient density a crude approximation for $`r>10\mathrm{kpc}`$ in A 2052 is (Blanton et al. 2001; Soker et al. 2002, their eq. 16)

$$\rho _c(r)10^{25}\left(\frac{r}{10\mathrm{kpc}}\right)^1\mathrm{g}\mathrm{cm}^3.$$
(12)
Because of the presence of bubbles, it is impossible to know the density prior to bubble inflation in A 2052, so this approximation is more accurate for $`r>30\mathrm{kpc}`$ and very crude for smaller radii. By neglecting the magnetic pressure inside the jet and relativistic effects, hence $`\dot{E}_j=\dot{M}_jv_j^2/2`$, we can take an equation similar to equation (2) in the previous section, and write for speed of the jet’s head through the cluster medium (e.g., Krause 2003)

$`v_h\left({\displaystyle \frac{\dot{E}_j}{2\pi \beta z^2v_j\rho _c}}\right)^{1/2}=1300\left({\displaystyle \frac{\dot{E}_j}{10^{45}\mathrm{erg}\mathrm{s}^1}}\right)^{1/2}\left({\displaystyle \frac{v_j}{10^4\mathrm{km}\mathrm{s}^1}}\right)^{1/2}`$
$`\times \left({\displaystyle \frac{\beta }{0.1}}\right)^{1/2}\left({\displaystyle \frac{\rho _c}{10^{25}\mathrm{gcm}^3}}\right)^{1/2}\left({\displaystyle \frac{z}{10\mathrm{kpc}}}\right)^1\mathrm{km}\mathrm{s}^1`$ (13)
where the symbols have their meaning as in section 2. For the parameters used here, at $`z10\mathrm{kpc}`$ the cluster density to jet density ratio is $`100`$. We scaled the jet velocity to the non-relativistic regime, as Omma et al. (2003) do. The velocities (which depend on the time and location) $`v_b`$ and $`v_h`$ given above, are mildly supersonic, or even subsonic. This means the ambient pressure cannot be neglected for late time. Still, there is justification for neglecting pressure (Soker et al. 2002), in particular in the approximate condition derived below. An accurate expression cannot be derived because of the poor knowledge of the exact density profile in the very inner regions of the cluster. For the scaling used above, the jet reaches a distance of $`r=10\mathrm{kpc}`$ in $`1.4\times 10^7\mathrm{yr}`$. Therefore, the time and distance scaling used above are compatible. Inserting equations (11) and (13) in the condition to inflate a fat bubble $`v_bv_h`$ (eq. 1), I derive a condition on the half opening angle of the jet (similar to that used in deriving eq. 8),

$`\alpha 65^{}\left({\displaystyle \frac{\dot{E}_j}{10^{45}\mathrm{erg}\mathrm{s}^1}}\right)^{3/10}\left({\displaystyle \frac{v_j}{10^4\mathrm{km}\mathrm{s}^1}}\right)^{1/2}\left({\displaystyle \frac{\rho _c}{10^{25}\mathrm{gcm}^3}}\right)^{3/10}\left({\displaystyle \frac{z}{10\mathrm{kpc}}}\right)^1\left({\displaystyle \frac{\tau }{10^7\mathrm{yr}}}\right)^{2/5}.`$ (14)
Because several factors were neglected, as mentioned earlier, the last condition is a crude one. Still, a strong constraint can be deduced from it, namely in order to inflate a fat bubble in cooling flow clusters the CFW should be wide open. Even for $`\dot{E}_j10^{44}\mathrm{erg}\mathrm{s}^1`$ and $`v_j3\times 10^4\mathrm{km}\mathrm{s}^1`$, equation (14) gives $`\alpha 20^{}`$. The jets which inflated the bubbles in the cluster A2052 were probably such jets. (Weak relativistic jets can be narrow though.) This may point to a significant difference between the AGN activity in normal large elliptical galaxies, and the AGN activity of the cD galaxies in cooling flow clusters. In any case, keeping in mind the uncertainties in the derivations of the last two sections, the conditions on the opening angle of the CFW (jets) in PNs and clusters to inflate fat bubbles are similar. This further strengthen the idea of a unified process.

## 4 COMPARISON WITH NUMERICAL SIMULATIONS

Most numerical simulations of jets aim at a specific environment. I could not find in the literature simulations of conical, rather than cylindrical, jets propagating into a denser medium in the parameter-space relevant to the present cases. This is the main motivation for the order-of-magnitude estimates presented in this paper: To encourage simulations with wide-angle jets.

The most relevant simulations in the literature are in a very recent paper by Omma et al. (2003), who conduct numerical simulations of non-relativistic jets propagating into a cooling flow cluster environment. They simulate two cases where they inject cylindrical jets, one case with a jet radius of $`r_{\mathrm{jet}}=2\mathrm{kpc}`$, and the other with $`r_{\mathrm{jet}}=3\mathrm{kpc}`$. Close to the source, where the distance along he jet axis is not much larger than the jet’s diameter $`z2r_{\mathrm{jet}}`$, they obtain an inflated fat-bubble. As the cylindrical jet expands to larger distances, $`z2r_{\mathrm{jet}}`$, the bubble loses its general spherical structure, and it is detached from the center. Such simulations can’t reproduce fat bubbles as in Perseus or A2052, unless the radius of the simulated jets at injection is very large. They do reproduce structures similar to the detached (from the center) bubbles in M 87 (Virgo A; See fig. 15 of Omma et al. 2003) and Hydra A. The wider jet in the simulations of Omma et al. (2003) forms fatter bubbles and proceeds at a lower velocity. This implies that at a specific location in the simulated cluster environment, when the jet radius is larger than some critical radius, a bubble is formed rather than a narrow cocoon. However, eventually the jet proceeds forward faster than the bubble expands (until it slows down when injection ends). To maintain a bubble attached to the center, or close to it, the jet radius should increase, i.e., a conical jet, as was shown in the previous section of the present paper.

A more quantitative comparison is possible. The ambient density, as they take from David et al. (2001), is $`\rho _c10^{25}\mathrm{gcm}^3`$, and their injection power and velocity are $`\dot{E}_j=6\times 10^{43}\mathrm{erg}\mathrm{s}^1`$ and $`V_j=10^4\mathrm{km}\mathrm{s}^1`$, respectively. Substituting these values in equation (14) gives $`\alpha 28^{}`$, for the same scaling of time and distance as in equation (14). For the $`r_{\mathrm{jet}}=2\mathrm{kpc}`$ jet this implies $`z4\mathrm{kpc}`$, while for $`r_{\mathrm{jet}}=3\mathrm{kpc}`$ this implies $`z6\mathrm{kpc}`$. This value of $`r_{\mathrm{jet}}/z0.5`$ explains why the bubbles in Omma et al. (2003) simulations lose their ‘fatty’ structure when the jet reaches a distance of about the jet diameter. Although it is hard to compare the cylindrical-jet simulations of Omma et al. (2003) with the present paper, the above discussion suggests that there is a solid ground for the order of magnitude analysis conducted here.

There are less relevant simulations, which non the less shed some light on the present analysis; for these, only qualitative comparison is meaningful. Krause & Camenzind (2003b) find that as they inject the jet a spherical bubble is formed around the jet. Later the bubble becomes elongated, first to an elliptical shape, then to a narrow extension along the jet. This is compatible with the results of Omma et al. (2003) cited above. This is along the claim of the present paper that a wide-angle jet will form a long-lasting large spherical (fat) bubble. Loken et al. (1995) simulate the formation of wide-angle tailed radio sources, which are though to form when a radio galaxy which blows two jets moves through the intracluster medium. This is quantitatively similar to precessing jets, in that the jets continuously encounters fresh ambient medium. Their three-dimensional simulations clearly show how jets are disrupted, and large bubbles are formed. These simulations support the claims that wide-angle jets or precessing jets can be disrupted more easily and form fat bubbles. Lee et al. (2001) conduct isothermal simulations of jets and winds from young stellar objects. Despite that isothermal flows are less favorable to inflate bubbles than the cases studied here where radiative cooling time is long, Lee et al. (2001) find that the overall width of the jet-driven shell is smaller than that of the wide-angle wind-driven shell. This supports the notion that wider flows form fatter bubbles (shells in their case).

Finally, it should be noted that recollimation of wide-angle jets, such that the opening angle $`\alpha `$ significantly decreases, will not change much the conclusions. To recollimate a continuous jet the ambient medium must impart transverse momentum to the jet’s material. For a wide-angle jet, the transverse speed is not much below the jet’s speed, a factor of $`\mathrm{sin}\alpha `$. Therefore, the ambient medium and the jet’s material will go through strong shocks, heat-up, and form hot low density bubbles, as assumed in the previous sections.

## 5 DISCUSSION AND SUMMARY

Paper I pointed out the similarities in morphology and some non-dimensional quantities between X-ray-deficient bubbles in clusters of galaxies and the optical-deficient bubbles in PNs. The comparison here, as in paper I, is for pairs of fat, almost spherical, bubbles in these two classes of objects. The main goal of the first paper was to point out these similarities; from that I argue that pairs of fat bubbles in PNs are formed by jets, most probably blown by a companion. In the present paper the main goal is to examine some similarities between the jets in these two classes of systems, assuming indeed that jets, or CFW (for collimated fast wind) blow the fat bubbles. The main conclusion is that for such bubbles to be formed, in both classes of systems the opening angle of the jets must be large. In equation (8) for PNs, and equation (14) for clusters, the values for the half opening angle (measured from the symmetry axis) for blowing fat bubbles are given. For such wide-opening angle jets, a CFW is the more appropriate term. Narrow jets will form elongated lobes rather than fat bubbles (Krause 2003; Krause & Camenzind 2003a); some examples of these are given in section 1. Presently, most simulations (e.g, Krause 2003; Krause & Camenzind 2003a, and references in these papers) simulate narrow jets, and examine the influence of the jet to ambient density ratio. The present study suggests that the opening angle is a crucial parameter.

Beside the similarities studied in paper I and here, other similar processes may exist. An example is heat conduction from the hot bubbles to the cooler environment. In the past this process was studied in more detail for PNs (Soker 1994; Zhekov & Myasnikov 2000; Soker & Kastner 2003). This process, among other things, was claimed to be able to heat somewhat the cooler ($`10^4\mathrm{K}`$) environment. The heating will not be isotropic because it is regulated by magnetic fields (Soker 1994; Zhekov & Myasnikov 2000). With the search for heating sources of the radiatively cooling ICM in cooling flow clusters, heat conduction from the hot bubble to the cooler ($`17\times 10^7\mathrm{K}`$) ICM brought this process to the attention of the cooling flow community. The hot spots found around the X-ray depression region in the cooling flow cluster in MKW3s (Mazzotta et al. 2002) may have been heated by such heat conduction, as suggested by some of these authors during the Cooling Flow meeting held in Charlottesville, Virginia, in June 2003.

I summarize by emphasizing again the need to include jets with large opening angles, i.e., half opening angles of $`\alpha 3070^{}`$, in simulating both PNs and cooling flow clusters (and galaxies). Such jets are more appropriately called CFW.

This research was supported in part by the Israel Science Foundation.