††thanks: Corresponding author. E-mail<EMAIL_ADDRESS>
# Observation of bosonic stimulation in light scattering

Yu-Kun Lu  Research Laboratory of Electronics, MIT-Harvard Center for
Ultracold Atoms, and Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Yair Margalit  Research
Laboratory of Electronics, MIT-Harvard Center for Ultracold Atoms, and
Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA  Wolfgang Ketterle  Research Laboratory of
Electronics, MIT-Harvard Center for Ultracold Atoms, and Department of
Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA

###### Abstract

For bosons, the transition rate into an already occupied quantum state is
enhanced by its occupation number: the effect of bosonic stimulation. Bosonic
enhancement of light scattering has been predicted more than 30 years ago but
has not been observed before. Here we theoretically investigate and
experimentally demonstrate this effect. We show that the bosonic enhancement
factor for a harmonically trapped gas is bounded by a universal constant
$\zeta(2)/\zeta(3)$ above the phase transition to a Bose-Einstein condensate
(BEC), and depends linearly on the BEC fraction just below the phase
transition. Bosonically enhanced light scattering is observed and
characterized above and below the phase transition, and the effect of
interactions is discussed. For a multi-level system, bosonic enhancement is
reduced because bosonic stimulation occurs only for Rayleigh scattering, but
not for Raman scattering.

Bosonic stimulation occurs for bosonic particles transitioning into a final
state with non-zero occupation number $n$: the transition rate is enhanced by
a factor of $n+1$. Since most elementary particles are fermions, most
observations of bosonic stimulations were related to photons, with the laser
as the most dramatic manifestation. The field of ultracold atoms has created
new opportunities to realize other paradigmatic phenomena for bosonic
stimulation: it was observed in the formation of BEC [1], superradiance [2],
matter-wave amplification [3, 4] and enhanced density fluctuations [5, 6].
However, one phenomenon has never been observed, although it was predicted
more than 30 years ago: bosonic enhancement of light scattering [7, 8, 9, 10,
11, 12, 13]. One reason why the experimental observation has been elusive [14,
15, 16, 17] is the requirement of high densities, ideally many atoms per
$\lambdabar^{3}$ ($\lambdabar=\lambda/2\pi$ is the reduced wavelength). At
such densities, ultracold atomic clouds become short-lived due to inelastic
collisions [18, 19]. For the same reason, the fermionic counterpart, Pauli
blocking of light scattering, was demonstrated only very recently [20, 17,
21].

In this work, we theoretically investigate and experimentally observe bosonic
enhancement of light scattering in an ultracold Bose gas of sodium atoms. We
obtain new analytic results of the enhancement factor for non-interacting
bosons in several regimes. Experimentally, we prepared a dense cloud of bosons
with a peak density up to $2\times 10^{15}\text{cm}^{-3}$ and a BEC transition
temperature of $T_{c}\approx 10\mu K$, achieving the highest $T_{c}$ in
ultracold atoms except for atomic hydrogen [22] and the highest critical
density. The bosonic enhancement of light scattering was observed in two
different regimes: for scattering within the thermal cloud and for scattering
between the BEC and thermally occupied states. We show how interactions modify
bosonic stimulation by modifying the overlap between condensate and thermal
cloud, and by changing the pair correlation function.

Figure 1:  Schematic for bosonic stimulation in light scattering. Shown are
momentum space pictures of the cloud above and below the phase transition.
When a photon is scattered, the photon recoil changes the momentum of the atom
by $\vec{q}$ (red arrow). (a) Above the phase transition ($T>T_{c}$), both the
initial and final momentum states are in the thermal cloud. (b) Below the
phase transition ($T<T_{c}$), a condensate emerges around zero momentum, and
scattering between the condensate and the thermal cloud is also possible. In
either case, light scattering is enhanced if the final state is occupied.
Figure 2:  Bosonic stimulation of light scattering in a trapped 23Na gas
across the phase transition. (a) Experimental data of light scattering versus
temperature for different (normalized) recoil momentum
$\kappa=\sqrt{E_{\text{rec}}/k_{B}T_{c}}$. The dots (diamonds) represents the
data above (below) the phase transition. Data points here are each averaged
over 12 to 16 samples, and the uncertainty is reflected by the scattering of
the data. (b) Theoretical prediction of the bosonic enhancement factor for a
harmonically trapped non-interacting gas. The red trajectory corresponds to
the traces in (a), while the black trajectory corresponds to the trace when an
almost pure condensate is prepared at low density and compressed.

The principle behind bosonic stimulation in light scattering is illustrated in
Fig. 1. Light scattering is a two-photon process where an atom absorbs a
photon (with momentum $\vec{k}_{\text{i}}$) from the probe beam and emits a
scattered photon (momentum $\vec{k}_{\text{f}}$). The momentum of the atom
changes due to the photon recoil by
$\vec{q}=\vec{k}_{\text{i}}-\vec{k}_{\text{f}}$ (Fig. 1a). In general, the
(normalized) light scattering rate is given by the structure factor which can
be obtained in the local density approximation by averaging over the phase-
space [23]:

$S(\vec{q})=1\pm\langle n(\vec{r},\vec{p}+\vec{q})\rangle,$ (1)

with the plus (minus) sign corresponding to bosons (fermions). The phase-space
average is weighted by the occupation number $n(\vec{r},\vec{p})$. In the
limit of zero recoil momentum, the second term reduces to the average phase-
space density (PSD). This term becomes $1$ for fermions at zero temperature,
predicting complete Pauli blocking [20]. For bosons, the physics is richer
because of the presence of the BEC phase transition. For a non-degenerate
thermal gas, the averaged PSD is $n\lambda_{t}^{3}/2^{3/2}$ ($n$ is the
density and $\lambda_{t}$ the thermal de Broglie wavelength). Extrapolating
this formula to the BEC phase transition with the criterion
$n\lambda_{t}^{3}=\zeta(3/2)\approx 2.612$ gives an enhancement factor of
$1.92$ (where $\zeta$ is the Riemann zeta function). However, we will show
that the enhancement factor at the transition temperature can diverge (in free
space) or is limited to 1.37 for harmonic confinement.

The phase-space occupation number $n(\vec{r},\vec{p})$ in Eq. 1 is given by
the Bose-Einstein occupation number in the thermal cloud
$n_{\text{th}}=\left[\exp\left[\left(\epsilon_{\vec{p}}+U(\vec{r})-\mu\right)/k_{B}T\right]-1\right]^{-1}$.
The presence of a non-interacting condensate adds a $\delta$-function in phase
space:

$n(\vec{r},\vec{p})=\begin{cases}n_{\text{th}}(\vec{r},\vec{p})&\text{$T\geq
T_{c}$}\\\
n_{\text{th}}(\vec{r},\vec{p})|_{\mu=0}+h^{3}N_{0}\delta(\vec{r},\vec{p})&\text{$T<T_{c}$}\end{cases}.$
(2)

Here $\epsilon_{\vec{p}}=p^{2}/2m$, $U(\vec{r})$ and $\mu$ represent the
kinetic energy, potential energy, and chemical potential, respectively.
$N_{0}=f(T/T_{c})N$ is the number of atoms in the condensate with $f$ being
the condensate fraction.

Figure 3:  Bosonic enhancement of light scattering above and below the BEC
phase transition. (a) Above the phase transition, the scattering happens
within the thermal cloud. The dashed curve represents the theoretical
prediction for ideal gas while the solid curve represents the interacting
theory. The data is identical to the data in Fig. 2 with $T>T_{c}$ but
averaged over nearby points. The error bars throughout the whole paper are
purely statistical and reflect one standard error of the mean. (b) Below the
phase transition, the scattering between BEC and thermal cloud dominates. The
dashed curves are predictions for ideal gas, while the solid curves used semi-
ideal gas approximation and the modified pair correlation function. The data
is identical to the data in Fig. 2 with $T<T_{c}$.

Above $T_{c}$, scattering occurs between two thermally occupied states (Fig.
1a). The scattering rate is enhanced as the temperature approaches $T_{c}$
because of the increase of $n_{\text{th}}(\vec{r},\vec{p})$ at low momenta.
Below the transition temperature, scattering into and out of the condensate
becomes possible (Fig. 1b). Combining the expression for $n(\vec{r},\vec{p})$
with Eq. 1, the total scattering rate below $T_{c}$ reads:

$S(\vec{q})=1+2fn_{\text{th}}(0,\vec{q})+\langle
n_{\text{th}}(\vec{r},\vec{p}+\vec{q})|_{\mu=0}\rangle.$ (3)

The first term represents the contribution of single-particle Rayleigh
scattering events, the second term is the condensate-to-thermal scattering
events, while the third term describes thermal-to-thermal scattering events.
Already for moderate condensate fractions, the scattering is dominated by
scattering events involving the condensate. As the temperature becomes lower,
the condensate fraction $f$ increases but the occupation number
$n_{\text{th}}(0,\vec{q})$ ($n_{\text{th}}(0,\vec{q})\approx
k_{B}T/E_{\text{rec}}$ for $k_{B}T\gg E_{\text{rec}}$) in the thermal cloud
decreases. As a result, the total scattering rate first increase and then
decrease as the temperature is further decreased below $T_{c}$. Ultimately at
zero temperature, all atoms are in the condensate and bosonic stimulation is
absent because of zero occupation in the final states. Fig. 2b shows the
theoretical calculation for a non-interacting Bose gas in a harmonic trap.

Bosonic enhancement monotonically increases when the (normalized) recoil
momentum $\kappa=\sqrt{E_{\text{rec}}/k_{B}T_{c}}$ decreases (where
$E_{\text{rec}}=\hbar^{2}k^{2}/m$ is the atomic recoil energy). In a harmonic
trap, for $\kappa\rightarrow 0$, the enhancement above the BEC phase
transition is bounded by $\zeta(2)/\zeta(3)\approx 1.37$, while the
enhancement factor below $T_{c}$ diverges as $1/\kappa^{2}$. In contrast, in
free space (3D box potential) the enhancement factor diverges already above
$T_{c}$ as $\pi^{3/2}\kappa^{-1}/\zeta(3/2)$. This major difference
illustrates that in a harmonic trap at $T_{c}$, only a small volume is at or
near the conditions for the phase transition. Derivation of these results, as
well as a generalization to arbitrary power-law potentials, can be found in
the Supplementary Materials [24].

Figure 4:  Bosonic enhancement for Rayleigh and Raman scattering. (a) For
large detuning, the atom's internal states are labelled by $\ket{J,m_{J}}$
where $J$ represents the total angular momentum of the electron and $m_{J}$
its projection along the quantization axis. The blue (red) straight arrow
represents the excitation paths for probe light with $\sigma_{+}$
($\sigma_{-}$) polarization, while the wavy arrows represent the emitted
photons with branching ratio labeled on the side. (b) Observed photon
scattering for $\sigma_{+}$ (blue) and $\sigma_{-}$ (red, scaled by $9/5$)
polarized probe light. The solid blue curve is a fit to the $\sigma_{+}$ data
(using an exponential function as a convenient trial function), the red dashed
light is obtained from the blue fit and Eq. 4. The solid red line is a fit to
Eq. 4 with adjustable $\gamma$. Each data point is the average over 12
samples.

