text:
# Λ_𝑐⁺/Λ_𝑐⁻ Asymmetry in Hadroproduction from Heavy-Quark Recombination

## Abstract

Asymmetries in the hadroproduction of charm particles directly probe power corrections to the QCD factorization theorems. In this paper, the heavy-quark recombination mechanism, a power correction that explains charm meson asymmetries, is extended to charm baryons. In this mechanism, a light quark participates in the hard scattering that creates a charm quark and they hadronize together into a charm baryon. This provides a natural and economical explanation for the $`\mathrm{\Lambda }_c^+/\mathrm{\Lambda }_c^{}`$ asymmetries measured in $`\pi N`$ and $`pN`$ collisions.

The production of charm particles in fixed-target hadroproduction experiments exhibits large asymmetries that are commonly referred to as the “leading particle effect” Aitala:1996hf ; Adamovich:1998mu ; Aitala:2000rd ; Garcia:2001xj . Charm hadrons that have a valence parton in common with the beam hadron are produced in greater numbers than other charm hadrons in the forward region of the beam. Asymmetries have also been observed in the production of light particles, such as pions and kaons. Asymmetries in charm particles are particularly interesting, because one can exploit the fact that the charm quark mass $`m_c`$ is much larger than $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ to make closer contact with fundamental aspects of Quantum Chromodynamics (QCD). The large mass guarantees that even at small transverse momentum the production process involves short-distance effects that can be treated using perturbative QCD. Furthermore, the nonperturbative long-distance effects of QCD can be organized as an expansion in $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c`$.

There have been many measurements of the asymmetries for charm mesons Aitala:1996hf . Several proposed charm production mechanisms are able to explain these asymmetries by tuning nonperturbative parameters models-M ; models-IC . Recent experiments have also measured the asymmetry for the charm baryon $`\mathrm{\Lambda }_c^+`$ Adamovich:1998mu ; Aitala:2000rd ; Garcia:2001xj , defined by

$`\alpha [\mathrm{\Lambda }_c]={\displaystyle \frac{\sigma [\mathrm{\Lambda }_c^+]\sigma [\mathrm{\Lambda }_c^{}]}{\sigma [\mathrm{\Lambda }_c^+]+\sigma [\mathrm{\Lambda }_c^{}]}}.`$ (1)
The WA89 Adamovich:1998mu  and SELEX Garcia:2001xj  experiments observe a large positive asymmetry for $`\mathrm{\Lambda }_c`$ produced in the forward direction of $`p`$ and $`\mathrm{\Sigma }^{}`$ beams. These asymmetries are consistent with the leading particle effect, but much larger than those observed for charm mesons. For $`\pi ^{}`$ beams, the leading particle effect predicts no $`\mathrm{\Lambda }_c`$ asymmetry, but a small positive asymmetry has been observed by the E791 Aitala:2000rd  and SELEX Garcia:2001xj  experiments. Explaining the $`\mathrm{\Lambda }_c`$ asymmetries is a severe challenge for any of the proposed mechanisms for generating charm asymmetries models-IC ; models-B .

The factorization theorems of perturbative QCD Collins:gx  imply that the cross section for $`\mathrm{\Lambda }_c^+`$ in a collision between two hadrons $`h,h^{}`$ is given by

$`d\sigma [hh^{}\mathrm{\Lambda }_c^++X]`$ (2)
$`={\displaystyle \underset{i,j}{}}f_{i/h}f_{j/h^{}}d\widehat{\sigma }[ijc\overline{c}+X]D_{c\mathrm{\Lambda }_c^+}+\mathrm{}.`$
Here $`f_{i/h}`$ is a parton distribution, $`d\widehat{\sigma }[ijc\overline{c}+X]`$ is the parton cross section, and $`D_{c\mathrm{\Lambda }_c^+}`$ is the fragmentation function for a $`c`$ quark hadronizing into a $`\mathrm{\Lambda }_c^+`$. The ellipsis represents corrections that are suppressed by powers of $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c`$ or $`\mathrm{\Lambda }_{\mathrm{QCD}}/p_{}`$. The leading order processes $`ggc\overline{c}`$ and $`q\overline{q}c\overline{c}`$ produce $`c`$ and $`\overline{c}`$ symmetrically. Charge conjugation invariance requires that $`D_{c\mathrm{\Lambda }_c^+}=D_{\overline{c}\mathrm{\Lambda }_c^{}}`$, so $`\alpha [\mathrm{\Lambda }_c]=0`$ at leading order in perturbation theory. Next-to-leading order perturbative corrections Nason:1989zy ; Frixione:1994nb  generate asymmetries that are an order of magnitude too small to explain the data. Therefore the observed asymmetries in charm production must come from the power corrections to Eq. (2).

Recent work has shown that $`D`$ meson asymmetries in hadroproduction and photoproduction can be explained by an $`O(\mathrm{\Lambda }_{\mathrm{QCD}}/m_c)`$ power correction called heavy-quark recombination Braaten:2001bf ; Braaten:2001uu ; Braaten:2002yt . In the $`c\overline{q}`$ recombination process, a light antiquark $`\overline{q}`$ that participates in the hard scattering emerges with momentum of $`O(\mathrm{\Lambda }_{\mathrm{QCD}})`$ in the rest frame of the charm quark $`c`$ and the $`c\overline{q}`$ pair then hadronizes into a $`D`$ meson. In this paper, we extend the heavy-quark recombination mechanism to charm baryons. The most important process is $`cq`$ recombination, which is like $`c\overline{q}`$ recombination except the $`\overline{q}`$ is replaced by a light quark $`q`$ and the $`cq`$ diquark hadronizes into a charm baryon. We will show that this mechanism can explain the observed $`\mathrm{\Lambda }_c`$ asymmetries.

The $`c\overline{q}`$ recombination cross section for a $`D`$ meson is

$`d\widehat{\sigma }[\overline{q}gD]={\displaystyle \underset{n}{}}d\widehat{\sigma }[\overline{q}gc\overline{q}(n)+\overline{c}]\rho [c\overline{q}(n)D].`$ (3)
The factor $`\rho [c\overline{q}(n)D]`$ takes into account the nonperturbative hadronization of a $`c\overline{q}`$ with color and spin quantum numbers $`n`$ into a final state that includes the $`D`$ meson. Since the process is inclusive, the quantum numbers of the $`c\overline{q}`$ produced in the short-distance process can be different from that of the $`D`$. The color and spin quantum numbers can both be changed by the emission of soft gluons in the hadronization process. The flavor of the recombining $`\overline{q}`$ can also be different from that of the valence antiquark in the $`D`$, but this requires creating a light quark-antiquark pair which is suppressed in the large $`N_c`$ limit. Neglecting such contributions, the heavy-quark recombination cross section for $`D^+`$ depends on four independent parameters:

$`\rho _1=\rho [c\overline{d}(^1S_0^{(1)})D^+],\stackrel{~}{\rho }_1=\rho [c\overline{d}(^3S_1^{(1)})D^+],`$ (4)
$`\rho _8=\rho [c\overline{d}(^1S_0^{(8)})D^+],\stackrel{~}{\rho }_8=\rho [c\overline{d}(^3S_1^{(8)})D^+].`$
Explicit expressions for these parameters in terms of nonperturbative QCD matrix elements can be found in Ref. Chang:2003ag . They scale with the heavy quark mass as $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c`$. Analogous parameters for $`D^0`$ and $`D^{}`$ mesons are obtained by using isospin symmetry and charge conjugation invariance, while parameters for $`D^+`$ states are related to those in Eq. (4) by heavy-quark spin symmetry. One might have expected the cross section in Eq. (3) to involve a convolution with a nonperturbative function that depends on the fraction of the light-cone momentum of the $`D`$ meson carried by the light antiquark $`\overline{q}`$. However, to lowest order in $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c`$, only the leading moment of such a distribution is relevant. Therefore, the $`c\overline{q}`$ recombination cross sections are calculable using perturbative QCD up to an overall multiplicative factor $`\rho `$.

The direct $`c\overline{q}`$ recombination process is not expected to be a significant source of charm baryons, since baryon production requires creating at least two light quark-antiquark pairs and is therefore suppressed by $`1/N_c^2`$ relative to Eq. (4). The leading recombination mechanism for charm baryon production is $`cq`$ recombination. A leading order Feynman diagram for this process is shown in Fig. 1. Creation of a light quark-antiquark pair is required for the $`cq`$ to hadronize into either a charm meson or a charm baryon, so there is a $`1/N_c`$ suppression in either case. This factor makes $`cq`$ recombination a subleading mechanism for charm mesons, but the leading mechanism for charm baryons. The $`cq`$ recombination cross section for $`\mathrm{\Lambda }_c^+`$ has the form

$`d\widehat{\sigma }[qg\mathrm{\Lambda }_c^+]={\displaystyle \underset{n}{}}d\widehat{\sigma }[qgcq(n)+\overline{c}]\eta [cq(n)\mathrm{\Lambda }_c^+].`$ (5)
The factor $`\eta [cq(n)\mathrm{\Lambda }_c^+]`$ takes into account the nonperturbative hadronization of a $`cq`$ with color and spin quantum numbers $`n`$ into a final state that includes the $`\mathrm{\Lambda }_c^+`$. Isospin symmetry implies that $`\eta [cq(n)\mathrm{\Lambda }_c^+]`$ is the same for $`q=u,d`$, while it is suppressed by $`1/N_c`$ for $`q=s`$. There are four possible color and spin states of the $`cq`$ and therefore four independent $`\eta `$ parameters:

