# Ultrafast, low-energy, all-optical switch in polariton waveguides

D. G. Suárez-Forero  Present address: Joint Quantum Institute, University of
Maryland, College Park, MD, USA CNR NANOTEC, Institute of Nanotechnology, Via
Monteroni, 73100 Lecce, Italy  Dipartimento di Fisica, Università del Salento,
Strada Provinciale Lecce-Monteroni, Campus Ecotekne, Lecce 73100, Italy  F.
Riminucci  Dipartimento di Fisica, Università del Salento, Strada Provinciale
Lecce-Monteroni, Campus Ecotekne, Lecce 73100, Italy  Molecular Foundry,
Lawrence Berkeley National Laboratory, One Cyclotron Road, Berkeley,
California, 94720, USA  V. Ardizzone<EMAIL_ADDRESS>CNR
NANOTEC, Institute of Nanotechnology, Via Monteroni, 73100 Lecce, Italy  A.
Gianfrate  CNR NANOTEC, Institute of Nanotechnology, Via Monteroni, 73100
Lecce, Italy  F. Todisco  CNR NANOTEC, Institute of Nanotechnology, Via
Monteroni, 73100 Lecce, Italy  M. De Giorgi  CNR NANOTEC, Institute of
Nanotechnology, Via Monteroni, 73100 Lecce, Italy  D. Ballarini  CNR NANOTEC,
Institute of Nanotechnology, Via Monteroni, 73100 Lecce, Italy  G. Gigli  CNR
NANOTEC, Institute of Nanotechnology, Via Monteroni, 73100 Lecce, Italy  K.
Baldwin  PRISM, Princeton Institute for the Science and Technology of
Materials, Princeton Unviversity, Princeton, NJ 08540  L. Pfeiffer  PRISM,
Princeton Institute for the Science and Technology of Materials, Princeton
Unviversity, Princeton, NJ 08540  D. Sanvitto  CNR NANOTEC, Institute of
Nanotechnology, Via Monteroni, 73100 Lecce, Italy  INFN, Sez. Lecce, Via per
Monteroni, Lecce, Italy

###### Abstract

The requirement for optical-electrical-optical conversion of signals in
optical technologies is often one of the majors bottleneck in terms of speed
and energy consumption. The use of dressed photons (also called polaritons),
that allows for intrinsic sizable interactions, could significantly improve
the performances of optical integrated elements such as switches or optical
gates. In this work we demonstrate the ultrafast switch of a laser coupled
into a polaritonic waveguide triggered by an optical pulse resonant with the
same dispersion but at a lower energy. Our experiments show two effects
capable to interrupt the transmission of the laser in two different time
ranges: a sub-picosecond time range due to the optical Stark effect, and a
picosecond range governed by the creation of a charge reservoir. In the latter
regime we found that at certain power of excitation the activation of dark
states allows for a long persistence of the switching much beyond the bright
exciton lifetime.

A fundamental component in the development of optical data processing is the
ultrafast all-optical switch. Indeed, there is a growing interest in the
development of fully optical devices as promising alternative for the latency
of electro-optical systems, especially for applications in which electronics
is only needed for amplification or modification of optical processes.
Prototypes of systems capable to switch an optical signal in ultra-short times
have been demonstrated in microring resonatorsSethi2014 , photonic
crystalsNozaki2010 , III-V semiconductor quantum wellsTakahashi2000 ; Li2007 ,
grapheneLi2014 , combined systems of gold and grapheneono2020  and exciton-
polaritons systems Ma2020 . Exciton-polaritons, or polaritons, are hybrid
light-matter quasi-particles arising when an electromagnetic mode confined in
an optical resonator is coupled to an excitonic transition. Polaritons have
demonstrated great potentiality for the development of all-optical
technologies such as transistorsLeyder2007 ; Ballarini2013 ,
routersMarsault2015 , all-optical switchesArdizzone2013  and polarization
rotatorsGnusov2021 . Among all the platforms sustaining polaritons,
semiconductor waveguidesWalker2013  (WG) could be the cornerstone for building
integrated polaritonic circuits. They are attracting for their capability to
sustain fast-propagating modes with high quality factors, and importantly,
with simple fabrication processes. For example, optical confinement is
achieved through total internal reflection with no need for growing thick
distributed Bragg reflectors (DBRs) as it is the case for standard
semiconductor microcavities. Waveguide geometry also allows the manipulation
of the polariton-polariton interaction (and hence the optical nonlinearity)
through external electric fieldsRosenberg2016 ; Rosenberg2018 ; Suarez-
Forero2021 . Integrated devices such as polariton lasersJamadi2018 ; Suarez-
Forero2020-2  have also been demonstrated using polariton waveguides. Here we
demonstrate an ultrafast switching off of a polariton flow by means of another
optical pulse that uses the same excitation path.

As it is well known any excitonic system undergoes the so called optical Stark
effect when hit by a strong optical pulse off resonance with the semiconductor
optical transition. In case of polaritons this effect results in an abrupt
loss of the strong coupling due to the change in the exciton
resonanceHayat2012 ; Lange2018 . As a consequence the polariton mode shifts in
energy, regaining its original position only after the pulse is gone. In this
work we will exploit this effect to turn off a cw laser for as long as the
pulse duration. As we will see, a second effect can also take part in this
process, depending on the pulse power. While for low power the Stark effect
dominates the dynamics, for higher power an exciton reservoir is created due
to two photons absorption. This adds a longer recovery time caused by the
blueshift induced by the exciton-exciton interactions and similarly observed
in Ref. Giorgi2012 . At the highest powers not only a significant reservoir is
created, but also several dark excitonic states which lead to very long
recovery times beyond our laser repetition rate. Such dark states have also
been discussed in other works Walker .

Figure 1: a) Sketch of the switch-off mechanism. A CW laser (orange) is
injected resonantly inside the WG, creating a continuum flux of polaritons.
When the fs pulse (red) couples into the system the transmission is
interrupted via two different effects: the Optical Stark effect that modulates
the Rabi splitting of the polariton mode, and the creation of a charge-
carriers reservoir via two-photons absorption that quenches the transmission
of the CW by drastically shifting the exciton energy. b) Angular and energy
resolved photoluminescence emission of the system collected from the
outcoupling grating. The theoretical fits of the $TE_{0}$ and $TM_{0}$ modes
are indicated by dashed orange and red lines, respectively. c) Far field
emission of the outcoupling grating when both the CW and fs lasers are
resonantly injected into the system.

We use a WG with 12 GaAs quantum wells (QWs) 20 nm thick, separated by
Al0.4Ga0.6As barriers of 20 nm, covered with a 10 nm GaAs capping layer. The
structure has a 500 nm cladding of Al0.8Ga0.2As. To reduce the mode volume, a
microwire structure 40 nm deep, 1 $\mu$m wide and 140 $\mu$m long is etched on
the surface. A pair of focusing gratings is fabricated at both ends of the
microwire for the injection and extraction of the signal and control beams. A
sketch of the experimental realization is displayed in Fig. 1a and a detailed
description of the full experimental setup can be found in the SM. The
polaritonic flux inside the WG is interrupted when a laser pulse (control) at
lower energy is coupled inside the structure. Figures 1b-c show the guided
modes dispersion and the energies and in-plane momenta at which the signal and
control lasers are coupled into the system. As it can be observed in the
figure, the signal laser is tuned to inject population in the $TE_{o}$ mode at
1525.14 meV, 4.9 meV below the excitonic resonance. The control laser is
injected into the $TM_{o}$ mode and tuned at around $1500$ meV; far from the
exciton resonance to avoid the presence of population created by
photoluminescence and to exploit the Optical Stark effect (Fig. 1c).

We start by using a CW laser as signal. The transmitted light is directed into
a streak camera and we observe the modification in the transmission of the
signal as a function of time for different powers of the control beam and for
different energies of the signal beam. In each case the transmission is
normalized to the intensity of the continuous signal during the first 200 ps
($\sim-400$ ps to $-200$ ps). The need for this normalization comes from the
fact that the total CW signal at the initial time depends on the pulse power,
suggesting that the system does not fully empty between consecutive pulses (12
ns), as we will discuss later (see also SM). The results are shown in Fig. 2.
Panel a corresponds to the case in which the CW laser has an energy of 1525.14
meV; 4.9 meV below the excitonic resonance. As the fs laser power increases,
the effect becomes evident: at $\tau=0$, when the pulse reaches the system,
the transmission is reduced to values lower than 0.5 in less than 20 ps (see
inset). Notice from Fig. 2a how a power of 66.5 mW in the pulsed laser
corresponds to a switching time of less than 20 ps with a switching energy of
$\approx 0.8$ nJ. The increasing recovery time with pump power, indicates the
creation of a population of charge carriers through two-photons absorption.

Figure 2: a) Time resolved transmission of the CW laser when injected
resonantly into the WG. At $\tau=0$ the fs pulse reaches the system, modifying
the excitonic resonance and hence, reducing the transmission of the CW laser
by a factor $\approx 2$ in less than 20 ps. The effect has a strong dependence
on the fs laser power. 5 different powers are displayed in the figure. The
inset shows the dynamics for a pump power of 280 mW, but measured with a
higher temporal resolution, allowing to appreciate the switch-off speed. b)
Analysis of the switching effect at different energies of the CW laser. The
energies closer to the excitonic resonance are more affected than the photonic
ones, evidencing the excitonic origin of the effect.

The excitonic nature of the effect, implies a dependence on the energy of the
continuous flux of polaritons (respect to the excitonic resonance). In the
case of Fig. 2a, the CW population is injected at an excitonic fraction
$C_{X}=0.69$. After fixing the pulsed laser power to 350 mW, the experiment is
repeated for other excitonic fractions (Fig. 2b). These results confirm that
the switching-off has an excitonic origin: the higher the excitonic fraction
of the CW beam, the lower the minimum transmission reached after the pulse
arrival.

