# Generalized Mean Field Approach to a Resonant Bose-Fermi Mixture

D.C.E. Bortolotti  JILA, NIST and Department of Physics, University of
Colorado, Boulder, CO 80309-0440  LENS and Dipartimento di Fisica, Universitá
di Firenze, and INFM, Sesto Fiorentino, Italy  A.V. Avdeenkov  Institute of
Nuclear Physics, Moscow State University, Vorob’evy Gory, 119992, Moscow,
Russia  J.L. Bohn  JILA, NIST and Department of Physics, University of
Colorado, Boulder, CO 80309-0440

###### Abstract

We formulate a generalized mean-field theory of a mixture of fermionic and
bosonic atoms, in which the fermion-boson interaction can be controlled by a
Feshbach resonance. The theory correctly accounts for molecular binding
energies of the molecules in the two-body limit, in contrast to the most
straightforward mean-field theory. Using this theory, we discuss the
equilibrium properties of fermionic molecules created from atom pairs in the
gas. We also address the formation of molecules when the magnetic field is
ramped across the resonance, and present a simple Landau-Zener result for this
process.

##  I introduction

The use of magnetic Feshbach resonances to manipulate the interactions in
ultracold quantum gases has greatly enriched the study of many-body physics.
Notable examples include the crossover between BCS (Bardeen-Cooper-Schrieffer
BCS1957 ) and BEC (Bose-Einstein Condensate bose1924 ; einstein1924 )
superfluidity in ultracold Fermi gases regal2003nature_sup ;
bartenstein2004cmb ; zwierlein2004cpf , and the “Bose Nova” collapse in Bose
gases donley2002nature . Recent experimental developments, olsen2004prl_rbk ;
zaccanti2006pra_rbk ; stan2004ofr ; ospelkaus:120402 ; jin2008 ; ketterle2008
have enabled the creation of an ultracold mixture of bosons and fermions,
where an interspecies Feshbach resonance may introduce a rich source of new
phenomena. From the theoretical point of view, studies of Bose-Fermi mixtures
to date have been mostly limited to non resonant physics, focusing mainly on
mean field effects in trapped systems roth2002mfi ; roth2002sas ;
modugno2002cdf ; hu2003ttb ; liu2003fte ; salasnich2007sbd ; pedri2005tdb ;
adhikari:043617 ; buchler2004psa ; kanamoto2006psf ; adhikari2005fbs , phases
in optical lattices albus2003mba ; lewenstein2004abf ; roth2004qpa ;
sanpera2004afb ; gunter2006bfm ; tit2008 ,or equilibrium studies of
homogeneous gases, focusing mainly on phonon induced superfluidity or beyond-
mean-field effects bijlsma2000ped ; heiselberg2000iii ; efremov2002pwc ;
viverit2002bis ; matera2003fpb ; albus2002qft ; wang2006sct . Pioneering
theoretical work on the resonant gas include Ref. powell2005prb , in which a
mean-field equilibrium study of the gas is supplemented with a beyond-mean-
field analysis of the bosonic depletion; and Ref. storozhenko2005pra , where
an equilibrium theory is developed using a separable-potentials model.

The aim of this article is to develop and solve a mean-field theory describing
an ultracold atomic Bose-Fermi mixture in the presence of an interspecies
Feshbach resonance. This goal appears innocuous enough at first sight, since
mean-field theories for resonant Bose-Bose Timmermans1999prl ;
holland2002prl_ramsey  and Fermi-Fermi milstein2002pra ; griffin2002prl ;
stoof2003joptb  gases exist, and have been studied extensively. In both of
these theories, the mean-field approximation consists of considering the
bosonic Feshbach molecules as being fully condensed, and this greatly
simplifies the treatment, since the Hamiltonian reduces to a standard
Bogoliubov-like integrable form landau9 .

The fundamental difference between these examples and the Bose-Fermi mixtures
is that in the latter the Feshbach molecules are fermions, and therefore their
center of mass-momentum must be included explicitly. The most obvious mean-
field approach consists in considering the atomic Bose gas to be fully
condensed. However, as we will show below, resonant molecules are really
composed of two bound atoms, which spend their time together vibrating around
their center of mass. It follows that outright omission of the bosonic
fluctuations of the atoms, disallows the bosonic constituents to oscillate
(i.e. fluctuate) at all, and therefore this leads to an improper description
of the physics of atom pairs.

This article is organized as follows. In section II we introduce the field
theory model used to study Feshbach resonances, describing briefly the
parametrization used, and outlining the exact solution of this model in the
two-body limit. The section ends with a test of this two-body theory, by
comparing the binding and resonance energies predicted by the model and the
virtually-exact analogues obtained from two-body close coupling calculations.
Section III introduces the simplest mean-field many-body theory of the gas,
obtained by disregarding all bosonic fluctuations. The solution of the theory
is outlined, and its limitations highlighted. In spite of these limitations,
mean-field theory provides a useful language for dealing with the problem, the
utility of which will persist even beyond the limits of applicability of the
theory itself.

Finally in section IV we introduce our generalized mean-field theory, which
is, in short, similar to the mean-field theory described in section III, but
with the notable improvement of using properly renormalized molecules as
building blocks, instead of their bare counterparts. This approach is not
trivially described in the Hamiltonian formalism, where substituting dressed
molecules for free ones would lead to double counting of diagrams. In this
section we therefore shift to the Green-function/path-integral language, where
this double-counting can be avoided quite easily. Finally we proceed to the
numerical solution of this theory, and note that for narrow resonances the
results are consistent with their mean-field equivalents. This encourages us
to develop a simple theory to study the molecular formation via magnetic field
ramps, and, using an approach based on the Landau-Zener formalism lz1 ; lz2 ,
we derive analytic expressions in section V.

Throughout this article, we work with zero temperature gases, in the free
space thermodynamic limit. These are limitations which render the results
obtained here hard to directly compare with experimental results. One of the
main possible future directions of this work should include solving the same
problem in a trap, and generalization to higher temperatures.

##  II The Model

We are interested primarily in the effects of resonant behavior on the
otherwise reasonably understood properties of the system. To this end we use a
model which has become standard in the last few years. This model has been
useful in studying the effect of resonant scattering in Bose
holland2002prl_ramsey ; stoof2004pr ; burnett2003pra  and Fermi holland2002pra
; stoof2004prl ; gurarie2004prl ; griffin2003pra ; sasha2004jpb ; sasha2005pra
gases. In the case of the bose-fermi mixture, this model has been used in
Ref.powell2005prb ,and bortolotti2006jpb ; bortolotti2006sfm . We refer to
these works for further details about the origin and justification of the
Hamiltonian we use here, and for details on the solution in the two body
regime.

The Hamiltonian for the system reads:

$H=H_{0}+H_{I},$ (1)

where

$\displaystyle H_{0}$ $\displaystyle=$ $\displaystyle\sum_{p}\epsilon_{p}^{F}\
\hat{a}_{p}^{\dagger}\hat{a}_{p}+\sum_{p}\epsilon_{p}^{B}\
\hat{b}_{p}^{\dagger}\hat{b}_{p}+\sum_{p}\left(\epsilon_{p}^{M}+\nu\right)\
\hat{c}_{p}^{\dagger}\hat{c}_{p}$ $\displaystyle+$ $\displaystyle{\gamma\over
2V}\sum_{p,p^{\prime},q}\hat{b}_{p-q}^{\dagger}\hat{b}_{p^{\prime}+q}^{\dagger}\hat{b}_{p}\hat{b}_{p^{\prime}}$
$\displaystyle H_{I}$ $\displaystyle=$ $\displaystyle{V_{bg}\over
V}\sum_{p,p^{\prime},q}\hat{a}_{p-q}^{\dagger}\hat{b}_{p^{\prime}+q}^{\dagger}\hat{a}_{p}\hat{b}_{p^{\prime}}$
$\displaystyle+$
$\displaystyle{g\over\sqrt{V}}\sum_{q,p}(\hat{c}_{q}^{\dagger}\hat{a}_{-p+q/2}\hat{b}_{p+q/2}+h.c.).$

Here $\hat{a}_{p},\hat{b}_{p}$, are the annihilator operators for,
respectively, fermions and bosons, $\hat{c}_{p}$ is the annihilator operator
for the molecular field holland2002pra ; sasha2004jpb ; sasha2005pra ;
$\gamma=4\pi a_{b}/m_{b}$ is the interaction term for bosons, where $a_{b}$ is
the boson-boson scattering length; $V_{bg},\nu$, and $g$ are parameters
related to the Bose Fermi interaction, yet to be determined; the single
particle energies are $\epsilon^{\alpha}=p^{2}/2m_{\alpha}$, where
$m_{\alpha}$ indicates the mass of bosons, fermions, or pairs; and $V$ is the
volume of a quantization box with periodic boundary conditions.

The first step is to find the values for $V_{bg},\nu,g$, in terms of
measurable parameters. We will, for this purpose, calculate the 2-body
T-matrix resulting from the Hamiltonian in Eq. LABEL:act-bfm. Integrating the
molecular field out of the real time path integral, leads to the following
Bose-Fermi interaction Hamiltonian:

$H_{I}^{2body}={1\over V}\left(V_{bg}+{g^{2}\over
E-\nu}\right)\sum_{p}\hat{a}_{p}^{\dagger}\hat{b}_{-p}^{\dagger}\hat{a}_{p}\hat{b}_{-p}.$
(3)

This expression is represented in center of mass coordinates, and $E$ is the
collision energy of the system. From the above equation we read trivially the
zero energy scattering amplitude in the saddle point approximation:

$T=(V_{bg}-{g^{2}\over\nu}),$ (4)

which corresponds to the Born approximation (This is akin to identifying the
scattering amplitude $f=a_{sc}$ in the Gross-Pitaevskii equation, where the
interaction term would be ${2\pi\over m_{bf}}a_{b}$ ). We emphasize that this
approximation is only valid in the zero energy limit, and it does not,
therefore, describe the correct binding energy as a function of detuning.
However, with this approach we obtain an adequate description of the behavior
of scattering length as a function of detuning, which allows us to relate the
parameters of our theory to experimental observables via the conventional
parametrization holland2002pra ; stoof2004pr

$T(B)={2\pi\over m_{bf}}a_{bg}\left(1-{\Delta_{B}\over(B-B_{0})}\right),$ (5)

where $a_{bg}$ is the value of the scattering length far from resonance,
$\Delta_{B}$ is the width, in magnetic field, of the resonance, $m_{bf}$ is
the reduced mass, and $B_{0}$ is the field at which the resonance is centered.