$`\eta _3=\eta [cu(^1S_0^{(\overline{3})})\mathrm{\Lambda }_c^+],\stackrel{~}{\eta }_3=\eta [cu(^3S_1^{(\overline{3})})\mathrm{\Lambda }_c^+],`$ (6)
$`\eta _6=\eta [cu(^1S_0^{(6)})\mathrm{\Lambda }_c^+],\stackrel{~}{\eta }_6=\eta [cu(^3S_1^{(6)})\mathrm{\Lambda }_c^+].`$
These parameters scale as $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c`$, so $`cq`$ recombination gives a power-suppressed contribution to the cross section.

The parton cross sections for $`cq`$ recombination can be calculated using a straightforward generalization of the method described in Ref. Braaten:2001uu  for $`c\overline{q}`$ recombination. Charge conjugation is applied to the line of spinors and Dirac matrices associated with the recombining light quark in Fig. 1. Angular momentum states can then be projected out using the operators of Ref. Braaten:2001uu . The amplitude is projected onto the appropriate color representation, which is either the $`\overline{3}`$ or $`6`$ of $`SU(3)`$. A simple prescription for projecting onto the leading moment of the light-cone momentum fraction of the $`q`$ can be found in Ref. Braaten:2001uu . The parton cross sections for $`cq`$ recombination are

$`{\displaystyle \frac{d\widehat{\sigma }}{d\widehat{t}}}[qgcq(n)+\overline{c}]`$ $`=`$ $`{\displaystyle \frac{2\pi ^2\alpha _s^3}{27}}{\displaystyle \frac{m_c^2}{\widehat{s}^3}}F(n|\widehat{s},\widehat{t}),`$ (7)
where $`\widehat{s}`$, $`\widehat{t}`$, and $`\widehat{u}`$ are the standard parton Mandelstam variables for $`c+qcq(n)+\overline{c}`$. The functions $`F(n|\widehat{s},\widehat{t})`$ for the four color and spin channels are

$`F(^1S_0^{(\overline{3})}|\widehat{s},\widehat{t})`$ $`=`$ $`{\displaystyle \frac{16U}{S}}\left(1{\displaystyle \frac{ST}{U^2}}\right){\displaystyle \frac{m_c^2}{T}}\left(3+{\displaystyle \frac{28U}{T}}+{\displaystyle \frac{16U^2}{T^2}}{\displaystyle \frac{16T^2}{U^2}}\right)+{\displaystyle \frac{4m_c^4S}{UT^2}}\left(3+{\displaystyle \frac{4T}{U}}+{\displaystyle \frac{8U}{T}}\right),`$ (8)
$`F(^3S_1^{(\overline{3})}|\widehat{s},\widehat{t})`$ $`=`$ $`3F(^1S_0^{(\overline{3})}|\widehat{s},\widehat{t})32({\displaystyle \frac{T}{U}}{\displaystyle \frac{U^2}{T^2}}){\displaystyle \frac{4m_c^2}{T}}(8{\displaystyle \frac{6U}{T}}{\displaystyle \frac{16U^2}{T^2}}+{\displaystyle \frac{13T}{U}}+{\displaystyle \frac{15T^2}{U^2}}),`$ (9)
$`F(^1S_0^{(6)}|\widehat{s},\widehat{t})`$ $`=`$ $`{\displaystyle \frac{4U}{S}}\left(2{\displaystyle \frac{5ST}{U^2}}\right){\displaystyle \frac{m_c^2}{T}}\left(27+{\displaystyle \frac{14U}{T}}+{\displaystyle \frac{8U^2}{T^2}}{\displaystyle \frac{20T^2}{U^2}}\right)+{\displaystyle \frac{2m_c^4S}{UT^2}}\left(9+{\displaystyle \frac{10T}{U}}+{\displaystyle \frac{8U}{T}}\right),`$ (10)
$`F(^3S_1^{(6)}|\widehat{s},\widehat{t})`$ $`=`$ $`3F(^1S_0^{(6)}|\widehat{s},\widehat{t}){\displaystyle \frac{8U}{S}}(3+{\displaystyle \frac{5ST}{U^2}}+{\displaystyle \frac{5U}{T}}+{\displaystyle \frac{2U^2}{T^2}})+{\displaystyle \frac{4m_c^2S}{U^2}}(27{\displaystyle \frac{U}{T}}{\displaystyle \frac{U^2}{T^2}}{\displaystyle \frac{8U^3}{T^3}}).`$ (11)
where $`S=\widehat{s}`$, $`T=\widehat{t}m_c^2`$, and $`U=\widehat{u}m_c^2`$.

The parton cross sections for both $`cq`$ and $`c\overline{q}`$ recombination are strongly peaked in the forward direction of the incoming $`q`$ or $`\overline{q}`$. For example, consider the ratio of the parton cross sections for $`cq`$ recombination to that for $`ggc\overline{c}`$, which dominates the fragmentation term in the cross section. We define $`\theta `$ to be the angle between the incoming $`q`$ and outgoing $`cq`$ in the parton center-of-momentum frame. In the backward direction $`\theta =\pi `$, the ratio is suppressed by $`m_c^2/S`$ in both $`{}_{}{}^{3}S_{1}^{}`$ channels and by $`m_c^6/S^3`$ in both $`{}_{}{}^{1}S_{0}^{}`$ channels. At $`\theta =\pi /2`$, the ratio is suppressed by $`m_c^2/S`$ in all four channels. In the forward direction $`\theta =0`$, there is no kinematic suppression of this ratio. The forward enhancement of the $`c\overline{q}`$ cross section gives charm meson asymmetries which are largest near $`x_F1`$. For $`\mathrm{\Lambda }_c`$ produced in $`pN`$ collisions, the fragmentation cross section is smaller relative to $`cq`$ recombination, so the asymmetry is large even for $`x_F=0.2`$.

In addition to direct recombination of $`cq`$ into $`\mathrm{\Lambda }_c^+`$, we need to include two additional effects: $`cq`$ recombination into a heavier charm baryon that subsequently decays into $`\mathrm{\Lambda }_c^+`$ and the fragmentation into $`\mathrm{\Lambda }_c^+`$ of a $`c`$ that is produced in a $`\overline{c}\overline{q}`$ or $`\overline{c}q`$ recombination process. The cross sections for $`\mathrm{\Lambda }_c^+`$ from the latter process, which we will call “opposite-side recombination”, are

$`d\widehat{\sigma }[qg\mathrm{\Lambda }_c^++X]`$ $`=`$ $`{\displaystyle \underset{n}{}}d\widehat{\sigma }[qg\overline{c}q(n)+c]`$ (12)
$`\times {\displaystyle \underset{\overline{D}}{}}\rho [\overline{c}q(n)\overline{D}]D_{c\mathrm{\Lambda }_c^+},`$
$`d\widehat{\sigma }[\overline{q}g\mathrm{\Lambda }_c^++X]`$ $`=`$ $`{\displaystyle \underset{n}{}}d\widehat{\sigma }[\overline{q}g\overline{c}\overline{q}(n)+c]`$ (13)
$`\times {\displaystyle \underset{\overline{B}}{}}\eta [\overline{c}\overline{q}(n)\overline{B}]D_{c\mathrm{\Lambda }_c^+}.`$
The recombination factors in Eq. (12) and Eq. (13) are summed over $`\overline{D}`$ mesons whose valence partons are $`\overline{c}q`$ and over antibaryons $`\overline{B}`$ whose valence partons inlcude $`\overline{c}\overline{q}`$.

The feeddown from heavier charm baryons that decay into $`\mathrm{\Lambda }_c`$ can be taken into account through inclusive $`\eta `$ parameters defined by

$`\eta _{\mathrm{inc}}[cq(n)\mathrm{\Lambda }_c^+]`$ $`=`$ $`\eta [cq(n)\mathrm{\Lambda }_c^+]`$
$`+`$ $`{\displaystyle \underset{B}{}}\eta [cq(n)B]B[B\mathrm{\Lambda }_c^++X].`$
The sum over $`B`$ includes all charm baryons that decay into $`\mathrm{\Lambda }_c^+`$. They include $`\mathrm{\Sigma }_c^+`$, $`\mathrm{\Sigma }_c^+`$, and the negative-parity excitations $`\mathrm{\Lambda }_c^+(2593)`$ and $`\mathrm{\Lambda }_c^+(2625)`$ states, all of which have branching fractions into $`\mathrm{\Lambda }_c^+`$ of nearly 100%. Charm baryons with strangeness do not contribute to $`\eta _{\mathrm{inc}}`$.

In our analysis, we choose $`m_c=1.5`$ GeV, use the one-loop running $`\alpha _s`$ with $`\mathrm{\Lambda }_{\mathrm{QCD}}=200`$ MeV, and set the renormalization and factorization scales equal to $`\sqrt{p_{}^2+m_c^2}`$. We use the parton distributions GRV 98 LO Gluck:1998xa  for the proton and GRV-P LO Gluck:1991ey  for the pion. For the fragmentation function for $`c\mathrm{\Lambda }_c^+`$, we use

$`D_{c\mathrm{\Lambda }_c^+}(z)=f_{\mathrm{\Lambda }_c^+}\delta (1z)`$ (15)
where $`f_{\mathrm{\Lambda }_c^+}=0.076`$ is the inclusive fragmentation probability Gladilin:1999pj . We also use delta-function fragmentation functions for the charm mesons, since this reproduces the shapes of their momentum distributions more accurately than Peterson fragmentation functions Frixione:1994nb . In the opposite-side $`\overline{c}q`$ recombination cross section, Eq. (12), we include the $`D`$ and $`D^{}`$ multiplets, but neglect the excited charm mesons. The best 1-parameter fit to all the $`D`$ meson asymmetries measured by E791 gives $`\rho _1=0.15`$ with $`\stackrel{~}{\rho }_1=\rho _8=\stackrel{~}{\rho }_8=0`$. The fit can be improved by using additional parameters, but not dramatically. This value of $`\rho _1`$ is larger than the value $`\rho _1=0.06`$ obtained in Ref. Braaten:2002yt  using Peterson fragmentation functions. Note that the sum of recombination parameters appearing in the opposite-side $`\overline{c}\overline{q}`$ recombination cross section, Eq. (13), differs from the inclusive parameter in Eq. ($`\mathrm{\Lambda }_c^+/\mathrm{\Lambda }_c^{}`$ Asymmetry in Hadroproduction from Heavy-Quark Recombination). The two are related if the the sum in Eq. (13) is dominated by the lowest mass $`J^P=\frac{1}{2}^+`$ and $`\frac{3}{2}^+`$ SU(3) multiplets. Then, using charge conjugation and a simple quark counting argument, we estimate $`_{\overline{B}}\eta [\overline{c}\overline{q}(n)\overline{B}]3/2\eta _{\mathrm{inc}}[cu(n)\mathrm{\Lambda }_c^+]`$ for $`q=u,d,s`$.

Because of the large uncertainties associated with the parton distributions of the $`\mathrm{\Sigma }^{}`$, we focus on the $`\mathrm{\Lambda }_c`$ asymmetry data from $`\pi ^{}`$ and $`p`$ beams measured by the E791 and SELEX experiments. The best 1-parameter fit to all this data yields $`\stackrel{~}{\eta }_{\mathrm{\hspace{0.17em}3},\mathrm{inc}}=0.058`$, with the other three inclusive $`\eta `$ parameters set to 0. The asymmetry variable $`\alpha [\mathrm{\Lambda }_c]`$ is shown as a function of $`x_F`$ for the pion beam in Fig. 2 and for the proton beam in Fig. 3. The 1-parameter fit agrees well with both the pion beam and proton beam data. The same fits also yield good agreement with observed $`p_{}`$ dependence of the asymmetries. The fits can be improved by using additional parameters. In Figs. 2 and 3, we also show the predictions for $`\alpha [\mathrm{\Lambda }_c]`$ if all the $`\eta `$ parameters are set to 0. Note that opposite-side recombination into $`D`$ mesons generates a positive asymmetry even if all the $`\eta `$ parameters vanish. It gives a reasonable fit to the pion beam data, but it severely underpredicts the asymmetry for the proton beam. Therefore, the large asymmetry from the proton beam is convincing evidence for the $`cq`$ recombination mechanism.

We have shown that heavy quark recombination provides a natural and economical explanation of the production asymmetries for charm baryons as well as charm mesons. Further work includes a more systematic analysis of all the charm asymmetry data from hadroproduction experiments and the prediction of charm and bottom asymmetries in present and future experiments. Previous analyses of $`D`$ meson asymmetries Braaten:2001uu ; Braaten:2002yt  do not include contributions from opposite-side $`cq`$ recombination into charm baryons. This is particularly important for $`D_s`$ mesons since any asymmetry is generated by opposite side recombination.

E.B. and M.K. are supported in part by DOE grant DE-FG02-91-ER4069. Y.J. is supported by NSF grant PHY-0100677. T.M. is supported in part by DOE grants DE-FG02-96ER40945 and DE-AC05-84ER40150.


format: text

text:
# Mixtures of Bosonic and Fermionic Atoms in Optical Lattices

## I Introduction

Recent spectacular progress in the manipulation of neutral atoms in optical lattices Greiner1 ; Greiner2 ; OrzTuc01  has opened the way to the simulation of complex quantum systems of condensed matter physics, such as high–$`T_c`$ superconductors, Hall systems, and superfluid 4He, by means of atomic systems with perfectly controllable physical parameters Anglin . Optical lattices are stable periodic arrays of microscopic potentials created by the interference patterns of intersecting laser beams Jessen1 . Atoms can be confined to different lattice sites, and by varying the strength of the periodic potential it is possible to tune the interatomic interactions with great precision. They can be enhanced well into regimes of strong correlation, even in the dilute limit. The transition to a strong coupling regime can be realized by increasing the depth of the lattice potential wells, a quantity that is directly proportional to the intensity of the laser light. This is an experimental parameter that can be controlled with great accuracy. For this reason, besides the fundamental interest for the investigation of quantum phase transitions Sachdev  and other basic quantum phenomena JacBru98 ; OosStr02 ; RuoDun02 ; Hofstetter02 ; Pared02 ; Recati03 ; Zwerger03 , optical lattices have become an important practical tool for applications, ranging from laser cooling Kerman  to quantum control and information processing Jessen2 ; DuaDemLuk03 , and quantum computation Deutsch ; JakBriCirGarZol98 ; GarCir03 ; DorFed02 ; PacKni03 .

The theory of neutral Bosonic atoms in optical lattices has been developed JacBru98  by assuming that the atoms are confined to the lowest Bloch band of the periodic potential. It can then be shown that the system is effectively described by a single–band Bose–Hubbard model Hamiltonian Fisher . In such a model the superfluid–insulator transition is predicted to occur when the on site Boson–Boson interaction energy becomes comparable to the hopping energy between adjacent lattice sites. This situation can be experimentally achieved by increasing the strength of the lattice potential, which results in a strong suppression of the kinetic (hopping) energy term. In this way, the superfluid–Mott-insulator quantum phase transition has been realized by loading an ultracold atomic Bose–Einstein condensate in an optical lattice Greiner1 .

The present paper is concerned with the study of dilute mixtures of interacting Bosonic and Fermionic neutral atoms subject to an optical lattice and a superimposed trapping harmonic potential at zero temperature. We assume the fermions to be identical (for instance spin–polarized in a magnetic trap), so that there are only $`s`$–wave Boson–Boson and Fermion–Boson contact interactions present. We construct an effective single–band Bose–Fermi Hubbard Hamiltonian, and we determine the ground state energy and on site density distributions in different mean field approximations of increasing complexity. Our main aim at this level of description is to determine the basic ground state properties of the mixture and to study how the Bosonic superfluid–insulator transition is influenced by the presence of the Fermions. Besides the study of the latter issue, and the assessment of the properties of linear stability of the system against Boson–Fermion demixing, a remarkable finding of our analysis is that a quantum binary mixture loaded in a very deep optical lattice allows for a disordered phase of very many degenerate or quasi–degenerate ground states separated by very high potential energy barriers. In the limit of very large lattice potential strengths the basic mirror symmetry of the optical lattice is broken.

The plan of the paper is as follows: In Section II we set the notation and derive the Bose–Fermi Hubbard model Hamiltonian. We then discuss the range of validity of the approximations and the assumptions used in the derivation. In Section III we introduce some basic mean field descriptions to study the properties of stability of the mixture against phase separation, and we provide a simple analytical criterion for the onset of a superfluid phase for the Bosons starting from a Mott–insulating ground state.

In Section IV we present numerical simulations for a small number of particles in the framework of the Gutzwiller variational ansatz. Two important cases have to be distinguished: Boson–Fermion repulsive or attractive interaction (the Boson–Boson interaction is taken to be always repulsive). In the first instance, one can observe a continuous transition to a complete demixing of the Fermions from the Bosons, and to a Mott–insulating phase of the Bosons, as the strength of the lattice potential is increased. In the case of attractive Boson–Fermion interactions there is no transition to demixing or to collapse of the mixture, while one still observes a Mott–insulating transition for the Bosons. By studying the behavior of the superfluid order parameter we show that in both cases the transition takes place at the same critical value of the lattice potential strength. Moreover, due to the strong attraction, the Fermions and the Bosons tend to form together ordered block–crystalline structures at the center of the trap.

In Section V we present a numerical analysis that, although constrained to a small number of particles (five bosons and five fermions), seems to indicate the existence of a rich structure of degenerate energy minima for large enough values of the lattice potential strength. Such an energy landscape suggests the possible existence of disordered phases of the mixture due to the delicate interplay between the different physical parameters (Boson hopping, Fermion hopping, Boson on site energy, Boson–Fermion on site energy), the lattice depth, and the symmetries of the problem. Finally, a summary and an outlook to future research are shortly discussed in Section VI.

## II Model Hamiltonian

We start by introducing the Hamiltonian for a Bose–Fermi mixture loaded into optical lattice potentials and confined by additional, slowly varying, external (harmonic) trapping potentials. It is given by

$$\widehat{H}=\widehat{T}_B+\widehat{T}_F+\widehat{V}_B+\widehat{V}_F+\widehat{W}_{BB}+\widehat{W}_{BF},$$
(1)
where

$`\widehat{T}_B`$ $`=`$ $`{\displaystyle d^3𝐫\widehat{\mathrm{\Phi }}^{}(𝐫)\frac{\mathrm{}^2^2}{2m_B}\widehat{\mathrm{\Phi }}(𝐫)},`$ (2)
$`\widehat{T}_F`$ $`=`$ $`{\displaystyle d^3𝐫\widehat{\mathrm{\Psi }}^{}(𝐫)\frac{\mathrm{}^2^2}{2m_F}\widehat{\mathrm{\Psi }}(𝐫)},`$ (3)
represent the Boson and Fermion kinetic energies, respectively, while

$`\widehat{W}_{BB}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{4\pi ^2\mathrm{}^2a_{BB}}{m_B}}{\displaystyle d^3𝐫\widehat{\mathrm{\Phi }}^{}(𝐫)\widehat{\mathrm{\Phi }}^{}(𝐫)\widehat{\mathrm{\Phi }}(𝐫)\widehat{\mathrm{\Phi }}(𝐫)},`$ (4)
$`\widehat{W}_{BF}`$ $`=`$ $`{\displaystyle \frac{2\pi ^2\mathrm{}^2a_{BF}}{m_R}}{\displaystyle d^3𝐫\widehat{\mathrm{\Phi }}^{}(𝐫)\widehat{\mathrm{\Psi }}^{}(𝐫)\widehat{\mathrm{\Psi }}(𝐫)\widehat{\mathrm{\Phi }}(𝐫)},`$ (5)
denote the Boson–Boson and the Fermion–Boson contact interaction energies. They are parametrized by the Boson–Boson and the Fermion–Boson $`s`$–wave scattering lengths $`a_{BB}`$ and $`a_{BF}`$, respectively, and by the Boson mass $`m_B`$ and the reduced mass $`m_R=m_Bm_F/(m_B+m_F)`$, where $`m_F`$ denotes the Fermion mass. The potential energies

$`\widehat{V}_B`$ $`=`$ $`{\displaystyle d^3𝐫\widehat{\mathrm{\Phi }}^{}(𝐫)\left(V_B(𝐫)+P_B(𝐫)\right)\widehat{\mathrm{\Phi }}(𝐫)},`$ (6)
$`\widehat{V}_F`$ $`=`$ $`{\displaystyle d^3𝐫\widehat{\mathrm{\Psi }}^{}(𝐫)\left(V_F(𝐫)+P_F(𝐫)\right)\widehat{\mathrm{\Psi }}(𝐫)},`$ (7)
are due to the trapping and lattice potentials. We consider pure magnetic trapping for Bosons and Fermions, so that Fermions are spin–polarized and their $`s`$–wave interaction energy $`\widehat{W}_{FF}`$ can be neglected: $`\widehat{W}_{FF}=0`$. In the subsequent analysis we will consider the harmonic approximation of a typical quadrupolar magnetic field with strong anisotropy in the transverse directions $`y`$ and $`z`$, i.e.,

$`V_B(𝐫)m_B\omega _B^2(x^2+\lambda ^2y^2+\lambda ^2z^2)/2,`$ (8)
and

$$V_F(𝐫)m_F\omega _F^2(x^2+\lambda ^2y^2+\lambda ^2z^2)/2,$$
(9)
where $`\lambda 1`$ is the anisotropy parameter. Moreover, if we assume trapping in the same magnetic state for the Bosons and the Fermions, then the trapping frequencies are related according to $`\omega _F/\omega _B=(m_B/m_F)^{1/2}`$, so that the two potentials coincide: $`V_B(𝐫)=V_F(𝐫)`$. The ground–state harmonic oscillator lengths, however, are different due to the different masses, and also differ for the $`x`$-direction on the one hand and the $`y`$ and $`z`$-directions on the other hand:

$$\mathrm{}_{B/F}^{}=\sqrt{\mathrm{}/(m_{B/F}\omega _{B/F})}$$
(10)
in the $`x`$ direction, and $`\mathrm{}_{B/F}^{}=\mathrm{}_{B/F}^{}/\sqrt{\lambda }`$ in the $`y`$ and $`z`$ directions. We next consider a lattice structure for the Bosons and the Fermions in the $`x`$–direction, associated to the corresponding Bosonic and Fermionic one–dimensional optical lattice potentials $`P_B(x)`$ and $`P_F(x)`$ are

$`P_B(x)`$ $`=`$ $`V_B^0\mathrm{sin}^2(\pi x/a),`$
$`P_F(x)`$ $`=`$ $`V_F^0\mathrm{sin}^2(\pi x/a),`$ (11)
where $`a`$ is the lattice spacing associated to the wave vector $`k=\pi /a`$ of the standing laser light. If the lattice potentials are produced by a far off–resonant laser for both species, the lattice potential strengths are equal for both Fermions and Bosons: $`V_F^0=V_B^0=V_0`$, and the two optical lattices coincide exactly. This is the situation we will always consider in the following.

In the presence of a strong optical lattice and a sufficiently shallow external confinement in the $`x`$ direction, the field operators can be expanded in terms of the single–particle Wannier functions localized at each lattice site $`x_i`$. Further, the typical interaction energies involved are normally not strong enough in order to excite higher vibrational states, and we can retain only the lowest vibrational state in each lattice potential well both for Bosons and Fermions (single–band approximation). In case of stronger external confinements, or interactions, one should include higher Bloch bands as well in the expansion of the field operators, a case we do not consider in the present context. In the harmonic approximation, the Wannier functions $`w(𝐫)`$ factorize in the product of harmonic oscillator states in each direction, with the trapping potential almost constant between adjacent lattice sites. We then have

$`\widehat{\mathrm{\Phi }}(𝐫)`$ $`=`$ $`{\displaystyle \underset{i}{}}\widehat{a}_iw_x^B(xx_i)w_y^B(y)w_z^B(z),`$ (12)
$`\widehat{\mathrm{\Psi }}(𝐫)`$ $`=`$ $`{\displaystyle \underset{i}{}}\widehat{b}_iw_x^F(xx_i)w_y^F(y)w_z^F(z),`$ (13)
where $`\widehat{a}_i`$ and $`\widehat{b}_i`$ are respectively the Bosonic and Fermionic annihilation operators at the $`i`$–th lattice site, $`x_i=ia`$, and the index $`i`$ runs on positive and negative integers, the origin of the lattice being fixed at $`i=0`$ so that it coincides with the center (the minimum) of the external trapping potential. In each lattice potential well the Wannier local ground states for Bosons and Fermions are Gaussians in the harmonic approximation:

$`w_y^{B/F}(y)`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}\left[y^2/2(\mathrm{}_{B/F}^{})^2\right]}{\pi ^{1/4}(\mathrm{}_{B/F}^{})^{1/2}}},`$ (14)
$`w_z^{B/F}(z)`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}\left[z^2/2(\mathrm{}_{B/F}^{})^2\right]}{\pi ^{1/4}(\mathrm{}_{B/F}^{})^{1/2}}},`$ (15)
and

$$w_x^{B/F}(xx_i)=\frac{\mathrm{exp}\left[(xx_i)^2/2(\mathrm{}_{B/F}^0)^2\right]}{\pi ^{1/4}(\mathrm{}_{B/F}^0)^{1/2}},$$
(16)
where

$$\mathrm{}_{B/F}^0=a/[\pi (V_0/E_{B/F}^R)^{1/4}],$$
(17)
is the width of the harmonic oscillator potential wells at each lattice site, with $`E_B^R=(\pi \mathrm{})^2/2a^2m_B`$ and $`E_F^R=(\pi \mathrm{})^2/2a^2m_F`$ being the Boson and Fermion recoil energies, respectively.

In this paper we will consider the physical situation of very shallow trapping potentials, such that $`\mathrm{}_{B/F}^{}aN_{B/F}`$ and consequently local density approximation (LDA) can be applied in the study of the ground–state properties of the system. Therefore, when exploiting the Wannier function expansions (12) and (13) to map the full Hamiltonian (1) into its lattice version, we discard all terms that are of order $`(aN_{B/F}/\mathrm{}_{B/F}^{})^2`$ or of higher powers of it. Otherwise, nonlocal effects caused by the trapping potential, like site–dependent hopping terms, have to be considered Inpreparation . Finally, this approximation scheme leads to the following Hubbard–type Hamiltonian:

$`\widehat{H}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left(J_B\widehat{a}_{i+1}^{}\widehat{a}_i+J_F\widehat{b}_{i+1}^{}\widehat{b}_i+\text{H. c.}\right)`$
$`+`$ $`{\displaystyle \frac{U_{BB}}{2}}{\displaystyle \underset{i}{}}\widehat{n}_B^{(i)}(\widehat{n}_B^{(i)}1)+U_{BF}{\displaystyle \underset{i}{}}\widehat{n}_B^{(i)}\widehat{n}_F^{(i)}`$
$`+`$ $`{\displaystyle \underset{i}{}}V_B^{(i)}\widehat{n}_B^{(i)}+{\displaystyle \underset{i}{}}V_F^{(i)}\widehat{n}_F^{(i)}`$
$`+`$ $`\mathrm{}\left(\lambda \omega _B+\omega _B^0/2\right)\widehat{N}_B+\mathrm{}\left(\lambda \omega _F+\omega _F^0/2\right)\widehat{N}_F.`$
The first line in the above Bose–Fermi Hubbard Hamiltonian describes independent nearest–neighbor hopping of Bosons and Fermions, with amplitudes $`J_B`$ and $`J_F`$, respectively. The terms in the second line describe Boson–Boson on site repulsion (with $`U_{BB}>0`$) and Boson–Fermion on site interaction. This interaction can be repulsive or attractive, depending on the sign of $`U_{BF}`$. The third line describes the energy offset at each lattice site due to the $`x`$ component of the external trapping potentials $`V_{B/F}(𝐫)`$, and the last line contains the overall constant zero–point energy terms due to the $`y`$ and $`z`$ components of $`V_{B/F}(𝐫)`$ and to the lattice potential $`P(x)`$. The on site interaction and offset energy terms are simple functions of the on site Boson and Fermion occupation number operators $`\widehat{n}_B^{(i)}=\widehat{a}_i^{}\widehat{a}_i`$ and $`\widehat{n}_F^{(i)}=\widehat{b}_i^{}\widehat{b}_i`$, while the zero–point energy terms are proportional to the total particle number operators $`\widehat{N}_B=_i\widehat{a}_i^{}\widehat{a}_i`$ and $`\widehat{N}_F=_i\widehat{b}_i^{}\widehat{b}_i`$. The frequency

$$\omega _{B/F}^0=\mathrm{}/[(\mathrm{}_{B/F}^0)^2m_{B/F}]$$
(19)
fixes the Bosonic and Fermionic harmonic oscillations in each lattice well. The relevant parameters entering in the Hamiltonian are the on site values of the trapping harmonic potential

$`V_{B/F}^{(i)}`$ $`=`$ $`{\displaystyle \frac{m_{B/F}}{2}}\omega _{B/F}^2x_i^2,`$ (20)
the nearest–neighbor hopping amplitudes between adjacent sites $`x_i`$ and $`x_{i+1}`$ for Bosons and Fermions

$`J_{B/F}`$ $`=`$ $`{\displaystyle }dxw_x^{B/F}(xx_i)[{\displaystyle \frac{\mathrm{}^2}{2m_{B/F}}}{\displaystyle \frac{d^2}{dx^2}}`$ (21)
$`+`$ $`V_0\mathrm{sin}^2\left(\pi {\displaystyle \frac{x}{a}}\right)]w^{B/F}_x(xx_{i+1}),`$
the strength of the on site repulsion energy between two Bosonic atoms at the same lattice site

$`U_{BB}`$ $`=`$ $`{\displaystyle \frac{4\pi \mathrm{}^2a_{BB}}{m_B}}{\displaystyle 𝑑x(w_x^B(xx_i))^4}`$ (22)
$`\times `$ $`{\displaystyle 𝑑y(w_y^B(y))^4𝑑z(w_z^B(z))^4},`$
and the strength of the on site interaction energy (either repulsive or attractive) between a Bosonic and a Fermionic atom at the same lattice site

$`U_{BF}`$ $`=`$ $`{\displaystyle \frac{2\pi \mathrm{}^2a_{BF}}{m_R}}{\displaystyle 𝑑x\left[w_x^B(xx_i)w_x^F(xx_i)\right]^2}`$
$`\times `$ $`{\displaystyle 𝑑y\left[w_y^B(y)w_y^F(y)\right]^2𝑑z\left[w_z^B(z)w_z^F(z)\right]^2}.`$
In typical situations we may neglect next–to–nearest neighbor hopping amplitudes and nearest–neighbor interaction couplings that are usually some orders of magnitude smaller, so that the Hamiltonian (II) provides a rather accurate model for the dynamics of a Bose–Fermi mixture with three–dimensional scattering in a one–dimensional periodic potential. Terms involving nearest–neighbor interaction strengths and/or next–to–nearest neighbor hopping amplitudes can become relevant and need to be included, e.g., when considering phonon exchange between Fermions, and this would lead to a Bose–Fermi analog of the so–called extended Hubbard models. To evaluate estimates for the parameters entering the Bose–Fermi Hubbard Hamiltonian (II) using Eqns. (14),(15) and (16), we will set the Boson recoil energy $`E_B^R=\mathrm{}^2\pi ^2/(2m_Ba^2)`$ as the unit of energy. We then introduce the dimensionless quantity $`\stackrel{~}{V}_0=V_0/E_B^R`$, and, analogously, the dimensionless quantities $`\stackrel{~}{U}_{BB}`$, $`\stackrel{~}{U}_{BF}`$, $`\stackrel{~}{V}_B^{(i)}`$, $`\stackrel{~}{V}_F^{(i)}`$, $`\stackrel{~}{J}_B`$, and $`\stackrel{~}{J}_F`$. We then have

$`\stackrel{~}{U}_{BB}`$ $`=`$ $`\sqrt{{\displaystyle \frac{8}{\pi ^3}}}{\displaystyle \frac{a_{BB}a}{(\mathrm{}_B^{})^2}}\stackrel{~}{V}_0^{1/4},`$ (24)
$`\stackrel{~}{U}_{BF}`$ $`=`$ $`\sqrt{{\displaystyle \frac{8}{\pi ^3}}}\left(1+{\displaystyle \frac{m_B}{m_F}}\right){\displaystyle \frac{a_{BF}a}{(\mathrm{}_B^{})^2+(\mathrm{}_F^{})^2}}\stackrel{~}{V}_0^{1/4},`$ (25)
$`\stackrel{~}{V}_B^{(i)}={\displaystyle \frac{i^2}{\pi ^2(\mathrm{}_B^{}/a)^4}},\stackrel{~}{V}_F^{(i)}={\displaystyle \frac{m_B}{m_F}}{\displaystyle \frac{i^2}{\pi ^2(\mathrm{}_F^{}/a)^4}},`$ (26)
$`\stackrel{~}{J}_B`$ $`=`$ $`\left({\displaystyle \frac{\pi ^2}{4}}1\right)\stackrel{~}{V}_0\mathrm{exp}\left[{\displaystyle \frac{\pi ^2}{4}}\sqrt{\stackrel{~}{V}_0}\right],`$ (27)
$`\stackrel{~}{J}_F`$ $`=`$ $`\left({\displaystyle \frac{\pi ^2}{4}}1\right)\stackrel{~}{V}_0\mathrm{exp}\left[{\displaystyle \frac{\pi ^2}{4}}\sqrt{{\displaystyle \frac{m_F}{m_B}}\stackrel{~}{V}_0}\right].`$ (28)
In FIG. 1 we show the dependencies of the these parameters on the potential strength $`\stackrel{~}{V}_0`$ (compare also Ref. DuaDemLuk03 ). For reference we have included as well the overlap integral
$`w(xx_i)|w(xx_{i+1})`$ of adjacent Wannier functions. The overlap is negligible but for very small values of the potential strength, confirming that terms of the order of the overlap integral can be neglected in the Hamiltonian. The Gaussian approximation holds rather well as can be seen by comparing the associated Bosonic hopping amplitude $`J_B`$ with the one obtained by using the exact 1–$`D`$ Mathieu equation Zwerger03 .

Besides the conditions mentioned earlier, all the expressions derived in the present section are justified under the following circumstances. First of all, we must require that the two–body scattering processes are not influenced by the confinements, a condition that is guaranteed if the lengths of the confining and lattice potentials in all directions are much larger than the Boson–Boson and Fermion–Boson scattering lengths. Next, the single–band structure of the lattice Hamiltonian is assured if the lattice spacing $`a`$ is much greater than the harmonic confinements in each direction at all lattice sites. On the other hand, in this limit the harmonic approximation for the Wannier functions at each lattice well is automatically satisfied. Finally, as mentioned earlier, the assumption of a slowly varying confining potential such that LDA is applicable leads to the condition $`\mathrm{}_{B/F}^{}aN_{B/F}`$. We can summarize all the above conditions with the following chain of inequalities:

$$\{|a_{BF}|,a_{BB}\}\{\mathrm{}_{B/F}^0,\mathrm{}_{B/F}^{}\}a\mathrm{}_{B/F}^{}/N_{B/F}.$$
(29)
Our model is for some aspects unrealistic, since in present experimental situations the transverse confinements cannot be made very strong. Therefore a multi–band structure can appear with several radial states being occupied, as reported in a recent experiment by the Florence group on Bose–Fermi mixtures in a $`1`$$`D`$ optical lattice GiovanniModugno .

## III Phase stability and the superfluid transition

In this section we investigate the zero temperature ground state properties of the system in a mean field approximation. In the following we will adopt a grand-canonical description through the Hamiltonian

$$\widehat{K}=\widehat{H}\mu _B\widehat{N}_B\mu _F\widehat{N}_F,$$
(30)
where $`\mu _B`$ and $`\mu _F`$ are the Bosonic and Fermionic chemical potentials. According to the Hohenberg–Kohn theorem, the ground state energy

$$E=\mathrm{\Psi }_0|\widehat{K}|\mathrm{\Psi }_0$$
(31)
is a functional of the on site Bosonic and Fermionic densities $`n_B^{(i)}=\widehat{a}_i^{}\widehat{a}_i`$ and $`n_F^{(i)}=\widehat{b}_i^{}\widehat{b}_i`$, where the expectation values are taken with respect to the ground state with state vector $`|\mathrm{\Psi }_0`$. We decompose the functional $`E`$ according to

$`E=E_B+E_F+E_{BF}\mu _B{\displaystyle \underset{i}{}}n_B^{(i)}\mu _F{\displaystyle \underset{i}{}}n_F^{(i)},`$ (32)
where $`E_B`$ is the energy contribution depending only on the Boson parameters $`J_B`$, $`U_{BB}`$, $`V_B^{(i)}`$; $`E_F`$ is the energy depending only on the Fermion parameters; and $`E_{BF}`$ is the term due to Boson–Fermion interactions. We treat this latter term in mean field approximation: neglecting exchange correlation effects:

$$E_{BF}=U_{BF}\underset{i}{}n_B^{(i)}n_F^{(i)}.$$
(33)
Exchange correlation effects have been recently studied for the case of homogeneous mixtures in the continuum Uns ; Giorgini . For the Fermion energy $`E_F`$, we take the energy of the noninteracting homogeneous system and exploit local density approximation (LDA) on it,

$`E_F={\displaystyle \frac{2J_F}{\pi }}{\displaystyle \underset{i}{}}\mathrm{sin}(\pi n_F^{(i)})+{\displaystyle \underset{i}{}}V_F^{(i)}n_F^{(i)}.`$ (34)
This approximate description of the Fermions is well justified in the presence of a slowly varying trapping potential (so that LDA can be applied), when there are no direct interactions among the Fermions (as in our case), and moreover when one can neglect induced phonon–mediated self–interactions due to the presence of the Bosons. Therefore, in this situation, the nontrivial features of different quantum phases will regard only the Bosonic sector and not the Fermionic one. However, the presence of the Fermions will indirectly contribute to the properties of the different Bosonic phases, and this is the subject that we will study in the following.

In order to find an expression for the Boson energy $`E_B`$ we will proceed in steps of increasing accuracy. First we perform a very simple mean field analysis in two extreme limits: a completely superfluid Boson ground state and a totally Mott–insulating Boson ground state. In the latter case we will provide a simple criterion for stability of the mixture against demixing. Next, we will perform a perturbation expansion around the Mott–insulating Boson ground state to recover perturbatively the phase boundary against transition to superfluidity. Finally, in the next section, we will study the ground state properties of the mixture using a Gutzwiller ansatz for the Bosons capable of describing the intermediate regimes between the insulating and superfluid Bosonic phases.

We first consider the Bosons to be superfluid. In this regime the chemical potential and the number of particles in a homogeneous system are related, to lowest order in $`U_{BB}`$, via OosStr02 :

$`\mu _B=U_{BB}n_02J_B,`$ (35)
where $`n_0`$ is the density of condensed Bosons. Additionally, for very weak interaction $`n_0n_B`$. Exploiting this result in LDA and using the mean field expression for the Bose–Fermi interaction energy we can then write for the inhomogeneous Bose–Fermi mixture at a given lattice site:

$`U_{BB}n_B^{(i)}=\mu _B+2J_BV_B^{(i)}U_{BF}n_F^{(i)}.`$ (36)
Next, we consider the case of a Mott–insulating Bosonic phase. To lowest order in $`J_B`$ we neglect the kinetic term altogether. Then it is easily shown that the relation between the Bosonic chemical potential and the Bosonic density for a homogeneous system is given by

$`\mu _B=U_{BB}n_BU_{BB}/2.`$ (37)
Exploiting LDA as before, we have in the inhomogeneous case at a given lattice site:

$`U_{BB}n_B^{(i)}=\mu _B+U_{BB}/2V_B^{(i)}U_{BF}n_F^{(i)}.`$ (38)
Comparing Eqns. (36) and (38), we observe the same behavior of the on site density profiles but for a constant correction to the Boson chemical potential depending whether the Bosons are in a superfluid or in a Mott–insulating state. Finally, differentiating the energy functional with respect to the on site populations of the Fermions, we determine the associated density field and the set of coupled equations describing the ground state of the mixture at any lattice site,

$`U_{BB}n_B^{(i)}`$ $`=`$ $`\mu _B^{}V_B^{(i)}U_{BF}n_F^{(i)},`$ (39)
$`2J_F\mathrm{cos}(\pi n_F^{(i)})`$ $`=`$ $`\mu _FV_F^{(i)}U_{BF}n_B^{(i)},`$ (40)
where $`\mu _B^{}`$ is the proper expression of the Boson chemical potential according to whether the Bosons are in the Mott–insulating or superfluid regime. These equations are valid at a given lattice site $`i`$ if $`\mu _B^{}V_B^{(i)}U_{BF}n_F^{(i)}>0`$, otherwise one must set $`n_B^{(i)}=0`$. On the other hand, if

$$(\mu _FV_F^{(i)}U_{BF}n_B^{(i)})/2J_F<0$$
(41)
we must impose $`n_F^{(i)}=0`$ at the given lattice site, while $`n_F^{(i)}=1`$ must be imposed when $`(\mu _FV_F^{(i)}U_{BF}n_B^{(i)})/2J_F>1`$. These expressions are the lattice analogs of the Thomas–Fermi description of Boson–Fermion mixtures in the continuum. We remark that in the Mott–insulating regime the Boson on site populations $`n_B^{(i)}`$ must be rounded off to the integer closest to the solutions of Eqns. (39)–(40).

In the Mott–insulating regime we can determine a criterion of linear stability against phase separation of the two species if we expand the energy functional $`E`$ to second order in the small density variations $`\delta n_{B/F}^{(i)}`$ around the minimum provided by the solution of Eqns. (39)–(40):

$`\delta ^2E`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left(\begin{array}{c}\delta n_B^{(i)}\\ \delta n_F^{(i)}\end{array}\right)[\left(\begin{array}{cc}U_{BB}& U_{BF}\\ U_{BF}& 2\pi J_F\mathrm{sin}(\pi n_F^{(i)})\end{array}\right)`$ (49)
$`\times \left(\begin{array}{c}\delta n_B^{(i)}\\ \delta n_F^{(i)}\end{array}\right)].`$
This quadratic form is positive at a given site $`i`$ if and only if