The data collected with the streak camera allow to clearly reconstruct the
recovery of the cw transmission after the fs control has passed. To
reconstruct the much faster switch off of the transmission we use a homodyne
detection technique. The fs control is used to generate a fs signal through
the super-continuum (SC) generation effect (Fig. S1). This technique allows us
to achieve sub-ps temporal resolution. We reconstruct the dynamics by using a
delay line that changes the arrival time of the fs control pulse respect to
the probe signal (see also SM). In this regime, the Optical Stark effect
governs the dynamics Hayat2012 ; Lange2018 . It modulates the light-matter
coupling during the pulse temporal width. A delay line provides the temporal
control needed to obtain a precision on the order of tens of fs. By varying
the arrival time of the fs control pulse, and using the SC pulse to probe the
system, we reconstruct these ultrafast dynamics ($\tau\equiv t_{SC}-t_{fs}$).

Figure 3: Characterization of the optical Stark effect in the microwire
structure. a) Differential transmission of the SC laser at a fixed in-plane
momentum and for different arrival times of the fs pulse. b) vertical profile
at a fixed energies (dashed rectangle in panel a). With the arrival of the
control pulse, the light-matter coupling is modulated, generating an
instantaneous energy shift, that causes a reduction of the transmitted
intensity.

The results, presented in Fig. 3 show the differential transmission, defined
as: $T_{total}-T_{fs}-T_{SC}$, where $T_{total}$ is the transmission of the SC
pulse in presence of the fs pulse (for different time delays), $T_{fs}$ is the
system transmission when only the fs pulse is present and $T_{SC}$ is the
transmission when only the SC pulse reaches the WG. The instantaneous shifting
at $\tau\sim 0$ (Fig. 3a) is the clear signature of the optical Stark effect,
that modulates the light-matter coupling, shifting the lower polariton branch.
Panel b displays a profile taken at the original polariton energy for a given
k, showing the instantaneous modification of the transmission thanks to the
Stark modulation at $\tau=0$. For longer times, the control beam arrives
before the SC pulse, generating a bath of charge carriers via two-photon
absorption. This explains the lower differential transmission at $\tau>0$.
Although a streak camera cannot resolve the Stark ultrafast dynamics, it
allows to observe the effect of the reservoir on the transmission of a
continuum polariton flux.

It is worth to mention that the additional transverse confinement provided by
the etched microwire, is essential to observe the switching effect. The
experiment was repeated in identical conditions in wider structures of 20
$\mu$m, but no reduction in the transmission was observed after the injection
of the fs pulse, indicating that not only the pulse must be injected into the
structure, but its modal volume has to be reduced. Moreover, by pumping with
the same energy of the fs laser, but with an in-plane momentum that does not
match the guided mode (and hence it does not couple into the WG), the
continuous polariton flux is not interrupted after the pulse arrival,
indicating that the effect is not present at the grating coupler for the two
lasers. In a similar way, by tuning the energy of the pulsed laser out of the
resonant condition to lower energies, the flux is not interrupted either. The
population of charge carriers could also be created by tuning the fs pulse at
energies above the material’s band-gap. The dynamics obtained in this
configuration (Fig. S2) show the desired switching effect, but the unavoidable
PL signal constitutes an undesirable noise in the detection system, hindering
a practical application of this effect. Finally, no dependence on the CW power
was detected.

As mentioned above, the initial transmission for each fs pulse power varies
(see Fig. S3), which suggests that the cavity does not fully empty between
consecutive pulses. This indicates that the reservoir has an important
component of dark-long-lived excitons whose lifetime is, at least, comparable
with the delay between fs pulses (12 ns). Importantly, this particular
behavior does not affect whatsoever the switching condition, since the
important quantity in this case is the relative transmission, moreover it
could be exploited to constrain the system to remain in its “off” state during
very long times compared to the switching time. This point is illustrated and
discussed in the Supplementary Material. For a theoretical modelization of
this effect we rely on a Gross-Pitaevskii (GP) formalism for a population of
polaritons in presence of a reservoir of excitons and charge
carriersGargoubi2016 . As demonstrated in previous works, the nonlinearity has
an apparent increase due to the presence of such particlesWalker , which is
the origin of the energy shift that affects the polariton transmission through
the WG. The 1D confinement provided by the $1$ $\mu$m wire, justifies the one-
dimensional model used here. In this conditions the system behavior can be
described by the set of equations 1.

$\displaystyle i\hbar\frac{\partial\Psi(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x^{2}}+g_{XX}\cdot
C_{X}\left[|\Psi(x,t)|^{2}+n_{R}(x,t)\right]-i\left(C_{X}\gamma_{e}+\left[1-C_{X}\right]\gamma_{p}\right)\right]\Psi(x,t)+P_{CW}(x,t)$
(1a) $\displaystyle\hbar\frac{\partial n_{R}(x,t)}{\partial t}$
$\displaystyle=$ $\displaystyle-2\gamma_{R}n_{R}(x,t)+P_{fs}(x,t)$ (1b)

In this set of coupled equations $\Psi(x,t)$ represents the wavefunction of
the polariton injected resonantly by the CW laser, $g_{XX}$ is the strengh of
the exciton exchange interaction, $C_{X}$ the excitonic fraction,
$\gamma_{e}$, $\gamma_{p}$ and $\gamma_{R}$ the decay rates of the (coupled)
exciton, photonic mode and long-lived reservoir population, respectively.
Finally, $n_{R}(x,t)$ represents the particles density of the reservoir. The
pumping terms $P_{CW}$ and $P_{fs}$ stand for the external injection of
resonant polaritonic population and exciton reservoir by two photon
absorption. Their functional form will be determined by the experimental pump
configuration. In our case we focus the CW (that resonantly injects population
of polaritons) and fs (that indirectly injects population into the long-lived
charge-carriers reservoir) lasers into the input grating at $x=0$, so (in the
rotating frame) their functional form is:

$\displaystyle P_{CW}(x,t)$ $\displaystyle=$ $\displaystyle
p_{CW}e^{-\left(x/\sigma_{CW}\right)^{2}}e^{i\Delta_{\omega}t}$ (2a)
$\displaystyle P_{fs}(x,t)$ $\displaystyle=$ $\displaystyle
N(1-e^{-\epsilon\cdot
p_{fs}^{2}})e^{-(x/\sigma_{fs})^{2}}e^{-\left(\frac{t-t_{o}}{\tau}\right)^{2}}$
(2b)

Where $p_{CW}$ and $p_{fs}$ are proportional to the continuum and pulsed
lasers intensity, respectively, $\sigma_{CW}$ and $\sigma_{fs}$ determine the
waist size of the CW and pulsed lasers, respectively, $\hbar\Delta_{\omega}$
is the energy difference between the excitonic resonance and guided mode at
the in-plane momentum of the resonant pump, $t_{o}$ is the arrival time of the
fs pulse and $\tau$ depends on its temporal width. The dependence of $P_{fs}$
with $p_{fs}$ takes into account the saturable nature of an excitonic
reservoir, while the quadratic dependence of $P_{fs}$ takes into account the
effect that creates the reservoir: two photons absorption with an efficiency
$\epsilon$. The complete explanation of the model and a table with the values
of the parameters employed in the simulations (and their experimental
deduction) can be found in the SM.

The transmission calculated with this method and the experimental data of the
time-resolved transmission for the pump power 66.5 mW are presented in Fig.
4a. The good correspondence between theory and experiment confirms the
hyphotesis of the reduction in the transmission due to the renormalization of
the energy dispersion induced by the charge-carriers reservoir.

Figure 4: a) Experimental time-resolved transmission for the case of a pulsed
laser power of 66.5 mW (continuous line) with the theoretical fitting from the
GPE system of equations 1 (dashed line). b) Power dependence of the CW laser
transmission minimum. The continuous red line corresponds to the reproduction
of this behavior by the theoretical model.

In summary, we demonstrated an ultrafast all-optical switching of a CW laser
in a polariton WG system. The ultrafast switching behaviour comprises two
distinct physical effects spanning different temporal scales and powers: an
ultrafast modulation of the light-matter coupling that occurs in a time range
of hundreds of fs due to the optical Stark effect, and a renormalization of
the polariton energy due to the creation of a charge carriers reservoir. While
the former effect takes place within the duration of the pulse beam, the
latter effect occurs in a time range of tens of ps. A comparison of our system
with state-of-the-art all-optical switchesRutckaia2020  reveals that the
system studied here is among the fastest all-optical switches demonstrated so
far. Both switching effects require energies of approximately 0.8 nJ, a
quantity that can be reduced by either further reducing the photonic modal
volume or by the optimization of the pulse injection through the grating.
These results provide an additional step in the road towards fully optical
technologies, capable to operate without the need of an optical-electronic-
optical conversion.

Acknowledgments
We thank Paolo Cazzato for technical support.
We are grateful to Prof. Ronen Rapaport for inspiring discussions and for
sharing information about the sample design.
The authors acknowledge the Italian Ministry of University (MUR) for funding
through the PRIN project “Interacting Photons in Polariton Circuits” –
INPhoPOL (2017P9FJBS_001).
Work at the Molecular Foundry was supported by the Office of Science, Office
of Basic Energy Sciences, of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231.
We thank Scott Dhuey at the Molecular Foundry for assistance with the electron
beam lithography.
We acknowledge the project FISR - C.N.R. “Tecnopolo di nanotecnologia e
fotonica per la medicina di precisione” - CUP B83B17000010001 and ”Progetto
Tecnopolo per la Medicina di precisione, Deliberazione della Giunta Regionale
n. 2117 del 21/11/2018.
This research is partly funded by the Gordon and Betty Moore Foundation’s
EPiQS Initiative, Grant GBMF9615 to L. N. Pfeiffer, and by the National
Science Foundation MRSEC grant DMR 1420541.