In the experiment, around $4\times 10^{5}$ partially condensed 23Na atoms in
hyperfine state $\ket{F=2,m_{F}=2}$ are prepared in an optical dipole trap
[24]. The cloud is then heated up to different temperatures by modulating the
trap intensity while keeping the atom number almost constant. Light is
scattered by exciting the cycling transition with $\sigma_{+}$ polarized
light. The optical setup is described in our previous work [20]. For the
measurements to be perturbative, the scattering rate was kept below 0.14
photon per atom during 4 ms. We observed light scattering for $4$ different
densities over a range of temperatures (Fig. 2a). The density is controlled by
adjusting the depth of the optical dipole trap and determines $T_{c}$ and
$\kappa$. The red line in Fig. 2 shows the parameters of these measurements
(the traces for the different values of $T_{c}$ overlap).

The data in Fig. 2, even without quantitative analysis, show the salient
features of light scattering for bosons: above the phase transition, there is
a modest enhancement in scattering due to the increase in phase space density,
as the temperature is reduced. Once the condensate forms, the enhancement is
more substantial due to scattering into and out of the condensate. We also
observed the absence of enhanced light scattering for an almost pure Bose-
Einstein condensate. However, due to the strong increase of three-body loss at
high density, we could prepare such condensates only at lower densities in
more weakly confining traps where the maximum enhancement is small due to the
larger value of $\kappa$. When the cloud was compressed to higher densities
(shown as the black trajectory in Fig. 2b), an increase of light scattering
was observed, qualitatively consistent with the prediction [24].

For temperatures above $T_{c}$, light scattering occurs within the thermal
distribution. In Fig. 3a, we compare the results with the predictions for an
ideal gas and an interacting gas, without any free parameters except the
overall scaling. Data with different $T_{c}$ are binned together because of
their small difference in $\kappa$. Usually, thermal Bose gases away from
Feshbach resonances are treated as non-interacting. However, at our high
densities, the mean-field repulsion is relevant even for thermal clouds; it
softens the trapping potential and increases the bosonic enhancement by up to
$50\%$ because of the higher density of states at low energy [24]. Bosonic
stimulation is reduced by a second effect of interactions: the reduction of
the pair correlation function suppresses the enhancement by a factor
$1-8\sqrt{2}a/\lambda_{t}\approx 0.6$ at high temperatures [25, 24]. Close to
the phase transition, the thermal de Broglie wavelength $\lambda_{t}$ is
effectively replaced by the correlation length $\xi$. The divergence of $\xi$
at the phase transition will diminish the suppression caused by pair
correlation. To the best of our knowledge, the suppression of bosonic
stimulation by repulsive interactions has not been discussed before. In the
supplement, we discuss some evidence for the effect of pair correlations and
its reduction near criticality. Fig. 3a demonstrates that the two effects of
interactions almost cancel each other.

For temperatures below the BEC phase transition, the scattering between the
condensate and the thermal cloud dominates the bosonic enhancement. For non-
interacting Bose gas, the enhancement factor increases almost linearly with
condensate fraction close to the critical point: $S(\vec{q})\approx
1+2fk_{B}T/E_{\text{rec}}$ (dashed curves in Fig. 3b). However, because of the
mean-field repulsion between the thermal cloud and the condensate, the density
of thermal atoms inside the condensate volume is strongly reduced [26]
diminishing the bosonic enhancement. The enhancement $S(\vec{q})-1$ is thus
reduced by the factor $(\int
d\vec{r}n_{0}(\vec{r})n_{\text{th}}^{\text{(int)}}(\vec{r},\vec{q}))/(N_{0}n_{\text{th}}(0,\vec{q}))$,
with $n_{0}(\vec{r})$ being the density distribution of the interacting BEC,
and $n_{\text{th}}^{\text{(int)}}(\vec{r},\vec{q})$ the phase-space occupation
number for the interacting thermal cloud. The experimentally observed
enhancement is reduced by up to a factor of 3 compared to the non-interacting
theory (Fig. 3b). We model the effect of interactions using the semi-ideal gas
approximation [26] which assumes an interacting condensate and repulsion of
the thermal cloud by the condensate, and include the modified pair correlation
function [25] (Fig. 3b). Note that the theory curves in Fig. 3b has no free
parameters (the overall scaling was the same as Fig. 3a). Our model explains
most of the reduction in light scattering compared to the predictions for a
non-interacting gas. The remaining difference could be due to the backaction
from the thermal cloud to the BEC which will further reduce their overlap
which can be approximately treated using a Hartree-Fock approximation [24].
However, a fully consistent treatment of the interacting bose gas is
challenging both theoretically and numerically.

Until now, we have assumed a two-level system, realized by exciting the
initial $\ket{F=2,m_{F}=2}$ hyperfine state on a cycling transition with
$\sigma_{+}$ polarized light (see Fig. 4a, blue). We now discuss the situation
when the atoms are excited with light of different polarization. In the limit
of detunings much larger than the hyperfine interaction in the excited state,
the nuclear hyperfine structure can be neglected since the nuclear spin state
does not change during the light scattering. Therefore, we realize the
simplified level diagram shown in Fig. 4a.

When the atoms are initially prepared in a single internal state, only
Rayleigh scattering shows bosonic enhancement. Raman scattering (or optical
pumping) populates initially empty states and is not enhanced. We therefore
expect the enhancement factor $\eta$ to be reduced to $\gamma\eta+(1-\gamma)$
where $\gamma$ represents the branching ratio (ratio of Rayleigh scattering to
Raman scattering). To describe our experimental situation, we express the
observed scattering rate for $\sigma_{+}$ light as $R_{+}=\eta R$, with $R$
representing the Rayleigh scattering rate without any enhancement. The
expected observed scattering rate for $\sigma_{-}$ light is then predicted as

$R_{-}=\frac{1}{3}(\gamma\eta R+2(1-\gamma)R).$ (4)

Here the prefactor $1/3$ comes from the reduced matrix element, the branching
ratio $\gamma=1/3$, and the extra factor of two comes from the different
angular distributions for the scattered $\sigma$ and $\pi$ light. By comparing
Eq. 4 in the high temperature and low-temperature limit, we find the
enhancement factor for $\sigma_{-}$ polarized light to be $(4+\eta)/5$. The
results presented in Fig. 4b confirm that the bosonic enhancement for
$\sigma_{-}$ probe light is always below the enhancement for $\sigma_{+}$. We
can use Eq. 4 and regard $\gamma$ as a fitting parameter which characterizes
the bosonic enhancement. The best fit to the data gives $\gamma=0.5(2)$, and
is consistent with the theoretical prediction ($\gamma=1/3$) within the
experimental uncertainty (which is larger for $\sigma_{-}$ scattering due to
the small number of photons observed). If the light scattering is non-
perturbative, the population can build up in the initially empty state via
optical pumping which will then lead to bosonic enhancement of the Raman
transition which has been observed in the form of Raman superradiance [27,
28].

In conclusion, we have shown how quantum statistics and interaction modifies
the optical properties of a Bose gas. We have worked in a regime where the
effect of quantum statistics dominate over other effects such as dipole-dipole
interactions [29], superradiance [2], light-assisted collisions [30], multiple
scattering and atomic lensing. Besides the fundamental interest, understanding
light scattering is crucial for quantitative diagnostics of bosonic systems.
For future work, it is promising to study bosonic enhancement of light
scattering in a box potential [31] because of the absence of density
inhomogeneity. The large extension of the critical region will give much
stronger enhancement and enable studies of the $1/\kappa$ divergence of the
enhancement factor. In addition, in a box potential, interactions can affect
light scattering only through the modification of the pair correlation
function. It should be possible to clearly observe how bosonic enhancement is
reduced above the phase transition, but not right at the critical point.
Furthermore, for strong repulsive interactions near a Feshbach resonance,
light scattering in bosons could possibly be even suppressed below the single-
particle Rayleigh rate, in analogy with Pauli blocking for fermions. Both
Pauli blocking and strong interactions reduce the pair correlation function.
More generally, light scattering provides a way to measure the static
structure factor which could be used to characterize strongly interacting
systems, including systems with strong dipolar interactions which are
anisotropic and long-range.

## References

  * Miesner  _et al._ [1998] H.-J. Miesner, D. Stamper-Kurn, M. Andrews, D. Durfee, S. Inouye, and W. Ketterle, Science 279, 1005 (1998).
  * Inouye  _et al._ [1999a] S. Inouye, A. Chikkatur, D. Stamper-Kurn, J. Stenger, D. Pritchard, and W. Ketterle, Science 285, 571 (1999a).
  * Inouye  _et al._ [1999b] S. Inouye, T. Pfau, S. Gupta, A. Chikkatur, A. Görlitz, D. Pritchard, and W. Ketterle, Nature 402, 641 (1999b).
  * Kozuma  _et al._ [1999] M. Kozuma, Y. Suzuki, Y. Torii, T. Sugiura, T. Kuga, E. Hagley, and L. Deng, Science 286, 2309 (1999).
  * Esteve  _et al._ [2006] J. Esteve, J.-B. Trebbia, T. Schumm, A. Aspect, C. I. Westbrook, and I. Bouchoule, Phys. Rev. Lett. 96, 130403 (2006).
  * Blumkin  _et al._ [2013] A. Blumkin, S. Rinott, R. Schley, A. Berkovitz, I. Shammass, and J. Steinhauer, Physical Review Letters 110, 265301 (2013).
  * Svistunov and Shlyapnikov [1990] B. Svistunov and G. Shlyapnikov, Sov. Phys. JETP 70, 460 (1990).
  * Politzer [1991] H. D. Politzer, Phys. Rev. A 43, 6444 (1991).
  * Zhang and Walls [1994] W. Zhang and D. F. Walls, Phys. Rev. A 49, 3799 (1994).
  * Morice  _et al._ [1995] O. Morice, Y. Castin, and J. Dalibard, Phys. Rev. A 51, 3896 (1995).
  * Javanainen and Ruostekoski [1995] J. Javanainen and J. Ruostekoski, Phys. Rev. A 52, 3033 (1995).
  * Javanainen [1994] J. Javanainen, Phys. Rev. Lett. 72, 2375 (1994).
  * Politzer [1997] H. D. Politzer, Phys. Rev. A 55, 1140 (1997).
  * foo [a] Reference 15 claimed the observation of bosonic enhancement of light scattering. However, their observable was atom loss caused by optical pumping. As we show in the last part of the paper, optical pumping (i.e. Raman scatteing) into an unoccupied state does not show any bosonic enhancement. The larger observed loss is possibly caused by the higher atomic density in a Bose-Einstein condensate.
  * Bons  _et al._ [2016] P. C. Bons, R. de Haas, D. de Jong, A. Groot, and P. van der Straten, Phys. Rev. Lett. 116, 173602 (2016).
  * foo [b] After the experimental work has been completed and presented at the BEC 2021 conference, possible evidence of bosonic stimulation effects was observed in optical absorption measurements which showed an enhancement of 10 - 15 % across the BEC phase transition (Ref. 17). However, this effect was explained without considering light scattering involving the condensate, and no full analysis was given (since the topic of the paper was Pauli blocking for fermions).
  * Deb and Kjærgaard [2021] A. B. Deb and N. Kjærgaard, Science 374, 972 (2021).
  * Inguscio  _et al._ [1999] M. Inguscio, S. Stringari, and C. Wieman,  _Bose-Einstein condensation in atomic gases_ , Vol. 140 (IOS Press, 1999).
  * Top  _et al._ [2021] F. Ç. Top, Y. Margalit, and W. Ketterle, Physical Review A 104, 043311 (2021).
  * Margalit  _et al._ [2021] Y. Margalit, Y.-K. Lu, F. Çağrı Top, and W. Ketterle, Science 374, 976 (2021).
  * Sanner  _et al._ [2021] C. Sanner, L. Sonderhouse, R. B. Hutson, L. Yan, W. R. Milner, and J. Ye, Science 374, 979 (2021).
  * Fried  _et al._ [1998] D. G. Fried, T. C. Killian, L. Willmann, D. Landhuis, S. C. Moss, D. Kleppner, and T. J. Greytak, Phys. Rev. Lett. 81, 3811 (1998).
  * Shuve and Thywissen [2009] B. Shuve and J. Thywissen, Journal of Physics B: Atomic, Molecular and Optical Physics 43, 015301 (2009).
  * [24] See Supplementary Materials.
  * Naraschewski and Glauber [1999] M. Naraschewski and R. J. Glauber, Phys. Rev. A 59, 4595 (1999).
  * Naraschewski and Stamper-Kurn [1998] M. Naraschewski and D. M. Stamper-Kurn, Phys. Rev. A 58, 2423 (1998).
  * Schneble  _et al._ [2004] D. Schneble, G. K. Campbell, E. W. Streed, M. Boyd, D. E. Pritchard, and W. Ketterle, Phys. Rev. A 69, 041601 (2004).
  * Yoshikawa  _et al._ [2004] Y. Yoshikawa, T. Sugiura, Y. Torii, and T. Kuga, Phys. Rev. A 69, 041603 (2004).
  * Zhu  _et al._ [2016] B. Zhu, J. Cooper, J. Ye, and A. M. Rey, Phys. Rev. A 94, 023612 (2016).
  * Weiner  _et al._ [1999] J. Weiner, V. S. Bagnato, S. Zilio, and P. S. Julienne, Rev. Mod. Phys. 71, 1 (1999).
  * Navon  _et al._ [2021] N. Navon, R. P. Smith, and Z. Hadzibabic, Nature Physics 17, 1334–1341 (2021).