$`2\pi J_F\mathrm{sin}(\pi n_F^{(i)})U_{BB}>U_{BF}^2`$ (50)
and $`2\pi J_F\mathrm{sin}(\pi n_F^{(i)})+U_{BB}0`$. This last condition is always satisfied for $`U_{BB}>0`$ and identical Fermions. If this is not the case for every site $`i`$, then the ground state is not stable against demixing. This result is similar to that recently obtained for a mixture of two different Boson species on a lattice CheWu03 , which states that the mixture is stable if $`U_1U_2>U_{12}^2`$, where $`U_1`$ and $`U_2`$ are the Boson–Boson interaction strengths of species $`1`$ and $`2`$ respectively, and $`U_{12}`$ is the interspecies coupling. The form of expression (50) then suggests that the Pauli on site energy $`2\pi J_F\mathrm{sin}(\pi n_F^{(i)})`$ has the meaning of a density–dependent interaction strength. A similar correspondence was previously pointed out for homogeneous Bose–Fermi mixtures in the continuum Viverit .

Introducing a perturbation expansion with respect to $`J_B`$ around the Mott–insulating ground state we can recover the zero–temperature phase transition to the superfluid phase. The reverse, i.e., to build a perturbative expansion in powers of $`U_{BB}`$ around the superfluid ground state fails to describe the transition to a Mott insulator, as pointed out in Ref. OosStr02  for the pure Bose case. We follow the procedure adopted in Ref. CheWu03  for the two–component Boson mixture, with the due modifications for the present case of a Boson–Fermion mixture, by treating the Bosonic kinetic (hopping) term as the perturbation with respect to the Bosonic Mott–insulating ground state. This scheme was first introduced for one–component Bose systems in Refs. OosStr02 ; Freericks ; SheKri93 . We proceed by expanding the ground state energy with respect to the (local) Bosonic superfluid parameter $`\psi ^{(i)}=a_i`$. At the phase boundary between a Mott insulator (MI) and a superfluid (SF) the expansion coefficients must vanish, yielding the following criterion for the onset of the transition to the (local) SF state:

$`U_{BB}(2n_B^{(i)}1)2J_B`$
$`\left(U_{BB}^24U_{BB}^2(2n_B^{(i)}+1)+4J_B^2\right)^{1/2}`$
$`<\mu _BV_B^{(i)}U_{BF}n_F^{(i)}`$
$`<U_{BB}(2n_B^{(i)}1)2J_B`$
$`+\left(U_{BB}^24U_{BB}^2(2n_B^{(i)}+1)+4J_B^2\right)^{1/2}.`$ (51)
The minimum value of $`U_{BB}/J_B`$, where a MI phase can exist, is given by the condition

$$U_{BB}/J_B=4n_B^{(i)}+2+2\sqrt{(2n_B^{(i)}+1)^21},$$
(52)
and it involves the Fermionic sector indirectly through the dependence of $`n_B^{(i)}`$ on the Fermionic parameters and density distributions provided by Eqns. (39)–(40). Apart from this important modification, the phase diagram, at this level of approximation, is analogous to that of a one–component Bose system.

## IV Number-Conserving Gutzwiller Ansatz and Numerical Analysis

The simplest ansatz for the Boson ground state being capable of describing both the SF and the MI phases is the Gutzwiller Ansatz, which contains the mean field approximations previously discussed as special cases. It consists of factorizing the amplitudes of superpositions of all possible Fock states consistent with a fixed number of Bosons $`N_B`$, in the following way KrauCaf91 :