Competing Interest
The authors declare no competing interest.

## References

  * [1] Purnima Sethi and Sukhdev Roy.  All-optical ultrafast XOR/XNOR logic gates, binary counter, and double-bit comparator with silicon microring resonators.  Applied Optics, 53(28):6527, oct 2014.
  * [2] Kengo Nozaki, Takasumi Tanabe, Akihiko Shinya, Shinji Matsuo, Tomonari Sato, Hideaki Taniyama, and Masaya Notomi.  Sub-femtojoule all-optical switching using a photonic-crystal nanocavity.  Nature Photonics, 4(7):477–483, jul 2010.
  * [3] Ryo Takahashi, Hiroki Itoh, and Hidetoshi Iwamura.  Ultrafast high-contrast all-optical switching using spin polarization in low-temperature-grown multiple quantum wells.  Applied Physics Letters, 77(19):2958–2960, nov 2000.
  * [4] Yan Li, Anirban Bhattacharyya, Christos Thomidis, Theodore D. Moustakas, and Roberto Paiella.  Ultrafast all-optical switching with low saturation energy via intersubband transitions in GaN/AlN quantum-well waveguides.  Optics Express, 15(26):17922, dec 2007.
  * [5] Wei Li, Bigeng Chen, Chao Meng, Wei Fang, Yao Xiao, Xiyuan Li, Zhifang Hu, Yingxin Xu, Limin Tong, Hongqing Wang, Weitao Liu, Jiming Bao, and Y. Ron Shen.  Ultrafast all-optical graphene modulator.  Nano Letters, 14(2):955–959, feb 2014.
  * [6] Masaaki Ono, Masanori Hata, Masato Tsunekawa, Kengo Nozaki, Hisashi Sumikura, Hisashi Chiba, and Masaya Notomi.  Ultrafast and energy-efficient all-optical switching with graphene-loaded deep-subwavelength plasmonic waveguides.  Nature Photonics, 2020.
  * [7] Xuekai Ma, Bernd Berger, Marc Aßmann, Rodislav Driben, Torsten Meier, Christian Schneider, Sven Höfling, and Stefan Schumacher.  Realization of all-optical vortex switching in exciton-polariton condensates.  Nature Communications, 11(1):1–7, dec 2020.
  * [8] C. Leyder, T. C.H. Liew, A. V. Kavokin, I. A. Shelykh, M. Romanelli, J. Ph Karr, E. Giacobino, and A. Bramati.  Interference of coherent polariton beams in microcavities: Polarization-controlled optical gates.  Physical Review Letters, 99(19):196402, nov 2007.
  * [9] D. Ballarini, M. De Giorgi, E. Cancellieri, R. Houdré, E. Giacobino, R. Cingolani, A. Bramati, G. Gigli, and D. Sanvitto.  All-optical polariton transistor.  Nature Communications, 4(1):1–8, apr 2013.
  * [10] Félix Marsault, Hai Son Nguyen, Dimitrii Tanese, Aristide Lemaître, Elisabeth Galopin, Isabelle Sagnes, Alberto Amo, and Jacqueline Bloch.  Realization of an all optical exciton-polariton router.  Applied Physics Letters, 107(20):201115, nov 2015.
  * [11] Vincenzo Ardizzone, Przemyslaw Lewandowski, M. H. Luk, Y. C. Tse, N. H. Kwong, Andreas Lücke, Marco Abbarchi, Emmanuel Baudin, Elisabeth Galopin, Jacqueline Bloch, Aristide Lemaitre, P. T. Leung, Philippe Roussignol, Rolf Binder, Jerome Tignon, and Stefan Schumacher.  Formation and control of turing patterns in a coherent quantum fluid.  Scientific Reports, 3(1), oct 2013.
  * [12] I Gnusov, H Sigurdsson, J D Töpfer, S Baryshev, and P G Lagoudakis.  All-optical exciton-polariton polarization rotator.
  * [13] P. M. Walker, L. Tinkler, M. Durska, D. M. Whittaker, I. J. Luxmoore, B. Royall, D. N. Krizhanovskii, M. S. Skolnick, I. Farrer, and D. A. Ritchie.  Exciton polaritons in semiconductor waveguides.  Applied Physics Letters, 102(1):012109, jan 2013.
  * [14] Itamar Rosenberg, Yotam Mazuz-Harpaz, Ronen Rapaport, Kenneth West, and Loren Pfeiffer.  Electrically controlled mutual interactions of flying waveguide dipolaritons.  Physical Review B, 93(19):195151, may 2016.
  * [15] Itamar Rosenberg, Dror Liran, Yotam Mazuz-Harpaz, Kenneth West, Loren Pfeiffer, and Ronen Rapaport.  Strongly interacting dipolar-polaritons.  Science Advances, 4(10):eaat8880, oct 2018.
  * [16] D G Suárez-Forero, F Riminucci, V Ardizzone, N Karpowicz, E Maggiolini, G Macorini, G Lerario, F Todisco, M De Giorgi, L Dominici, D Ballarini, G Gigli, A S Lanotte, K West, K Baldwin, L Pfeiffer, and D Sanvitto.  Enhancement of Parametric Effects in Polariton Waveguides Induced by Dipolar Interactions.  Physical Review Letters, 126, 2021.
  * [17] O. Jamadi, F. Reveret, P. Disseix, F. Medard, J. Leymarie, A. Moreau, D. Solnyshkov, C. Deparis, M. Leroux, E. Cambril, S. Bouchoule, J. Zuniga-Perez, and G. Malpuech.  Edge-emitting polariton laser and amplifier based on a ZnO waveguide.  Light: Science and Applications, 7(1):2047–7538, dec 2018.
  * [18] D. G. Suárez-Forero, F. Riminucci, V. Ardizzone, M. de Giorgi, L. Dominici, F. Todisco, G. Lerario, L. N. Pfeiffer, G. Gigli, D. Ballarini, and D. Sanvitto.  Electrically controlled waveguide polariton laser.  Optica, 7(11):1579, mar 2020.
  * [19] Alex Hayat, Christoph Lange, Lee A. Rozema, Ardavan Darabi, Henry M. Van Driel, Aephraim M. Steinberg, Bryan Nelsen, David W. Snoke, Loren N. Pfeiffer, and Kenneth W. West.  Dynamic stark effect in strongly coupled microcavity exciton polaritons.  Physical Review Letters, 109(3):033605, jul 2012.
  * [20] Christoph Lange, Emiliano Cancellieri, Dmitry Panna, David M. Whittaker, Mark Steger, David W. Snoke, Loren N. Pfeiffer, Kenneth W. West, and Alex Hayat.  Ultrafast control of strong light-matter coupling.  New Journal of Physics, 20(1):013032, jan 2018.
  * [21] M. De Giorgi, D. Ballarini, E. Cancellieri, F. M. Marchetti, M. H. Szymanska, C. Tejedor, R. Cingolani, E. Giacobino, A. Bramati, G. Gigli, and D. Sanvitto.  Control and ultrafast dynamics of a two-fluid polariton switch.  Physical Review Letters, 109(26), dec 2012.
  * [22] P. M. Walker, L. Tinkler, B. Royall, D. V. Skryabin, I. Farrer, D. A. Ritchie, M. S. Skolnick, and D. N. Krizhanovskii.  Dark solitons in high velocity waveguide polariton fluids.  Phys. Rev. Lett., 119:097403, Aug 2017.
  * [23] Hamis Gargoubi, Thierry Guillet, Sihem Jaziri, Jalloul Balti, and Brahim Guizal.  Polariton condensation threshold investigation through the numerical resolution of the generalized Gross-Pitaevskii equation.  Physical Review E, 94(4):043310, oct 2016.
  * [24] Viktoriia Rutckaia and Joerg Schilling.  Ultrafast low-energy all-optical switching, 2020.
  * [25] Alexey Kavokin, Jeremy J. Baumberg, Guillaume Malpuech, and Fabrice P. Laussy.  Microcavities.  Oxford Science, second edition, 2008.
  * [26] E. Estrecho, T. Gao, N. Bobrovska, D. Comber-Todd, M. D. Fraser, M. Steger, K. West, L. N. Pfeiffer, J. Levinsen, M. M. Parish, T. C.H. Liew, M. Matuszewski, D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya.  Direct measurement of polariton-polariton interaction strength in the Thomas-Fermi regime of exciton-polariton condensation.  Physical Review B, 100(3):035306, jul 2019.
  * [27] F. Tassone and Y. Yamamoto.  Exciton-exciton scattering dynamics in a semiconductor microcavity and stimulated scattering into polaritons.  Physical Review B - Condensed Matter and Materials Physics, 59(16):10830–10842, apr 1999.
  * [28] Iacopo Carusotto and Cristiano Ciuti.  Quantum fluids of light.  Reviews of Modern Physics, 85(1):299–366, feb 2013.

Supporting Information

##  I Experimental setup

Two different experimental configurations are used along this work, depending
on the studied effect. The sub-picosecond modulation of the energy dispersion
by the Optical Stark effect cannot be resolved with a streak camera, due to
its typical time range, which depends on the pump pulse temporal extension. In
our case, we use a laser with a pulse duration of 100 fs, a time window below
the camera’s resolution ($>$250 fs). Therefore, to characterize this effect in
our system we used an ultrafast pulse, generated with a super-continuum (SC)
source to probe the polariton dispersion. This configuration is illustrated in
Fig. S1a. A delay line was implemented for the fs pulse, in order to vary its
arrival time respect to the SC pulse. While the SC pulse was coupled
resonantly into the TE mode, the fs pulse is injected resonantly into the TM
mode. The transmission of the system when both lasers are injected is
displayed in Fig. S1b. The effect was then characterized by registering the
transmission of the SC pulse as the arrival time of the fs pulse was
systematically modified.