Acknowledgments
We thank Julius de Hond and Pierre Barral for comments on the manuscript.
Funding: This work is supported by NSF through the Center for Ultracold Atoms
and Grant No. 1506369, and from a Vannevar-Bush Faculty Fellowship. Author
Contributions: All authors contributed to the concepts of the experiment,
Y.-K.L. and Y.M performed the experiment, Y.-K.L. analyzed the data and
performed the theoretical calculations, Y.-K.L., Y.M and W.K. wrote the paper.
Competing interests: The authors declare no competing interests. Data
Availability: The data that support the plots within this paper and other
findings of this study are available from the corresponding authors upon
reasonable request.

Supplementary Materials
Materials and Methods
Supplementary Text
Figs. S1-S6
Table S1-S2
References (S1-S9)

Supplemental Materials: Observation of bosonic stimulation in light scattering

##  S1 Experimental details

Sample preparation Around 1 billion sodium atoms are loaded into a magneto-
optical trap (MOT) and transferred into a quadrupole magnetic trap with an
optical plug where they are cooled by microwave evaporative cooling.
Subsequently, around $10^{7}$ atoms in hyperfine state $\ket{F=2,m_{F}=2}$ are
transferred into a single-beam 1064 nm optical dipole trap. The beam waist and
intensity of the trapping beam are controlled by a variable-size iris and an
acousto-optic modulator, respectively. The beam waist was initially chosen to
be big to increase the loading efficiency into a large trap volume. Because
the bosonic enhancement effect requires $k_{B}T_{c}\gg E_{\text{rec}}$, the
atomic gas is compressed to a higher density with a higher $T_{c}$ by fully
opening up the iris. Decreasing the intensity of the optical trapping beam
induces evaporative cooling to a Bose-Einstein condensed cloud with around
$4\times 10^{5}$ atoms and around 30% condensate fraction. The final $T_{c}$
can be adjusted from $6.8\mu K$ to $11\mu K$ by changing the final power of
the optical trap. Subsequent heating of the cloud was done by sinusoidal
modulating the power of the optical trapping beam by 50% at 1 kHz frequency
with variable duration (from 1 ms to 49 ms).

Light scattering and collection The probe light was generated from an
independent 1178 nm laser with frequency doubling. The frequency of the probe
light was red detuned from the D2 transition of sodium by $\Delta=25.7\text{
GHz}$. The large detuning was chosen to avoid multiple scattering of photons,
to minimize the lensing effect and photo-association loss. In addition, since
the detuning is much larger than the hyperfine interaction in the
electronically excited state ($\Delta\gg A_{\text{HFS}}=18.5\text{ MHz}$), the
nuclear spin of the atom is preserved during the light scattering. The
$1/e^{2}$ diameter of the probe beam is $3.3$ mm in the atoms' plane. The
intensity of the beam was stabilized to avoid shot-to-shot fluctuations.

Fluorescence light is collected at a 90-degree angle relative to the probe
beam with a collection efficiency of around 0.3%. To detect the small number
of scattered photons, we used an EMCCD camera with quantum efficiency of 97%
and readout noise of 3 - 5$e^{-}$ per pixel. Hardware binning of pixels was
implemented to make sure the photon shot noise dominates over the readout
noise. Background noise was suppressed by choosing a region of interest only
slightly larger than the image of the cloud.

For the light scattering in Figs. 2 and 3 at most 0.14 photons were scattered
per atom. The pulse duration was 4 ms. For the data in Fig. 4, since the
$\sigma_{-}$ polarized probe light will cause optical pumping, we used weaker
pulses with at most 0.05 (0.02) photon per atom for Rayleigh (Raman)
scattering and a pulse duration of 3 ms.

We verified that the atoms were probed in the perturbative regime where
optical pumping and nonlinear effects such as superradiance are absent, by
measuring the scattered light as a function of probe beam power (Fig. S1). The
measurement was done for both a cloud with and without BEC. The linear
behavior shows that the light scattering is indeed perturbative.

Thermometry The temperature of the cloud was determined from time-of-flight
absorption measurement. The fitting of absorption images was done using the
Bose function for the thermal cloud and Thomas-Fermi distribution for the BEC
[1]. Since the atoms were released from the trap 8 ms after the light pulse,
we had to account for a possible change in temperature during this hold time.
The hold time was chosen to allow the atoms to thermalize, and was also
necessary for a mechanical camera shutter to open which was necessary to block
the probe beam during the light scattering.

Because of the high density of the cloud, the lifetime of the cloud due to
three-body recombination was around 15 ms. Heating due to particle loss
(``anti-evaporation'') during the light scattering and the extra hold time was
therefore non-negligible. We accounted for this effect by measuring the
temperature for identical cycles, but without the light scattering pulse and
the extra hold time. This temperature (averaged over multiple cycles) was
measured for each sample heating time and used for the temperature axes in
Figs. 2 and 3.

Figure S1:  Linearity check for light scattering. Scattered photons are
measured for different probe beam powers. The blue (red) data points
represents a cloud with (without) BEC. The absence of nonlinearity shows that
we are working in the perturbative regime. Data points here are each averaged
over 3 samples.

Properties of the high-density cloud Significant Bose enhancement of light
scattering at large angles requires the critical temperature to be higher than
the recoil energy. Therefore, the density $n$ of our sample needed to be
around $10^{15}\text{ cm}^{-3}$. Since the three-body loss rate coefficient
for thermal sodium atom in $F=2$ states at zero magnetic field is $9.6\times
10^{-29}\text{ cm}^{6}\text{/s}$, the lifetime of the cloud is limited to
around 15 ms. Another consequence of the high density is the short mean free
path ($l=1/n\sigma\approx 1\mu m$ where $\sigma$ is the cross-section for
elastic collisions). Since the typical size of the cloud is around $1\mu m$ in
the radial direction and $30\mu m$ in the axial direction, this implies slow
hydrodynamic transport in the axial direction. The typical trap frequencies in
the radial (axial) direction are 9.5 kHz (290 Hz). Due to hydrodynamic
collisional effects, the equilibrium time scale in the axial direction is much
slower than the trap period and can potentially lead to non-fully equilibrated
clouds in combination with fast three-body loss. More precise measurement of
bosonic light scattering should therefore be done at lower densities and
smaller scattering angles, or by choosing atoms with even lower three-body
loss rate (e.g., ${}^{23}\text{Na}$ in the $F=1$ state or
${}^{164}\text{Dy}$).

##  S2 Light scattering for a pure BEC and adiabatic compression

A pure Bose-Einstein condensate shows no bosonic enhancement of light
scattering. We therefore studied light scattering using an almost pure
condensate, and then compressed the condensate adiabatically by ramping up the
power of the optical trap. An almost pure BEC could be prepared only in a
shallow trap. Subsequently, the trap power was ramped up to increase the
density and therefore $T_{c}$. Light scattering was measured for different
final trap depths. Without compression, the photon count rate is identical to
the signal for a non-degenerate thermal gas (Fig. S2), demonstrating the
absence of bosonic enhancement in a pure condensate. However, for the low
values of $\kappa=1.5$, the maximum predicted enhancement for a mixed cloud is
only 9%. For higher densities, bosonic enhancement of light scattering was
observed, in qualitative agreement with the theory. During compression, the
cloud approximately followed the black trajectory in Fig. 2. Due to losses and
heating during compression, the condensate fraction decreased, and light
scattering was enhanced. Due to the heating, we were not able to map out the
region of low $T/T_{c}$ and low $\kappa$, where the enhancement would increase
with temperature.

Figure S2:  Light scattering after compressing the cloud. The photon
scattering signal was measured for different final trap depths. The dashed
line represents the signal observed for a non-degenerate cloud with no bosonic
enhancement. Each data point was averaged over 6 samples.

##  S3 Theory for the bosonic enhancement factor

The bosonic enhancement factor approaches zero in the limit of low density or
large recoil, $\kappa\rightarrow\infty$, and is maximized in the high-density
limit $\kappa=0$ where the normalized light scattering is given by the
$S(\vec{q}=0)$ at zero momentum. We now rigorously calculate $S(0)$ for power-
law densities of state $g(\epsilon)=A\epsilon^{x}$. For a generic power-law
trapping potential $V(r)\propto r^{\alpha}$ in $d$ dimensional space, the
exponent $x$ in the density of state reads [2]:

$x=\frac{d}{2}+\frac{d}{\alpha}-1.$ (S5)

The box potential has $\alpha=\infty$, and in three dimensions $x=1/2$.

Rewriting the integrals in Eq. 1 in terms of energy $\epsilon$, we have:

$S(\vec{q}=0)=1+\frac{\int_{0}^{\infty}d\epsilon
n(\epsilon)^{2}g(\epsilon)}{\int_{0}^{\infty}d\epsilon
n(\epsilon)g(\epsilon)}.$ (S6)

Here $n(\epsilon)=(z^{-1}\exp(\beta\epsilon)-1)^{-1}$ is the occupation number
with $z=\exp(\beta\mu)$ being the fugacity and $\beta=1/k_{B}T$.

The integrals have analytic solutions:

$\int_{0}^{\infty}d\epsilon
n(\epsilon)g(\epsilon)=A\Gamma(1+x)g_{1+x}(z)/\beta^{1+x},$ (S7)
$\int_{0}^{\infty}d\epsilon
n(\epsilon)^{2}g(\epsilon)=A\Gamma(1+x)(-g_{1+x}(z)+g_{x}(z))/\beta^{1+x}.$
(S8)

Here $g_{n}(z)$ represents the $n$th order polylogarithm function. Finally, we
get a simple expression for the bosonic enhancement factor:

$S(0)=\frac{g_{x}(z)}{g_{1+x}(z)}.$ (S9)