$`|\mathrm{\Psi }_B`$ $``$ $`{\displaystyle \underset{_jn_j=N_B}{}}{\displaystyle \underset{i}{}}f_{n_i}^{(i)}{\displaystyle \frac{(\widehat{a}_i^{})^{n_i}}{\sqrt{n_i!}}}|0.`$ (53)
In this way correlations between the Bosons under particle exchange are included in a natural way. Using the number–conserving Gutzwiller ansatz in the determination of the energy functional, while keeping the same approximations previously introduced for the Boson–Fermion interaction and the Fermion energy, the total ground state energy reads

$`E=E_B+E_F+U_{BF}{\displaystyle \underset{i}{}}n_B^{(i)}n_F^{(i)},`$ (54)
where the subsidiary conditions ensuring particle number conservation are

$`{\displaystyle \underset{i}{}}n_B^{(i)}`$ $`=`$ $`{\displaystyle \underset{i}{}}\widehat{a}_i^{}\widehat{a}_i=N_B,`$ (55)
$`{\displaystyle \underset{i}{}}n_F^{(i)}`$ $`=`$ $`{\displaystyle \underset{i}{}}\widehat{b}_i^{}\widehat{b}_i=N_F.`$ (56)
The Boson energy contribution is now

$`E_B`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_B\left({\displaystyle \underset{i}{}}\psi _{}^{(i+1)}{}_{}{}^{}\psi ^{(i)}+\text{C. c.}\right)`$ (57)
$`+`$ $`{\displaystyle \frac{U_{BB}}{2}}(\sigma _B^{(i)}n_B^{(i)})+V_B^{(i)}n_B^{(i)},`$
and the Bosonic observables are related to the Gutzwiller amplitudes by

$`n_B^{(i)}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}n(f_n^{(i)})^2,`$ (58)
$`\sigma _B^{(i)}`$ $`=`$ $`\widehat{a}_i^{}\widehat{a}_i\widehat{a}_i^{}\widehat{a}_i={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}n^2(f_n^{(i)})^2,`$ (59)
$`\psi ^{(i)}`$ $`=`$ $`\widehat{a}_i={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\sqrt{n+1}f_n^{(i)}f_{n+1}^{(i)}.`$ (60)
Moreover, we must impose the natural constraints that

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(f_n^{(i)})^2`$ $`=`$ $`1,`$ (61)
$`0n_F^{(i)}`$ $``$ $`1,`$ (62)
for each lattice site $`i`$, reflecting the fact that the Gutzwiller amplitudes form a probability distribution for each lattice site, and that the on site Fermion occupation number cannot exceed one.

To identify the ground state amounts to solving a constrained optimization problem: one has to minimize the energy functional (57) subject to the constraints given by Eqns. (55) and (56), together with Eqns. (61) and (62).

We have solved the problem numerically for a small system of ten particles (five Bosons and five Fermions). The first observation is that the optimization problem is not a convex optimization problem. Hence, one has to expect several local, “poorer” extrema in addition to the (not necessarily unique) global one. The numerical solution of this optimization problem has been performed first using a simulated annealing method KirGelVec83  with an appropriate logarithmic annealing schedule. The quadratic constraints (61) and (62) have been incorporated in a dynamical penalty formulation (see, e.g., Ref. WahWan99 ). Finally, for the local refinement the Nelder–Mead downhill simplex method NelMea65  has been applied.

In FIG. 2 we show the change of the on site Bosonic densities with increasing lattice potential strength $`\stackrel{~}{V}_0`$ for a system of five Bosons and five Fermions with moderate repulsive Boson–Fermion interaction. We note from FIG. 2 that as the strength of the lattice potential increases the Bosons go in a complete Mott–insulating phase, forming a block crystalline configuration around the center of the trap (which coincides with the origin of the optical lattice) with exactly one Boson per lattice site. The corresponding on site Fermionic densities are plotted in FIG. 3. From both figures we can see that, if $`U_{BF}>0`$, by increasing the lattice potential strength the system eventually undergoes simultaneously a Boson MI transition and complete phase separation, in accordance with Eqn. (50) along with Eqns. (24) – (28).

The local Bosonic superfluid parameter $`\psi ^{(i)}`$ for the same physical situation is shown in FIG. 4. Although we are dealing with a finite system, we can already see a rather clear signature of the onset of a phase transition to a Mott–insulator regime when the superfluid parameter suddenly drops to very low values at an approximate critical lattice potential strength $`\stackrel{~}{V}_0^c7`$.

We next consider the ground–state properties in the case of an attractive Boson–Fermion interaction. Because of the strong attraction with growing lattice depth, the Fermions follow the Bosons in building a sharp crystalline block around the center of the trap, as can be seen from FIGs. 5 and 6. We cannot expect in this case to observe a simultaneous mean field collapse like the one predicted for a trapped Bose–Fermi mixture in the continuum Roth ; Cazalilla  (for the effects beyond mean field see Uns2 ), as this possibility is forbidden in a single–band approximation. Finally, we consider the behavior of the Bosonic superfluid on site parameter in the case of a Boson–Fermion attractive interaction.

Comparing FIG. 7 with FIG. 4, we see that the transition to a Mott insulating phase for the Bosons takes place at the same lattice potential strength, irrespectively of the repulsive or attractive nature of the Boson–Fermion interaction. This finding confirms the results of the mean field analysis presented in the previous Section.

## V Mirror symmetry breaking and transition to degeneracy

The above optimization problem associated with the constrained minimization of the energy is not convex, hence there can be many local minima in addition to the global one. However, even the ground state may be approximately or exactly degenerate. In fact, this is what happens in the case of Boson–Boson and Boson–Fermion repulsion for large values of the lattice potential strength $`\stackrel{~}{V}_0`$. As $`\stackrel{~}{V}_0`$ grows, it becomes eventually energetically more favorable for the bosons to be arranged in single–particle occupancy of the available sites around the center of the external trap. The Bosonic and Fermionic on–site occupation numbers can only assume the values $`0`$ or $`1`$, and a definite Boson–Fermion symmetry is established in the Bose–Fermi Hubbard Hamiltonian assuming that the on site Fermionic and Bosonic trapping potentials coincide.

A similar transmutation of Bosons into Fermions in strong optical lattices has been pointed out by Paredes and Cirac in a recent paper Pared02 . They consider a model of pure Bosons in an optical lattice and show that in the limit of very strong Boson–Boson on site interaction, the Bosonic operators can be mapped into Fermionic operators by means of the well known Jordan–Wigner transformation. Let us consider what happens in the case of a Boson–Fermion mixture. As the lattice strength grows, configurations of lowest energy that are mirror–symmetric with respect to the center of the lattice, like e.g. those of FIGs. 2 and 3, become approximately energetically equivalent to other symmetric configurations (e.g., a checkerboard of alternating Bosons and Fermions with one particle per lattice site), as well as to nonsymmetric configurations (like a succession of four Fermions followed by five Bosons and then a last Fermion, again with one particle per lattice site), and mirror symmetry breaking takes place. We may thus consider sequences of energy functionals with increasing lattice potential strengths $`\stackrel{~}{V}_0`$. For each value of $`\stackrel{~}{V}_0`$, one may identify a ground state. Then, the difference in energy of this ground state to those states that can be obtained by interchanging the role of Fermions and Bosons will converge to zero as $`\stackrel{~}{V}_0`$ grows. The Boson hopping contribution will become negligible, whereas the behavior of $`\stackrel{~}{U}_{BB}`$ will enforce the mean Bosonic on site occupation number to be at most one. Hence, for each lattice site, the constraints on the Boson and Fermion occupation numbers become identical (at most one Boson or one Fermion per lattice site). Notice that the suppression of the hopping terms is exponential. Moreover, since $`\stackrel{~}{V}_B^{(i)}=\stackrel{~}{V}_F^{(i)}`$ for all lattice sites $`i`$, the larger the value of $`\stackrel{~}{V}_0`$, the more symmetric is the role of Bosons and Fermions. There are then many ground states that are degenerate in energy with respect to any permutation of lattice sites – as long as all particles are located around the minimum of the confining external potential $`\stackrel{~}{V}_B^{(0)}=\stackrel{~}{V}_F^{(0)}=0`$. These degenerate configurations will be given by all possible symmetric and nonsymmetric Fermion and Boson distributions in a region around the center of the lattice, with every site of the region occupied by one and only one particle. Such possible configurations are for example checkerboard alternating patterns of Bosons and Fermions, or Mott Bosonic central configurations with Fermionic wings on the sides, or consecutive block crystalline arrangements of variable length of Bosons and Fermions. In brief, while the Hamiltonian formally retains its mirror symmetry under reflection of the lattice around its center, the degenerate ground states need not, and spontaneous mirror symmetry breaking occurs. At the same time complete Boson–Fermion exchange symmetry sets on. No ground state is a priori favored compared to any other: any random pattern of consecutive Bosons and Fermions located around the minimum of the external trapping potential is a legitimate ground state. FIG. 8 shows representative on site Bosonic densities in the regime of large values of $`\stackrel{~}{V}_0`$ around $`\stackrel{~}{V}_0=50`$ for the case of Boson–Boson and Boson–Fermion repulsion in a system composed of five Bosons and five Fermions: at each value of the lattice potential strength, a particular state is selected from the set of those with same energy. Each vanishing value of the on site Bosonic density means that exactly one Fermion has filled that particular lattice site. The large value chosen for $`\stackrel{~}{V}_0`$ allows to clearly stress the random nature of the configuration patterns even for very small changes of the lattice potential strength, whereas degeneracy and disorder can set in already at lower values of the lattice depth, depending on the tuning of the harmonic oscillator and scattering lengths (see below).

The degenerate states are separated by energy barriers. The system is non–ergodic, and hysteresis should be observed: what particular state is chosen, depends on the exact mechanism of preparation of the system and of loading of the mixture into the optical lattice.

The criterion for the onset of degeneracy and nonperiodic ground states in the bulk region around the center of the lattice and of the trapping potential is easily identified, by looking at the relative importance of the trapping on site energy with respect to the on site Boson or Fermion interaction energy. For instance, to allow for the Fermionic behavior of the Boson on site occupation numbers (either $`0`$ or $`1`$) one must require that the energy is lower having one Boson at the edge of the bulk central region rather than having it sitting on top of another Boson at the center of the lattice:

$$\stackrel{~}{U}_{BB}>\stackrel{~}{V}_B\left(i=(N_B+N_F)/2\right).$$
(63)
The analogous condition for the Bose–Fermi on site interaction is:

$$\stackrel{~}{U}_{BF}>\stackrel{~}{V}_{B/F}\left(i=(N_B+N_F)/2\right).$$
(64)
For smaller values of $`\stackrel{~}{V}_0`$, the Boson hopping contribution will become more and more important. A representative situation of this intermediate regime is depicted in FIGs. 2 and 3: here, the repulsion between Bosons and Fermions is strong enough to allow for phase separation, while the non–negligible hopping terms still favor configurations where Bosons have Bosons as nearest neighbors. The transition to degeneracy and disorder, exact in the limit of infinite lattice depth, is a novel and peculiar feature of Bose–Fermi mixtures and it should hold in general for any multicomponent Bose and/or Fermi dilute atomic system loaded in a deep optical lattice at zero temperature, provided that intercomponent interactions are repulsive and the on site confining potentials coincide for the different components. It clearly cannot take place in a single–component system, say a pure single–component Bose gas, where only a SF–MI transition occurs JacBru98 . The rather complex and rich interplay between ordered and disordered configurations of Bose–Fermi mixtures in very deep optical lattices will be considered in more detail elsewhere.

## VI Summary and Outlook

In conclusion, we have studied the zero–temperature properties of a mixture of weakly interacting gases of neutral Bosonic and Fermionic atoms loaded in one–dimensional optical lattices and confined by harmonic trapping potentials. We have derived a single–band Bose–Fermi Hubbard Hamiltonian, and performed some mean field studies of the zero–temperature phase diagram. We have considered the case of a quasi–free Fermion sea acting on the Bosons, which have been treated in their full dynamical range. We have always worked in the approximation of $`s`$–wave Boson–Boson and Boson–Fermion contact interactions. According to the possible different combinations of intraspecies and interspecies attractive and repulsive interactions, the system displays a rich phase structure, including the onset of a SF–MI transition in the Boson sector, and a simultaneous transition to demixing in the Boson–Fermion sector. The optical lattice potential plays a crucial role, allowing to tune the system into regimes of strong Boson–Boson and Boson–Fermion couplings as the lattice depth is increased. For very deep lattices the system displays a remarkable transition to a multiply degenerate phase in which all possible permutations of configurations with one Bosonic or Fermionic atom per site are legitimate ground states. The transition is related with breaking of the lattice mirror symmetry for very large values of the lattice depth. This peculiar disordered pattern of degenerate ground–state configurations separated by very large barriers is somehow reminiscent of the behavior of classical disordered systems like glasses and spin glasses, but it takes place in a quantum system at zero temperature.

The setting that has been investigated in detail in the present paper can be extended in various ways. Certainly a larger number of Bosons and Fermions have to be considered in order to obtain a more realistic description of the system. While our previous analytical findings are applicable to any numbers of atoms, in order to extend the numerical calculations to larger numbers more powerful numerical methods have to be introduced. So far, Monte Carlo simulations with a fairly large number of particles have been carried out only for an inhomogeneous Bose–Hubbard model Batrouni . The authors have also speculated that the qualitative phase diagram does not depend on the dimensionality of the system.

In order to extend the work presented in this paper to the case of interacting Fermions, one may either allow for different Fermionic species in magnetic traps or for spin–unpolarized identical Fermions in optical traps. Exchange–correlation effects, that are already included in the Gutzwiller Ansatz for the Bosons, will then become important for the Fermions as well Lieb . In this way the Bose–Fermi interaction has to be consistently incorporated beyond the mean–field level. In the present paper the mean–field treatment of the Bose–Fermi interaction is consistent, since we only included the Fermions in a mean–field and LDA prescription.

At the opposite extreme, one can consider the Fermions to be static impurities for the Bosons, as very recently discussed by Vignolo and collaborators Vignolo . Finally, after the first version of the present paper had been submitted, a preprint appeared on the induced Boson–Boson interaction due to the Fermions for mixtures loaded in a $`2`$$`D`$ lattice Blatter .

Besides these fundamental theoretical aspects related to the theory of quantum phase transitions and the statistical mechanics of complex systems, ultracold Bose–Fermi mixtures in an optical lattice qualify for potential applications in the physics of quantum information. As with systems involving either Bosons or Fermions that have been studied so far JakBriCirGarZol98 ; Deutsch ; PacKni03 ; GarCir03 ; DuaDemLuk03 , mixtures could be used for the preparation of multi-particle entangled states DuaDemLuk03  such as cluster states or certain instances of graph states Briegel , as well as for the implementation of quantum gates. With Bosons and Fermions serving two different purposes, Bose-Fermi mixtures could in fact allow for new possibilities of quantum information processing in optical lattices. The Fermions would be suitable for storage of quantum information due to their non–interacting behavior, whereas the Bosons could be used to let the systems interact and perform operations.

In turn, the study of such coupled quantum systems exhibiting collective phenomena with the methods of the theory of multi-particle entanglement is an attractive investigation in its own right. The aim here is not to use the coupled system to prepare strongly entangled systems that form the starting point for applications in quantum information processing. Instead, the motivation of such investigations is to go beyond conventional methods to characterize the natural correlations present in the distributed system at zero temperature Osborne ; Vidal ; Chain ; Zanardi .

The quantitative theory of entanglement that abstracts from the actual physical realization of a Bose–Fermi Hubbard model could provide the tools to understand how global properties emerge here from quantum correlations between the elementary constituents.

## VII Acknowledgments

We thank Kai Bongs and Giovanni Modugno for very useful discussions. A. A. and J. E. thank the DFG and the ESF for financial support. F. I. thanks the INFM for financial support as well as COSLAB and BEC2000+ ESF programs.


text:
# Anomalous fluctuations in phases with a broken continuous symmetry

## Abstract

It is shown that the Goldstone modes associated with a broken continuous symmetry lead to anomalously large fluctuations of the zero field order parameter at any temperature below $`T_c`$. In dimensions $`2<d<4`$, the variance of the extensive spontaneous magnetization scales as $`L^4`$ with the system size $`L`$, independent of the order parameter dynamics. The anomalous scaling is a consequence of the $`1/q^{4d}`$ divergence of the longitudinal susceptibility. For ground states in two dimensions with Goldstone modes vanishing linearly with momentum, the dynamical susceptibility contains a singular contribution $`(q^2\omega ^2/c^2)^{1/2}`$. The dynamic structure factor thus exhibits a critical continuum above the undamped spin wave pole, which may be detected by neutron scattering in the Néel-phase of 2D quantum antiferromagnets.

It is one of the basic properties of any thermodynamic system that the fluctuations Var$`\widehat{A}=\widehat{A}^2\widehat{A}^2V`$ of an extensive variable $`\widehat{A}`$ scale linearly with the system volume $`V`$. This property guarantees that the intensive variables $`\widehat{A}/V`$ are self-averaging, with rms fluctuations vanishing proportional to $`V^{1/2}`$. Physically, the extensive nature of Var$`\widehat{A}V`$ is due to the existence of a finite, microscopic correlation length $`\xi `$. A large system can thus be partitioned into an extensive number of $`V/\xi ^d`$ subvolumes which are statistically independent. Since the variance in each subvolume is finite, the central limit theorem then quite generally implies a Gaussian distribution for $`\widehat{A}`$ with a variance of order $`V`$, in agreement with the standard Einstein theory of fluctuations in macroscopic thermodynamics.

The above argument indicates that the linear scaling Var$`\widehat{A}V`$ may break down only at a critical point of a continuous phase transition where the correlation length $`\xi `$ diverges. In fact, from the standard relation Var$`\widehat{M}=VT\chi `$ between the fluctuations of the total ’magnetization’ $`\widehat{M}`$ and the corresponding linear susceptibility, the finite size scaling $`\chi (T_c,L)L^{2\eta }`$ of the susceptibilty right at $`T_c`$ barber  implies a nontrivial dependence Var$`\widehat{M}L^{d+2\eta }`$ on system size $`L`$ at the critical point. Below, it will be shown that anomalous fluctuations of the order parameter are present not only at $`T_c`$ but in fact at any temperature below $`T_c`$ provided the broken symmetry is continuous. This is a result of the presence of Goldstone modes, which imply that correlations decay algebraically below $`T_c`$ with exponents that are independent of temperature. A specific example of recent interest are fluctuations of the condensate number $`\widehat{N}_0`$ in a Bose Einstein condensate. Using a weak coupling Bogoliubov approach, Giorgini et.al have shown giorgini  that Var$`\widehat{N}_0T^2L^4`$ scale anomalously at low temperatures. In fact, this result also applies to strongly interacting superfluids, being essentially a consequence of the long wavelength phase fluctuations meier . In the present Letter, we will show, that the anomalous fluctuations in a Bose-Einstein condensate are just one particular example of a rather general phenomenon which appears in any phase with a broken continuous symmetry. Consider, for example, an isotropic ferromagnet in zero field below $`T_c`$. Because of the invariance under spin rotations, it costs no energy to rotate the direction of the magnetization vector. Provided that the microscopic interactions are short ranged, this implies that the transverse susceptibility $`\chi _{}(q)`$ diverges precisely as $`1/q^2`$ at small wavevectors $`q0`$. As a result, there are strong fluctuations of the direction of the magnetization, leading to a complete destruction of long range order at finite temperature in dimensions $`d2`$, the well known Mermin-Wagner-Hohenberg theorem. Regarding the magnitude of the magnetization, the naive expectation is, that its fluctuations are just like that of a standard thermodynamic variable because there is a finite restoring force for deviations from the equilibrium value. However, as was noted a long time ago by Patashinski and Pokrovski patashinski , the inevitable coupling between longitudinal and transverse order parameter fluctuations entails that the longitudinal susceptibility is also singular at small wave vectors, diverging as $`\chi _{}(q)1/q^{4d}`$ in dimensions $`2<d<4`$. Our aim in the following is to show that
a) the $`1/q^{4d}`$ divergence of the longitudinal susceptibility leads to anomalous fluctuations Var$`\widehat{M}_sL^4`$ of the zero field order parameter in $`2<d<4`$ and $`T0`$ for an arbitrary broken continuous symmetry phase. For superfluids, the corresponding relative condensate fluctuations are universal at low $`T`$.
b) there is an analog of the singular nature of $`\chi _{}`$ for zero temperature phases with a broken continuous symmetry. In two dimensions and with Goldstone modes whose frequencies $`\omega =cq`$ vanish linearly with momentum $`q=|\stackrel{}{q}|`$, the dynamical susceptibility has a singular contribution $`(q^2\omega ^2/c^2)^{1/2}`$ as noted first by Sachdev on the basis of a $`1/N`$ expansion sachdev . The dynamic structure factor of 2D quantum antiferromagnets thus exhibits a critical continuum above the standard $`\delta `$-function spin wave peak. This effect may be observed in high resolution neutron scattering experiments.