The identification of parameters between Eqns., (4) and (5) proceeds as
follows: far from resonance, $|B-B_{0}|>>\Delta_{B}$, the interaction is
defined by a background scattering length, via $V_{bg}={2\pi a_{bg}\over
m_{bf}}$. To relate magnetic field dependent quantity $B-B_{0}$ to its energy
dependent analog $\nu$, requires defining a parameter
$\delta_{B}=\partial\nu/\partial B$, which may be thought of as a kind of
magnetic moment for the molecules. It is worth noting that $\nu$ does not
represent the position of the resonance nor the binding energy of the
molecules, and that, in general $\delta_{B}$ is a field-dependent quantity,
since the thresholds move quadratically with field, because of nonlinear
corrections to the Zeeman effect. For current purposes we identify
$\delta_{B}$ by its behavior far from resonance, where it is approximately
constant. Careful calculations of scattering properties using the model in Eq.
(LABEL:act-bfm), however, leads to the correct Breit-Wigner behavior of the
2-body T-matrix newton .

Finally we get the following identifications:

$\displaystyle V_{bg}={2\pi a_{bg}\over m_{bf}}$ $\displaystyle
g=\sqrt{V_{bg}\delta_{B}\Delta_{B}}$ $\displaystyle\nu=\delta_{B}(B-B_{0}).$

Synthesizing the approach descibed in bortolotti2006jpb , and diagrammatically
represented in Fig. 1, we can obtain the exact two-body T-matrix of the system
by solving the Dyson equation

$\displaystyle T$ $\displaystyle=$ $\displaystyle gD^{0}g+gD^{0}g\ \Pi\
gD^{0}g+gD^{0}g\ \Pi\ gD^{0}g\ \Pi\ gD^{0}g+\cdots$ (7) $\displaystyle=$
$\displaystyle gDg$

where $T$ is the T-matrix for the collision, and which has formal solution

$T=gDg={g^{2}\over(D^{0})^{-1}-g^{2}\Pi}.$ (8)

These quantities take the explicit form

$\displaystyle D^{0}(E)$ $\displaystyle=\left({V_{bg}\over g^{2}}+{1\over
E-\nu}\right)$ (10) $\displaystyle\Pi(E)$ $\displaystyle=-i\int\ {d\omega\over
2\pi}{d{\bf p}\over(2\pi)^{3}}{1\over(\hbar\omega-{p^{2}\over
2m_{b}}+i0^{+})}\times$ $\displaystyle{1\over(E-\hbar\omega-{p^{2}\over
2m_{f}}+i0^{+})}$ $\displaystyle\approx
i{m_{bf}^{3/2}\over\sqrt{2}\pi}\sqrt{E}+{m_{bf}\Lambda\over\pi^{2}},$

where $m_{bf}$ is the boson-fermion reduced mass, and $\Lambda$ is an
ultraviolet momentum cutoff needed to hide the unphysical nature of the
contact interactions. Note that $D^{0}$ represents an effective molecular
field, accounting for the fermion-boson background interaction, and, as
described in detail in bortolotti2006jpb , it is obtained by integrating the
original molecular field, and performing a Hubbard-Stratonovich transformation
negele  to eliminate the direct boson-fermion interaction in favour of the
effective molecular field $D^{0}$.

Figure 1: Feynman diagrams representing the resonant collision of a fermion
and a boson. Solid lines represent fermions, dashed lines bosons, and double
solid-dashed lines represent the effective composite fermions.

Regularization of the theory bortolotti2006jpb ; milstein2002pra  is obtained
by the substitutions

$\displaystyle\bar{V_{bg}}=V_{bg}\left(1\over 1-{m_{bf}\Lambda
V_{bg}\over\pi^{2}}\right)$ $\displaystyle\bar{g}=g\left(1\over
1-{m_{bf}\Lambda V_{bg}\over\pi^{2}}\right)$
$\displaystyle\bar{\nu}=\nu+\bar{g}g{m_{bf}\Lambda V_{bg}\over\pi^{2}}.$ (11)

Finally the two-body T-matrix takes the form

$T(E)=\left[{1\over V_{bg}+{g^{2}\over
E-\nu}}+i{m_{bf}^{3/2}\over\sqrt{2}\pi}\sqrt{E}\right]^{-1}.$ (12)

###  II.1 Poles of the T-Matrix: Testing the Model

Bound states and resonances of the two-body system are identified in the
structure of poles of the T-matrix (Eq.12). This is illustrated in Figure 2,
where real and imaginary parts of the poles’ energies are plotted as a
function of magnetic field. The resonance portrayed in the figure is the
$544.7G$ resonance present in the $|9/2,-9/2\rangle|1,1\rangle$ state of
40K-87Rb. For $B<544.7G$ (corresponding to detunings $\nu<0$) the two-body
system possesses a true bound state, whose binding energy is denoted by the
solid line. In this case, the pole occurs for real energies. This bound state
vanishes as the detuning goes to zero, where the resonance occurs.

For positive detunings, $\nu>0$, on the other hand, the poles are complex, and
the inverse of the imaginary part is proportional to the lifetime of the
metastable resonant state.

Figure 2: The top panel shows the scattering length versus magnetic field for
the 544.7 G resonance, present in the $|9/2,-9/2\rangle|1,1\rangle$ states of
the 40K-87Rb collision. The bottom panel shows the poles of the model two-body
T-matrix (Eq. 12) parameterized for the same resonance, as a function of
magnetic field. Thick solid and dashed lines denote the real parts of relevant
poles, representing bound and resonance states respectively. The thin dashed
lines are the real parts of unphysical poles. The empty circles and squares,
represent the position of the resonance and the bound state, obtained via a
virtually exact close coupling calculation, and are presented to show the
level of accuracy of the model.

In this regime, there is no longer a true bound state, but there may be a
scattering resonance, indicated in Fig. 2 by a thick dashed line. This
resonance appears for magnetic fields $B>544.7$ for this particular resonance,
well before the disappearance of the bound state. This value is highly
dependent upon the value of the background potential. We will see in section
IV that for $V_{bg}=0$, the resonance actually appears at positive detunings.
In the case of 40K-87Rb, $V_{bg}<0$, implying that there is a weak potential
resonance in the open channel which interferes with the closed-channel
resonance, and causes it to cross the axis at negative detunings. For
$V_{BG}>0$ (unpublished ) the positive background scattering length is set by
a bound state in the open channel, which does not affect the resonance states,
but which interferes with the bound state at negative detunings.

The thin dashed lines in Fig. 2 are physically meaningless solutions to the
Schrödinger equation, in which the amplitude in the resonant state would grow
exponentially in time, rather than decay. These poles do not therefore
identify any particular features in the energy-dependent cross section of the
atoms, and will not modify the physics of the system. Finally Fig. 2 contains
data obtained from virtually exact close-coupling calculations which show the
extent of validity of the model. For the purposes at hand this agreement is
sufficient.

It should be noted that the agreement is not as good for positive background
scattering length systems, since the open-channel bound state determining this
scattering length is not adequately described by the model, which treats the
background physics as an essentially zero range interaction. This implies that
the relation between the background scattering length and open channel bound-
state energy is exactly $E_{b}=1/2\mu a_{bg}^{2}$, while in the physical
system this relation depends on the details of the interaction potential. This
problem has been addressed in the literature kokkelmans2004pra , but no
treatable field theory has yet been proposed.

$\hskip 14.22636ptB_{0}(G)\hskip 14.22636pt$ | $\hskip 14.22636pt\Delta_{B}(G)$ |  $\delta_{B}(K/G)$ |  $a_{bg}(a.u.)$
---|---|---|---
$492.49$ | $0.134$ | $3.624\times 10^{-5}$ | $-176.5$
$544.7$ | $3.13$ | $1.576\times 10^{-4}$ | $-176.5$
$659.2$ | $1.0$ | $2.017\times 10^{-4}$ | $-176.5$
Table 1: Parametrization of the three main Feshbach resonances used in this
thesis. All Three resonances are in the $|9/2,-9/2\rangle|1,1\rangle$ states
of the ${}^{40}K-^{87}Rb$ collision

##  III Mean-Field Theory

In this section we introduce the many-body physics of the system, by first
analyzing it in a mean-field approach. Because of the statistical properties
of the system, we will see right away that mean field theory does not recover
the correct two body physics in the low density limits. In spite of this
substantial weakness, however, the approach has several qualitative features
which persist even in the improved theory that we introduce below.
Furthermore, since the model is exactly solvable, it will allow us to develop
a language which will help us to understand the problem in simpler terms, and
to identify some small physical effects, which, when ignored, can greatly
simplify the beyond mean-field approach presented in the next section.

###  III.1 The Formalism

Starting with the Hamiltonian described by equation (1), we obtain the mean-
field Hamiltonian by substituting the boson annihilator $\hat{b}$ by its
expectation value $\phi=\langle\hat{b}\rangle$, a complex number. The number
operator $\hat{b}_{p}^{\dagger}\hat{b}_{p}$ therefore becomes
$|\phi|^{2}=N_{b}$, where $N_{b}$ is the number of condensed bosons. The grand
canonical Hamiltonian therefore becomes

$\displaystyle H$ $\displaystyle=$ $\displaystyle
E_{b}+\sum_{p}\left(\epsilon_{p}^{F}-\mu_{f}+V_{bg}n_{b}\right)\hat{a}_{p}^{\dagger}\hat{a}_{p}$
$\displaystyle+$
$\displaystyle\sum_{p}\left(\epsilon_{p}^{M}+\nu-\mu_{m}\right)\
\hat{c}_{p}^{\dagger}\hat{c}_{p}$ $\displaystyle+$ $\displaystyle
g\sqrt{n_{b}}\sum_{p}\left(\hat{c}_{p}^{\dagger}\hat{a}_{p}+h.c.\right),$

where $n_{b}$ is the density of condensed bosons, $E_{b}/V=\gamma
n_{b}^{2}-\mu_{b}n_{b}$ is the energy per unit volume of the (free) condensed
bosons, a constant contribution to the total energy of the system, and
$\mu_{(b,f,m)}$ are the chemical potentials. These are Lagrange multipliers
that serve to keep the densities constant as we minimize the energy to find
the ground state. In the following we will drop the volume term, absorbing it
in the definition of the creator/annihilator operators, such that the expected
value of the number operator represents a density, instead of a number.

Before proceeding with the analysis of this Hamiltonian, we introduce the set
of self-consistent equations we wish to solve. To this end we define the
quantities $n_{b(f)}^{0}$, representing the total density of bosons (fermions)
in the system, at detuning $\nu\rightarrow\infty$. At finite detunings some of
these atoms will combine into molecules, and the densities will be denoted as
$n_{(b,f,m)}$ for bosons, fermions and molecules respectively.

The system, therefore, is described by six quantities, namely three densities
and three chemical potentials, which require six equations to determine. These
equations, which can be derived by number-conservation constraints and energy
minimization arguments, are:

$\displaystyle n_{f}+n_{m}-n_{f}^{0}=0$ (14.a) $\displaystyle
n_{b}+n_{m}-n_{b}^{0}=0$ (14.b) $\displaystyle n_{f}={d\ \Omega\over
d\mu_{f}}$ (14.c) $\displaystyle n_{m}={d\ \Omega\over d\mu_{m}}$ (14.d)
$\displaystyle{d\Omega\over d\phi}=0$ (14.e)
$\displaystyle\mu_{b}+\mu_{f}=\mu_{m},$ (14.f)

where $\Omega=\langle H\rangle/V$ is the Gibbs free energy. Equations (14.a)
and (14.b) follow from the simple counting argument that for every molecule
created, there is one less free boson and one less free fermion in the gas.
Equations 14.c, and 14.d are simply the Lagrange multiplier constraint
equations, equation (14.e) follows from the mean-field approximation, whereby
the bosonic field is simply a complex number, and minimization of the energy
can therefore be done directly. Finally equation (14.f) is the law of mass
action, which follows from the fact that to make a molecule it takes one free
atom of each kind.

The next step is to write down $\Omega$ for the system, by taking the
expectation value of the Hamiltonian in equation( III.1), obtaining

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle
E_{b}/V+\sum_{p}\left(\epsilon_{p}^{F}-\mu_{f}+V_{bg}n_{b}\right)\eta_{f}(p)$
$\displaystyle+$
$\displaystyle\sum_{p}\left(\epsilon_{p}^{M}+\nu-\mu_{m}\right)\eta_{m}(p)$
$\displaystyle+$ $\displaystyle 2g\sqrt{n_{b}}\sum_{p}\eta_{mf}(p),$

where $\eta_{f}(p)=\langle\hat{a}_{p}^{\dagger}\hat{a}_{p}\rangle$ and
$\eta_{m}(p)=\langle\hat{c}_{p}^{\dagger}\hat{c}_{p}\rangle$ are the fermionic
and molecular momentum distributions,
$\eta_{mf}(p)=\langle\hat{c}_{p}^{\dagger}\hat{a}_{p}\rangle$ is an off-
diagonal correlation term arising from the interactions in the system, and the
densities are given by $n_{f,m,mf}=\int{dp\over 2\pi^{2}}\eta_{f,m,mf}(p)$.
Equations 14.a-14.f then read

$\displaystyle n_{f}+n_{m}-n_{f}^{0}=0$ (16.a) $\displaystyle
n_{b}+n_{m}-n_{b}^{0}=0$ (16.b) $\displaystyle n_{f,m,mf}=\int{dp\over
2\pi^{2}}\eta_{f,m,mf}(p)$ (16.c) $\displaystyle g\
n_{mf}-\mu_{b}\sqrt{n_{b}}+\gamma n_{b}^{3/2}=0$ (16.d)
$\displaystyle\mu_{b}+\mu_{f}=\mu_{m}$ (16.e)

The remaining task is now to find expressions to calculate the expected values
$\eta_{f,m,mf}(p)$. To this end we follow a Bogoliubov-like approach, similar
to that described in powell2005prb . The mean-field Hamiltonian is bilinear in
all creation/annihilation operators, which means that it can be diagonalized
via a change of basis, whereby introducing the operators

$\displaystyle\hat{\alpha}_{p}=A_{\alpha}\hat{a}_{p}+C_{\alpha}\hat{c}_{p}$
(17) $\displaystyle\hat{\beta}_{p}=A_{\beta}\hat{a}_{p}+C_{\beta}\hat{c}_{p},$

for some appropriately chosen coefficients $A_{\alpha,\beta}$ and
$B_{\alpha,\beta}$, the Hamiltonian will read

$H^{\prime}=E_{0}+\sum_{p}\lambda_{\alpha}(p)\hat{\alpha}_{p}^{\dagger}\hat{\alpha}_{p}+\sum_{p}\lambda_{\beta}(p)\hat{\beta}_{p}^{\dagger}\hat{\beta}_{p}.$
(18)

At this point we note that the Hamiltonian is just a separable sum of free-
particle Hamiltonians, where the free particles are fermions, with dispersion
relations $\lambda_{\alpha,\beta}(p)$. We can readily write down the
distribution

$\eta_{\alpha,\beta}(p)=\Theta(-\lambda_{\beta,\alpha}(p))$ (19)

where $\Theta$ is the step function, and calculate the densities
$n_{\alpha,\beta}$. The step function could be replaced by the free Fermi
distribution for non-zero temperatures, but the mean-field assumption that all
bosons are condensed would no longer hold. If these were ordinary free
fermions with dispersion $p^{2}/2m-\mu$, equation (19) would reduce to the
standard zero-temperature Fermi distribution. We will see below that
$\lambda_{\alpha,\beta}(p)$ are dispersion relations of quasi-particles that
are a mixture of atoms and molecules.

Below we show how these ideas, together with equation (14.a) - 14.f, give us
the tools we require to calculate the observable atomic and molecular
densities as a function of the chemical potentials.

To illustrate more explicitly the diagonalization procedure we define the
vectors

$A=\left(\begin{array}[]{r}\hat{a}_{p}\\\ \hat{c}_{p}\end{array}\right)\ \ \
A^{\dagger}=(\ \hat{a}_{p}^{\dagger}\ \ \hat{c}_{p}^{\dagger}\ ),$ (20)

and

$B=\left(\begin{array}[]{r}\hat{\alpha}_{p}\\\
\hat{\beta}_{p}\end{array}\right)\ \ \ B^{\dagger}=(\
\hat{\alpha}_{p}^{\dagger}\ \ \hat{\beta}_{p}^{\dagger}\ ),$ (21)

whereby the Hamiltonian can be written as $A^{\dagger}\hat{H}A$, and
$B^{\dagger}\hat{H^{\prime}}B$, where

$\hat{H}=\left(\begin{array}[]{cc}\left(\epsilon_{p}^{F}-\mu_{f}+V_{bg}n_{b}\right)&g\sqrt{n_{b}}\\\
g\sqrt{n_{b}}&\left(\epsilon_{p}^{M}+\nu-\mu_{m}\right)\end{array}\right),$
(22)

and

$\hat{H^{\prime}}=\left(\begin{array}[]{cc}\lambda_{\alpha}(p)&0\\\
0&\lambda_{\beta}(p)\end{array}\right).$ (23)

Diagonalizing $\hat{H}$, we get the two eigenvalues

$\displaystyle\lambda_{\alpha,\beta}(p)$ $\displaystyle=$
$\displaystyle{h_{f}(p)+h_{m}(p)\over 2}$ $\displaystyle\pm$
$\displaystyle{1\over 2}\sqrt{4g^{2}n_{b}+(h_{m}(p)-h_{f}(p))^{2}},$

where we have defined
$h_{f}(p)=\left(\epsilon_{p}^{F}-\mu_{f}+V_{bg}n_{b}\right)$, and
$h_{m}(p)=\left(\epsilon_{p}^{M}+\nu-\mu_{m}\right)$, and the unitary
eigenvector matrix

$U=\left(\begin{array}[]{cc}A_{\alpha}&B_{\alpha}\\\
A_{\beta}&B_{\beta}\end{array}\right).$ (25)

The transformation in eq. 17 can then be written as $A=U^{\dagger}B$, and its
inverse $B=UA$.

Our goal now is to write the densities $\eta_{m,f,mf}(p)$ in terms of the
known densities $\eta_{\alpha,\beta}(p)$. In component notation, (where
$A_{i}=\hat{a}_{p}$, etc.), we can write

$\langle A_{i}^{\dagger}A_{l}\rangle=\langle
B_{j}^{\dagger}U_{ji}(U^{\dagger})_{lj}B_{k}\rangle=U_{lj}U^{*}_{ij}\langle
B_{j}^{\dagger}B_{j}\rangle,$ (26)

where we have used the fact that since the Hamiltonian is diagonal in the $B$
basis, then $\langle B_{j}^{\dagger}B_{k}\rangle=\langle
B_{j}^{\dagger}B_{j}\rangle\delta_{jk}$.

Using this formalism we obtain the relations:

$\displaystyle\eta_{f}(p)$ $\displaystyle=$
$\displaystyle|A_{\alpha}|^{2}\eta_{\alpha}(p)+|B_{\alpha}|^{2}\eta_{\beta}(p)$
(27) $\displaystyle\eta_{f}(p)$ $\displaystyle=$
$\displaystyle|A_{\beta}|^{2}\eta_{\alpha}(p)+|B_{\beta}|^{2}\eta_{\beta}(p)$
$\displaystyle\eta_{fm}(p)$ $\displaystyle=$ $\displaystyle
A_{\alpha}^{*}A_{\beta}\eta_{\alpha}(p)+B_{\alpha}^{*}B_{\beta}\eta_{\beta}(p).$

Using these expressions in conjunction with eqs. 16.a-16.e will then allow us
to compute the equilibrium properties of the system.

This simplified, mean-field version of the solution can only approximately
reproduce the energies of atomic and molecular states, as is shown in figure
3. In this example, we have assumed a uniform mixture of 40K and 87Rb atoms
with densities $8.2\times 10^{14}$ cm-3 and $4.9\times 10^{15}$ cm-1,
respectively. These densities correspond to the central density of each
species assuming it is confined to a $100$ Hz spherical trap. Far from the
resonance, these energies asymptote to zero (representing the atomic state)
and to the detuning (representing the bare molecular state that went into the
theory). Near zero detuning, these levels cross owing to the coupling term
$g\sqrt{n_{b}}$ in Eqn. (22). The size of this crossing is therefore larger
the larger the bosonic density is. In the crossing region the eigenstates do
not clearly represent either atoms or molecules, but linear combinations of
the two.

Figure 3:  The thick solid lines represent the “renormalized” mean-field
energy levels $\lambda_{\alpha,\beta}(k_{f})$ , while the thin dashed lines
represent their bare counterparts.

###  III.2 Mean field: noninteracting case

To better understand the structure of the mean-field theory, in this section
we detail its results for a noninteracting gas, by setting
$g=V_{bg}=\gamma=0$. We contrast two different physical regimes, based on the
ratio of bosons to that of fermions, $r_{bf}\equiv n_{b}/n_{f}$.

Figure 4: Equilibrium chemical potentials (top panel) and populations (bottom
panel) as a function of detuning for a non interacting gas with $r_{bf}=.6$.
The solid lines represent fermions, dashed lines molecules, and dashed-dotted
lines bosons. The dotted lines in the top panel represent the bare molecular
and fermionic internal energies, respectively $\nu$ and $0$. The vertical
lines labeled a) and b) are discussed in the text, and represent the detuning
at which molecular formation begins and ends, respectively

Figure 5:  Equilibrium chemical potentials (top panel) and populations (bottom
panel) as a function of detuning for a non interacting gas with $r_{bf}=6$.
The solid lines represent fermions, dashed lines molecules, and dashed-dotted
lines bosons. The dotted lines in the top panel represent the bare molecular
and fermionic internal energies, respectively $\nu$ and $0$.