Figure S1: Characterization of the optical Stark effect. a) scheme of the
pumping setup. The fs pulse feeds a supercontinuum source, generating another
pulse with a broad spectrum. The desired energy is selected afterwards with a
bandpass filter, and then injected resonantly into the system. b) System
transmission when both pulses are injected resonantly. The control pulse is
resonant with the $TM_{o}$ mode, while the SC laser is resonant with the
$TE_{o}$ mode.

The picosecond regime, product of the creation of a charge-carriers reservoir
via two photons absorption is displayed in Fig. S2. We use a tunable
Ti:sapphire CW laser. For the different configurations, the laser is set at
variable energies in the range of the NIR, from 1525.14 meV to 1521.55 meV.
Since the aim was to pump resonantly the $TE_{o}$ polariton branch, after the
linear polarizer the HWP rotates the polarization to match such guided mode.
The pulsed laser is tuned as far to the exciton as the numerical aperture of
the pumping setup (including the patterned gratings) allows. It is a
Ti:sapphire with 80 MHz of repetition rate and 100 fs of pulse duration. It is
tuned (centered) at around 1500 meV. Both lasers have independent optical
lines, in order to have fine control on the injection angle and polarization
of each one. The sample is kept at a cryogenic temperature of 4K during the
full experiment, and to avoid any thermal effect (specially due to the CW
pump), we use a shutter that opens 10 ms every 3 seconds and we trigger the
streak camera to open in synchrony.
After the injection, propagation and extraction of the optical signals, the
real space plane is reconstructed and the reflection of the input lasers is
suppressed with a slit, so only the light from the outcoupling grating is
collected in detection. After the real space filter, the Fourier plane is
focused on a monochromator+streak camera system. The signal from the pulsed
laser is focused on the detection system to synchronize the time of arrival of
the pulse, after which it is filtered in energy (by the monochromator),
polarization (with the linear polarizer) and angle (since CW and pulsed
signals are resonant with the guided modes at different in-plane momenta, they
are injected and extracted with different angle). The CW signal is then
analyzed in the streak camera with a temporal resolution that varies from 0.25
ps to 2.3 ps, depending on the selected time range.

Figure S2: Depiction of the experimental setup. LP stands for ”Linear
Polarizer” and ”HWP” stands for Half-wave plate.

##  II Pulse energy above the bandgap

As mentioned in the main text, the switch effect has its origin in an
instantaneous modification of the polaritonic dispersion, which takes the
coupled CW signal out of resonance, reducing its transmission. Since such
modification is due to the creation of a charge carrier reservoir, one might
think that it is more efficient to generate directly the reservoir by pumping
with a pulse whose energy overcomes the QW bandgap. This indeed creates a
larger reservoir, however, it is not the best option, and the reason can be
understood by looking at Fig. S3.

Figure S3: Time resolved output signal for the case in which the fs pulsed
laser is tuned at an energy above the QW bandgap and with a power of $100$ mW.
The PL emission composes a great part of the detected light, adding undesired
noise to the CW signal. The inset corresponds to the energy selection
indicated by the red dashed lines, which is the CW signal affected by the
light generated by PL emission.

The fs pulse generates a PL emission that overlaps with the CW signal, adding
considering noise to the output. There is still a reminiscence of the
switching effect (as displayed in the inset), but it is critically affected by
the added optical noise. Although this configuration does not represent the
best option for the implementation of the ultrafast switch of the CW, it is in
good agreement with a transmission reduction originated by an excitonic shift
and consequent polariton dispersion renormalization.

##  III Formation and effects of Dark Excitons

As it has been previously discussed, the origin of the switching effect is the
creation of a charge reservoir by two-photons absorption of the pulse after
injected inside the structure. Such reservoir is composed, among others, by
dark excitons, whose optical recombination is prevented by the spin selection
rules. The effects of such particles in the optical properties of the material
is rather critical, though. As mentioned in the main text, a large population
of dark excitons could carry a renormalization of the polaritonic dispersion,
very hard to quantify due to the lack of emission from these particles. They
are characterized by very long lifetimes [25]. In our case, after a critical
power of the pulsed laser, such states become relevant, and affect the
transmission in time ranges longer than the delay between pulses. The initial
transmission of the CW laser is then affected by the presence of dark-excitons
accumulated in consecutive pulses. It is worth to mention that the “switching”
condition is not affected whatsoever by this effect, since, independent from
the initial transmission, the CW laser is reduced up to half its initial
intensity. The reservoir of dark excitons implies then a renormalization of
the CW flux transmission, changing, precisely, the CW flux inital condition.
The details of such behavior are displayed in Fig. S4.

Figure S4: Time-resolved unnormalized transmission of CW laser for different
powers of the pulsed laser. The reduction of the initial transmitivity for
increasing powers witnesses the creation of a reservoir of long-lived dark
excitons.

Notice how the initial condition of the CW laser is modified by the pumping
power but the switching events still takes place. This suggests that the
reservoir does not totally empty between consecutive pulses, indicating that
the created dark excitons have lifetimes at least comparable with the delay
between laser pulses (12 ns).

##  IV Theoretical model

For the theoretical modelization of the polariton scattering by the excitonic
reservoir we implement a Gross-Pitaevskii model for a density of polaritons in
presence of a reservoir of long-lived, uncoupled excitons. The system of
equations describing the system is as follows[23]:

$\displaystyle i\hbar\frac{\partial\Psi(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x^{2}}+g_{XX}\cdot
C_{X}\left[|\Psi(x,t)|^{2}+n_{R}(x,t)\right]-i\left(C_{X}\gamma_{e}+\left[1-C_{X}\right]\gamma_{p}\right)\right]\Psi(x,t)+P_{CW}(x,t)$
(S1a) $\displaystyle\hbar\frac{\partial n_{R}(x,t)}{\partial t}$
$\displaystyle=$ $\displaystyle-2\gamma_{R}n_{R}(x,t)+P_{fs}(x,t)$ (S1b)

Where $\Psi(x,t)$ is the wavefunction of the polariton injected resonantly.
The horizontal (y) dependence can be neglected due to the additional
confinement provided by the microwire, that restrains the mode to 1D. The
polariton mass $m$ has tricky details. If it is calculated as
$\frac{\partial^{2}E}{\partial k^{2}}$ the parameter will take negative
values. We decided to employ a common value for a DBR microcavity: $3\times
10^{-5}m_{e}$, where $m_{e}$ is the free electron mass. $g_{XX}$ is the
strength of the exchange exciton-exciton interaction, whose value has been
widely debated and very different values have been reported in different
experiments[26]. In this work we took as reference the commonly accepted value
of $6$ $\mu eV\cdot\mu m$[27]. The excitonic fraction $C_{X}$ determines how
close in energy the injected CW polaritons are from the excitonic resonance,
and hence, how strong the interaction is [25]. The decay rates of the
excitonic ($\gamma_{e}$) and photonic ($\gamma_{p}$) modes are determined
experimentally by pumping the WG system out of resonance while systematically
vary the distance from the pumping spot to the outcoupling grating. The
reduction of the PL intensity with the pumping distance, allows to determine
the decay rate of each polariton component. $\gamma_{R}$, the decay rate of
the exciton reservoir is much more challenging to determine, given the dark
nature of a large portion of the composing particles, therefore, this
parameter is used as a fitting degree of freedom. In equation S1b $n_{R}(x,t)$
is the reservoir particles density. The pumping terms $P_{CW}(x,t)$ and
$P_{fs}(x,t)$ depend on the details of each laser. In particular they have the
form:

$\displaystyle P_{CW}(x,t)$ $\displaystyle=$ $\displaystyle
p_{CW}e^{-\left(x/2\sigma_{CW}\right)^{2}}e^{i\Delta_{\omega}t}$ (S2a)
$\displaystyle P_{fs}(x,t)$ $\displaystyle=$ $\displaystyle
N(1-e^{-\epsilon\cdot
p_{fs}^{2}})e^{-(x/2\sigma_{fs})^{2}}e^{-\left(\frac{t-t_{o}}{\tau}\right)^{2}}$
(S2b)

In the rotating frame, $\Delta_{\omega}$ represents the difference between the
exction and the photonic mode energies at the resonant in-plane momentum. As
it can be observed, the temporal dependence of each expression varies due to
the nature of the pump: the CW case is represented by an harmonic function,
while the pulsed laser has a temporal gaussian profile whose amplitude depends
on the pulse duration. Consider that this time must be much longer than the
pulse temporal width, because it takes into account the formation of the
reservoir. $\sigma_{CW}$ and $\sigma_{fs}$ depend on the spatial profile of
the CW and pulsed lasers, respectively. Both values can be easily determined
from an imaging of the real space plane. The pump powers fix the values
$p_{CW}$ and $p_{fs}$, however, their determination is not straightforward for
the following reasons: on the one hand the in-coupling and out-coupling
efficiency of the etched gratings cannot be precisely measured; it is not
correct to assume that the coupling efficiency is equal to the extraction one,
since the spatial profiles in each case are fundamentally different (gaussian
for the input and exponential for the output). On the other hand, the exact
amount of the pulsed energy that becomes dark excitons is difficult to obtain
due to the lack of an optical signal. The pump parameters $p_{fs}$ and
$p_{CW}$ are then treated as fitting parameters. As mentioned in the main
text, the dynamics has a negligible dependence on $p_{CW}$, making this a
minor issue. Notice the particular dependence of $P_{fs}$ with $p_{fs}$; this
expression takes into account the saturability of the bath of excitons. Pauli
exclusion principle prevents the formation of two excitons separated by a
distance shorter than the Bohr radius[28], precluding the formation of a
reservoir with a population higher than $N$. Since the process by which the
bath of charge carriers is created is two-photons absorption, the process have
a quadratic dependence with $p_{fs}$ and it has an efficiency determined by
$\epsilon$. The values of the parameters employed for the numerical simulation
are summarized in Table 1.