As an effective description of an arbitrary phase with a broken continuous symmetry at finite temperature we use the nonlinear $`\sigma `$-model parisi . It describes directional fluctuations of an order parameter $`\stackrel{}{\psi }(x)=m_s^{(0)}\stackrel{}{\mathrm{\Omega }}(x)`$ with a fixed magnitude $`m_s^{(0)}`$ in terms of an $`N`$-component unit vector $`\stackrel{}{\mathrm{\Omega }}(x)`$. At zero external field, the effective action for the fluctuations of $`\stackrel{}{\mathrm{\Omega }}`$ is

$$S\left[\mathrm{\Omega }\right]=\frac{\rho _s}{2T}d^dx\left(\stackrel{}{\mathrm{\Omega }}(x)\right)^2|\stackrel{}{\mathrm{\Omega }}(x)|^2=1,$$
(1)
with the spin stiffness (or helicity modulus) $`\rho _s`$ as the single phenomenological parameter. In a finite system of volume $`V=L^d`$, there is of course no spontaneous magnetization at zero field. The breaking of a continuous symmetry below $`T_c`$ however shows up in the integrated $`O(N)`$-symmetric correlation function at finite volume

$$_Vd^dx\stackrel{}{\psi }(x)\stackrel{}{\psi }(0)_L=Vm_L^2Vm_s^2,$$
(2)
defining an intensive nonzero order parameter $`m_L^2`$ which approaches the spontaneous magnetization $`m_s^2`$ of the infinite system in the thermodynamic limit $`L\mathrm{}`$. For superfluids this is just the number of particles in the condensate. In dimensions $`d>2`$, where $`m_s`$ is nonzero at finite $`T`$, the leading long distance behaviour of the two point function $`G(r)=\stackrel{}{\psi }(x)\stackrel{}{\psi }(0)`$ may be obtained from a simple Gaussian spin wave calculation. Following a standard procedure parisi , the unit vector $`\stackrel{}{\mathrm{\Omega }}(x)=(\stackrel{}{\mathrm{\Pi }}(x),\sqrt{1\stackrel{}{\mathrm{\Pi }}(x)^2})`$ is decomposed into $`N1`$ ’transverse’ Goldstone fields $`\mathrm{\Pi }_a(x),a=1,\mathrm{}N1`$ and a longitudinal component $`\mathrm{\Omega }_N(x)=\sqrt{1\stackrel{}{\mathrm{\Pi }}(x)^2}`$. At low enough temperatures, where spin wave interactions can be neglected, the $`\mathrm{\Pi }`$ fields are Gaussian random variables with variance $`\mathrm{\Pi }_a(q)\mathrm{\Pi }_b(q^{})=\delta _{a,b}\delta _{q,q^{}}T/\rho _sq^2`$. The zero field correlation function below $`T_c`$