In the case of high fermion density, we set $r_{bf}=0.6$, and plot chemical
potentials and populations of the various states in Figure 4. Consider what
happens in an infinitely slow ramp from positive detuning (no molecules) to
negative detuning (introducing bound molecular states). For large positive
detuing, the fermionic and molecular chemical potentials are the same, since
the chemical potential of the condensed bosons vanishes. For large enough
detuning, the molecular chemical potential remains below the detuning, so it
is energetically unfavorable to make molecules. When the detuning dips below
the chemical potential (detuning a in the figure), fermions begin to pair with
bosons to make molecules (see populations in lower panel). This process
continues until all the bosons are consumed (detuning b), at which point the
populations stabilize. For detunings less that this, there remain both
fermions and fermionic molecules in the gas, and there are two Fermi surfaces
present. Because the internal energy of the molecules continues to diminish at
lower detuning, so does the molecular chemical potential. The two Fermi
surfaces therefore split from one another, although the relative population of
the two fermions is fixed.

By contrast, the case where the bosons outnumber the fermions is shown in
Figure 5, where we have set the density ratio to $r_{bf}=6$. As in the
previous case, no molecules are generated until the detuning drops lower than
the chemical potential of the atomic Fermi gas. Since there are enough bosons
to turn all the fermions into molecules, there there are no fermionic atoms at
sufficiently negative detuning, and the gas possesses only a single Fermi
surface. The chemical potential of the remaining bosons is still zero, since
these bosons are condensed. Formally, then, the chemical potential for atomic
fermions is negative, meaning that their formation at negative detuning is
energetically forbidden.

###  III.3 Mean field: interacting case

In most experimental circumstances, the density of bosons is larger than that
of fermions, since condensed bosons cluster to the center of the trap, whereas
fermions are kept away by Pauli blocking. We therefore focus on this case
hereafter, setting the bose and fermi densities to $n_{b}=4.9\times 10^{15}$
cm-3 and $n_{f}=8.2\times 10^{14}$ cm-3. The coupling term $g\sqrt{n_{b}}$ in
equation (22) is the perturbative expansion parameter for the problem, and
since it has units of energy, it must be compared with the characteristic non-
perturbed energy of the gas, which in this case is $E_{f}$. Also since in the
perturbative expansions it always appears squared (see eq. 7), we can define
the unitless small parameter for the system as
$\epsilon_{SM}=g^{2}n_{b}/E_{f}^{2}=g^{2}n_{f}/E_{f}^{2}r_{bf}$, For the
$492G$ resonance in table 1, we have $\epsilon_{SM}=6.35\times 10^{-2}r_{bf}$,
and since $r_{bf}=6$, the small parameter is of order $0.1$, appropriate for
perturbative treatment.

Figure 6 shows the equilibrium chemical potentials for the system obtained via
a self-consistent solution of equations (III.1) and 16.a-16.e. This figure is
qualitatively similar to the corresponding non-interacting result in Fig. 5,
but contains important differences. First, the nonzero boson-boson interaction
generates a nonzero bosonic chemical potential $\mu_{b}$ that breaks the
degeneracy between the molecular and fermionic chemical potentials. In
physical terms, this means that there is an energy cost in maintaining bosons
unpaired, and therefore we need to take this into account in the kinematic
analysis. Namely, to make molecules energetically favorable no longer requires
a detuning $\nu$ such that $\nu=\mu_{f}$, but now requires
$\nu=\mu_{f}+\mu_{b}$. The net result is to shift the chemical potential up
and to the left by an amount $\mu_{b}$, and to shift the molecular population
curve to the left by this amount.

Figure 6:  Mean-field equilibrium chemical potentials (top panel) and
populations (bottom panel) as a function of detuning. The solid lines
represent fermions, dashed lines molecules, and dashed-dotted lines bosons.

A second difference is that molecule creation takes place more gradually as a
function of detuning in the interacting case. This is simply the result of the
avoided crossing smearing out the molecular energy.