$m$ | $g_{XX}$ | $C_{X}$ | $\gamma_{e}$ | $\gamma_{p}$ | $\gamma_{R}$ | $\sigma_{CW}$ | $\sigma{fs}$ | $p_{CW}$ | $N$ | $\epsilon$ | $\tau$
---|---|---|---|---|---|---|---|---|---|---|---
$3\cdot 10^{-5}m_{e}$ |  $6$ $\mu eV\cdot\mu m$ | $0.69$ |  $223$ $\mu eV$ |  $30$ $\mu eV$ |  $1.8$ $\mu eV$ |  $35$ $\mu m$ |  $35$ $\mu m$ |  $0.3$ $eV\cdot\mu m^{-1/2}$ |  $1.8$ $eV\cdot\mu m^{-1}$ | $7\cdot 10^{-4}mW^{-2}$ |  $80$ ps
Table 1: Parameters employed in the numerical simulations.


# Diagrammatic Monte Carlo study of the Fermi polaron in two dimensions

Jonas Vlietinck  Department of Physics and Astronomy, Ghent University,
Proeftuinstraat 86, 9000 Gent, Belgium  Jan Ryckebusch  Department of Physics
and Astronomy, Ghent University, Proeftuinstraat 86, 9000 Gent, Belgium  Kris
Van Houcke  Department of Physics and Astronomy, Ghent University,
Proeftuinstraat 86, 9000 Gent, Belgium  Laboratoire de Physique Statistique,
Ecole Normale Supérieure, UPMC, Université Paris Diderot, CNRS, 24 rue
Lhomond, 75231 Paris Cedex 05, France

###### Abstract

We study the properties of the two-dimensional Fermi polaron model in which an
impurity attractively interacts with a Fermi sea of particles in the zero-
range limit. We use a diagrammatic Monte Carlo (DiagMC) method which allows us
to sample a Feynman diagrammatic series to very high order. The convergence
properties of the series and the role of multiple particle-hole excitations
are discussed. We study the polaron and molecule energy as a function of the
coupling strength, revealing a transition from a polaron to a molecule in the
ground state. We find a value for the critical interaction strength which
complies with the experimentally measured one and predictions from variational
methods. For all considered interaction strengths, the polaron $Z$ factor from
the full diagrammatic series almost coincides with the one-particle-hole
result. We also formally link the DiagMC and the variational approaches for
the polaron problem at hand.

###### pacs:

05.30.Fk, 03.75.Ss, 02.70.Ss

##  I Introduction

Experiments with ultracold gases are a powerful tool to investigate the
(thermo)dynamics of quantum many-body systems under controlled circumstances.
With Feshbach resonances feshbach , for example, one has the ability to tune
the interaction strength. Optical potentials potential  can be exploited to
modify the dimensionality of the studied systems. The properties of a single
impurity that interacts strongly with a background gas, for example, can be
addressed with ultracold atoms.

The so-called Fermi polaron problem refers to a single spin-down impurity that
is coupled to a non-interacting spin-up Fermi sea (FS). This problem
corresponds to the extreme limit of spin imbalance in a two-component Fermi
gas Partridge06 ; Shin08 ; Nascim09  and has implications on the phase diagram
of the strongly spin-polarized Fermi gas Pilati08 ; pietro ; Bertaina . At
weak attraction, one expects a “polaron” state Chevy06 , in which the impurity
is dressed with density fluctuations of the spin-up Fermi gas. Recent
experiments have observed indications of a transition from this polaronic
state to a molecular state (a two-body bound state of the impurity and an atom
of the sea) upon increasing the attraction strength in three dimensions (3D)
polaronMIT  and in two dimensions (2D) kos . Experimentally, the 2D regime can
be accomplished by means of a transverse trapping potential
$V(z)=\frac{1}{2}m\omega_{z}^{2}z^{2}$ (here, $\omega_{z}$ is the frequency
and $z$ is the transverse direction) that fulfills the condition
$k_{B}T\ll\epsilon_{F}\ll\hbar\omega_{z}$ ($T$ is the temperature and
$\epsilon_{F}$ is the Fermi energy of the FS). When excitations in the $z$
dimension are possible, one reaches the so-called quasi-2D regime jesper2 ;
dyke . The purely 2D limit is reached for
$\epsilon_{F}/\hbar\omega_{z}\rightarrow 0$ and will be the subject of this
paper.

The existence of a polaron-molecule transition in 3D has been predicted with
the aid of the diagrammatic Monte Carlo (DiagMC) method polaron1 ; polaron2 ;
polaron  and of variational methods Chevy06 ; Punk09 ; Combescot09 ; Mora .
For the latter, the maximum number of particle-hole (p-h) excitations of the
FS is limited to one or two Chevy06 ; Punk09 ; Combescot09 ; Mora . One might
naively expect that the role of quantum fluctuations increases in importance
with decreasing dimensionality and that high-order p-h excitations could
become more important in one and two dimensions. For the one-dimensional (1D)
Fermi polaron the known analytical solution displays no polaron-molecule
transition oneDexact . Like for the 3D polaron, the approximate method in
which the truncated Hilbert space contains one p-h and two p-h excitations of
the FS gives results for the 1D polaron approaching the exact solution oneD1 ;
oneD2 . In 2D, the Fermi polaron properties have been studied with variational
wave functions Zollner11 ; Parish11 ; Levinsen13 . To observe a polaron-
molecule transition in 2D it is crucial to include particle-hole excitations
in both the polaron and molecule wave functions Parish11 . In the limit of
weak interactions, the 1p-h and 2p-h variational Ansätze for the polaron
branch provide similar results. Surprisingly, this is also the situation for
strong correlations Levinsen13 .

In this work we focus on the 2D Fermi polaron for attractive interactions and
study the role of multiple particle-hole (mp-h) excitations for the ground-
state properties of the system. The quasiparticle properties of the polaron
are computed with the DiagMC method. This technique evaluates stochastically
to high order a series of Feynman diagrams for the one-particle and two-
particle self-energies. For the details of the DiagMC method and the adopted
method for determining the ground-state energies from the computed self-
energies, we refer to Refs. polaron ; polaron2 . In this work we present
DiagMC predictions for the interaction-strength dependence of the polaronic
and molecular ground-state properties in 2D. We first briefly discuss the
model and the diagrammatic method. We then discuss the results of the
simulations, with particular emphasis on the role of the mp-h excitations. We
also discuss how variational results for the polaron problem can be obtained
within the DiagMC formalism.

##  II Formalism

We consider a two-component Fermi gas confined to 2D at temperature $T=0$.
Even though we will consider the zero-range interaction in continuous space,
we start from a lattice model to avoid ultraviolet divergences from the onset.
The corresponding Hamiltonian reads

$\hat{H}=\sum_{\mathbf{k}\in\mathcal{B},\sigma=\uparrow\downarrow}\epsilon_{\mathbf{k}\sigma}~{}\hat{c}^{\dagger}_{\mathbf{k}\sigma}\hat{c}^{\phantom{\dagger}}_{\mathbf{k}\sigma}\\\
+g_{0}\sum_{\mathbf{r}}b^{2}~{}\hat{\Psi}^{\dagger}_{\uparrow}(\mathbf{r})\hat{\Psi}^{\dagger}_{\downarrow}(\mathbf{r})\hat{\Psi}^{\phantom{\dagger}}_{\downarrow}(\mathbf{r})\hat{\Psi}^{\phantom{\dagger}}_{\uparrow}(\mathbf{r})\;,$
(1)

with $\hat{\Psi}^{\phantom{\dagger}}_{\sigma}(\mathbf{r})$ and
$\hat{c}^{\phantom{\dagger}}_{\mathbf{k},\sigma}$ being the operators for
annihilating a spin-$\sigma$ fermion with mass $m_{\sigma}$ and dispersion
$\epsilon_{\mathbf{k}\sigma}=k^{2}/2m_{\sigma}$ in position and momentum
space. The components of the position vector $\mathbf{r}$ are integer
multiples of the finite lattice spacing $b$. Further, $g_{0}$ is the bare
interaction strength. The wave vectors $\mathbf{k}$ are in the first Brillouin
zone $\mathcal{B}=]-\pi/b,\pi/b]$. The continuum limit is reached for
$b\rightarrow 0$. We adopt the convention $\hbar=1$ and consider the mass-
balanced case $m_{\uparrow}=m_{\downarrow}=m$. lattice with spacing $b$. We
make use of the $T$ matrix landau  for a single spin-$\uparrow$ and
spin-$\downarrow$ fermion in vacuum,

$-\frac{1}{g_{0}}=\frac{1}{\mathcal{V}}\sum_{\mathbf{k}\in\mathcal{B}}\frac{1}{\varepsilon_{B}+\epsilon_{\mathbf{k}\uparrow}+\epsilon_{\mathbf{k}\downarrow}}\;,$
(2)

where $\mathcal{V}$ is the area of the system and $\varepsilon_{B}$ is the
two-body binding energy [which depends on $m$, $g_{0}$, and $b$ and
$\varepsilon_{B}(m,g,b)>0$] of a weakly bound state. Such a state always
exists for an attractive interaction in 2D. With the above relation we
eliminate the bare interaction strength $g_{0}$ in favor of the quantity
$\varepsilon_{B}$. Moreover, the diagrammatic approach allows us to take the
continuum limit $b\to 0$ and $g_{0}\to 0^{-}$ while keeping $\varepsilon_{B}$
fixed. Summing all ladder diagrams gives a partially dressed interaction
vertex $\Gamma^{0}$:

$\mbox{\raisebox{0.0pt}{\psfig{scale=0.3,clip={true}}}~{}}\;,$ (3)