$$G(r)=m_s^2\left[1+C_{}(r)+(N1)C_{}(r)\right]$$
(3)
is thus split into longitudinal and transverse parts, with $`m_s^2=G(\mathrm{})`$ the renormalized value of the spontaneous magnetization. To lowest nontrivial order in the small fluctuations $`\mathrm{\Pi }`$, the transverse correlation function decays proportional to $`T/\rho _sr^{d2}`$, as expected from the standard $`1/q^2`$ divergence of the transverse susceptibility

$$\chi _{}(q)=m_s^2C_{}(q)/T=\frac{m_s^2}{\rho _sq^2}$$
(4)
in the symmetry broken phase. In leading order, the longitudinal function is given by ($`c`$ stands for ’connected’)

$$C_{}(r)\frac{1}{4}\stackrel{}{\mathrm{\Pi }}^2(x)\stackrel{}{\mathrm{\Pi }}^2(0)_c=\frac{N1}{2}C_{}^2(r)$$
(5)
and thus again decays algebraically as $`1/r^{2(d2)}`$. Contrary to the naive mean field picture, where the londitudinal susceptibility $`\chi _{}(q)`$ below $`T_c`$ is finite as $`q0`$, the slow decay of $`C_{}(r)`$ implies that $`\chi _{}(q0)T/\rho _s^2q^{4d}`$ is divergent in $`d<4`$ patashinski . While the relation (5) is valid only to leading order in an expansion in powers of $`\mathrm{\Pi }`$, the behavior $`C_{}(r)r^{2(d2)}`$ is in fact expected to be exact at arbitrary temperatures below $`T_c`$ parisi , consistent with a rigorous correlation inequality for $`N=2`$ dunlop .

To discuss the fluctuations of the spontaneous magnetization in a finite system at zero field, we define a fluctuating, extensive variable $`(d1=d^dx_1)`$

$$\widehat{M}_s=\frac{1}{V}d1d2\stackrel{}{\psi }(1)\stackrel{}{\psi }(2)$$
(6)
with thermal average $`\widehat{M}_s=Vm_L^2`$. Its fluctuations are then determined by the connected four point function

$$u_4(1234)=\stackrel{}{\psi }(1)\stackrel{}{\psi }(2)\stackrel{}{\psi }(3)\stackrel{}{\psi }(4)_c.$$
(7)
Expanding this consistently up to order $`\mathrm{\Pi }^4`$, one finds that in an infinite system this function may be expressed simply as

$`u_4(1234)=m_s^4{\displaystyle \frac{N1}{2}}`$ (8)
$`\left[C_{}(r_{13})C_{}(r_{14})C_{}(r_{23})+C_{}(r_{24})\right]^2.`$
In order to obtain the scaling of Var $`\widehat{M}_s`$ in a finite system, we switch to a momentum representation and replace integrals $`_q`$ by discrete sums $`V^1_q^{}`$ over wavevectors. The $`q=0`$ contribution is excluded since it describes an irrelevant global rotation of the order parameter $`\stackrel{}{\psi }(x)`$. It is then straightforward to show that