Finally, we exploit the simplicity of the mean field approach to test some
approximations that will simplify the beyond-mean-field approach in the next
section. These approximations have been tested numerically, and they give
corrections of the order of $.1\%$ or less in calculated molecular populations
for all regimes of interest here. The approximations are: i) incorporate the
boson-boson interaction $\gamma N_{b}^{3}$ by shifting the detuning and
chemical potential as discussed above; ii) disregard the background scattering
between bosons and fermions, i.e., set $V_{bg}=0$, since this interaction is
dominated by its resonance part; and iii) disregard the correlation function
$\langle\eta_{b}f(p)\rangle$ ((analogous to the boson polarization operator in
the Green function formalism), since its contribution to $\mu_{b}$ is much
smaller than that of $\gamma n_{b}^{3}$. It is difficult to directly verify
the validity of these approximations in the beyond-mean-field approach.
Nevertheless, we expect that these approximations remain valid, since the
generalized mean field theory is, after all, a mean-field theory at heart.

##  IV Generalized Mean-Field Theory

In section III we reached the conclusion that the mean-field approach to the
resonant Bose-Fermi system does not properly account for the correct two body
physics of the system. In this section we wish to improve on this, by
introducing a generalization to mean-field theory, via an appropriate
renormalization of the molecular propagator, which is able to reproduce the
correct two-body physics in the low-density limit. To accomplish this, we will
have to abandon the Hamiltonian treatment of the previous section, in favor of
a perturbative approach based on the Green’s function formalism, much as was
done for two bodies in II. Throughout this section we use the approximations
made above, namely, $\gamma n_{b}^{3}=V_{bg}=\langle\eta_{b}f(p)\rangle=0$.

We begin by recasting the rn field result from Sec. III in the language of
Green functions. The self-consistent Dyson equations that describe this system
are

$\displaystyle G^{MF}_{F}(E,P)$ $\displaystyle=$ $\displaystyle
G^{0}_{F}(E,P)+g^{2}n_{b}\ D^{0}(E,P)\ G^{MF}_{F}(E,P)$ $\displaystyle
D^{MF}(E,P)$ $\displaystyle=$ $\displaystyle D^{0}(E,P)+g^{2}n_{b}\
G^{0}_{F}(E,P)\ D^{MF}(E,P),$ (28)

where the free propagators are simply

$\displaystyle D^{0}={1\over\omega-\xi^{M}(p)+i\eta\ {\rm sign}(\xi(p))}$
$\displaystyle G^{0}_{F}={1\over\omega-\xi^{F}(p)+i\eta\ {\rm sign}(\xi(p))},$
(29)

and $\xi^{M,F}(p)=(\epsilon_{p}^{M,F}-\mu_{m,f})$ These propagators are
described diagrammatically in Fig. 7. They represent the fact that a free
fermion may encounter a condensed boson and associate with it, temporarily
creating a molecule; or that a free molecule may temporarily split into a
fermion and a condensed boson. Self-consistency ensures that these processes
may be repeated coherently an infinite number of times. We neglect the bosonic
renormalization equation $\phi^{MF}=n_{b}^{0}+g^{2}n_{b}\ G^{0}_{F}(E,P)\
D^{0}(E,P)$, whereby a condensed boson may pick-up a fermion to create a
molecule; this is equivalent to the condition $\langle\eta_{b}f(p)\rangle=0$.

Figure 7:  Feynman diagrams included in the mean-field theory. Thin (thick)
solid lines represent free (renormalized) fermions, thin double dashed-solid
lines represent free molecules, and thick double dashed lines represent
renormalized molecules. The little lightning bolts represent condensed bosons,
whereby the arrow indicates whether they are taken from or released into the
condensate.

Solutions to these equations take the form

$\displaystyle G^{MF}_{F}(E,P)$ $\displaystyle=$ $\displaystyle{1\over
G^{0}_{F}(E,P)^{-1}-g^{2}n_{b}\ D^{0}(E,P)}$ $\displaystyle D^{MF}(E,P)$
$\displaystyle=$ $\displaystyle{1\over D^{0}(E,P)^{-1}-g^{2}n_{b}\
G^{0}(E,P)}.$ (30)

Using the definitions of $G^{0}_{F}(E,P)$ and $D^{0}(E,P)$ from Eq. 29, we can
find the poles corresponding to many-body bound states. In this case, it can
be shown that the poles are exactly the mean-field eigenvalues
$\lambda_{\alpha,\beta}(p)$ from section II. Moreover, the equations (IV) are
symmetric with respect to interchange of $G_{F}$ and $D$, which implies that
both renormalized green functions have the same poles, and the same residues.
We can therefore study the properties of the fermions by only looking at the
molecules. This is not completely surprising, since, given that the condensed
bosons are relatively inert, every molecule corresponds exactly to a missing
fermion, and vice-versa.

The most important deficiency of this mean-field approach is that it only
allows molecules to decay into a free fermion and a condensed boson,
disregarding the possibility that the bosonic byproduct may be noncondensed.
We must allow noncondensed bosons somehow, and yet these bosons make a
perturbation to the result, as seen by the following argument. The fundamental
mean field assumption is that the gas is at zero temperature, and therefore
the noncondensed population should be negligible at equilibrium. Furthermore,
if a molecule is composed of a zero-momentum boson and a fermion from the
Fermi sea, dissociating into a noncondensed boson implies that the outgoing
fermion would have momentum lower than the Fermi momentum, an event which
Pauli blocking makes quite unlikely. Therefore, if a molecule does indeed
decay yielding a non-condensed boson, it should immediately recapture the
boson in a virtual process such as that described in Fig. 1. It is only
convenient that these events are exactly the kind of events which will
correctly renormalize the binding energy of the molecules, leading to a theory
which will reproduce the exact two-body resonant physics.

The Dyson equation describing this generalized mean-field theory are:

$\displaystyle G^{GMF}_{F}(E,P)$ $\displaystyle=$ $\displaystyle
G^{0}_{F}(E,P)+g^{2}n_{b}\ D(E,P)\ G^{GMF}_{F}(E,P)$ $\displaystyle
D^{GMF}(E,P)$ $\displaystyle=$ $\displaystyle D(E,P)+g^{2}n_{b}\
G^{0}_{F}(E,P)\ D^{GMF}(E,P).$

Here we have replaced the free propagator $D^{0}$ by the renormalized
molecular propagator $D$ from equation (8). A diagrammatic representation of
this theory appears in Fig. 8. By analogy with the mean field version, the
solution to these equations are:

$\displaystyle G^{GMF}_{F}(E,P)$ $\displaystyle=$ $\displaystyle{1\over
G^{0}_{F}(E,P)^{-1}-g^{2}n_{b}\ D(E,P)}$ $\displaystyle D^{GMF}(E,P)$
$\displaystyle=$ $\displaystyle{1\over D(E,P)^{-1}-g^{2}n_{b}\ G^{0}(E,P)}.$
(32)

These equations preserve the symmetrical nature of the mean-field theory
described above, and also the avoided crossing of atomic and molecular levels.
This is demonstrated in Figure 9, where we reproduce the $P=k_{f}$ pole of $D$
from Fig. 2 as dashed lines. We also present in this figure the corresponding
poles for the generalized mean-field theory, as solid lines. As in figure 3 we
note the splitting in two energy levels, avoiding each-other around $\nu=0$.
The fundamental difference in this case is that the molecular curve does not
asymptote to the bare detuning, but rather to the correct molecular binding
and resonance energies.

Figure 8:  Feynman diagrams included in the generalized mean-field theory.
Like in the mean field case (Fig. 7), thin (thick) solid lines represent free
(renormalized) fermions, thin double dashed-solid lines represent free
molecules, and thick double dashed lines represent renormalized molecules. The
little lightning bolts represent condensed bosons, whereby the arrow indicates
whether they are taken from or released into the condensate. The novelty here
in the inclusion of the 2-body dressed molecules from Fig. 1.

Figure 9: The thick lines represent the poles of $D^{GMF}(k_{f})$, and
correspond to molecular and fermionic energies. The dashed lines represent the
molecular two body poles, corresponding to the molecular binding and resonant
energies obtained disregarding the background interaction bortolotti2006sfm .
As in the mean-field theory case, the effect of the interaction with the
condensate is to create an avoided crossing between the atomic and molecular
states. The “bulge” in the upper solid curve is a consequence of Pauli
blocking which has the consequence of favouring molecular stability. See
bortolotti2006sfm  for more details.

Studying the equilibrium properties of the system is now a matter solving the
self-consistent set of equations (16.a-16.e), while setting $\eta_{mf}$ and
$\lambda$ equal to zero. To do this we first need to extract the distributions
$\eta_{f,m}$ from the Green functions $D^{GMF}$ and $G^{GMF}_{F}$. To avoid
taking a distracting detour here, we refer the reader to appendix A for
details. As in the previous section, we will consider a mixture composed of a
free gas of fermionic ${}^{40}K$ atoms, with a density of $8.2\times 10^{14}$
cm-3, and a gas of condensed ${}^{87}Rb$ bosons with density $4.9\times
10^{15}$ cm-3 (corresponding to the respective Thomas-Fermi densities of
$10^{6}$ atoms of either species in the center of a $100Hz$ spherical trap).

Figure 10 shows the equilibrium molecular population as a function of
detuning, for the $492.5G$ resonance. For the densities assumed, the mean-
field parameter $\epsilon_{SM}=g^{2}n_{b}/E_{f}^{2}\approx 0.4$ is indeed
perturbative. For this narrow resonance, the agreement between mean-field and
generalized mean-field is quite good, and we could as easily have used the
bare molecular positions to calculate this quantity. However, the situation is
completely different for the wide resonance at 544.7G, for which the
equilibrium molecular populations are shown in Fig. 11. Here the
“perturbative” parameter has the value $\epsilon_{SM}=38.7$, and is not
perturbative at all. For a given small detuning, the simple mean-field
approximation would greatly overestimate the number of molecules in the gas at
equilibrium. The more realistic generalized mean field theory accounts for the
fact that the actual molecular bound-state energy is higher than the bare
detuning. This fact in turn hinders molecular formation, according to chemical
potential arguments analogous to those in Sec. III.2.

Figure 10:  Equilibrium molecular population as a function of detuning for the
narrow $492.49G$ resonance. The solid line represents results obtained via the
generalized mean field theory presented in the text, while the dashed-dotted
line represents the mean field results.

Figure 11:  Equilibrium molecular population as a function of detuning for the
wider $544.7G$ resonance. The solid line represents results obtained via the
generalized mean field theory presented in the text, while the dashed-dotted
line represents the mean field results.

##  V Molecule Formation

Moving beyond equilibrium properties, we are also interested in the prospects
for molecule creation upon ramping a magnetic field across the resonance. In
reference bortolotti2006jpb , the mean-field equations of motion for the
system at hand were derived, as follows:

$\displaystyle i\hbar{\partial\over\partial t}\phi$ $\displaystyle=$
$\displaystyle(V_{bg}\rho_{F}+\gamma|\phi|^{2})\phi+g\rho_{MF}^{*}$ (33.a)
$\displaystyle\hbar{\partial\over\partial t}\eta_{F}(p)$ $\displaystyle=$
$\displaystyle-2g\ \Im m(\phi\eta_{MF}(p))$ (33.b)
$\displaystyle\hbar{\partial\over\partial t}\eta_{M}(p)$ $\displaystyle=$
$\displaystyle 2g\ \Im m(\phi\eta_{MF}(p))$ (33.c) $\displaystyle
i\hbar{\partial\over\partial t}\eta_{MF}(p)$ $\displaystyle=$
$\displaystyle\left[\epsilon_{p}^{F}-\epsilon_{p}^{M}-\nu+V_{bg}|\phi|^{2}\right]\eta_{MF}(p)-$
(33.d) $\displaystyle\ \ \ \ g\phi^{*}\left(\eta_{F}(p)-\eta_{M}(p)\right),$

where $\eta_{F}(p)=\langle\hat{a}_{p}^{\dagger}\hat{a}_{p}\rangle$ is the the
fermionic distribution, $\eta_{M}(p)$ its molecular counterpart, and
$\rho_{M,F}=\int{dp\over 2\pi^{2}}p^{2}\eta_{M,F}(p)$ the fermionic and
molecular densities. Similarly
$\eta_{MF}(p)=\langle\hat{c}_{p}^{\dagger}\hat{a}_{p}\rangle$ and is the
distribution for molecule-fermion correlation, with the associated density
$\rho_{MF}$.

Reference bortolotti2006jpb  outlined the limitations of the non-equilibrium
theory, by claiming that to obtain the correct two-body physics in the low-
density limit it would be necessary to include three point and possibly higher
correlations. While this fact is indeed true, we have amended it in the
previous section. For experimentally reasonable parameters, the mean-field
theory can be complemented by the correct renormalized propagator to
accurately describe the equilibrium properties of the gas. Encouraged by this
argument, we now apply it to the problem of a field ramp as well.

In the following we wish to study molecular formation via a time dependent
ramp of the magnetic field across the resonance. To this end we use two
approaches: the first consists of propagating equations (33.a-33.d), ramping
the detuning linearly in time from a large positive value to a large negative
one, and plotting the final molecular population as a function of detuning
ramping rate $R$. The second approach consists in noticing that if $\nu(t)$ is
a linear function of time, then the mean-field Hamiltonian (eq. 22) is ideally
suited to a Landau-Zener treatment, whereby the final molecular population as
a function of detuning can be readily written as

$n_{m}/min(n_{b},n_{f})=1-e^{-R\over\tau}$ (34)

Here $n_{m}/max(n_{b},n_{f})$ is the fraction of possible molecules formed,
and $R=1/{\partial B\over\partial t}$ is the inverse ramp rate, and the
exponential time constant is given by $\tau={\hbar\delta_{B}\over
g^{2}n_{b}}={m_{bf}\over ha_{bg}n_{b}\Delta_{B}}$, where $a_{bg}$ is the
background scattering length, $h$ is Plank’s constant, and $\Delta_{B}$ is the
magnetic field width of the resonance.

Remarkably, the characteristic sweep rate $\tau$ does not depend on the
fermionic density. This arises from the fact that in mean-field theory the
momentum-states of the fermionic gas are uncoupled, except via the depletion
of the condensate. Since in the Landau-Zener approach the depletion is assumed
small, it follows that the various fermionic momentum states are considered
independently, and thus the probability of transition of the gas is equal to
the probability of transition of each individual momentum state. This
approximation is only valid for narrow resonances, such as the $492.5G$
resonance in table 1.

Molecular formation rate versus ramp rate is shown in Figure 12, for the cases
$r_{bf}=6$ (more fermions than bosons) and $r_{bf}=0.6$ (more fermions than
bosons. In both cases the Landau-Zener result agrees nearly perfectly with
direct numerical integration. This agreement is surprising in the case of more
fermions, since the width of the crossing is proportional to density of
leftover bosons, and we expect that this number will change substantially as
the bosonic population is depleted via the formation of molecules. This type
of time dependent crossing should not be properly described by the Landau-
Zener formula. However, monitoring the time evolution of the molecular
population as a function of time shows that the majority of the transfer takes
place quite abruptly somewhat after crossing the zero-detuning region, whereby
the change in bosonic density does not modify the energy levels substantially.

Figure 12:  Transition probability into molecular state via a magnetic field
ramp across the $492.5G$ Feshbach resonance. The solid lines are obtained via
numerical solutions of equations (33.a-33.d), while the dashed lines represent
the Landau-Zener equivalent. In the top panel the gas is composed of more
bosons than fermions $r_{bf}=6$, while in the bottom panel the opposite is
true $r_{bf}=.6$

##  VI Conclusion

In this article, we developed and solved a generalized mean-field theory
describing an ultracold atomic Bose-Fermi mixture in the presence of an
interspecies Feshbach resonance. The theory is “generalized” in the sense that
it correctly incorporates bosonic fluctuations, at least to the level that it
reproduces the correct two-body physics in the extreme dilute limit. This
theory, like any mean-field theory, presents undeniable limitations.
Nevertheless, any useful many-body treatment must start from a well conceived
mean-field theory. Future directions of this work should include the
generalization to finite temperature, and the inclusion of a trap, initially
in a local-density approximation. These advances would be essential to check
for empirical confirmation of the theory.

###### Acknowledgements.

DCEB and JLB acknowledge support from the DOE and the Keck Foundation.

##  Appendix A Green Function Methods for Fermions

In this Appendix we briefly introduce some of the Green function techniques
that we found useful in our calculations.

###  A.1 Free Green functions

We start from the Green function for a gas of free fermions, which is given,
in the frequency-momentum representation by

$G^{0}(w,{\bf q})={1\over\omega-\xi({\bf p})+i\eta\ {\rm sign}(\xi({\bf
p}))},$ (A-35)

where $\xi({\bf p})={\bf p}^{2}/2m-\mu$. The momentum distribution, at
equilibrium, is given by

$n({\bf p})=-i\lim_{\eta\to 0^{+}}\int{{d\omega\over
2\pi}e^{i\omega\eta}G^{0}(w,{\bf q})},$ (A-36)

Here the limit comes from the equilibrium condition; the frequency, in the
Green function definition, is the fourier space equivalent of time, whereby
the real time green function represents the evolution of the system from from
time $t$ to $t^{\prime}$, and the observables obtained this way, represent
expected values of the kind
$\langle\psi(t)|\mathcal{O}|\psi(t^{\prime})\rangle$. However, since we want
equilibrium conditions, we need to take the limit $t\to t^{\prime}$, which is
non trivial, since $G^{0}$ is defined by a green function equation of the form
$\mathcal{L}\ G^{0}(t-t^{\prime})\propto\delta(t-t^{\prime})$, where
$\mathcal{L}$ is some linear operator, and which highlights a peculiar
behavior in the limit we desire. However, since we know on physical grounds
that observables, such as the momentum distribution, must be defined and well
behaved at equilibrium, then by first taking the expectation value integral,
and then the limit, we can circumvent the problem. In eq. A-36 this implies we
cannot quite get rid of the fourier transform exponent
$e^{i\omega(t-t^{\prime})}$ until after the $\omega$ integral.

To perform the integral in A-36, we exploit the fourier exponent, by noting
that since $\eta$ is positive, $e^{i\omega\eta}\to 0$ as $\omega\to+i\infty$,
so that the integral is identical to a contour integral over the path defined
by the real $\omega$ axis, closed in the upper complex $\omega$ plain by an
infinite radius semicircle, which, as we have just seen, gives no contribution
to the integral. We can now integrate using the residue theorem.

We note that the integrand in A-36 has a simple pole at $\omega=\xi({\bf
p})-i\eta\ sign(\xi({\bf p}))$. Thus, if $\xi({\bf p})>0$, then the pole is in
the lower complex plane, and the integral vanishes, and if $\xi({\bf p})>0$,
then the pole is in the upper complex plane, with residue 1. Using the residue
theorem, and summarizing these results we finally get

$n({\bf p})=\Theta(-\xi({\bf p})),$ (A-37)

which we recognize as the zero temperature fermi distribution.

###  A.2 Interacting Green functions

According to Dyson’s equation, the green function for an interacting system
has the form

$G(w,{\bf q})={1\over\omega-\xi({\bf p})-\Sigma(\omega,{\bf p})},$ (A-38)

where $\Sigma(\omega,{\bf p})$ is an arbitrarily complicated function
summarizing all the interactions in the system, which is known as self energy.

The prescription to find $\Sigma$ is quite straightforward, and it consists of
adding all amputated connected feynman diagrams for the system. The fact that,
in general, the number of such diagrams is infinite, makes this task virtually
impossible. Nonetheless, eq. A-38 is very powerful, since it allows one to
include the effect of infinite subsets of the total number of diagrams in the
system, by only having to explicitly calculate a few representative ones.

An alternative standard approach, leads to the exact result (NOTE: Abrikosov
measures energy from $\mu$, here we measure from 0, which is more standard.)

$G(\omega,{\bf p})=\int_{0}^{\infty}{d\omega^{\prime}\left[{A(\omega,{\bf
p})\over\omega-\omega^{\prime}+i\eta}+{B(\omega,{\bf
p})\over\omega+\omega^{\prime}-i\eta}\right]},$ (A-39)

Where $A$ and $B$ are, again, arbitrary complicated functions, though they are
known to be finite.

To understand $A$ and $B$ more closely, we need to introduce the following
well known identity:

$\lim_{\nu\to 0}{1\over x\pm i\nu}=\mathcal{P}{1\over x}\mp i\pi\delta(x),$
(A-40)

where $\mathcal{P}$ is a Cauchy principal value, which represents the
contribution due to a discontinuity in a Riemann sheet (branch cut), and the
delta function represents the contribution due to the pole.

Applying A-40 to A-39 we get

$\displaystyle Re\ G(\omega,{\bf
p})=\mathcal{P}\int_{0}^{\infty}{d\omega^{\prime}\left[{A(\omega,{\bf
p})\over\omega-\omega^{\prime}+i\eta}+{B(\omega,{\bf
p})\over\omega+\omega^{\prime}-i\eta}\right]}$ (A-41) $\displaystyle Im\
G(\omega,{\bf p})=\left\\{\begin{array}[]{r c c}-\pi A(\omega,{\bf
p})&if&\omega>0\\\ \pi B(-\omega,{\bf p})&if&\omega<0\end{array}\right.$
(A-44)

Finally, eq. A-36 represents a fundamental property of green functions, and it
can be generalized to interacting systems simply substituting $G^{0}$ with
$G$. Applying it to eq. A-39, and performing the $\omega$ integral first, we
get

$n({\bf p})=\int_{0}^{\infty}{d\omega^{\prime}B(\omega,{\bf
p})}=\int_{-\infty}^{0}{d\omega{-1\over\pi}Im\ G(\omega,{\bf p}).}$ (A-45)

Introducing the function $\rho(\omega,{\bf p})=-2ImG(\omega,{\bf p})$,
generally called spectral function, the above equation can be written as

$n({\bf p})=\int{{d\omega\over 2\pi}\rho(\omega,{\bf p})\Theta(-\omega)}.$
(A-46)

An important property of the spectral function is that for all ${\bf p}$,

$\int{{d\omega\over 2\pi}\rho(\omega,{\bf p})}=1.$ (A-47)

This can be understood as a sum rule in the following sense: if we wish to
calculate the number of holes in the system, we would take the $\eta\to 0^{-}$
limit in equation A-36. The distribution would then have been $n_{holes}({\bf
p})=1-n({\bf p})=\int_{0}^{\infty}{d\omega^{\prime}A(\omega^{\prime},{\bf
p})}$, so that
$1=\int_{0}^{\infty}{d\omega^{\prime}\left[A(\omega^{\prime},{\bf
p})+B(\omega^{\prime},{\bf p})\right]}=\int{{d\omega\over
2\pi}\rho(\omega,{\bf p})}$.

Using eq.A-38, together with the definition of $\rho$, we can write

$\rho(\omega,{\bf p})={-2Im\Sigma(\omega,{\bf p})\over\left[\omega-\xi({\bf
p})-Re\Sigma(\omega,{\bf p})\right]^{2}+\left[Im\Sigma(\omega,{\bf
p})\right]^{2}}.$ (A-48)

Furthermore, if $\Sigma$ were to be real, or if, equivalently, the pole of the
green function were to be real, for some momentum ${\bf p}$, then taking the
limit $Im\Sigma\to 0$ of A-38, and using eq,A-40, we get

$\rho(\omega,{\bf p})=2\pi\delta(\omega-\xi({\bf p})-Re\Sigma(\omega,{\bf
p})),$ (A-49)

which can be simplified, using the properties of the delta function, to

$\rho(\omega,{\bf p})=2\pi Z({\bf p})\delta(\omega-\omega_{0}({\bf p})),$
(A-50)

where Z, known as spectral weight is given by

$Z({\bf p})={1\over\left|1-{\partial\over\partial\omega}Re\Sigma(\omega,{\bf
p})\right|_{w=\omega_{0}({\bf p})}},$ (A-51)

and $\omega_{0}({\bf p})$ is the pole of the green function, defined by

$\omega_{0}({\bf p})-\xi({\bf p})-\Sigma(\omega_{0}({\bf p}),{\bf p})=0.$
(A-52)

The momentum distribution in this case is thus given by

$n({\bf p})=Z({\bf p})\int{\delta(\omega-\omega_{0}({\bf
p}))\Theta(-\omega)}=Z({\bf p})\Theta(-\omega_{0}({\bf p})).$ (A-53)

## References

  * [1] J. Bardeen, L. N. Cooper, and J. R. Schrieffer.  Phys. Rev., 108(5):1175–1204, Dec 1957.
  * [2] S.N. Bose.  Z. Phys, 26:178–181, 1924.
  * [3] A. Einstein.  Physik- Mathematik, pages 261–267, 1924.
  * [4] M. Greiner, C.A. Regal, and D.S. Jin.  Nature, 426(6966):537–540, 2003.
  * [5] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J.H. Denschlag, and R. Grimm.  Physical Review Letters, 92(12):120401, 2004.
  * [6] MW Zwierlein, CA Stan, CH Schunck, SMF Raupach, AJ Kerman, and W. Ketterle.  Physical Review Letters, 92(12):120403, 2004.
  * [7] E.A. Donley, N.R. Claussen, S.T. Thompson, and C.E. Wieman.  Arxiv preprint cond-mat/0204436, 2002.
  * [8] S. Inouye, J. Goldwin, M. L. Olsen, C. Ticknor, J. L. Bohn, and D. S. Jin.  Phys. Rev. Lett., 93(18):183201, Oct 2004.
  * [9] Francesca Ferlaino, Chiara D’Errico, Giacomo Roati, Matteo Zaccanti, Massimo Inguscio, Giovanni Modugno, and Andrea Simoni.  Phys. Rev. A, 73(4):040702, Apr 2006.
  * [10] CA Stan, MW Zwierlein, CH Schunck, SMF Raupach, and W. Ketterle.  Physical Review Letters, 93(14):143001, 2004.
  * [11] J.J.Zirbel, K.-K.Ni, S.Ospelkaus, J.P D Incao, C.E.Wieman, J.Ye, and D.S.Jin.  Physical Review Letters, 100(11):143201, 2008.
  * [12] Yong-il Shin, A.Schirotzek, C.H.Schunck, and W.Ketterle  Physical Review Letters, 101(15):070404, 2008.
  * [13] C. Ospelkaus, S. Ospelkaus, L. Humbert, P. Ernst, K. Sengstock, and K. Bongs.  Physical Review Letters, 97(12):120402, 2006.
  * [14] R. Roth and H. Feldmeier.  Physical Review A, 65(2):21603, 2002.
  * [15] R. Roth.  Physical Review A, 66(1):13614, 2002.
  * [16] G. Modugno, G. Roati, F. Riboli, F. Ferlaino, R.J. Brecha, and M. Inguscio.  2002\.
  * [17] H. Hu and X.J. Liu.  Physical Review A, 68(2):23608, 2003.
  * [18] X.J. Liu, M. Modugno, and H. Hu.  Physical Review A, 68(5):53605, 2003.
  * [19] L. Salasnich, S.K. Adhikari, and F. Toigo.  Physical Review A, 75(2):23616, 2007.
  * [20] P. Pedri and L. Santos.  Two-Dimensional Bright Solitons in Dipolar Bose-Einstein Condensates.  Physical Review Letters, 95(20):200404, 2005.
  * [21] Sadhan K. Adhikari.  Physical Review A (Atomic, Molecular, and Optical Physics), 70(4):043617, 2004.
  * [22] HP Büchler and G. Blatter.  Physical Review A, 69(6):63603, 2004.
  * [23] R. Kanamoto and M. Tsubota.  Physical Review Letters, 96(20):200405, 2006.
  * [24] S.K. Adhikari.  Physical Review A, 72(5):53608, 2005.
  * [25] A. Albus, F. Illuminati, and J. Eisert.  Physical Review A, 68(2):23606, 2003.
  * [26] M. Lewenstein, L. Santos, MA Baranov, and H. Fehrmann.  Physical Review Letters, 92(5):50401, 2004.
  * [27] R. Roth and K. Burnett.  Physical Review A, 69(2):21601, 2004.
  * [28] A. Sanpera, A. Kantian, L. Sanchez-Palencia, J. Zakrzewski, and M. Lewenstein.  Physical Review Letters, 93(4):40401, 2004.
  * [29] K. Günter, T. Stöferle, H. Moritz, M. Köhl, and T. Esslinger.  Physical Review Letters, 96(18):180402, 2006.
  * [30] I.Titvinidze, M.Snoek, and W.Hofstetter.  Physical Review Letters, 100(14):100401, 2008.
  * [31] MJ Bijlsma, BA Heringa, and HTC Stoof.  Physical Review A, 61(5):53601, 2000.
  * [32] H. Heiselberg, CJ Pethick, H. Smith, and L. Viverit.  Physical Review Letters, 85(12):2418–2421, 2000.
  * [33] DV Efremov and L. Viverit.  Physical Review B, 65(13):134519, 2002.
  * [34] L. Viverit.  Physical Review A, 66(2):23605, 2002.
  * [35] F. Matera.  Physical Review A, 68(4):43624, 2003.
  * [36] AP Albus, SA Gardiner, F. Illuminati, and M. Wilkens.  Physical Review A, 65(5):53607, 2002.
  * [37] D.W. Wang.  Physical Review Letters, 96(14):140404, 2006.
  * [38] S. Powell, S. Sachdev, and H.P. Büchler.  Physical Review B, 72(2), 2005.
  * [39] A. Storozhenko, P. Schuck, T. Suzuki, H. Yabu, and J. Dukelsky.  Physical Review A (Atomic, Molecular, and Optical Physics), 71(6):063617, 2005.
  * [40] Eddy Timmermans, Paolo Tommasini, Mahir Hussein, and Arthur Kerman.  Phys. Rep., 315:199, 1999.
  * [41] S. J. J. M. F. Kokkelmans and M. J. Holland.  Phys. Rev. Lett., 89(18):180401, Oct 2002.
  * [42] S. J. J. M. F. Kokkelmans, J. N. Milstein, M. L. Chiofalo, R. Walser, and M. J. Holland.  Phys. Rev. A, 65(5):053617, May 2002.
  * [43] Y. Ohashi and A. Griffin.  Phys. Rev. Lett., 89(13):130402, Sep 2002.
  * [44] RA Duine and HTC Stoof.  Journal of Optics B: Quantum and Semiclassical Optics, 5(2):S212–S218, 2003.
  * [45] E.M. Lifshitz and L.P. Pitaevskii.  Statistical Physics part 2.  Butterworth-Heinemann, 1996.
  * [46] C Zener.  C. Proc. R. Soc. London A, 137(696), 1932.
  * [47] L.D. Landau.  Phys. Z. Sowjetunion, 1(89), 1932.
  * [48] RA Duine and HTC Stoof.  Physics Reports, 396(3):115–195, 2004.
  * [49] Thorsten Köhler, Thomas Gasenzer, and Keith Burnett.  Phys. Rev. A, 67(1):013601, Jan 2003.
  * [50] M. Holland, S. J. J. M. F. Kokkelmans, M. L. Chiofalo, and R. Walser.  Phys. Rev. Lett., 87(12):120406, Aug 2001.
  * [51] GM Falco and HTC Stoof.  Physical Review Letters, 92(13):130401, 2004.
  * [52] AV Andreev, V. Gurarie, and L. Radzihovsky.  Physical Review Letters, 93(13):130402, 2004.
  * [53] Y. Ohashi and A. Griffin.  Physical Review A, 67(3):33603, 2003.
  * [54] A.V. Avdeenkov.  Journal of Physics B Atomic Molecular and Optical Physics, 37(1):237–246, 2004.
  * [55] AV Avdeenkov and JL Bohn.  Physical Review A, 71(2):23609, 2005.
  * [56] D.C.E. Bortolotti, A.V. Avdeenkov, C. Ticknor, and J.L. Bohn.  Journal of Physics B: Atomic, Molecular, and Optical Physics, 39(1):189–203, 2006.
  * [57] A.V. Avdeenkov, D.C.E. Bortolotti, and J.L. Bohn.  Physical Review A, 74(1):12709, 2006.
  * [58] R.G. Newton.  Scattering Theory of Waves and Particles.  Courier Dover Publications, 2002.
  * [59] J.W. Negele and H. Orland.  Quantum Many-Particle Systems.  Addison-Wesley, 1988.
  * [60] D.C.E. Bortolotti.  unpublished, 2006.
  * [61] B. Marcelis, E. G. M. van Kempen, B. J. Verhaar, and S. J. J. M. F. Kokkelmans.  Phys. Rev. A, 70(1):012701, Jul 2004.


††thanks: corresponding author’s Email<EMAIL_ADDRESS>
#  Existence of an optimized stellarator with simple coils

Guodong Yu  Zhichen Feng  Peiyou Jiang  GuoYong Fu  Institute for Fusion
Theory and Simulation and Department of physics, Zhejiang University, Hangzhou
310027, China

###### Abstract

An optimized compact stellarator with four simple coils is obtained from
direct optimization via coil shape. The new stellarator consists of two
interlocking coils and two vertical field coils similar to those of the
Columbia Non-neutral Torus (CNT)[Pedersen et al. Phys. Rev. Lett. 88, 205002
(2002)]. The optimized configuration has global magnetic well and a low
helical ripple level comparable to that of Wendelstein 7-X (W7-X)[Wolf et al.
Nucl. Fusion 57, 102020 (2017)]. The two interlocking coils have a smooth
three-dimensional shape much simpler than those of advanced stellarators such
as W7-X. This result opens up possibilities of future stellarator reactors
with simplified coils.

The two main approaches of magnetic fusion energy (MCF) are tokamak’s and
stellarator’s. Tokamak is currently the dominant approach with advantages of
axisymmetric geometry and achieved plasma parameters significantly better than
those of other MCF devices. However, stellarators have recently enjoyed a
renaissance as recent results of the advanced stellarator Wendelstein 7-X
(W7-X)[1] demonstrated the reduced neoclassical energy transport[2][3]. It is
expected that plasma performance of W7-X can reach a level comparable to that
of an equivalent tokamak in next few years. Stellarators have advantages of
naturally steady state operations because their magnetic fields are mainly
generated by external coils and plasma current is not needed for confinement
and is usually quite small. The harmful current-driven instabilities such as
disruptions in tokamaks are absent in stellarators. The confinement properties
of stellarator plasmas are largely determined by external coils. The recent
success of W7-X demonstrates that it is possible to build the optimized three
dimensional (3D) coils to such high precisions that the designed plasma
confinement performance can be achieved as predicted in actual experiments. It
has been argued recently that the stellarator approach provides the fastest
track to the realization of fusion energy because favorable plasma confinement
properties can be designed, realized and controlled almost fully by external
3D coils[4].

The advanced stellarators such as W7-X are usually designed for improved
neoclassical confinement and MHD stability by two stage approaches. First, the
shape of plasma boundary is varied to obtain optimized properties. Second, the
optimized boundary shape is realized by designing the shapes of 3D coils.
Unfortunately, the resulting 3D coils are usually quite complex and are
difficult and costly to build. This is evidenced by the long delay of W7-X
project [5] and the cancellation of the NCSX project [6] due to difficulties
in engineering and building of 3D coils. Therefore it is necessary to simplify
the 3D coils for stellarator to become an economical platform for fusion
reactors. For this reason, recently, coil simplification of stellarator has
been a subject of intensive studies including reducing coil complexity in the
optimization[7] [8][9] as well as using permanent magnets to replace 3D
coils[10][11][12].

In this work, we explore the possibilities of optimized stellarators with
simplified coils. We use a direct optimization method of varying coil shape to
optimize both plasma confinement and MHD stability. In this way we can
effectively control the complexity of coils while optimizing the plasma
properties. The details of optimization method will be given later. Following
our recent work[13], we choose the coil topology of the Columbia Non-neutral
Torus (CNT)[14] for our optimization space. CNT consists of two circular
interlocking(IL) coils and two circular vertical field (VF) coils. It is
arguably the simplest compact stellarator ever built and successfully
operated[14]. Here we carry out a global optimization of both neoclassical
confinement and MHD stability by varying the 3D shape of the two IL coils. An
initial phase of this work has been done recently by targeting only the
neoclassical transport. The results showed that the effective helical ripple
level was reduced by an order of magnitude as compared to that of CNT
indicating significant improvement in neoclassical confinement[13]. This was
realized by two planar IL coils with elliptical shape. In present work, a much
broader parameter space of coil shape is explored allowing 3D shape of the two
IL coils. As a result, an optimized Compact Stellarator with Simple
Coils(CSSC), to be called CSSC, has been discovered. The CSSC has both global
magnetic well and a low level of helical ripple comparable to that of W7-X.
This breakthrough opens up possibilities of stellarator reactors with
simplified coils. The design of CSSC can be used as basis for a low-cost
experiment to study the physics of compact stellarators. The design could also
be used as a candidate stellarator for electron-positron plasma
experiments[14].

We now describe the configuration specification and physics properties of CSSC
in details. In our optimization, the radius of the two circular VF coils is
varied while the distance between them is fixed. The shapes of the two IL
coils are varied and are described by following Fourier representation[15]

$\displaystyle x=x_{c,0}+\sum_{n=1}^{n_{f}}x_{c,n}\cos(nt),$ (1)
$\displaystyle y=\sum_{n=1}^{n_{f}}y_{s,n}\sin(nt),$ (2) $\displaystyle
z=\sum_{n=1}^{n_{f}}z_{s,n}\sin(nt),$ (3)

where $(x,y,z)$ are Cartesian coordinates of line coils, $t$ is an angle
variable ranges $[0,2\pi]$ and $n$ is the harmonic number. The total number of
Fourier coefficients for each IL coil is $3\times n_{f}+1$ with $n_{f}$ being
the maximum harmonic number. In order to maintain two periods of flux
surfaces, the two IL coils are constrained to have same shape. Furthermore,
$n_{f}=3$ is chosen in this work. This choice of $n_{f}=3$ is a compromise
between allowing enough degree of freedom while constraining the coil
complexity. Finally, taking into account the current ratio between IL coils
and VF coils, the optimization parameter space has a total of 12 degree of
freedom.

Table 1 shows the Fourier coefficients of the optimized configuration CSSC
(4th column) along with the optimization ranges for each Fourier coefficient
(3rd column).

| Parameter | Range | Results |  |  |
---|---|---|---|---|---|---
IL coil | $x_{c,0}$ | 0.3$\sim$0.6 | 0.509 |  |  |
$x_{c,1}$ | 0.6$\sim$1.2 | 1.051 |  |  |
$x_{c,2}$ | -0.3$\sim$0.3 | 0.227 |  |  |
$x_{c,3}$ | -0.1$\sim$0.1 | -2.622e-2 |  |  |
$y_{s,1}$ | 0.6$\sim$1.2 | 0.806 |  |  |
$y_{s,2}$ | -0.3$\sim$0.3 | 0.162 |  |  |
$y_{s,3}$ | -0.1$\sim$0.1 | -6.257e-2 |  |  |
$z_{s,1}$ | -0.4$\sim$0.8 | 0.667 |  |  |
$z_{s,2}$ | -0.2$\sim$0.2 | -3.87e-2 |  |  |
$z_{s,3}$ | -0.1$\sim$0.1 | -6.052e-2 |  |  |
PF coil | $y_{s,1}=x_{c,1}$ | 1.0$\sim$2.0 | 1.850 |  |  |
$z_{c,0}$ | 0.77 | 0.770 |  |  |
Current ratio | $I_{IL}/I_{PF}$ | 1$\sim$3 | 1.323 |  |  |
Table 1: The optimization ranges of coil Fourier coefficients and coil current
ratio are in the third column. The optimized results are in the fourth column.

Fig. 1 shows the 3D schematic of coils and the last closed flux surface (LCFS)
of CSSC. We observe that the two IL coils have a smooth shape with modest 3D
variation. The averaged radius of the IL coils is $1.11m$. The radius of
vertical field coils have a radius of $1.85m$. The distance between the two
vertical field coils is $1.54m$. The current ratio of IL coils and VF coils is
$1.323$. The minor and major radius of LCFS is $0.15m$ and $0.58m$
respectively, corresponding to a low aspect ratio of $R/a=3.87$. The volume of
LCFS is $0.27m^{3}$. The Poincare plots of vacuum magnetic surfaces are shown
in Fig. 2 at toroidal angels $0^{\circ}$ and $90^{\circ}$ respectively. The
configuration has good magnetic flux surfaces except for a small magnetic
island at the $\iota=0.25$ surface. The rotational transform profile is shown
in Fig. 3 and is seen to vary from $0.2$ at the magnetic axis to $0.28$ at
LCFS with modest magnetic shear.

Figure 1: The optimized stellarator with four simple coils: the gold color
denotes the two inner interlocked 3D coils located between two circular
vertical field coils; the purple color denotes the last closed flux surface.

(a) $\phi=0^{\circ}$

(b) $\phi=90^{\circ}$

Figure 2: Poincare plot of vacuum magnetic surfaces at toroidal angels
$0^{\circ}$ (left) and $90^{\circ}$ (right).  Figure 3:  Rotational transform
iota versus normalized toroidal flux.

One of our main optimization targets is neoclassical confinement, particularly
the so-called $1/\nu$ neoclassical transport with $\nu$ being the collision
frequency of plasmas. This $1/\nu$ scaling is very unfavorable for stellarator
confinement for high temperature regime of fusion plasmas because the
transport increases strongly with plasma temperature[16]. Therefore, the
minimization of neoclassical transport in the $1/\nu$ regime is necessary for
stellarator reactors. It has been shown that the $1/\nu$ transport is
proportional the effective helical ripple coefficient
$\epsilon_{eff}^{3/2}$[17]. We use this coefficient as a target. Fig. 4 shows
the radial profile of $\epsilon_{eff}^{3/2}$ of the CSSC as well as W7-X and
the Large Helical Device (LHD)[18]. We observe that $\epsilon_{eff}^{3/2}$ of
CSSC is comparable to that of W7-X and is much smaller than that of LHD.

Figure 4:  Effective helical ripple $\epsilon_{eff}^{3/2}$ versus square root
of the normalized toroidal flux $\sqrt{\psi/\psi_{edge}}$ for the optimized
configuration (red line), W7-X (blue line)[3] and LHD (green line)[3].

The other important target of optimization is magnetic well which is a key
parameter for MHD stability. The magnetic well is defined via the surface-
averaged magnetic pressure[19] as following:

$W=2\frac{V}{\left\langle
B^{2}\right\rangle}\frac{d}{dV}\left\langle\frac{B^{2}}{2}\right\rangle$ (4)

where $B$ is magnetic field strength, $V(\psi)$ is the volume within flux
surface $\psi$. A positive gradient of the surface-averaged magnetic field
corresponds to a magnetic well. The magnetic well profile of CSSC is
calculated with free boundary condition using the VMEC code for several plasma
beta values. The results are plotted in Fig. 5. We see that CSSC possess a
global magnetic well profile with well depth increasing with plasma beta. In
addition to magnetic well, we have also carried out initial MHD stability
calculation using the 3D MHD stability code TERPSICHORE [20] in order to
confirm a more complete MHD stability of the CSSC. The results show that low-n
global MHD modes are stable up to the volume-averaged beta of $\beta=1\%$. It
should be noted that our optimization is done for vacuum magnetic field.
Future work will consider optimization at finite plasma beta.

Figure 5: Magnetic well depth $W$ versus the normalized toroidal flux
$\psi/\psi_{edge}$ for three values of plasma $\beta$.

So far, we have shown that the optimized configuration CSSC has both magnetic
well and low level of helical ripple. Now we show that these favorable
properties are achieved with relatively simple coils, i.e., the shape of the
two IL coils is relatively smooth with modest 3D variation as compared to that
of W7-X. Fig. 6 compares the curvature of the optimized configuration with
that of one of W7-X’s module coils at the same coil length of $1m$. We observe
that the coil curvature of CSSC is significantly smaller than that of the W7-X
coil. Also the variation of coil curvature along the coil loop is much less
complex indicating that the coil complexity of CSSC is much less in
comparison.

Figure 6: Coil curvature $\kappa$ variation along the loop of coil for the
optimized configuration (black line) and W7-X (red).

We now describe some details of the optimization method used to obtain above
results. In addition to neoclassical transport and magnetic well, we also
target plasma volume, rotational transform and flux surface quality. The
magnetic field generated by current-carrying-coils is calculated from the
Biot-Savart law. The flux surfaces and rotational transform profile are
calculated by following the magnetic field lines. The effective ripple of each
surface is calculated by integrating along a magnetic field line [17] as done
in our recent work[13]. We vary the shape of the two IL coils via the ten
Fourier coefficients as described above to search for desirable configurations
that meet our targets. Since the total degree of freedom is not a small
number, a multi-stage random search algorithm is adopted in the optimization
process. Specifically, at the first stage, the initial ranges of each degree
of freedom are selected as given in the TABLE 1. It should be noted that the
range of each Fourier coefficient is constrained to be smaller for higher $n$
in order to keep the coil complexity low. Given the ranges of parameters, a
large collection of stellarator configurations are generated in the
12-dimension parameter space with value of each parameter randomly selected
within its range. For the first stage, a total of $80$ million cases are
evaluated. A few good configurations are found that meet the criteria of
$\epsilon_{eff}^{3/2}<0.01$, $\iota>0.2$, plasma volume $V>0.25m^{3}$, and
having a global magnetic well profile. This concludes the first stage of
search. At the second stage, this random search process is repeated starting
from each of the good configurations found in the first stage. The ranges of
parameters are updated based on each new configuration and are chosen to be
narrower than those of the first stage. A total of four such stages of search
have been carried out in arriving the final optimized configuration CSSC as
defined in TABLE 1. It should be noted that there are other optimized
configurations with comparable quality. The details of other configurations
will be given in another paper.

In conclusion, a direct optimization from coils demonstrates the existence of
a new optimized stellarator with four simple coils. The optimized
configuration has favorable properties of magnetic well and a low effective
helical ripple level comparable to that of W7-X. The two 3D interlocking coils
have a smooth shape much simpler as compared to those of advanced stellarators
such as W7-X. This work opens up possibilities of future stellarator reactors
with simplified coils. In near term, the optimized configuration can be used
as a basis for a low-cost experiment to study the MHD stability and plasma
confinement of compact stellarators. The design could also be used as a
candidate stellarator for electron-positron plasma experiments.

## Acknowledgement

We thank Dr. Steve Hirshman for the use of the 3D equilibrium code VMEC code.
We also thank Dr. W. A. Cooper for use of the 3D MHD stability code
TERPSICHORE, and Dr. David Gates and Dr. Caoxiang Zhu for use of the
stellarator optimization code STELLOPT. This work was funded by the start-up
funding of Zhejiang University for one of the authors (Prof. Guoyong Fu).

## References

  * [1] R. C. Wolf, A. Ali, A. Alonso, J. Baldzuhn et al., Nucl. Fusion 57, 102020 (2017).
  * [2] A. Dinklage, C. D. Beidler, P. Helander et al., Nat. Phys. 14, 855-860 (2018).
  * [3] C. D. Beidler et al., Nature 596, 221 (2021).
  * [4] A. H. Boozer, Nucl. Fusion 61, 096024 (2021).
  * [5] K. Riße et al., Fusion Engng Des. 84, 1619 (2009).
  * [6] G. H. NEILSON, C. O. GRUBER, J. H. HARRIS et al., IEEE Trans. Plasma Sci. 38, 320 (2010).
  * [7] J. F. Lobsien, M. Drevak, T. S. Pedersen, Nucl. Fusion 58, 106013 (2018).
  * [8] T. G. Kruger, C. Zhu, A. Bader et al., J. Plasma Phys 87, 175870201 (2021).
  * [9] S. A. Henneberg, S. R. Hudson, D. Pfefferle et al. J. Plasma Phys. 87, 905870226 (2021).
  * [10] P. Helander, M. Drevak, M. Zarnstorff, S. C. Cowley, Phys. Rev. Lett. 124, 095001 (2020).
  * [11] C. X. Zhu, M. Zarnstorff, D. Gates and A. Brooks, Nucl. Fusion 60, 076016 (2020).
  * [12] G. S. Xu, Z. Y. Lu, D. H. Chen. Nucl. Fusion 61, 026025 (2021).
  * [13] G. D. Yu, Z. C. Feng, P. Y. Jiang and G. Y. Fu, Phys. Plasmas 28, 092501 (2021).
  * [14] T. S. Pedersen, J. R. Danielson, C. Hugenschmidt et al. New Journal of Physics 14 (2012) 035010.
  * [15] C. X. Zhu, S. R. Hudson, Y. T. Song et al. Nucl. Fusion 58, 016008 (2017).
  * [16] D. D.-M. Ho and R. M. Kulsrud. Physics of Fluids, 30, 442-461 (1987).
  * [17] V. V. Nemov, S. V. Kasilov, W. Kernbichler, and M. F. Heyn. Physics of Plasmas, 6, 4622-4632 (1999).
  * [18] A. Iiyoshi, M. Fujiwara, O. Motojima et al. Fusion Technol. 17, 169-187 (1990).
  * [19] [1] J. P. Freidberg, Ideal MHD[M], Cambridge University Press (2014).
  * [20] D.V. Anderson, W.A. Cooper, R. Gruber, S. Merazzi, and U. Schwenn, The International Journal of Supercomputing Applications, 4, 34-47 (1990).