It can be shown that ${g_{x}(z)}/{g_{1+x}(z)}$ has its maximum at the BEC
phase transition $z=1$. The behavior of the polylogarithm functions at $z=1$
divides the physics into three regimes depending on the value of $x$. The
results are summarized in table S1.

When $x\leq 0$, both $g_{x}(z)$ and $g_{1+x}(z)$ diverge at $z=1$. The
divergence of $g_{1+x}(z)$ implies the absence of the BEC phase transition.
Therefore, $z<1$ for all temperatures, and the enhancement factor
$g_{x}(z)/g_{1+x}(z)$ is always finite.

When $0<x\leq 1$, $g_{x}(z)$ diverges at $z=1$ but $g_{1+x}(z)$ remains
finite. This implies the presence of a BEC phase transition and (for
$\kappa=0$) an unbounded light scattering rate with a divergent enhancement
factor at the BEC phase transition point.

When $x>1$, both $g_{x}(z)$ and $g_{1+x}(z)$ remains finite at $z=1$. This
implies the presence of a BEC phase transition and an bounded light scattering
rate. The maximum possible enhancement factor at the tranistion temperature is
${\zeta(x)}/{\zeta(1+x)}$.

Exponent $x$ | $x\leq 0$ | $0<x\leq 1$ | $x>1$
---|---|---|---
BEC phase transition | No | Yes | Yes
Bosonic enhancement | Bounded | Diverges as $\kappa^{2x-2}$ when $T\rightarrow T_{c}$ | Bounded by ${\zeta(x)}/{\zeta(1+x)}$ for $T\geq T_{c}$
Table S1:  Bosonic enhancement factor for different exponent $x$ in the
density of states. In the limit of zero recoil momentum, the enhancement
factor is bounded for $x\leq 0$ and $x>1$ while diverging for $0<x\leq 1$.

For the relevant case of the 3D box potential ($x=1/2$), light scattering can
diverge, whereas for the 3-dimensional harmonic trapping potential ($x=2$) the
maximum enhancement factor is $\zeta(2)/\zeta(3)\approx 1.37$. Figs. S3 and S4
illustrate those cases.

Figure S3:  Bosonic enhancement factor for an ideal gas in a 3D box potential
for different recoil momenta $\kappa$. At the phase transition point, the
bosonic enhancement factor diverges as $1/\kappa$.

To study the divergence of the enhancement factor near the phase transition
for $0<x\leq 1$, we now determine the behavior for small, but non-zero recoil
momentum $\kappa$.

We start with the exact expression for the enhancement factor at finite recoil
$\kappa$ at the phase transition:

$S(q)=1+\frac{1}{\sqrt{\pi}\Gamma(x+1/2)\zeta(x+1)}\int_{0}^{\infty}da\int_{-\infty}^{\infty}dy\frac{1}{\exp(a+y^{2})-1}\frac{a^{x-1/2}}{\exp(a+(y+\kappa)^{2})-1}.$

Here $y$ represents the dimensionless momentum in the recoil direction, and
$a$ represents the total energy from all other dimensions in phase space.
Since the enhancement factor diverges as $\kappa^{2-2x}$ for $0<x\leq 1$, we
need to calculate:

$I=\lim_{\kappa\rightarrow
0}\int_{0}^{\infty}da\int_{-\infty}^{\infty}dy\frac{\kappa^{2-2x}}{\exp(a+y^{2})-1}\frac{a^{x-1/2}}{\exp(a+(y+\kappa)^{2})-1}.$
(S10)

After transforming to polar coordinates $a=r^{2}\sin(\theta)^{2}$ and
$y=-r\cos(\theta)$, one obtains $I=\lim_{\kappa\rightarrow
0}\int_{0}^{\infty}dr\int_{0}^{\pi}d\theta f(r,\theta)$ with

$f(r,\theta)=\frac{\kappa^{2-2x}}{\exp(r^{2})-1}\frac{2r(r\sin(\theta))^{2x}}{\exp(r^{2}-2r\kappa\cos(\theta)+\kappa^{2})-1}.$
(S11)

We first perform the integral over $r$. The function $f(r,\theta)$ has a
branch cut from $r=0$ to $r=+\infty$. It also has poles at
$r_{\pm,n}=\kappa\cos(\theta)\pm\sqrt{-\kappa^{2}\sin(\theta)^{2}+2\pi ni}$
with $n$ being an integer. The residues at $r=r_{\pm.n}$ are:

$\text{Res}(f(r_{\pm,n}))=\pm\frac{2\kappa^{2-2x}\sin(\theta)^{2x}}{\exp((r_{\pm})^{2})-1}\frac{(r_{\pm})^{2x+1}}{r_{+}-r_{-}}.$
(S12)

In the limit of $\kappa\rightarrow 0$, the residues vanishes if $n\neq 0$. For
$n=0$, the residues remain finite:

$\text{Res}(f(r_{\pm,0}))=\pm\frac{2\kappa^{2-2x}\sin(\theta)^{2x}}{\exp(\kappa^{2}\exp(\pm
2i\theta))-1}\frac{(\kappa\exp(\pm
i\theta))^{2x+1}}{2i\kappa\sin(\theta)}\approx-\pm i(\exp(\pm
i\theta)\sin(\theta))^{2x-1}.$ (S13)

By using the Hankel contour in the complex $r$ plane, the integral $I$ can be
related to the sum of all residues in the complex plane by:

$I(1-\exp(2i\pi x))=2\pi
i\sum_{n}\int_{0}^{\pi}d\theta\text{Res}(f(r_{\pm,n},\theta)),$ (S14)

with the result:

$I=\frac{2\pi
i\int_{0}^{\pi}d\theta(-i(\exp(i\theta)\sin(\theta))^{2x-1}+i(\exp(-i\theta)\sin(\theta))^{2x-1})}{1-\exp(2i\pi
x)}=\frac{2\pi^{2}}{4^{x}\sin(\pi x)}.$ (S15)

Putting all pre-factors in, one obtains the divergent behavior of the bosonic
enhancement factor as:

$S(\kappa\rightarrow 0)\approx\frac{2\pi^{3/2}\kappa^{2x-2}}{4^{x}\sin(\pi
x)\Gamma(x+1/2)\zeta(x+1)}.$ (S16)

For $0<x<1$, the divergence in Eq. S16 can be understood in an intuitive way.
In $d$ spatial dimensions, the first-order correlation function of the non-
interacting Bose gas approaches $g^{(1)}(r)\sim 1/r^{d-2}$ at the critical
point. Because the structure factor of a homogeneous sample
$S_{\text{box}}(k)$ is related to the Fourier transform of the pair
correlation function $g^{(2)}(r)-1=g^{(1)}(r)^{2}\sim 1/r^{2d-4}$,
$S_{\text{box}}(k)$ diverges as $\kappa^{d-4}$ when $\kappa\rightarrow 0$. For
a Bose gas trapped in a power-law potential $V(r)=V_{0}(r/L)^{\alpha}$ ($L$ is
the characteristic length of the system), only a fraction of the gas is close
to local criticality at $T=T_{c}$. This fraction can be estimated from the
characteristic length scale $r_{0}$ where the trapping potential energy is
similar to the recoil energy: $V_{0}(r_{0}/L)^{\alpha}\sim k^{2}/2m$.
Therefore, the fraction $f\sim(r_{0}/L)^{d}\sim\kappa^{2d/\alpha}$.
Qualitatively, we can treat the sample within $r_{0}$ as homogeneous and apply
the result for $S_{\text{box}}(k)$ at criticality to this small region. In the
end, we get the structure factor for the trapped gas as
$S(k)\sim\kappa^{2d/\alpha}\kappa^{d-4}\sim\kappa^{2x-2}$.

For $x>1$, the bosonic enhancement mainly comes from regions outside the
critical region which is very small, therefore the overall enhancement factor
is finite.

Figure S4:  Bosonic enhancement factor for an ideal gas in a 3D harmonic
trapping potential for different recoil momenta $\kappa$. The enhancement
factor is bounded at and above the phase transition, but will diverge as
$1/\kappa^{2}$ below the phase transition.

##  S4 Interaction effects on light scattering

Density distribution The repulsive interaction between atoms modifies the
density distribution in the trapping potential and greatly reduces the overlap
between the thermal cloud and the condensate. The density distribution of the
cloud was obtained using the local density approximation. Locally, the Bose
gas is treated as a homogeneous gas with a Hartree-Fock (HF) mean-field
Hamiltonian [3]:

$H_{\text{HF}}-\mu
N=-gnN-\frac{gn_{0}N_{0}}{2}+\sum_{k}\left(\epsilon_{k}-\mu+2gn\right)a^{\dagger}_{k}a_{k}.$
(S17)

Here $n_{0}$ represents the density of the BEC while $n$ is the total density.
$\mu$ is the chemical potential, and $g=4\pi\hbar^{2}a/m$ represents the
interaction strength. $N$ and $N_{0}$ are the total atom number and atom
number in the condensate, respectively. $a_{k}^{\dagger}$ and $a_{k}$ are the
creation and annihilation operators for a free particle with momentum $k$ and
kinetic energy $\epsilon_{k}=k^{2}/2m$. For a given chemical potential, the
local density is obtained by solving the following pair of equations self-
consistently:

$\left\\{\begin{aligned}
n&=n_{0}+\frac{1}{\lambda_{t}^{3}}\,g_{3/2}(e^{\beta(\mu-2gn)})\\\ \mu&=2gn-
gn_{0}\\\ \end{aligned}\right.$ (S18)

We used an open-source package to solve these equations numerically for a
harmonic trapping potential [4]. Note that here all calculations were done by
assuming an isotropic harmonic trap with the mean trap frequency
${\omega}=(\omega_{x}\omega_{y}\omega_{z})^{1/3}$. Unlike our previous work
with fermions [5], the anharmonicity corrections for bosons here are much
smaller and will be neglected. Fig. S5 shows the density profiles of a Bose
gas for an ideal gas, in the semi-ideal gas approximation and HF
approximation. Ref. [6] discusses the differences between the three models and
compares them to experimental data. The ideal gas approximation assumes no
interaction between atoms. The semi-ideal gas approximation considers the
interaction within the BEC and between BEC and the thermal cloud, but not the
interaction within the thermal cloud and the backaction of the thermal cloud
onto the condensate. The mean-field repulsion of the thermal cloud by the
condensate significantly reduces their overlap, causing the reduction of Bose
enhancement. As a comparison, the HF approximation takes all interactions into
account at the mean-field level. The calculations in Fig. S5 show that the
backaction from the thermal cloud compresses the BEC to higher densities and
further reduces their overlap by around 20%. However, the HF calculation
predicts an unphysical jump of the BEC density (Fig. S5). For modeling the
density distribution in the presence of a condensate, we used the semi-ideal
gas approximation since it is numerically a much simpler approach. For a fully
consistent calculation of the density profiles, one would have to go beyond
the HF approximation.

To understand the importance of the different interaction terms, it is useful
to compare the peak density of the thermal cloud at the phase transition
($n_{\text{th}}$) to the peak density of the BEC at $T=0$ ($n_{0}$). By using
$n_{\text{th}}\propto N/(k_{B}T_{c}/m\omega^{2})^{3/2}$ and
$k_{B}T_{c}=\hbar\omega(N/\zeta(3))^{1/3}$, we get $n_{\text{th}}\propto
N^{1/2}/a_{\text{ho}}^{3}$ ($a_{\text{ho}}=\sqrt{\hbar/m\omega}$ is the
oscillator length). For the condensate, we can use $n_{0}=\mu_{0}/g$ with
$\mu_{0}=\hbar\omega(15Na/a_{\text{ho}})^{2/5}/2$ and $g=4\pi\hbar^{2}a/m$ to
obtain $n_{0}\propto N^{2/5}/(a^{3/5}a_{\text{ho}}^{12/5})$. Thus, the ratio
of the densities of the thermal cloud and the BEC becomes
$n_{\text{th}}/n_{0}\propto N^{1/6}a/a_{\text{ho}}$. Although the conventional
wisdom is that the density of the BEC is usually much higher than the thermal
cloud, $n_{\text{th}}$ and $n_{0}$ differ by only a factor of two in our
experiment because of the tight confinement (Fig. S5). Since the mean-field
interaction of the thermal cloud has an extra factor of two due to the
exchange term, we have a situation where the interactions in the condensate
and the thermal cloud are equally important.

Figure S5:  The density distribution of the condensate and thermal cloud for
different models. In the calculations we used the following values: scattering
length of $a=85a_{0}$ ($a_{0}$ is the Bohr radius), trap frequency
$\omega=2\pi\times 2.7\text{kHz}$, atom number $N=4\times 10^{5}$ and
condensate fraction of 30%.

Pair correlation function Atomic interactions will affect light scattering
also through the pair correlation function. The exact expression of the pair
correlation function for the non-interacting Bose gas reads [7]:

$G^{(1)}(r)=\sum_{l=1}^{\infty}\frac{e^{l\beta\mu}}{\lambda_{t}^{3}l^{3/2}}e^{\frac{-\pi
r^{2}}{l\lambda_{t}^{2}}}.$ (S19)

At high temperatures, the $l=1$ term dominates, so the correlation function is
a Gaussian with a correlation length $\lambda_{t}/\sqrt{2\pi}$. Close to the
critical point, the correlation function decays exponentially with a
correlation length $\xi$ [8]:

$g^{(1)}{(r)}\sim\frac{\lambda_{t}}{\zeta(3/2)r}e^{-r/\xi},\quad\xi=\frac{1}{\sqrt{-2m\mu}}=\frac{\lambda_{t}}{\sqrt{-4\pi\beta\mu}}.$
(S20)

Here $g^{(1)}(r)=G^{(1)}(r)/G^{(1)}(0)$ is the normalized first-order
correlation function. The correlation length $\xi$ diverges as
$\left|T-T_{c}\right|^{-1}$ near the critical point for the non-interacting
Bose gas. Using the correlation function for the non-interacting gas, we can
calculate the structure factor $S(k)$ for the ideal gas and the first-order
correction of the structure factor with interactions [7], $S_{\text{int}}(k)$.
Near the critical point, we obtain:

$\begin{array}[]{l}S(k=0)=\int_{0}^{\infty}\left|g^{(1)}{(r)}\right|^{2}4\pi
r^{2}dr\propto\lambda_{t}^{2}\int_{0}^{\infty}\frac{e^{-2r/\xi}}{r^{2}}4\pi
r^{2}dr=2\pi\lambda_{t}^{2}\xi,\\\ S_{\text{int}}(k=0)\approx
S(k=0)-\int_{\lambda_{t}}^{\infty}\left|g^{(1)}{(r)}\right|^{2}\frac{4a}{r}4\pi
r^{2}dr\propto
2\pi\lambda_{t}^{2}\xi-16\pi\lambda_{t}^{2}a\Gamma\left(0,\frac{2\lambda_{t}}{\xi}\right)\sim
2\pi\lambda_{t}^{2}\xi-16\pi\lambda_{t}^{2}a\ln\frac{\xi}{2\lambda_{t}}.\\\
\end{array}$ (S21)

This expression shows that the relative correction due to interactions
$S_{\text{int}}/S-1\sim-8a\ln(\xi/\lambda_{t})/\xi$ vanishes as we approach
the critical point. Note that the calculations were done by assuming $k=0$
because the correction for finite momentum is quadratic in $k$ and can be
neglected.

At high temperatures, $g^{(1)}(r)\propto e^{{-\pi r^{2}}/{\lambda_{t}^{2}}}$
and the structure factor and its interacting correction are given by [7]:

$\begin{array}[]{l}S(k=0)=\int_{0}^{\infty}\left|g^{(1)}{(r)}\right|^{2}4\pi
r^{2}dr\propto\int_{0}^{\infty}e^{{-2\pi r^{2}}/{\lambda_{t}^{2}}}4\pi
r^{2}dr=\frac{\lambda_{t}^{3}}{2\sqrt{2}},\\\
S_{\text{int}}(k=0)=S(k=0)-\int_{0}^{\infty}\left|g^{(1)}{(r)}\right|^{2}\frac{4a}{r}4\pi
r^{2}dr\propto\frac{\lambda_{t}^{3}}{2\sqrt{2}}\left(1-8\sqrt{2}\frac{a}{\lambda_{t}}\right).\\\
\end{array}$ (S22)

Therefore, the correction factor for interaction effects approaches
$1-8\sqrt{2}{a}/{\lambda_{t}}$ at high temperatures. Note that the interacting
pair correlation function describes the true correlation in the position of
the atoms, whereas the non-interacting second-order correlation function
$g^{(2)}_{0}(r=0)$ which is 1 for a pure condensate and 2 for a thermal cloud
appears in the interaction parameters for the condensate and the thermal
cloud, respectively [9, 7].

Figure S6:  Comparison between different models above the phase transition.
(a) Predictions for the Bose enhancement factor for different models. (b)
Fitting of the models to the experimental data. The only free parameter is the
overall scaling. In the calculations we used the following parameters:
scattering length $a=85a_{0}$ ($a_{0}$ is the Bohr radius), atom number
$N=4\times 10^{5}$ and dimensionless recoil momentum $\kappa=0.51$.

Comparing interacting models to experimental data above the phase transition
The bosonic enhancement of light scattering is about 30% just above the phase
transition temperature. The mean-field modification of the density
distribution increases the light scattering by about 10%, whereas the
modification of the pair correlation function leads to a decrease by about
10%. This is why the predictions for the full theory for interactions are
close to the ideal gas prediction. This is shown in Fig. S6 where three
different models are compared to the ideal gas: The first mean-field (MF)
model includes only the modified density distribution of the cloud, but not
the pair correlation function, while the second model (MF + overall
suppression) approximates the effect of the modified pair correlation as an
overall reduction factor of $1-8\sqrt{2}a/\lambda_{t}$. The third model (full
interacting theory) includes the dependence of the pair correlation function
on the local chemical potential. So we locally calculated the pair correlation
function (without assuming a Gaussian or exponential form) and the structure
factor $S_{\text{int}}(k)$ which gives the local scattering rate. The total
light scattering rate is obtained by integrating the local rates across the
cloud.

We fit the data using these three different models by only adjusting the
overall scaling. A comparison of the $\chi^{2}$ values for the fits shows that
the best fit of the data is provided by the full interacting model (table S2).
Therefore, our observations provide some evidence for the modification of pair
correlations by interactions, and even for the different effects of the pair
correlations near criticality. However, the statistical significance is
modest, as illustrated in Fig. S6b. Note that some of the differences in Fig.
S6a disappeared because the high-temperature asymptotic value is a fit
parameter in Fig. S6b. Unfortunately, we couldn't measure the asymptotic value
at higher temperatures because of the limitation by the trap depth. We suggest
that future experimental studies use a box potential instead of a harmonic
potential to avoid the partial cancellation of mean-field and pair correlation
effects.

Details of the interacting model below the phase transition The full
interacting theory below the phase transition becomes very complex, and we did
not aim for a fully quantitative description. Therefore, we used the following
approximations for the theory curves in Fig. 3b: (1) the density distribution
was obtained using the semi-ideal gas approximation. (2) The effect of the
pair correlation function was introduced by a reduction factor of
$1-8\sqrt{2}{a}/{\lambda_{t}}$ (Eq. S22). Note that the theory curves in Fig.
3b has no free parameters (the overall scaling was the same as Fig. 3a).

| Ideal gas | Mean-field (MF) only | MF + overall suppression | Full interacting theory
---|---|---|---|---
$\chi^{2}$ value | 20.0 | 18.5 | 21.4 | 16.1
Normalized $\chi^{2}$ value | 16.1 | 14.9 | 17.3 | 13
Degrees of freedom (dof) | 13 | 13 | 13 | 13
Table S2:  $\chi^{2}$ values for fits using different models. The $\chi^{2}$
values are normalized in such a way that it becomes $\text{dof}=13$ for the
full interacting model. As a comparison, the probability for $\chi^{2}$
(dof=13) to be larger than 15.1 is 30% while the probability for $\chi^{2}$ to
be larger than 17.0 is 20%.

## References

  * Ketterle  _et al._ [1999] W. Ketterle, D. S. Durfee, and D. Stamper-Kurn, arXiv preprint cond-mat/9904034 (1999).
  * Li  _et al._ [1999] M. Li, L. Chen, and C. Chen, Phys. Rev. A 59, 3109 (1999).
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * Mordini [2019] C. Mordini, Github (2019).
  * Margalit  _et al._ [2021] Y. Margalit, Y.-K. Lu, F. Çağrı Top, and W. Ketterle, Science 374, 976 (2021).
  * Meppelink  _et al._ [2010] R. Meppelink, R. A. Rozendaal, S. B. Koller, J. M. Vogels, and P. van der Straten, Phys. Rev. A 81, 053632 (2010).
  * Naraschewski and Glauber [1999] M. Naraschewski and R. J. Glauber, Phys. Rev. A 59, 4595 (1999).
  * Proukakis  _et al._ [2017] N. P. Proukakis, D. W. Snoke, and P. B. Littlewood,  _Universal themes of Bose-Einstein condensation_ (Cambridge university press, 2017).
  * Ketterle and Miesner [1997] W. Ketterle and H.-J. Miesner, Phys. Rev. A 56, 3291 (1997).


# Distributionally Robust Chance-Constrained Optimization for Hierarchical
UAV-based MEC

Can Cui†, Ziye Jia†, Chao Dong†, Zhuang Ling∗, Jiahao You† and Qihui Wu†
†The Key Laboratory of Dynamic Cognitive System of Electromagnetic Spectrum
Space, Ministry of Industry and Information Technology, Nanjing University of
Aeronautics and Astronautics, Nanjing, Jiangsu, 210016, China
∗College of Communication Engineering, Jilin University, Changchun, Jilin,
130012, China.
{cuican020619, jiaziye, dch, yjiahao<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Multi-access edge computing (MEC) is regarded as a promising technology in the
sixth-generation communication. However, the antenna gain is always affected
by the environment when unmanned aerial vehicles (UAVs) are served as MEC
platforms, resulting in unexpected channel errors. In order to deal with the
problem and reduce the power consumption in the UAV-based MEC, we jointly
optimize the access scheme and power allocation in the hierarchical UAV-based
MEC. Specifically, UAVs are deployed in the lower layer to collect data from
ground users. Moreover, a UAV with powerful computation ability is deployed in
the upper layer to assist with computing. The goal is to guarantee the quality
of service and minimize the total power consumption. We consider the errors
caused by various perturbations in realistic circumstances and formulate a
distributionally robust chance-constrained optimization problem with an
uncertainty set. The problem with chance constraints is intractable. To tackle
this issue, we utilize the conditional value-at-risk method to reformulate the
problem into a semidefinite programming form. Then, a joint algorithm for
access scheme and power allocation is designed. Finally, we conduct
simulations to demonstrate the efficiency of the proposed algorithm.

###### Index Terms:

Unmanned aerial vehicles, multi-access edge computing, distributionally robust
optimization, conditional value-at-risk.

##  I Introduction

Due to the rapid development of smart devices and artificial intelligence,
there exist a growing number of application data to be processed, which is an
outstanding characteristic in the sixth generation of communication system
[1]. Besides, an emergent architecture, multi-access edge computing (MEC),
becomes prevalent in recent years, which contributes to reducing the system
delay and improving the quality of service (QoS) [2]. Furthermore, unmanned
aerial vehicles (UAVs) are well known due to the easy deployment and flexible
movement, which attract tremendous attentions in both academics and
industries. UAVs can serve as communication relay nodes, mobile base stations
(BSs), as well as computing platforms for MEC [3], [4]. However, since there
exist uncertain environmental factors such as wind, temperature, and airflow,
the gains of UAVs are under uncertainty sets due to the environmental factors.
For instance, the effect of wind results in the fluctuation of transmission
gain within a certain range, and the extreme temperature leads to transmission
gain errors. The effects of such unpredictable errors cannot be ignored in
practical applications.

As illustrated that UAV-enabled MEC system can provide benefits for energy-
efficient task offloading, there exist a couple of related works [5]. For
example, [6] studies a UAV-based mobile cloud computing system to minimize
total energy consumption. In [7], Zeng et al. propose an algorithm to maximize
the throughput by jointly optimizing the trajectory of the UAV swarm along
with the transmit power. [8] presents a UAV-assisted MEC for task offloading
which is solved with a deep reinforcement learning method. In [9], a response
delay optimization algorithm is suggested for the deployment of UAVs as well
as MEC servers. However, the uncertain errors of transmission gains are
ignored or described as a certain distribution in the aforementioned works,
which are not realistic in the practical environment. As above, in this paper,
we consider the errors without distribution information to minimize the total
power consumption in the hierarchical UAV-based MEC. Since there exist
unpredictable parameters such as gain errors, which are in an uncertainty set,
the proposed optimization problem is still intractable. Hence, a
distibutionally robust optimization (DRO) problem considering the uncertainty
of gains is further proposed to deal with this issue.

Both the stochastic programming and robust optimization can be applied to
uncertain optimization problems, and the distribution information of random
parameters is necessary, which cannot be obtained in some scenarios[10], [11].
However, DRO can be implemented without such detailed statistical information
and provide more conservative solutions [12]. Besides, without distribution
information on the uncertainty, the DRO problems are usually computationally-
prohibitive. One approach to solving the DRO problem utilizes the conditional
value-at-risk (CVaR) with limited statistical information [13], [14].
Utilizing CVaR with historical data, we formulate the chance-constrained
problem into a semidefinite programming (SDP) form, which can be handled with
feasible solutions [15]. In this paper, we build the uncertainty set for the
error parameters with historical data, enabling more practical applications.
Accordingly, the system model can be formulated as a DRO problem and solved
via the CVaR mechanism.

In detail, we propose a hierarchical UAV-assisted MEC system to complete data
processing and transmit data back to the ground station. The system is
applicable to many real scenarios, such as real-time monitoring and disaster
rescue for remote areas short of BS coverage [16],[17]. Note that the UAVs are
always constrained by energy, the total power consumption should be optimized,
and the QoS of users also should be guaranteed. Taking into account the errors
from antenna gains, we propose a chance-constrained problem with the
uncertainty set and handle it by the CVaR mechanism. The main contributions of
our work are summarized as follows.

  * •

A hierarchical UAV-based MEC system is proposed, including the lower-layer
UAVs to collect data and the upper-layer UAV for relay. Besides, both layers
of UAVs are equipped with computation resources for MEC. Furthermore, the gain
uncertainty from transceivers is considered for practical scenarios, and we
present a corresponding system model to minimize the total power consumption
under such uncertainty.

  * •

The problem is formulated with an uncertainty set for transceiver gain errors,
and to handle this issue, we propose a DRO-based mechanism, which is non-
convex. Then, we approximate the chance constraints by the CVaR method and
reformulate the issue into a tractable SDP form.

  * •

To deal with the reformulated SDP problem, we present an algorithm to jointly
optimize both the access scheme and power allocation. Finally, we conduct
simulations to evaluate the performance of the proposed algorithm, and the
results verify the effectiveness.

The rest of this paper is organized as follows. Section II presents the system
model and corresponding problem formulation. Section III employs the CVaR-
based mechanism to reformulate the original problem with chance constraints
into a SDP form, and the corresponding algorithm is designed. Simulation
results and the analyses are provided in Section IV. Finally, conclusions are
drawn in Section V.

##  II System Model and Problem Fromulation

In this paper, we consider a hierarchical UAV-assisted MEC system, in which
lower-layer UAVs are deployed to collect and compute the data produced by
ground users and transmit them back to the BS. Due to the limited computation
resources of UAVs, a UAV with strong computing ability is deployed in the
upper layer to handle the compute-intensive tasks. Lower-layer UAVs can
complete the computation tasks locally and transmit the results to the BS, or
transmit the data to the upper-layer UAV for further computation to reduce
power consumption. Then, the upper-layer UAV transmits the processed results
to the ground BS. The hierarchical UAV-based MEC scenario is shown in Fig. 1.

###  II-A System Model

####  II-A1 Communication Model

In this subsection, the communication model of UAV-to-UAV and UAV-to-BS is
investigated. We utilize $\mathcal{N}=\\{1,2,\ldots,i,\ldots,N\\}$ to indicate
the set of UAVs of the lower layer. Let _$p^{d}_{i}$_ denote the transmit
power of i-th lower-layer UAV. _g i_ denotes the i-th UAV transceiver antenna
gain and $B_{i}$ indicates the bandwidth. Consequently, the data transmission
rate from the i-th UAV in the lower layer to the upper-layer UAV or the BS is
expressed as:

$r_{i}=B_{i}\log_{2}(1+\frac{{\gamma_{i}p^{d}_{i}{\left|g_{i}\right|}^{2}}}{\sigma_{s}^{2}}),\forall
i\in\\{1,2,\ldots,N\\},$ (1)

where $\sigma_{s}^{2}$ is the variance of white Gaussian noise with a mean
value of $0$. $\gamma_{i}$ is the path loss related to the distance between
i-th lower-layer UAV and upper-layer UAV, or between i-th lower-layer UAV and
BS.

We assume the data transmission adopts the orthogonal multiple access
technology in this system, so the interference between different data can be
ignored. Moreover, only one powerful UAV is considered in the upper layer.

Let $B_{h}$ denote the bandwidth of the upper-layer UAV. _$p^{d}_{h}$_ is the
transmit power, _g h_ denotes the transceiver antenna gain of the UAV hovering
in the upper layer, and $\gamma_{h}$ represents the path loss. According to
the Shannon formula, the data transmission rate of the upper layer UAV towards
the BS is calculated as:

$r_{h}=B_{h}\log_{2}(1+\frac{{\gamma_{h}p^{d}_{h}{\left|g_{h}\right|}^{2}}}{\sigma_{s}^{2}}).$
(2)

Figure 1: A hierarchical UAV-based MEC scenario.

We assume _$L_{i}$_ represents the data length transmitted by the i-th UAV.
_$x_{i}\in\\{0,1\\}$_ is a binary indicator for computing mode. _$x_{i}=1$_
represents the i-th UAV transmits the users’ data to the upper-layer UAV for
further computation, and _$x_{i}=0$_ denotes the i-th UAV computes collected
data locally and transmit them back to the BS.

Based on the above formulas, we can obtain the transmission delay as follows:

$t_{i}^{d}=\frac{L_{i}}{r_{i}},\forall i\in\\{1,2,\ldots,N\\},$ (3)

and

$t_{h}^{d}=\frac{\sum_{i}^{N}x_{i}L_{i}}{r_{h}},$ (4)

where $t_{i}^{d}$ denotes the transmission delay of the i-th UAV in the lower
layer, and $t_{h}^{d}$ is the transmission delay of the upper-layer UAV,
respectively. Further, the total transmission power $p_{r}^{d}$ of UAVs in the
lower layer is

$p_{r}^{d}=\sum_{i=1}^{N}p_{i}^{d}.$ (5)

####  II-A2 Uncertainty Model

Take into account the transceiver antenna gain errors influenced by the
environmental factors denoted as $\Delta g_{i}$ and $\Delta g_{h}$, i.e.,

$g_{i}=\overline{g}_{i}+\Delta g_{i},\forall i\in\\{1,2,\ldots,N\\},$ (6)

and

$g_{h}=\overline{g}_{h}+\Delta g_{h},$ (7)

where $\overline{g}_{i}$ and $\overline{g}_{h}$ indicate the theoretical
antenna transmission gain of the i-th lower-layer UAV and upper-layer UAV,
respectively. Similarly, $\Delta{g_{i}}$ and $\Delta{g_{h}}$ represent the
antenna gain errors produced by actual environmental impacts. The estimation
error $\Delta{g_{i}}$ and $\Delta{g_{h}}$ are assumed to follow an unknown
distribution $\mathbb{P}$ with mean $\mathbb{E_{P}}(\Delta g_{i})=\mu_{i}$ and
$\mathbb{E_{P}}(\Delta g_{h})=\mu_{h}$, as well as the variance
$\mathbb{D_{P}}(\Delta g_{i})=\sigma_{i}^{2}$ and $\mathbb{D_{P}}(\Delta
g_{h})=\sigma_{h}^{2}$, respectively. Uncertainty set $\mathcal{P}$
reprensents for all possible probability distributions of $\Delta{g_{i}}$ and
$\Delta{g_{h}}$, i.e., $\mathbb{P}\in\mathcal{P}$.

####  II-A3 Computation Model

Let variable _$c_{i}$_ represent the complexity of the computing tasks
processed by the i-th UAV. _$\eta_{i}$_ and _$\eta_{h}$_ denote the
computational power consumption coefficient of the i-th lower-layer UAV and
the upper-layer UAV, respectively, which represent the power required to
calculate the data of each cycle. Then we can formulate the computation power
of i-th lower-layer UAV and the upper-layer UAV as:

$p_{i}^{c}=\eta_{i}c_{i},\forall i\in\\{1,2,\ldots,N\\},$ (8)

and

$p_{h}^{c}=\eta_{h}\sum_{i=1}^{N}x_{i}c_{i}.$ (9)

With the above formulations we can obtain the total computation power
consumption of UAVs in the lower layers as:

$p_{r}^{c}=\sum_{i=1}^{N}(1-x_{i})p_{i}^{c}.$ (10)

####  II-A4 Power Consumption

According to the above discussion, ignoring the power required for UAV
hovering and flying, which is a constant value, the total power is expressed
as:

$P_{total}=p_{r}^{d}+p_{h}^{d}+p_{r}^{c}+p_{h}^{c}.$ (11)

###  II-B Problem Formulation

In this part, we focus on reducing the total power in the system for a better
service. The target is to minimize $P_{total}$ by jointly optimizing the UAV
access scheme and power allocation when the QoS meets users’ demands. Hence,
the problem is formulated as

$\displaystyle\textbf{P0}:$
$\displaystyle\min_{x_{i},p_{i}^{d},p_{h}^{d}}P_{total},$ (12)
$\displaystyle\textrm{s.t.}\quad\mathit{C}1:\mathbf{Pr}\\{t_{i}^{d}\leq
t_{i,max}^{d}\\}\geq\alpha_{1},\forall i,$ (13)
$\displaystyle\qquad\mathit{C}2:\mathbf{Pr}\\{t_{h}^{d}\leq
t_{h,max}^{d}\\}\geq\alpha_{2},$ (14)
$\displaystyle\qquad\mathit{C}3:0\leq\sum_{i=1}^{N}x_{i}\leq m,$ (15)
$\displaystyle\qquad\mathit{C}4:x_{i}\in\\{0,1\\},\forall i,$ (16)
$\displaystyle\qquad\mathit{C}5:p_{i}^{d}\geq 0,\forall i,$ (17)
$\displaystyle\qquad\mathit{C}6:p_{h}^{d}\geq 0,$ (18)

where $\mathit{C}1$ and $\mathit{C}2$ are the chance constraints, indicating
that the transmission delays of UAVs are constrained in a probabilistic
manner. In other words, the delays should not be larger than the tolerable
delay $t_{i,max}^{d}$ and $t_{h,max}^{d}$ with probabilities of $\alpha_{1}$
and $\alpha_{2}$ at a minimum, respectively. Futhermore, $\alpha_{1}$,
$\alpha_{2}$ $\in$ $(0,1)$. $\mathit{C}3$ denotes that the upper-layer UAV
with powerful computation ability can at most access $\mathit{m}$ UAVs at the
same time. $\mathit{C}4$ denotes $\mathit{x_{i}}$ is a binary variable.
$\mathit{C}5$ and $\mathit{C}6$ guarantee that $p_{i}^{d}$ and $p_{h}^{d}$ are
non-negative continuous variables.

Due to the chance contraints (13) and (14), the transmission delay $t_{i}^{d}$
and $t_{h}^{d}$ in P0 are difficult to obtain. However, there exists one
effective way to handle the problem with distributionally robust method.
Specifically, let $\underset{\mathbb{P}\in\mathcal{P}}{inf}$ represent the
lower bound of the probability under the probability distribution
$\mathbb{P}$, and (13) and (14) can be rewritten by:

$\underset{\mathbb{P}\in\mathcal{P}}{inf}\quad\mathbf{Pr}_{\mathbb{P}}\\{t_{i}^{d}\leq
t_{i,max}^{d}\\}\geq\alpha_{1},\forall i\in\\{1,2,\ldots,N\\},$ (19)

and

$\underset{\mathbb{P}\in\mathcal{P}}{inf}\quad\mathbf{Pr}_{\mathbb{P}}\\{t_{h}^{d}\leq
t_{h,max}^{d}\\}\geq\alpha_{2},$ (20)

respectively, which are distributionally robust chance constraints (DRCCs).

##  III Reformulation and Algorithm

Since the DRCC problems (19) and (20) are intractable with the uncertainty
set, in this section, we employ the CVaR mechanism to handle the issue and
then design an algorithm to obtain the final solution.

###  III-A CVaR Based Mechanism

Generally, the value-at-risk (VaR) of a variable $u$ with the safety factor
$\alpha$ is defined as the minimal value of $v$, and $u$ is no more than $v$
at a possibility of $\alpha$ [18], i.e.,

$VaR_{\alpha}(u)=\min\\{v|P(u\leq v)\geq\alpha\\}.$ (21)

It is noted that VaR is non-convex and discontinuous. Based on VaR, we can
consequently propose the definition of CVaR, which is defined as the
conditional expectation of $u$ when $u\geq VaR_{\alpha}(u)$, i.e.,

$CVaR_{\alpha}(u)=E[u|u\geq VaR_{\alpha}(u)].$ (22)

The relationship between VaR and CVaR is shown in Fig. 2. It is obvious that
$CVaR_{\alpha}(u)\geq VaR_{\alpha}(u)$ [12]. Furthermore, CVaR is a
conservative approximate estimation of loss, which is robust. As is proposed
in [12], for a given measurable loss function $\varphi$($\xi$):
$\mathbb{R}^{k}\rightarrow\mathbb{R}$, CVaR under the safety factor $\alpha$
concerning the probability distribution $\mathbb{P}$ on $\mathbb{R}^{k}$ is
expressed as:

Figure 2: Relationship between VaR and CVaR.

$\mathbb{P}-CVaR_{\alpha}(\varphi(\xi))=\underset{\beta\in\mathbb{R}}{inf}\\{\beta+\frac{1}{1-\alpha}\mathbb{E_{P}}\left[\max\left(0,\varphi(\xi)-\beta\right)\right]\\}.$
(23)

As is shown in Fig. 2, according to the definition, we can then obtain the
following formula:

$\mathbb{P}\\{\varphi(\xi)\leq\mathbb{P}-CVaR_{\alpha}(\varphi(\xi))\\}\geq\alpha.$
(24)

As above, the CVaR constraint can be formed as:

$\begin{split}\underset{\mathbb{P}\in\mathcal{P}}{sup}\quad\mathbb{P}-CVaR_{\alpha}(\varphi(\xi))\leq
0,\forall\mathbb{P}\in\mathcal{P}\Leftrightarrow\\\
\underset{\mathbb{P}\in\mathcal{P}}{inf}\quad\mathbb{P}\\{\varphi(\xi)\leq
0\\}\geq\alpha.\end{split}$ (25)

Formula (25) represents that the conservative approximation of DRCC on the
right side can be constituted by the CVaR constraint on the left side. In
other words, CVaR represents a convex function of the random parameter. Then,
we can reformulate the CVaR constraint into a tractable SDP form.

Lemma 1: Let
$\varphi(\xi)=\xi^{\top}\mathbf{\Theta}\xi+\bm{\theta}^{\top}\xi+\bm{\theta}^{0}$
for $\bm{\Theta}\in\mathbb{S}^{k}$, $\bm{\theta}\in\mathbb{R}^{k}$, and
$\bm{\theta}^{0}\in\mathbf{R}$, and then, the worst-case CVaR
$\underset{\mathbb{P}\in\mathcal{P}}{sup}\quad\mathbb{P}-CVaR_{\alpha}(\varphi(\xi))$
can be rewritten into a tractable SDP as follows:

$\begin{split}&inf\quad\beta+\frac{1}{1-\alpha}\langle{\bm{\Omega},\bm{H}}\rangle,\\\
\textrm{s.t.}\quad&\bm{H}\in\mathbb{S}^{k+1},\quad\beta\in\mathbb{R},\\\
&\bm{H}\succeq\mathbf{0},\\\
&\bm{H}-\left[\begin{matrix}\bm{\Theta}&\frac{1}{2}\bm{\theta}\\\
\frac{1}{2}\bm{\theta}^{\top}&\bm{\theta}^{0}-\beta\\\
\end{matrix}\right]\succeq\mathbf{0},\end{split}$ (26)

where $\mathbf{\Omega=\left[\begin{matrix}\Sigma+\mu\mu^{\top}&\mu\\\
\mu^{\top}&1\\\ \end{matrix}\right]}$, $\Sigma\in\mathbb{S}^{k}$ is the
covariance of the random vector $\xi$, and $\mu\in\mathbb{R}^{k}$ is the mean
matrix. Furthermore, $\bm{H}$ is an auxiliary matrix and $\beta$ is an
auxiliary variable.
$\langle{\bm{\Omega},\bm{H}}\rangle=\mathrm{Tr}(\bm{\Omega}\bm{H})$ represents
the trace scalar product. $\bm{H}\succeq\mathbf{0}$ represents that the matrix
$\bm{H}$ is positive semidefinite[15].

Therefore, according to Lemma 1, we transform the DRO into the form of SDP by
the CVaR mechanism, which is a tractable and computationally efficient
approximation technique. This form can be regarded as an extension of linear
programming, and the effectiveness of this solution is proved in previous
works [19].

###  III-B Problem Reformulation

As for (19), a DRCC problem, the first-order Taylor expansion of the loss
function is

$\varphi(\xi_{1})=L_{i}\sigma_{s}^{2}-B_{i}\gamma_{i}p^{d}_{i}t_{i,max}^{d}\xi_{1}^{2},$
(27)

where $\xi_{1}$ denotes the random parameter $\left|g_{i}\right|$.

CVaR can construct a convex approximation for the chance constraints.
According to Lemma 1, we transform (27) into a SDP problem, i.e.,

$\begin{split}&inf\quad\beta_{1}+\frac{1}{1-\alpha_{1}}\langle{\bm{\Omega_{1}},\bm{H_{1}}}\rangle\leq
0,\\\
\textrm{s.t.}\quad&\bm{H_{1}}\in\mathbb{S}^{2},\quad\beta_{1}\in\mathbb{R},\\\
&\bm{H_{1}}\succeq\mathbf{0},\\\
&\bm{H_{1}}-\left[\begin{matrix}-B_{i}\gamma_{i}p_{i}^{d}t_{i,max}^{d}&0\\\
0&L_{i}\sigma_{s}^{2}-\beta_{1}\\\
\end{matrix}\right]\succeq\mathbf{0},\end{split}$ (28)

where $\bm{H_{1}}$ and $\beta_{1}$ are both auxiliary variables. Let
$\mu_{g_{i}}=\overline{g}_{i}+\mu_{i}$ present the mean of $g_{i}$, we have

$\Omega_{1}=\left[\begin{matrix}\sigma_{i}^{2}+\mu_{g_{i}}\mu_{g_{i}}^{\top}&\mu_{g_{i}}\\\
\mu_{g_{i}}^{\top}&1\\\ \end{matrix}\right].$ (29)

Likewise, let $\xi_{2}$ denote the random parameter $\left|g_{h}\right|$, and
the loss function in (20) can be rewritten as:

$\varphi(\xi_{2})=\sum_{i=1}^{N}x_{i}L_{i}\sigma_{s}^{2}-B_{h}\gamma_{h}p^{d}_{h}t_{h,max}^{d}\xi_{2}^{2},$
(30)

which can be further presented in a SDP form, i.e.,

$\begin{split}&inf\quad\beta_{2}+\frac{1}{1-\alpha_{2}}\langle{\bm{\Omega_{2}},\bm{H_{2}}}\rangle\leq
0,\\\
\textrm{s.t.}\quad&\bm{H_{2}}\in\mathbb{S}^{2},\quad\beta_{2}\in\mathbb{R},\\\
&\bm{H_{2}}\succeq\mathbf{0},\\\
&\bm{H_{2}}-\left[\begin{matrix}-B_{h}\gamma_{h}p^{d}_{h}t_{h,max}^{d}&0\\\
0&\sum_{i=1}^{N}x_{i}L_{i}\sigma_{s}^{2}-\beta_{2}\\\
\end{matrix}\right]\succeq\mathbf{0}.\end{split}$ (31)

Similarly, $\bm{H_{2}}$ is the auxiliary matrix, and $\beta_{2}$ is the
auxiliary variable for the CVaR constraint.
$\mu_{g_{h}}=\overline{g}_{h}+\mu_{h}$ is the mean value of $g_{h}$. Then we
obtain

$\Omega_{2}=\left[\begin{matrix}\sigma_{h}^{2}+\mu_{g_{h}}\mu_{g_{h}}^{\top}&\mu_{g_{h}}\\\
\mu_{g_{h}}^{\top}&1\\\ \end{matrix}\right].$ (32)

Finally, based on above discussions, the chance-constrained problem (12) under
the uncertainty set can be reformulated into a SDP form, i.e.,

$\begin{split}\textbf{P1}:\quad&\min_{\begin{subarray}{c}x_{i},p^{d}_{i},p^{d}_{h},\\\
\bm{H_{1}},\bm{H_{2}},\beta_{1},\beta_{2}\end{subarray}}P_{total},\\\
&\mathbf{\textrm{s.t.}}\quad(\ref{C3})-(\ref{C6}),(\ref{SDP-1}),(\ref{SDP-2}),\\\
&\forall i\in\\{1,2,\ldots,N\\},\end{split}\quad$

which is still non-convex with the binary variable $x_{i}$, and continuous
variables $p^{d}_{i}$ and $p^{d}_{h}$. The problem is in the form of mixed
integer non-linear programming. To tackle with this problem, we design a joint
optimization algorithm on access scheme and power allocation, as shown in
Algorithm 1.

###  III-C Algorithm Design

It is observed that due to the two decision variables $x$ and $p^{d}$, problem
P1 is still non-convex. To solve this problem, we divide it into two
subproblems. The first subproblem is concerning the access scheme $x$, and the
second subproblem is related to the transmission power allocation plan
$p^{d}$.

As shown in step 2 of Algorithm 1, in the first subproblem, an initial value
$p^{d}$ is assumed to minimize the computation power consumption and ensure
the access scheme, i.e., $x$. In other words, the lower-layer UAVs make a
choice to complete the computation work locally or transform it to the upper-
layer UAV based on the task complexity. The problem with integers is solved by
MOSEK. Then, in the second subproblem, according to the access scheme $x^{*}$
that is obtained in step 2, we can further substitute it to solve the
transmission power allocation program. Specifically, the subproblem can be
tackled by fixing the access scheme, i.e., $x^{*}$. Then, in step 3 of
Algorithm 1, the transmission power is allocated and variables $p^{d}_{i}$ as
well as $p^{d}_{h}$ are obtained. Repeat the process until the solutions of
two subproblems converge. The final result, i.e., the total power consumption
$P_{total}$ in the hierarchical UAV-assisted MEC system, can be figured out as
$\bm{optval}$.

Algorithm 1 Joint Optimization Algorithm on Access Sheme and Power Allocation.

0: Network parameters, $\bm{c_{i}}$, $\bm{L_{i}}$.

0: $\bm{x_{i}}$, $\bm{p^{d}_{i}}$, $\bm{p^{d}_{h}}$ and minimum system power
consumption $\bm{optval}$.

1: repeat

2: Solve the subproblem of access-scheme on $\bm{x_{i}}$ to obtain
$\bm{x_{i}^{*}}$ and the power consumption $\bm{\tau_{1}}$. $\bm{x_{i}}$
$\rightarrow$ $\bm{x_{i}^{*}}$.

3: Solve the subproblem of power $\bm{p^{d}}$ to obtain $\bm{p^{d*}}$ and the
power consumption $\bm{\tau_{2}}$. $\bm{p^{d}}$ $\rightarrow$ $\bm{p^{d*}}$.

4: until Results converge.

5: Obtain the final total power $\bm{optval}$.

##  IV Simulation Results

Numeral simulations are conducted to evaluate the performance of the proposed
algorithm. Briefly speaking, the data lengths are assumed to range from
$40kbits$ to $60kbits$. $B_{i}$ and $B_{h}$ are assumed as $10MHz$ and
$10MHz$, respectively. The task complexities range from $40kcycles$ to
$100kcycles$. The CPU computational power consumption coefficient $\eta_{i}$
is $5\times 10^{-6}cycles/W$, and $\eta_{h}$ is $10^{-6}cycles/W$. Antenna
gains $\left|g_{i}\right|$ and $\left|g_{h}\right|$ are $5$. The mean value of
errors is $0$ and the variance is $0.01$ for all. Furthermore, noise
$\sigma_{s}^{2}$ is assumed as $10^{-12}W$. The safety factors are
$\alpha_{1}=0.95$ and $\alpha_{2}=0.95$.

To verify the effectiveness and feasibility, we compare the results of the
CVaR method and non-robust method without uncertainty errors $\Delta{g_{i}}$
and $\Delta{g_{h}}$. In detail, Fig. 3 shows the performance of two methods in
total power v.s. the number of lower-layer UAVs. As the number of lower-layer
UAVs increases, the total power consumption of the MEC system increase as
well. It can be explained that as the scale of the whole problem becomes
larger, the users’ data collected are larger to be computed and transmitted.
Comparing the results of the two methods, it is observed the power consumption
in the CVaR method is smaller than that in the non-robust method. It is
explained that there is a tolerance in the CVaR method, so that the
transmission power of each UAV can be reduced within the tolerant range.
Moreover, different from the non-robust method, which optimizes the problems
under the perfect assumption which is apparently inconsistent with the real
circumstance, the designed algorithm takes into account the uncertainties in
reality and is applicable to the real situation.

Figure 3: System power v.s. the number of lower-layer UAVs for different
methods. Figure 4: Individial transmission power of the i-th UAV in the lower
layer for different methods.

Fig. 4 shows the individual transmission power of each lower-layer UAV when
there exist six UAVs in the lower layer. We compare the result with the non-
robust method. For tasks varying in data length, it demonstrates that the CVaR
method performs better in transmission power consumption. The conclusion can
be drawn that the proposed method, which is considered in a more realistic
world, can obtain a better result than that in the ideal situation.
Specifically, considering the errors between the real situation and the ideal
situation, the designed algorithm only needs to obtain the mean and variance
of the gain errors from the historical data, instead of obtaining the true
distribution or probability function of the errors. The experimental results
verify the effectiveness of the proposed joint optimization algorithm on
access scheme and power allocation.

In Fig. 5, we investigate the performance of total power consumption, i,e.,
the corresponding system power v.s. lower-layer UAV numbers for different
maximum transmission delays. It depicts the influence of the maximum tolerance
on the total power. With the increment of maximum delay to be satisfied, the
transmission power consumption reduces. It is accounted for the fact that
transmission rates become lower for all UAVs when the tolerance transmission
delay increases.

Figure 5: System power v.s. the number of lower-layer UAVs for different
maximum transmission delays.

##  V Conclusions

In this paper, we study a hierarchical UAV-assisted MEC scenario to optimize
the power consumption. To tackle with the potential uncertainty impacted by
environmental factors, which is not ignorable in practical circumstances, a
DRO problem based on the uncertainty set is proposed. Then, with the CVaR
mechanism, the original problem is reformulated into a SDP form. To jointly
optimize the reformulation on the access scheme and power allocation, we
further design an algorithm to handle it and obtain the final solution. By
conducting simulation experiments, the robustness and feasibility of the
proposed algorithm compared with the non-robust method are verified.

## References

  * [1] S. Chen, Y.-C. Liang, S. Sun, S. Kang, W. Cheng, and M. Peng, “Vision, requirements, and technology trend of 6G: How to tackle the challenges of system coverage, capacity, user data-rate and movement speed,” _IEEE Wirel Commun_ , vol. 27, no. 2, pp. 218–228, Apr. 2020.
  * [2] N. Abbas, Y. Zhang, A. Taherkordi, and T. Skeie, “Mobile edge computing: A survey,” _IEEE Internet Things J._ , vol. 5, no. 1, pp. 450–465, Feb. 2018\.
  * [3] C. Dong, Y. Shen, Y. Qu, K. Wang, J. Zheng, Q. Wu, and F. Wu, “UAVs as an intelligent service: Boosting edge intelligence for air-ground integrated networks,” _IEEE Network_ , vol. 35, no. 4, pp. 167–175, Jul. 2021.
  * [4] Z. Jia, M. Sheng, J. Li, D. Niyato, and Z. Han, “LEO-satellite-assisted UAV: Joint trajectory and data collection for internet of remote things in 6G aerial access networks,” _IEEE Internet Things J._ , vol. 8, no. 12, pp. 9814–9826, Jun. 2021.
  * [5] J. You, Z. Jia, C. Dong, L. He, Y. Cao, and Q. Wu, “Computation offloading for uncertain marine tasks by cooperation of UAVs and vessels,” _arXiv e-prints_ , p. arXiv:2209.07228, Feb. 2023. [Online]. Available: https://doi.org/10.48550/arXiv.2302.06055
  * [6] S. Jeong, O. Simeone, and J. Kang, “Mobile edge computing via a UAV-mounted cloudlet: Optimization of bit allocation and path planning,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 3, pp. 2049–2063, Mar. 2018.
  * [7] Y. Zeng, R. Zhang, and T. J. Lim, “Throughput maximization for UAV-enabled mobile relaying systems,” _IEEE Trans Commun_ , vol. 64, no. 12, pp. 4983–4996, Dec. 2016.
  * [8] Y. M. Park, S. Salman Hassan, Y. Kyaw Tun, Z. Han, and C. S. Hong, “Joint trajectory and resource optimization of MEC-assisted UAVs in sub-THz networks: A resources-based multi-agent proximal policy optimization DRL with attention mechanism,” _arXiv e-prints_ , p. arXiv:2209.07228, Sep. 2022. [Online]. Available: https://arxiv.org/abs/2209.07228v1
  * [9] Q. Zhang, J. Chen, L. Ji, Z. Feng, Z. Han, and Z. Chen, “Response delay optimization in mobile edge computing enabled UAV swarm,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 3, pp. 3280–3295, Mar. 2020.
  * [10] J. R. Birge and F. Louveaux, _Introduction to Stochastic Programming_. Springer New York, 2011.
  * [11] B.-T. Aharon and N. Arkadi, _Robust Optimization_. Princeton University Press, 2009.
  * [12] R. T. Rockafellar and S. Uryasev, “Optimization of conditional value-at-risk,” _J RISK_ , vol. 2, no. 3, pp. 21–41, Sep. 2000.
  * [13] Z. Ling, F. Hu, Y. Zhang, L. Fan, F. Gao, and Z. Han, “Distributionally robust chance-constrained backscatter communication-assisted computation offloading in wbans,” _IEEE Transactions on Communications_ , vol. 69, no. 5, pp. 3395–3408, May 2021.
  * [14] Z. Ling, F. Hu, H. Zhang, and Z. Han, “Age of information minimization in healthcare IoT using distributionally robust optimization,” _IEEE Internet Things J._ , vol. 9, no. 7, pp. 16 154–16 167, Feb. 2022.
  * [15] Steve, Zymler, Daniel, Kuhn, Ber, and Rustem, “Distributionally robust joint chance constraints with second-order moment information,” _Math Program_ , vol. 137, no. 1-2, pp. 167–198, Feb. 2011.
  * [16] Z. Jia, Q. Wu, C. Dong, C. Yuen, and Z. Han, “Hierarchical aerial computing for internet of things via cooperation of HAPs and UAVs,” _IEEE Internet Things J._ , Feb. 2022, early access.
  * [17] Z. Jia, M. Sheng, J. Li, D. Zhou, and Z. Han, “Joint HAP access and LEO satellite backhaul in 6G: Matching game-based approaches,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 4, pp. 1147–1159, Apr. 2021.
  * [18] S. Sarykalin, G. Serraino, and S. Uryasev, _Value-at-Risk vs. Conditional Value-at-Risk in Risk Management and Optimization_ , ch. Chapter 13, pp. 270–294. [Online]. Available: https://pubsonline.informs.org/doi/abs/10.1287/educ.1080.0052
  * [19] L. Vandenberghe and S. Boyd, “Semidefinite programming,” _SIAM Review_ , vol. 38, no. 1, pp. 49–95, 1996. [Online]. Available: https://doi.org/10.1137/1038003