$$\mathrm{Var}\widehat{M}_s=2m_s^4(N1)\left(\frac{T}{\rho _s}\right)^2\underset{q}{}^{}\frac{1}{q^4}$$
(9)
which is proportional to $`L^4`$ in $`2<d<4`$ by simple dimensional analysis note1 . Defining a numerical coefficient $`B`$ by $`_q^{}q^4=BL^4/8\pi ^2`$, we find that $`B=8E_3(2)/\pi ^2=0.501`$ for a 3D box with Dirichlet boundary conditions. Here $`E_d(t)=_{n_1=1,\mathrm{}n_d=1}^{\mathrm{}}(n_1^2+\mathrm{}+n_d^2)^t`$ is the generalized Epstein zeta function, convergent for $`d<2t`$ (for periodic boundary conditions $`B=0.8375`$ meier ).

Although our derivation of the general result $`(9)`$ is based on an expansion in powers of $`\mathrm{\Pi }`$ and thus appears to be restricted to low temperature, the exponent in Var $`\widehat{M}_sL^4`$ is universal below $`T_c`$ note2  just like the $`q^{(4d)}`$ divergence of the longitudinal susceptibility. Similarly, the temperature dependence of Var $`\widehat{M}_s`$ will be $`\left(T/\rho _s(T)\right)^2`$ for arbitrary temperatures below $`T_c`$, vanishing proportional to $`T^2`$ at low temperatures since $`\rho _s(T=0)`$ is finite. This follows from the fact that at any temperature $`T`$ in the broken symmetry phase, the dominant finite size dependence is determined by the leading low energy constant in the effective field theory for fluctuations of the order parameter, which is precisely $`\rho _s`$ as defined in equation $`(1)`$. New effective constants only enter in higher order terms, as discussed for instance in the context of chiral perturbation theory in QCD hasenfratz . From eqn.$`(9)`$, it follows that the relative fluctuations Var $`\widehat{M}_s`$/$`\widehat{M}_s^2`$ scale as $`L^{2(d2)}`$. The spontaneous magnetization is therefore self-averaging in $`d>2`$ although weaker than expected naively, unless $`d>4`$. In this context, it should be mentioned that, in the case of a non-Abelian symmetry $`N3`$, the spin stiffness $`\rho _s`$ has anomalous fluctuations

$$\frac{\mathrm{Var}\rho _s}{\rho _s^2}(N2)\left(\frac{T}{\rho _sL^{(d2)}}\right)^2$$
(10)
on its own, as shown by Chakravarty chakravarty . They are very similar to the fluctuations of the spontaneous magnetization discussed here and indeed arise from the same kind of Goldstone anomalies. A special situation appears in homogeneous superfluids $`(N=2)`$ at low temperatures. Since translation invariance requires the superfluid density to be equal to the full density $`n`$, the associated stiffness $`\rho _s(T0)=\mathrm{}^2n/m`$ is independent of the interaction. Equation $`(9)`$ then leads to the remarkable result that the relative fluctuations of the condensate number at low temperature $`(d=3)`$

$$\underset{T0}{lim}\frac{\mathrm{Var}\widehat{N}_0}{\widehat{N}_0^2}=\frac{B}{\left(n\lambda _T^2L\right)^2}$$
(11)
are completely universal, depending only on density, system size and the thermal wavelength $`\lambda _T=h/\sqrt{2\pi mT}`$, but not on whether the superfluid is weakly or strongly interacting remark . Finally, it is important to realize that - as in the case of the Mermin-Wagner-Hohenberg theorem - the result $`(9)`$ is independent of the order parameter dynamics. It applies equally to say ferro- or antiferromagnets even though the temperature dependence of the average order parameter $`\widehat{M}_s`$ is, of course, very different in both cases. This may be shown either by a microscopic derivation of equation $`(9)`$ on the basis of noninteracting spin wave theory around a perfectly ordered ground state or, alternatively, by using a quantum mechanical generalization of the nonlinear $`\sigma `$ model which describes magnons with dispersion either $`\omega q^2`$ or $`\omega q`$.

As shown above, the anomalous fluctuations $`(9)`$ of the order parameter are a consequence of the $`q^{(4d)}`$ divergence of the longitudinal susceptibility in $`d<4`$. In the following we want to discuss the analog of this phenomenon in zero temperature phases where a continuous symmetry is broken. In this case the specific dynamics of the order parameter is important. Here, we assume that the Goldstone modes have a linear spectrum $`\omega =cq`$ as in superfluids and antiferromagnets. An effective description of the ordered phase can then be obtained from a quantum mechanical nonlinear $`\sigma `$model, with a unit vector $`\stackrel{}{\mathrm{\Omega }}(x,\tau )`$ which depends both on space $`x`$ and imaginary time $`\tau [0,\beta \mathrm{}]`$. The corresponding effective action

$$S\left[\mathrm{\Omega }\right]=\frac{\rho _s}{2\mathrm{}}_0^\beta \mathrm{}𝑑\tau d^dx\left[\left(\stackrel{}{\mathrm{\Omega }}(x)\right)^2+\left(\frac{1}{c}_\tau \stackrel{}{\mathrm{\Omega }}\right)^2\right]$$
(12)
has the spin wave velocity $`c`$ as the only additional parameter. Together with the renormalized value $`m_s^2`$ of the long range order, $`\rho _s`$ and $`c`$ completely determine the low energy properties of the ordered phase. As an example, this model applies both to the Néel ordered state of 2D quantum antiferromagnets discussed extensively in the context of high temperature superconductors chn  or to the superfluid phase of cold atoms in an optical lattice which has been realized recently in 3D greiner . Using again a lowest order expansion in powers of the small fluctuations $`\stackrel{}{\mathrm{\Pi }}(x,\tau )`$ in the standard decomposition $`\stackrel{}{\mathrm{\Omega }}=(\stackrel{}{\mathrm{\Pi }},\sqrt{1\stackrel{}{\mathrm{\Pi }}^2})`$, the transverse correlation function at $`T=0`$ in two dimensions is given by note3

$$C_{}(x,\tau )=\frac{\mathrm{}c}{4\pi \rho _s}\frac{1}{\sqrt{r^2+(c\tau )^2}}.$$
(13)
Analytic continuation to real time $`t=i\tau `$ and Fouriertransformation give the standard form of the transverse dynamical susceptibility in any dimension

$$\chi _{}(q,\omega )=\frac{m_s^2}{\rho _s}\frac{1}{q^2\omega ^2/c^2}.$$
(14)
It leads to the expected quasiparticle pole at $`\omega =cq`$, reflecting the presence of undamped antiferromagnetic spin waves. Considering the longitudinal correlations, the analog of the relation $`(5)`$ again applies to leading order. As a result, the longitudinal dynamical susceptibility in 2D turns out to be

$$\chi _{}(q,\omega )=\frac{(N1)m_s^2\mathrm{}c}{16\rho _s^2}\frac{1}{\sqrt{q^2\omega ^2/c^2}},$$
(15)
which has a branch cut rather than a simple pole. Formally the result $`(15)`$ is completely analogous to the $`1/q`$ divergence of the classical static susceptibility $`\chi _{}(q)`$ in $`d=3`$ discussed above. Indeed the effective dimensionality of the 2d quantum antiferromagnet is $`d+z=3`$ and the dependence on frequency $`\omega `$ is dictated by the formal Lorentz-invariance of the model $`(12)`$.

In order to relate these results to directly observable quantities, we consider the rotationally averaged staggered susceptibility (for antiferromagnets $`N=3`$)

$$\chi _s(q,\omega )=\frac{N1}{N}\chi _{}(q,\omega )+\frac{1}{N}\chi _{}(q,\omega ).$$
(16)
For frequencies $`\mathrm{}\omega T`$, its imaginary part is equal to the dynamic structure factor $`S(q,\omega )`$ at $`\omega >0`$ up to a factor $`2\mathrm{}`$, giving

$$S(q,\omega )=2m_s^2\xi _J\frac{N1}{N}\left[\frac{\pi }{2q}\delta (\omega cq)+\frac{\xi _J}{16}\frac{\theta (\omega cq)}{\sqrt{\omega ^2c^2q^2}}\right].$$
(17)
The longitudinal fluctuations of the Néel order thus lead to a critical continuum above the spin wave pole at $`\omega =cq`$, which decays only algebraically. The continuum results from the decay of a normally massive amplitude mode with momentum $`\stackrel{}{p}`$ into a pair of spin waves with momenta $`\stackrel{}{q}`$ and $`\stackrel{}{p}\stackrel{}{q}`$, which is possible for any $`\omega >cq`$, with a singular cross section because of the large phase space. The amplitude mode is thus completely overdamped in two dimensions sachdev . The relative weight of these fluctuations compared to the dominant transverse contribution is determined by the dimensionless parameter $`q\xi _J`$, with $`\xi _J=\mathrm{}c/\rho _s`$ the Josephson correlation length. This length controls the crossover from the Goldstone regime $`q\xi _J1`$, where the spin dynamics is well described by small fluctuations around a Néel ordered ground state to quantum critical fluctuations at $`q\xi _J1`$. The quantum critical regime has a dynamical susceptibility of the form csy

$$\chi _s(q,\omega )=\frac{m_s^2}{\rho _s}\left(\frac{N\xi _J}{2\pi }\right)^\eta \frac{A_Q}{\left(q^2\omega ^2/c^2\right)^{1\eta /2}}$$
(18)
with an amplitude $`A_Q`$ close to one and an exponent $`\eta `$ which is nearly zero. Now, in the regime where the ground state exhibits strong Néel order, the Josephson correlation length is only several lattice constants chn  and thus the quantum critical regime is hardly accessible. In turn, there is a rather wide Goldstone regime characterized by $`q\xi _J1`$ and $`q\xi (T)1`$, where $`\xi (T)`$ is the finite 2d correlation length over which Néel order is lost at finite temperature. Due to the exponential dependence chn

$$\xi (T)\xi _J\mathrm{exp}(2\pi \rho _s/T)$$
(19)
the relevant range of wavevectors $`\xi ^1(T)<q<\xi _J^1`$ is rather wide at low temperatures. Experimentally, spin waves in 2D quantum antiferromagnets have been observed by inelastic neutron scattering in the vicinity of the Néel ordering wave vector $`(\pi ,\pi )`$. In a constant $`\omega `$ scan, the first term in $`(17)`$ gives rise to a sharp peak at $`q=\omega /c`$ with an amplitude $`1/q`$ aeppli . The second contribution in $`(17)`$ due to the longitudinal fluctuations implies an additional algebraic tail towards smaller wave vectors, provided $`q`$ is in the Goldstone regime. At room temperature, where $`\xi (T)`$ is several hundred Å and with typical values of $`\xi _J`$ this gives a range $`0.005<q`$$`{}_{}{}^{1})<0.05`$. Given the resolution in ref. aeppli , detection of this algebraic tail appears very difficult, however high resolution measurements may be able to observe the additional contribution from the longitudinal spin fluctuations, which apparently behave like the critical fluctuations $`(18)`$, however with a rather large exponent $`\eta =1`$.

In summary, it has been shown that Goldstone modes associated with a broken continuous symmetry lead to fluctuations of the zero field order parameter which scale proportional to $`L^4`$ at any temperature below $`T_c`$. For 2D ordered ground states like the Néel phase of high temperature superconductors, the underlying longitudinal fluctuations lead to an additional critical continuum beyond the simple spin wave pole in the dynamic structure factor, whose detection requires high resolution neutron scattering experiments.

I would like to thank the theory department at the ETH in Zürich, where this work was started, for providing a very stimulating environment. Useful remarks by A. Aharony, R. Balian, D. Belitz, J. Fröhlich, M. Lüscher, C. Newman and S. Sachdev are gratefully acknowledged.