where the dot represents the bare interaction vertex $g_{0}$ and the lines
represent bare-particle propagators for the spin-down impurity (dashed lines)
and the spin-up Fermi sea (solid lines). In momentum-imaginary frequency this
graphical representation corresponds to

$[\Gamma^{0}(p,i\Omega)]^{-1}=g_{0}^{-1}-\Pi^{0}(p,i\Omega)\;,$ (4)

with

$\displaystyle\Pi^{0}(p,i\Omega)=\frac{1}{\mathcal{V}}\sum_{\mathbf{k}\in\mathcal{B}}\frac{H(|\frac{\mathbf{p}}{2}+\mathbf{k}|-k_{F})}{i\Omega-\epsilon_{\frac{\mathbf{p}}{2}-\mathbf{k}\downarrow}-\epsilon_{\frac{\mathbf{p}}{2}+\mathbf{k}\uparrow}+\mu+\varepsilon_{F}}\;,$
(5)

with $H(x)$ being the Heaviside step function and $\mu<0$ being a free
parameter representing an energy offset of the impurity dispersion. Further,
$k_{F}$ and $\varepsilon_{F}=\frac{k_{F}^{2}}{2m}$ are the Fermi momentum and
the Fermi energy of the spin-up sea. The combination of Eqs. (2) and (4) gives

$\displaystyle\frac{1}{\Gamma^{0}(p,i\Omega)}=$ $\displaystyle-$
$\displaystyle\frac{1}{\mathcal{V}}\sum_{\mathbf{k}\in\mathcal{B}}\left[\frac{1}{\varepsilon_{B}+\epsilon_{\mathbf{k}\uparrow}+\epsilon_{\mathbf{k}\downarrow}}\right.$
(6) $\displaystyle+$
$\displaystyle\left.\frac{H(|\frac{\mathbf{p}}{2}+\mathbf{k}|-k_{F})}{i\Omega-\epsilon_{\frac{\mathbf{p}}{2}-\mathbf{k}\downarrow}-\epsilon_{\frac{\mathbf{p}}{2}+\mathbf{k}\uparrow}+\mu+\varepsilon_{F}}\right]\;.$

The relevant parameter that characterizes the interaction in Equation. (6) is
$\epsilon_{B}$. Eq. (6) is well defined in the thermodynamic and $b\to 0$
limits. One finds

$\displaystyle\frac{1}{\Gamma^{0}(p,i\Omega)}=\frac{m}{4\pi}~{}{\rm
ln}\bigg{[}\frac{2\varepsilon_{B}}{-z+\sqrt{(z-\epsilon_{\mathbf{p}})^{2}-4\varepsilon_{F}\epsilon_{\mathbf{p}}}}\bigg{]}\;,$
(7)

with $z\equiv i\Omega+\mu-\varepsilon_{F}$. In deriving the above expression
for $\Gamma^{0}(p,i\Omega)$ we have taken $\mu<-\varepsilon_{F}$. Since
Feynman diagrams for the self-energy will be evaluated in the momentum-
imaginary-time representation $(p,\tau)$, we need to evaluate the Fourier
transform

$\Gamma^{0}(p,\tau)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}~{}d\Omega~{}e^{-i\Omega\tau}~{}\Gamma^{0}(p,i\Omega)\;.$
(8)

In order to determine the leading behavior of $\Gamma^{0}(p,\tau)$ for small
$\tau$, we introduce the vertex function $\tilde{\Gamma}^{0}$, which differs
from $\Gamma^{0}$ by ignoring the Fermi surface when integrating out the
internal momenta. This amounts to ignoring the Heaviside function in Eq. (5).
We obtain

$\displaystyle\frac{1}{\tilde{\Gamma}^{0}(p,i\Omega)}=\frac{m}{4\pi}~{}{\rm
ln}\bigg{[}-\frac{\varepsilon_{B}}{i\Omega+\mu+\varepsilon_{F}-\frac{\epsilon_{\mathbf{p}}}{2}}\bigg{]}\;.$
(9)

In the $(p,i\Omega)$-representation,

$\displaystyle\frac{1}{\Gamma^{0}}-\frac{1}{\tilde{\Gamma}^{0}}=\frac{m}{4\pi}{\rm
ln}\bigg{[}\frac{-2(z+2\varepsilon_{F})+\epsilon_{\mathbf{p}}}{-z+\sqrt{(z-\epsilon_{\mathbf{p}})^{2}-4\varepsilon_{F}\epsilon_{\mathbf{p}}}}\bigg{]}\;.$
(10)

The $(p,\tau)$ representation of $\tilde{\Gamma}^{0}$ is

$\displaystyle\tilde{\Gamma}^{0}(p,\tau)$ $\displaystyle=$
$\displaystyle-\frac{4\pi\varepsilon_{B}}{m}e^{-(\frac{\epsilon_{\mathbf{p}}}{2}-\varepsilon_{F}-\mu)\tau}~{}\bigg{[}\int_{0}^{+\infty}dx\frac{e^{-x\varepsilon_{B}\tau}}{\pi^{2}+{\rm
ln}^{2}(x)}$ (11)
$\displaystyle+~{}e^{\varepsilon_{B}\tau}H(\frac{\epsilon_{\mathbf{p}}}{2}-\varepsilon_{F}-\varepsilon_{B}-\mu)\bigg{]}~{}H(\tau)$

for $\mu<-\varepsilon_{F}-\varepsilon_{B}$, ensuring that only $\tau>0$
contributes for all momenta $p$. The integral in Eq. (11) can be computed
numerically, but converges poorly for $\tau\rightarrow 0^{+}$. Under those
conditions we make use of the asymptotic behavior:

$\int_{0}^{+\infty}dx\frac{e^{-x\varepsilon_{B}\tau}}{\pi^{2}+{\rm
ln}^{2}(x)}\underset{\tau\rightarrow
0}{\sim}\frac{1}{\varepsilon_{B}\tau}\frac{1}{{\rm
ln}^{2}(\varepsilon_{B}\tau)}\;.$ (12)

To obtain $\Gamma^{0}(p,\tau)$ we computed numerically the following Fourier
transform:

$\displaystyle\Gamma^{0}(p,\tau)-\tilde{\Gamma}^{0}(p,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{+\infty}d\Omega~{}e^{-i\Omega\tau}$
(13) $\displaystyle\times$
$\displaystyle\left[\Gamma^{0}(p,i\Omega)-\tilde{\Gamma}^{0}(p,i\Omega)\right]\;.$

The left-hand side of Eq. (13) can be computed more easily than
$\tilde{\Gamma}^{0}(p,\tau)$ as it contains no singularities. Next, the
function $\Gamma^{0}(p,\tau)$ is obtained as
$\tilde{\Gamma}^{0}(p,\tau)+\left[\Gamma^{0}(p,\tau)-\tilde{\Gamma}^{0}(p,\tau)\right]$.
Although the functions $\tilde{\Gamma}^{0}(p,\tau)$ and $\Gamma^{0}(p,\tau)$
are extremely sharp and divergent for $\tau\to 0$, they are integrable.
Special care should be taken when using these functions in the Monte Carlo
code. It is important to correctly sample very short times, and one needs to
make sure there is no loss of accuracy when keeping track of imaginary time
differences of the $\Gamma^{0}$ lines in the diagrams.

Figure 1:  Graphical representation of the Dyson equation. The free (dressed)
one-body impurity propagator is denoted by $G_{\downarrow}^{0}$
($G_{\downarrow}$). $\Sigma$ and $\Pi$ are the one-body and two-body self-
energies, respectively. $\Gamma$ is the fully dressed interaction, whereas
$\Gamma^{0}$ is the partially dressed interaction as shown in Eq. (3).

Just like for the 3D polaron problem polaron1 ; polaron2 ; polaron , we
consider a diagrammatic series for the self-energy built from the free one-
body propagators for the impurity and the spin-up Fermi sea and from the
renormalized interaction $\Gamma^{0}$. We refer to this series as the _bare
series_ , which we evaluate with the DiagMC method. The diagram topologies in
2D and 3D are exactly the same. The major differences between the
diagrammatic-series evaluations in 2D and 3D are the renormalized interaction
$\Gamma^{0}(p,\tau)$ and the phase-space volume elements. The one and two-body
self-energies are related to the one-particle propagator $G$ and the fully
dressed interaction $\Gamma$ by means of a Dyson equation, as illustrated
schematically in Fig. 1. From the poles of $G$ and $\Gamma$ we can extract the
polaron and the molecule energy, respectively. The fully dressed interaction
is closely related to the two-particle propagator polaron .

For the 3D Fermi polaron problem there are two dominant diagrams at each given
order that emerge next to many diagrams with a much smaller contribution
polaron . These dominant diagrams contribute almost equally but have opposite
sign. In 2D, however, the numerical calculations indicate that at a given
order the very same two diagrams dominate, but to a lesser extent; that is,
the nondominant diagrams have a larger weight in the final 2D result. By
_weight_ of a given diagram we mean the absolute value of its contribution to
the self-energy. We note that the sign of a single diagram at fixed internal
and external variables depends only on its topology and not on the values of
the internal and external variables. We stress that this is not true for a
Fermi system with two interacting components with finite density vanhoucke12 ;
vanhoucke13 . In 2D the total weight of a given order (i.e., the sum of the
absolute values of the contributions of diagrams) is distributed over more
diagrams than in 3D. Because the sign alternation occurs over a broader
distribution of the weights, we get more statistical noise in sampling the
self-energy in 2D compared to 3D. In 3D we can evaluate the diagrammatic
series for the one-body self-energy accurately up to order $12$, whereas in 2D
we can reach order $8$.

In principle, other choices for the propagators (“bare” versus “dressed”
propagators) are possible, and this was discussed in detail for the 3D Fermi
polaron in our previous paper polaron . Replacing the bare propagators by
dressed ones reduces the number of diagrams at each given order. One may
expect that this replacement could allow one to reach higher orders. For the
3D polaron, however, the most favorable conditions of cancellations between
the contributions from the various diagrams were met in the bare scheme
polaron . In the DiagMC framework a higher accuracy can be reached under
conditions of strong cancellations between the various contributions. From
numerical investigations with various propagators for the 2D Fermi polaron we
could draw similar conclusions as in the 3D studies. Accordingly, all
numerical results for the quasiparticle properties presented below are
obtained for a series expansion with bare propagators.

Figure 2: (Color online) Dependence of the polaron energy $E^{N}_{p}$ on the
cutoff diagram order $N$. $E_{p}$ is the value obtained after extrapolation to
$N\to+\infty$ (and resummation for $\eta=-0.25$). Results are shown for
$\eta=-0.25,\eta=0.5,\eta=1.5$. The lines represent an exponential fit.

To characterize the magnitude of the interaction strength we use the
dimensionless parameter
$\eta\equiv\text{ln}[k_{F}a_{2D}]=\text{ln}[2\varepsilon_{F}/\varepsilon_{B}]/2$.
Here, $a_{2D}>0$ is the 2D scattering length, related to the dimer binding
energy by $\varepsilon_{B}=1/(2m_{r}a_{2D}^{2})$ with
$m_{r}=m_{\uparrow}m_{\downarrow}/(m_{\uparrow}+m_{\downarrow})$ being the
reduced mass. The BCS regime corresponds to $\eta\gg 1$ while the Bose-
Einstein condensate (BEC) regime corresponds to $\eta\ll-1$. The system is
perturbative in the regimes $|\eta|\gg 1$, while the strongly correlated
regime corresponds to $|\eta|\lesssim 1$. bloom  In the weak-coupling regime
[small interaction strengths $g_{0}$ in the Hamiltonian of Eq. (1) or large
positive $\eta$ in the zero-range limit], we find that the one-body and the
two-body self-energy $\Sigma$ and $\Pi$ converge absolutely as a function of
the maximum diagram order. This is demonstrated in Fig. 2 for $\eta=1.5$,
where the polaron energy $E^{N}_{p}$ converges exponentially as a function of
the cutoff diagram order $N$. Similar convergence is also found for the
molecule energy. Under conditions of convergence with diagram order,
extrapolation to order infinity can be done in a trivial way. Similar
convergence is also seen for $\eta=0.5$. In the strongly correlated regime the
series starts oscillating with order when $\eta\lesssim 0$, and the
oscillations get stronger the deeper we go into the BEC regime. The
oscillations in the extracted polaron energy are illustrated in Fig. 2 for
$\eta=-0.25$. To obtain meaningful results we rely on Abelian resummation
techniques polaron ; vanhoucke12 . We evaluate the series
$\sigma_{\epsilon}=\sum_{N}\sigma^{(N)}e^{-\epsilon\lambda_{N}}$, with
$\sigma^{N}$ being the one-body self-energy for diagram order $N$ and
$\lambda_{N}$ being a function that depends on the diagram order $N$. For each
$\epsilon$ the polaron energy $E_{p}$ is calculated from $\sigma_{\epsilon}$
and an extrapolation is done by taking the limit $\epsilon\rightarrow 0$. The
whole procedure is illustrated in Fig. 3. To estimate the systematic error of
the extrapolation procedure, different resummation functions $\lambda_{N}$ are
used. As becomes clear from Fig. 3 the whole resummation procedure is a stable
one and induces uncertainties on the extracted energies of the order of a few
percent. All the results of Fig. 4 are obtained with the Abelian resummation
technique. The stronger the coupling constant is the larger the size of the
error attributed to the resummation. An accurate extrapolation to infinite
diagram order could be achieved for all values of $\eta$.

Figure 3: (Color online) Abelian resummation of the bare series for the one-
body self-energy diagrams at $\eta=-0.25$. We evaluate
$\sigma_{\epsilon}=\sum_{N}\sigma^{(N)}e^{-\epsilon\lambda_{N}}$, with
$\sigma^{N}$ being the one-body self-energy for diagram order $N$. We use the
following functions $\lambda_{N}$ : (i) Gauss 1: $\lambda_{N}=(N-1)^{2}$ for
$N>1$ and $\lambda_{N}=0$ for $N=1$, (ii) Lindelöf 1: $\lambda_{N}$ $=$
$(N-1)~{}{\log}(N-1)$ for $N>2$ and $\lambda_{N}=0$ for $N\leq 2$, and (iii)
Gauss 2: $\lambda_{N}=(N-3)^{2}$ for $N>3$ and $\lambda_{N}=0$ for $N\leq 3$.
The polaron energy $E_{p}/\epsilon_{F}$ is extracted in the limit $\epsilon$
$=$ $0^{+}$ for various choices of $\lambda_{N}$. Figure 4: (Color online)
Polaron and molecule ground-state energies $E$ in units of the Fermi energy
$\varepsilon_{F}$ as a function of $\eta$. The momentum of the impurity is
equal to zero. Energies are shifted by the two-body binding energy
$\varepsilon_{B}/\varepsilon_{F}=2e^{-2\eta}$ to magnify the details. The
solid line is the DiagMC result for $N=1$. The symbols are the result of the
full DiagMC calculations (including diagrams up to order 8).

##  III Results and discussion

In Fig. 4, polaron and molecule energies are displayed for a wide range of the
parameter $\eta$. DiagMC results include all diagrams up to order 8 and
extrapolation to the infinite diagram order. In the region $\eta\lesssim 0$ a
small discrepancy (of the order of $0.1\%$ of the ground-state energy) is
found with the variational results Levinsen13  of Parish and Levinsen based on
the wave-function Ansatz up to 2p-h excitations. Clearly, a phase transition
appears at the critical value $\eta_{c}=-0.95\pm 0.15$. A variational result
which includes 2p-h excitations for the polaron and 1p-h excitations for the
molecule, gives $\eta=-0.97$. Levinsen13  Both mentioned calculations are in
agreement with the experimental result $\eta=-0.88(0.20)$ kos .

The DiagMC method allows one to include a large number of particle-hole
excitations that dress the impurity. Truncation of the Hilbert space to a
maximum number of p-h pairs can nonetheless be achieved within the DiagMC
approach. This allows one to arrive at the variational formulation. Previous
variational studies using a wave function Ansatz up to 1p-h or 2p-h
excitations showed that these truncations give remarkably accurate results
Combescot .

To understand why the truncation is possible within a Feynman diagrammatic
approach for the self-energy, we first remark that a variational approach is
easily established within a path-integral formalism. Path integrals with
continuous imaginary time, for example, are based on an expansion of the
evolution operator,

$\displaystyle e^{-\beta\hat{K}}$ $\displaystyle=$ $\displaystyle
e^{-\beta\hat{K}_{0}}\big{(}1-\int_{0}^{\beta}d\tau_{1}K_{1}(\tau_{1})$ (14)
$\displaystyle+$
$\displaystyle\int_{0}^{\beta}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}~{}\hat{K}_{1}(\tau_{1})\hat{K}_{1}(\tau_{2})-\ldots\big{)}\;,$

where $\hat{K}=\hat{H}-\mu\hat{N}=\hat{K}_{0}+\hat{K}_{1}-\mu\hat{N}$, with
$\left[\hat{K}_{0},\hat{K}_{1}\right]\neq 0$. The operator
$\hat{K}_{1}(\tau)=e^{\hat{K}_{0}\tau}\hat{K}_{1}e^{-\hat{K}_{0}\tau}$, which
defines the series expansion, is expressed in the interaction picture.
Further, $\beta=1/k_{B}T$, with $k_{B}$ being Boltzmann’s constant and $T$
being temperature, $\hat{H}$ is the Hamiltonian, $\hat{N}$ the number
operator, and $\mu$ is the chemical potential. The imaginary-time evolution
operator in Eq. (14) can be used as a ground-state projection operator: for
sufficiently long imaginary time $\beta$ the excited-state components of a
trial state are exponentially suppressed. One typically evaluates all the
terms in the expansion equation (14) in the eigenbasis of $\hat{K}_{0}$. This
procedure forms the basis of path-integral Monte Carlo simulation of lattice
models, where $\hat{K}_{1}$ is usually the kinetic energy term worm . A
discretized time version is used in path-integral Monte Carlo methods in
continuous space ceperley ; boninsegni . Either way, the contributions to the
path integral have the direct physical interpretation of a time history of the
many-particle system. At each instant of time, one can constrain the
accessible states of the Hilbert space, in line with what is done in a
variational approach. Within the standard Feynman diagrammatic formalism for
Green’s functions, however, this truncation of the Hilbert space is not easy
to accomplish for an arbitrary system, as one expands in powers of the two-
body interaction term of the Hamiltonian. This will be explained in the next
paragraph.

It turns out to be formally easier to start from finite $T$ and to take the
$\beta\to\infty$ limit in the end. For a many-fermion system, the finite-
temperature Green’s function in position and imaginary-time representation
$(\mathbf{x},\tau)$ is defined as

$G_{\alpha\sigma}(\mathbf{x},\tau)=-\frac{{\rm
Tr}[e^{-\beta\hat{K}}T_{\tau}\hat{\psi}^{\phantom{\dagger}}_{H\alpha}(\mathbf{x},\tau)\hat{\psi}^{\dagger}_{H\sigma}(\mathbf{x},0)]}{{\rm
Tr}[e^{-\beta\hat{K}}]}\;,$ (15)

with $\alpha$ and $\sigma$ denoting an appropriate set of quantum numbers
(such as spin) and $T_{\tau}$ being the time-ordering operator. The field
operator in the Heisenberg picture
$\hat{\psi}^{\phantom{\dagger}}_{H\alpha}(\mathbf{x},\tau)=e^{\hat{K}\tau}\hat{\psi}^{\phantom{\dagger}}_{\alpha}(\mathbf{x})e^{-\hat{K}\tau}$
annihilates a fermion in state $\alpha$ at position $\mathbf{x}$ and time
$\tau$. To arrive at the Feynman diagrammatic expansion, one makes a
perturbative expansion for the evolution operator $e^{-\beta\hat{K}}$ in both
the numerator _and_ the denominator of Eq. (15) (the finite $T$ ensures that
both exist). The expansion of the partition function $Z$ in the denominator
can be represented graphically by the series of all fully closed diagrams
(connected and disconnected). When $\beta$ approaches $+\infty$, the
denominator is proportional to $\langle\Psi_{0}|\Psi_{0}\rangle$
($|\Psi_{0}\rangle$ is the ground state of the interacting many-body system),
and the disconnected diagrams correspond to all possible vacuum fluctuations
of the system at hand. The expansion in the numerator factorizes into an
expansion of connected diagrams for $G_{\alpha\sigma}$ and disconnected
diagrams for $Z$. So the sum of disconnected diagrams drops out, as expected
for an intensive quantity like $G_{\alpha\sigma}(\mathbf{x},\tau)$. It is
exactly this factorization that prevents one from truncating the Hilbert space
at any instant of time in the evolution. In other words, variational
calculations based on Feynman diagrams for the self-energy are generally not
feasible.

In the polaron problem vacuum fluctuations are absent since $|\Psi_{0}\rangle$
corresponds to the spin-down vacuum and a non-interacting spin-up Fermi sea.
In other words, the vacuum cannot be polarized in the absence of an impurity.
As a consequence, we face a situation similar to the path integral with a
direct physical interpretation of the time history of the impurity. This
peculiar feature allows one to restrict the Hilbert space at each given time.
If we allow at most 1p-h excitations at each instant of time, only one diagram
survives: the lowest-order self-energy diagram built from $\Gamma^{0}$ and the
free spin-up single-particle propagator $G^{0}_{\uparrow}$. The equivalence
between this diagram and the 1p-h variational approach had already been
pointed out in Ref. oneD1 . An $n$p-h variational approach is achieved by
allowing at most $n$ backward spin-up lines at each step in the imaginary-time
evolution.

For large $\eta$ it is obvious from Fig. 4 that the polaron energy from the
full series expansion becomes equal to the 1p-h result. Even for stronger
interactions (smaller $\eta$) the first-order results remain close to the full
DiagMC one. Within the statistical accuracy of the numerical calculations,
convergence for the one-body self-energy is already reached after inclusion of
2p-h excitations. Indeed, for all values of $\eta$, we find agreement between
our 2p-h variational DiagMC approach and the full DiagMC approach within
statistical error bars. For the molecular branch, we retrieve the result for
the two-body self-energy from the full series expansion after including 1p-h
excitations. For the 3D Fermi-polaron a similar conclusion was drawn. Also in
3D, the first-order result is a very good approximation polaron . Going up to
2p-h pairs gives a perfect agreement with full DiagMC results. From the above
considerations it follows, however, that the diagrammatic truncations which
provide good results for the polaron problem may not be appropriate for the
more complex many-body problem with comparable densities for both components.

The quasiparticle residue or $Z$ factor of the polaron gives the overlap of
the noninteracting wave function and the fully interacting one. This overlap
is very small for a molecular ground state of the fully interacting system
polaron . The $Z$ factor as a function of $\eta$ is shown in Fig. 5. Note that
the polaron $Z$ factor does not vanish in the region $\eta\lesssim-1$ where
the ground state is a molecule. The $Z$ factor is, however, still meaningful
since the polaron is a well-defined (metastable) excited state of the 2D
system. Again, the first-order result gives a good approximation to the full
result. The measured $Z$ factor for the 3D situation has been reported in Ref.
pietro ; polaronMIT . The 2D experimental data are reported in Ref. kos , and
the $\eta$ dependence of the quasiparticle weight $Z$ is presented in
arbitrary units. We reproduce the observation that $Z$ strongly increases
between $\eta_{c}\lesssim\eta\lesssim 1$ and saturates to a certain value for
$\eta>1$.

##  IV Conclusion

Summarizing, we have developed a framework to study with the DiagMC method the
ground-state properties of the 2D Fermi polaron for attractive interactions.
We have shown that the framework allows one to select an arbitrary number of
$n$p-h excitations of the FS, thereby making a connection with typical
variational approaches which are confined to $n$=1 and $n$=2. We have studied
the quasiparticle properties of the ground state for a wide range of
interaction strengths. A phase transition between the polaron and molecule
states is found at interaction strengths compatible with experimental values
and with variational predictions. To a remarkable degree, it is observed that
for all interaction strengths the full DiagMC results (which include all
$n$p-h excitations) for the ground-state properties can be reasonably
approximated by n=1 truncations. In a n=2 truncation scheme the full result is
already reached within the error bars. This lends support for variational
approaches to the low-dimensional polaron problem, for which one could have
naively expected a large sensitivity to quantum fluctuations.

Figure 5:  (Color online) The quasiparticle residue $Z$ of the polaron as a
function of $\eta$. The solid line represents the 1p-h result ($N=1$ diagram).

## Acknowledgments

This work is supported by the Fund for Scientific research - Flanders. We
would like to thank C. Lobo, N. Prokof’ev, B. Svistunov, and F. Werner for
helpful discussions and suggestions, and we thank J. Levinsen and M. Parish
for sending us their data.

## References

  * (1) C. Chin, R. Grimm, P. Julienne and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * (2) I. Bloch, Nat. Phys. 1, 23 (2005).
  * (3) G.B. Partridge, W. Li, Y.A. Liao, R.G. Hulet, M. Haque and H.T.C. Stoof, Phys. Rev. Lett. 97, 190407 (2006).
  * (4) Y. Shin, C.H. Schunck, A. Schirotzek and W. Ketterle, Nature (London)451, 689 (2008).
  * (5) S. Nascimbène, N. Navon, K.J. Jiang, L. Tarruell, M. Teichmann, J. McKeever, F. Chevy and C. Salomon, Phys. Rev. Lett. 103, 170402 (2009).
  * (6) S. Pilati and S. Giorgini, Phys. Rev. Lett. 100, 030401 (2008).
  * (7) P. Massignan, Z. Yu and G.M. Bruun, Phys. Rev. Lett. 110, 230401 (2013).
  * (8) G. Bertaina and S. Giorgini, Phys. Rev. Lett. 106, 110403 (2011).
  * (9) F. Chevy, Phys. Rev. A 74, 063628 (2006).
  * (10) A. Schirotzek, C.-H. Wu, A. Sommer, and M. W. Zwierlein, Phys. Rev. Lett. 102, 230402 (2009).
  * (11) M. Koschorreck, D. Pertot, E. Vogt, B. Fröhlich, M. Feld and M. Köhl, Nature (London) 485, 619 (2012).
  * (12) P. Dyke ,E. D. Kuhnle, S. Whitlock, H. Hu ,M. Mark, S. Hoinka, M. Lingham, P. Hannaford and C.J. Vale, Phys. Rev. Lett.106, 105304 (2011).
  * (13) J. Levinsen and S. K. Baur, Phys. Rev. A 86, 041602 (2012).
  * (14) N.V. Prokof’ev and B.V. Svistunov, Phys. Rev. B 77, 020408 (2008).
  * (15) N.V. Prokof’ev and B.V. Svistunov, Phys. Rev. B 77, 125101 (2008).
  * (16) J. Vlietinck, J. Ryckebusch and K. Van Houcke, Phys. Rev. B 87, 115133 (2013).
  * (17) M. Punk, P.T. Dumitrescu and W. Zwerger, Phys. Rev. A 80, 053605 (2009).
  * (18) R. Combescot, S. Giraud and X. Leyronas, EPL 88, 60007 (2009).
  * (19) C. Mora and F. Chevy, Phys. Rev. A 80, 033607 (2009).
  * (20) J.B. McGuire, J. Math. Phys. 7, 123 (1966).
  * (21) R. Combescot, A. Recati, C. Lobo and F. Chevy, Phys. Rev. Lett. 98, 180402 (2007).
  * (22) S. Giraud and R. Combescot, Phys. Rev. A 79, 043615 (2009).
  * (23) S. Zöllner, G.M. Bruun and C.J. Pethick, Phys. Rev. A 83, 021603(R) (2011).
  * (24) M.M. Parish, Phys. Rev. A 83, 051603(R) (2011).
  * (25) M.M. Parish and J. Levinsen, Phys. Rev. A 87, 033616 (2013).
  * (26) L.D. Landau and E.M. Lifshitz, _Statistical Physics_ (Pergamon, New York, 1980), Vol.2.
  * (27) K. Van Houcke, F. Werner, E. Kozik, N. Prokof’ev, B. Svistunov, M.J.H. Ku, A.T. Sommer, L.W. Cheuk, A. Schirotzek and M.W. Zwierlein, Nat. Phys. 8, 366 (2012).
  * (28) K. Van Houcke, F. Werner, N. Prokof’ev and B. Svistunov, arXiv:1305.3901.
  * (29) P. Bloom, Phys. Rev. B 12, 125 (1975).
  * (30) R. Combescot and S. Giraud, Phys. Rev. Lett. 101, 050404 (2008).
  * (31) N.V. Prokof’ev, B.V. Svistunov and I.S. Tupitsyn, J. Exp. Theor. Phys. 87, 310 (1998).
  * (32) D.M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
  * (33) M. Boninsegni, N. Prokof’ev and B. Svistunov, Phys. Rev. Lett. 96, 070601 (2006).
  * (34) C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G.M. Bruun, F. Schreck, and R. Grimm, Nature (London) 485, 615 (2012)