# WKB approximation to bosonic dark matter

Lauren Street<EMAIL_ADDRESS>Department of Physics, University of
Cincinnati  Peter Suranyi<EMAIL_ADDRESS>Department of Physics,
University of Cincinnati  L.C.R. Wijewardhana<EMAIL_ADDRESS>Department of Physics, University of Cincinnati

###### Abstract

Galactic dark matter halos may be composed of ultralight axions (ULAs)
($m_{a}\lesssim 1$ eV) with wave functions that satisfy nonlinear Schrödinger-
Poisson equations (SPA). We find eigenstates of SPA in WKB approximation. The
expansion parameter of the WKB approximation is $\delta=1/\sqrt{S}$, where
$S=2MRGm_{a}^{2}$, with $M$ being the total mass, $R$ the radius of the halo,
and $G$ the gravitational constant. $S\gg 1$ for almost all galaxies, even if
the ULA mass is as small as $m_{a}=10^{-22}$ eV, making the leading order WKB
approximation almost exact. As the level spacing of bound states is roughly
proportional to $\delta$, the number of states in the gravitational well is
huge. We do not see a reason why not all or most of them contribute to the
halo. Using an appropriate distribution function allows the summation of
states to construct the profile of the halo as a function of the gravitational
potential, which can be found solving the Poisson equation. Using various
energy distribution functions, we obtain results similar to those in
simulations. Future plans include investigations of collapse through time
dependent generalizations, and inclusion of self-interactions, which also
induce decay processes of the halo.

The structure of galaxies and the rotation curves of stars in galaxies can
potentially be explained with the assumption that most of galactic matter is
composed of presently unknown particles, termed dark matter (DM), which
interact very weakly with particles of the Standard Model. One of the most
popular variants of DM is the weakly interacting massive particle (WIMP),
consisting of massive, non-relativistic particles, heavier than neutrinos
1982ApJ…263L…1P ; Bond:1982uy ; Blumenthal:1982mv . Since no such particles,
in the appropriate mass range, have been discovered yet, other alternatives
for DM have also been considered. Among others, a prominent candidate is
ultralight axions (ULAs) with Compton wavelengths ranging from cosmic size
Baldeschi:1983mq ; Turner:1983he ; Press:1989id ; Sin:1992bg ; Hu:2000ke ;
Guzman:2004wj ; Boehmer:2007qa ; Weinberg:2013aya ; Marsh:2015wka ;
Hui:2016ltb ; Lee:2017qve  to that of masses $m_{a}\sim 1$ eV Ferreira:2020fam
.

There have been many simulations of the collapse of ULAs on galactic scales
without Schive:2014dra ; Schive:2014hza ; Schwabe:2016rze ; Veltmaat:2016rxo ;
Du:2016aik ; Mocz:2017wlg ; Levkov:2018kau  and with baryonic feedback
Mocz:2019pyf . In both cases, ULA systems were shown to collapse to a
condensed core, surrounded by a virialized halo of non-relativistic ULA.
Recently, there have also been simulations performed for ULAs with self-
interactions, Amin:2019ums ; Glennon:2020dxs ; Mocz:2023adf , for systems
composed of multiple flavors without self-interactions Huang:2022ffc , and for
systems composed of multiple flavors with self-interactions
PhysRevD.107.063520 . In order to study dynamic perturbations to the central
soliton, multiple eigenstates of a Schrodinger Poisson system where the
gravitational potential for higher modes is generated by the solitonic ground
state was performed in reference Zagorac:2021qxq . These states can be used to
construct halos of low mass galaxies where due the the low value of S only a
small number of states contribute to the halo.

In another work Lin:2018whl , self-consistent simulations of the halo
constructed from excited states of the Schrödinger-Poisson (SP) equations were
performed. The authors found that the collapse of the system consisted of a
condensed soliton core surrounded by a halo composed of excited eigenstates.
In a subsequent work PhysRevD.101.081302 , systems satisfying SP equations
were analyzed assuming composition of a small number of energy eigenstates,
including the stability and virialization of the system.

The purpose of this work is to construct DM from self-adjoint or complex ULAs
with self-interactions using WKB approximation. We ignore self-interactions in
solving the equations of motion, but in a subsequent work we will consider the
effect of $2\rightarrow 2$ interactions on the stability of excited
eigenstates. We emphasize that, because only this particular interaction is
relevant, our model can be used for real, or just as easily for complex,
scalar fields. For the sake of simplicity, we focus on a real scalar field
giving rise to a ULA subject to a $\Phi^{4}$ self-interaction. Such an
interaction is the leading order expansion term of an axion potential,
$V=m_{a}^{2}f^{2}[1-\cos(\Phi/f)]$. For the range of galactic sizes and ULA
masses considered here, the contribution of self-interaction terms to the
equations of motion is negligible compared to that of the gravitational
interaction. The ratio of self-interactions to gravitational interactions
scales as

$\frac{\text{S}I}{GI}\sim\frac{M_{P}^{2}}{f_{a}^{2}}\frac{1}{m_{a}^{2}R^{2}},$

where $M_{P}=G^{-1/2}$ is the Planck mass and $G$ is Newton’s constant,
$f_{a}$ is the axion decay constant, $R$ is the radial scale of the system and
$m_{a}$ is the mass of the axion. Using $m_{a}=10^{-22}eV$, $f/M_{p}=10^{-3}$
for the Milky Way with $R=10^{5}$ light year, we obtain $SI/GI\simeq 10^{-4}$.

Self-interactions of ULAs may be important for extremely small galaxies. In
fact, based on studies of axion stars PhysRevD.84.043531 ; Chavanis:2011zm ;
Eby:2015hyx  they can possibly generate a cutoff in the mass of small, stable
galaxies with very large densities. That possibility will be investigated in
future publications.

##  I Equations of motion

Our basic assumption is that radial eigenfunctions of the halo satisfy a
Schrodinger-Gross-Pitaevskii equation

$\displaystyle
E_{nl}\psi_{nl}=-\frac{1}{2m_{a}}\left[\psi_{nl}^{\prime\prime}+\frac{2}{r}\psi_{nl}^{\prime}\right]+\left[\frac{1}{2m_{a}}\frac{l(l+1)}{r^{2}}+V_{g}\right]\psi_{nl},$
(I.1)

where quantum numbers $n$ and $l$ characterize eigenstates, and $V_{g}$ is the
gravitational potential. The most general stationary state wave function of
the halo is

$\displaystyle\Psi({\bf
r},t)=\sum_{nlm}\psi_{nl}(r)Y_{lm}(\theta,\phi)e^{i(E_{nl}t+\,\alpha_{nlm})},$
(I.2)

where $\alpha_{nlm}$ are random phases. Then the gravitational potential is

$\displaystyle V_{g}$ $\displaystyle=-G\,m_{a}\int
d^{3}r^{\prime}\frac{\Psi(\bf{r}^{\prime})^{2}}{|\bf{r}-\bf{r}^{\prime}|}\simeq-G\,m_{a}\sum_{nl}(2l+1)\int
d^{3}r^{\prime}\frac{\psi_{nl}(r^{\prime})^{2}}{|\bf{r}-\bf{r}^{\prime}|},$
(I.3)

where we average over time and random phases, and assume spherical symmetry.

The main result of this paper is the analytic derivation of the density
distribution of a ULA halo. Suppose the gravitational potential of the system,
$V_{g}(r)$, is known. Then the density distribution satisfies the Poisson
equation

$\nabla^{2}V_{g}(r)=-4\pi\,G\,m_{a}\rho(r).$ (I.4)

If we can construct the density distribution as a function of $V_{g}$, then
Eq. (I.4) constitutes a second order differential equation for $V_{g}$ which
can be solved, providing the density distribution. This density distribution
can then be compared to simulations or observations.

Expanding the wave function in radial coordinates, assuming again that
interference terms are negligible when taking averages and the density
distribution is spherically symmetric, the number density can be written as

$\rho(r)=m_{a}\sum_{nl}(2l+1)|\psi_{nl}|^{2},$ (I.5)

where we normalize wave functions as

$\int d^{3}r|\psi_{nl}|^{2}=N_{nl},$ (I.6)

where $N_{nl}\equiv M_{nl}/m_{a}$ and $M_{nl}$ is the total mass of states
having quantum numbers $n$ and $l$.

In phenomenological models (Navarro:1995iw ; Burkert:1995yz ) the central
density and radial scale are undetermined free parameters. To compare various
radial scales we introduce the universal radial scale parameter $R$, the
rescaled dimensionless coordinate, $z=r/R$, and the gravitational scaling
function, $w(z)$, through the equation

$\displaystyle V_{g}(z)=-G\frac{m_{a}^{2}}{R}\int
d^{3}z^{\prime}\frac{\tilde{\rho}(z^{\prime})}{|z-z^{\prime}|}=-\frac{GM\,m_{a}}{R}w(z),$
(I.7)

where we rescale the density as
$\tilde{\rho}(z)=\left(R^{3}/m_{a}\right)\rho(r)$. Choosing $R$ as the
harmonic average of $r$ weighted over the density,

$\displaystyle\frac{1}{R}=\left\langle\frac{1}{r}\right\rangle,$ (I.8)

the gravitational potential at the center is $V_{g}(0)=-GM\,m_{a}/R$. Then Eq.
(I.7) implies that $1\geq w(z)>0$, with $w(0)=1$.

Phenomenological models of halos are of the form $\rho(r)=\rho(0)f(r/R_{s})$,
where $R_{s}$ is a radial scale, differing from the one defined in Eq. (I.8).
The concentration is defined as $c=r_{\rm vir}/R_{s}$, the virial radius,
$r_{\rm vir}$, is usually defined as the radius of a sphere containing the
total mass of the halo and densities are cut off at $r_{\rm vir}$.

Using our method, all densities vanish at some finite value of $z=z_{\rm vir}$
which we identify with the scaled virial radius, $z_{\rm vir}=r_{\rm vir}/R$ .
In phenomenological models, the scales of models, e.g. concentration, are not
universal between models. As Eq. (I.8) implies $\langle 1/z\rangle=1$, we can
give a universal definition to the concentration $c=z_{\rm vir}$.

Note that $R$ is of the same order of magnitude but larger than the scaling
parameter of phenomenological models (Navarro:1995iw ; Burkert:1995yz ) and
unlike those, model independent. That can be shown by calculating the harmonic
average, Eq. (I.8), in those models. The harmonic radius of a NFW halo is

$R=R_{\rm NFW}\left[\left(1+\frac{1}{c_{\rm NFW}}\right)\log(1+c_{\rm
NFW})-1\right]>R_{\rm NFW},$ (I.9)

where $R_{\rm NFW}$ is the scaling radius of the NFW halo, and $c_{\rm NFW}$
is the NFW concentration. Similarly, the harmonic radius of the Burkert halo
is

$R=R_{\rm B}\frac{\log(1+c_{\rm B}^{2})+2\log(1+c_{\rm B})-2\tan^{-1}(c_{\rm
B})}{\log(1+c_{\rm B}^{2})-2\log(1+c_{\rm B})+2\tan^{-1}(c_{\rm B})}>R_{\rm
B},$ (I.10)

where $R_{\rm B}$ is the scaling radius of the Burkert halo, and $c_{\rm B}$
is the Burkert concentration.

Another interesting property of the gravitational scaling function, $w(z)$ is
related to its behavior at the virial radius. Virial radius is defined by
$w(z_{\rm vir})=0$. Taking Eq. (I.7) at $z=z_{\rm vir}$ and using Gauss’s
theorem we obtain

$V_{g}(z_{\rm vir})=-\frac{GM\,m_{a}}{r_{\rm vir}}=-\frac{GM\,m_{a}}{R\,z_{\rm
vir}}=-\frac{GM\,m_{a}}{R}w(z_{\rm vir}),$ (I.11)

with the implication $w(z_{\rm vir})=1/z_{\rm vir}$.

We also define rescaled dimensionless wave functions, with unit normalization,
as

$\phi_{nl}(z)=R^{-3/2}N_{nl}^{1/2}\psi_{nl}(r).$ (I.12)

Writing Eq. (I.1) in terms of rescaled wave functions we obtain

$\epsilon_{nl}\phi_{\epsilon
l}=\frac{1}{S}\left(-\phi_{nl}^{\prime\prime}-\frac{2}{z}\phi_{nl}+\frac{l(l+1)}{z^{2}}\phi_{nl}\right)-w(z)\phi_{\epsilon
l},$ (I.13)

where

$\epsilon_{nl}=\frac{1}{S}2m_{a}R^{2}E_{nl}$ (I.14)

and where the dimensionless halo size parameter is

$S=2MGR\,m_{a}^{2}.$ (I.15)

As the expectation value of the kinetic term is positive, Eq. (I.13) implies
that the range of the scaled energy parameter is $-1\leq\epsilon\leq 0$.

Rough estimates show that $S\gg 1$ even for moderate size galaxies, even if
the ULA mass is as small as $m_{a}=10^{-22}$ eV. For the Milky Way
$S_{MW}\gtrsim 10^{5}$. As a contrast to the galactic halo, $S=O(1)$ for a
soliton, or axion star.

Using Eq. (I.13) we estimate the number of bound states. The depth of the
rescaled potential well is $w(0)=1$ and as we see later the spectrum of
$\epsilon_{nl}$ fills most of the interval $-1\leq\epsilon\leq 0$. As we will
see in the next section, WKB quantization implies that $\epsilon$ is quantized
as $\epsilon\sim-n/\sqrt{S}$, where $n$ is the principal quantum number. Then
the principal quantum number takes up to $O(\sqrt{S})$ different values.
Including states with nonzero angular momentum, we estimate that the number of
energy levels in the potential well is $O(\sqrt{S}$, a very large number. We
find no reason why most of those states would not be occupied by the
astronomical number of ULAs. In a previous study Bernal:2009zy ;
PhysRevD.101.081302 , only a small number of excited states were considered.

##  II WKB approximation

Eq. (I.13) lends itself to a perturbative WKB expansion in $\delta=1/\sqrt{S}$
book  which gives a solution to the differential equation,

$\displaystyle\phi_{nl}(z)={\text{Exp}}\left[\frac{1}{\delta}\sum_{k=0}^{\infty}\delta^{k}P_{k}\right],$
(II.1)

where for the two independent solutions in the allowed, oscillating region are

$\displaystyle P_{0}$ $\displaystyle=\pm\,i\,\int_{z_{\rm
min}}^{z}dz^{\prime}\sqrt{F_{nl}(z^{\prime})},$ (II.2) $\displaystyle P_{1}$
$\displaystyle=\log\left[\frac{{\cal{N}}_{nl}}{zF_{nl}(z)^{1/4}}\right]$
(II.3)

while terms $P_{2},...$ are of $O(S^{-1/2})$, and negligible. Eq.(I.13) is
linear, so any linear combination of solutions is admissible.The unique
solution satisfying boundary conditions, finiteness at the turning points,
which are zeros of $F_{nl}[z]$, is

$\displaystyle\phi_{nl}=\frac{{\cal{N}}_{nl}}{zF_{nl}(z)^{1/4}}\sin\left(\sqrt{S}\int_{z_{\rm
min}}^{z}dz^{\prime}\sqrt{F_{nl}(z^{\prime})}\right),$ (II.4)

where $F_{nl}(z)$ is

$F_{nl}(z)=w(z)+\epsilon_{nl}-\frac{l(l+1)}{z^{2}\,S}.$ (II.5)

Note that wave functions $\phi_{nl}(z)$ are normalized to 1. Multiplier ${\cal
N}_{nl}$ is introduced to ensure the correct normalization of wave functions,

$\displaystyle{\cal N}_{nl}^{-2}$ $\displaystyle=2\pi\int_{z_{\rm
min}}^{z_{\rm max}}\frac{dz}{F_{nl}(z)^{1/2}}.$ (II.6)

In the classically forbidden region the wave function decreases exponentially
as $\phi\propto\exp(-\sqrt{S}c)$, where $c$ is finite as $S\to\infty$.. In the
limit $S\to\infty$ the approximate solution Eq. (II.4) becomes exact. As
usual, $z_{\rm min}$ and $z_{\rm max}$ are the turning points where $F_{nl}=0$
(II.4) is also known as the Wentzel ansätz to the WKB solution of Eq. (I.13).

Noting that in Eq. (II.4), factors other then the exponential function vary
slowly as function of $n$ and $l$, the quantization condition for energy
eigenvalues $\epsilon_{nl}$ can be read off from Eq. (II.4) PhysRev.51.669 ;
PhysRevA.53.3798 :

$\int_{z_{\rm min}}^{z_{\rm max}}\sqrt{F_{nl}}=\frac{n}{\sqrt{S}}\pi.$ (II.7)

The principal quantum number, $n$, and the orbital quantum, $l$, only appear
in combinations $\nu=n/\sqrt{S}$ and $\lambda=(l+1/2)/\sqrt{S}$.111We use
Langer’s method to define $\lambda$ PhysRev.51.669 PhysRevA.53.3798 . As the
separation of subsequent values of quantum numbers $\nu$ and $\lambda$
vanishes as $S$ increases, they will be replaced by continuous variables. That
replacement will facilitate the construction of the density as a function of
$w(z)$ in the next section. Using the continuous quantum numbers, we can write
the wave function as

$\displaystyle\phi_{\epsilon\lambda}=\frac{{\cal{N}}_{\epsilon\lambda}}{zF_{\epsilon\lambda}(z)^{1/4}}\sin\left(\sqrt{S}\int_{z_{\rm
min}}^{z}dz^{\prime}\sqrt{F_{\epsilon\lambda}(z^{\prime})}\right),$ (II.8)

In Eq. (II.8) we use an unique connection between quantum numbers $\nu$ and
$\epsilon$ through quantization condition (II.7), as it will be explained int
the next section.

##  III Construction of the halo

The rescaled wave function of the halo is obtained from Eq. (I.2) as

$\Phi(z,t)=\sum_{nlm}\phi_{nl}(z)Y_{lm}(\theta,\phi)e^{i(\alpha_{nlm}+\epsilon_{nl}t)}.$
(III.1)

Then the total density, averaged over rotations and random phases, and
normalized to 1, is

$\rho(z)=|\Phi(z,t)|^{2}\simeq\sum_{nl}(2l+1)\frac{M_{nl}}{M}|\phi_{nl}|^{2}.$
(III.2)

Consider now that wave functions depend on quantum numbers $n$ and $l$ only
through combinations $\nu=n/\sqrt{S}$, in Eq. (II.7) and
$\lambda=\sqrt{l(l+1)}/\sqrt{S}$, in Eq. (II.5). The rescaled quantum numbers,
$\nu$ and $\lambda$, have finite ranges and become dense in those ranges as
$S\to\infty.$ Therefore, the error of turning summations over $n$ and $l$ into
integrations over $\nu$ and $\lambda$ is only of $O(1/\sqrt{S})$ and
negligible for almost all galaxies. Then we obtain

$\displaystyle\rho(z)=S^{3/2}\int d\nu\int
d\lambda^{2}\frac{{\cal{N}}_{\epsilon\lambda}^{2}}{z^{2}F_{\epsilon\lambda}^{1/2}}b(\epsilon,\lambda),$
(III.3)

where dimensionless, continuous distribution function $b(\epsilon,\lambda)$
interpolates distribution function $M_{nl}/M$. It is normalized as

$\displaystyle S^{3/2}\int d\nu\int d\lambda^{2}b(\epsilon,\lambda)=1.$
(III.4)

The total variation of $\epsilon$ is 1, or $1+\epsilon_{c}$, in case there is
a gap of size $|\epsilon_{c}|$. There are $O(\sqrt{S})$ energy levels,
implying that the discrete values of $\epsilon$ are dense on interval
$(-1,0)$. Wave functions depend explicitly on $\epsilon$ only, therefore we
change integration variable $\nu$ to $\epsilon$. We use the quantization
condition (Eq. (II.7)), written in terms of continuous quantum numbers,

$\nu\,\pi=\int_{z_{\rm min}}^{z_{\rm
max}}dz\sqrt{w[z]+\epsilon-\frac{\lambda^{2}}{z^{2}}},$ (III.5)

to find the appropriate Jacobian. Taking the derivative of Eq. (III.5) with
respect to $\epsilon$ at constant $\lambda$, we note that terms coming from
differentiating the integral with respect to boundary values $z_{\rm min}$ and
$z_{\rm max}$ vanish. Then using Eq. (II.6) we obtain

$\frac{d\nu}{d\epsilon}\pi=\frac{1}{2}\int_{z_{\rm min}}^{z_{\rm
max}}\frac{dz}{\sqrt{w(z)+\epsilon-\frac{\lambda^{2}}{z^{2}}}}=\frac{1}{4\pi}{\cal
N}_{nl}^{-2},$ (III.6)

Substituting Eq. (III.6) into Eq. (III.3) we arrive at our final result for
the density:

$\displaystyle\rho(z)=\frac{S^{3/2}}{4\pi^{2}}\int_{-w(z)}^{\epsilon_{\rm
max}}d\epsilon\int_{0}^{z^{2}(w(z)+\epsilon)}\frac{d\lambda^{2}b(\epsilon,\lambda)}{z^{2}\sqrt{w[z]+\epsilon-\frac{\lambda^{2}}{z^{2}}}},$
(III.7)

where $\epsilon_{\rm max}\leq 0$ is an admissible integration constant
allowing for the existence of an energy gap. If the distribution function
$b(\epsilon,\lambda)$ is independent of the angular momentum quantum number,
$\lambda$, integration over $\lambda$ yields

$\displaystyle\rho(z)$
$\displaystyle=\frac{S^{3/2}}{2\pi^{2}}\int_{-w(z)}^{\epsilon_{\rm
max}}d\epsilon\sqrt{w[z]+\epsilon}\,b(\epsilon)=\rho(0)\frac{\int_{-w(z)}^{\epsilon_{\rm
max}}d\epsilon\sqrt{w(z)+\epsilon}\,b(\epsilon)}{\int_{-1}^{\epsilon_{\rm
max}}d\epsilon\sqrt{1+\epsilon}\,b(\epsilon)}.$ (III.8)

Even if $b(\epsilon,\lambda)$ depends on $\lambda$ it is expected that that an
expansion with respect to $\lambda^{2}$ converges rapidly. Then, using the
expansion

$\displaystyle b(\epsilon,\lambda)=\sum_{k=0}b(\epsilon)^{(k)}\lambda^{2k}$
(III.9)

we can integrate the series with respect to $\lambda^{2}$, term by term to
obtain

$\displaystyle\rho(z)=\frac{\rho(0)}{2}\frac{\sum_{k}z^{2k}\frac{\sqrt{\pi}\Gamma(1+k)}{\Gamma(3/2+k)}\int_{-w(z)}^{\epsilon_{\rm
max}}d\epsilon\,b^{(k)}(\epsilon)(w(z)+\epsilon)^{k+1/2}}{\int_{-1}^{\epsilon_{\rm
max}}d\epsilon\sqrt{1+\epsilon}\,b(\epsilon)}$ (III.10)

Using Eqs. (III.8) or (III.10), the Poisson equation,

$w^{\prime\prime}(z)+\frac{2}{z}w^{\prime}(z)=-4\pi\rho(z)$ (III.11)

becomes a differential equation for $w(z)$. Since $w(0)=1$ and
$w^{\prime}(0)=0$, to insure that the density is regular at the center (the
exception is the case when density is replaced by the NFW profile), after
providing a distribution function there are no undetermined integration
constants, other than $\epsilon_{\rm max}$ and $\rho(0)$. As we will see in
Sec. V, for some distribution functions $b(\epsilon,\lambda)$ the density (Eq.
(III.7)) or (III.8), as appropriate, can be analytically integrated.

The strategy for solving Eq. (III.11) is simpler for Eq. (III.8) then for Eq.
(III.10), because after rescaling $z$ as $x=z\,\sqrt{\rho(0)}$, eliminating
the central density from the equation, we find a unique solution at fixed
$\epsilon_{\rm max}$. Among others, we find $\langle 1/x\rangle$. Then we can
restore the original scaling variable $z=x\langle 1/x\rangle$, which is true
because $\langle 1/z\rangle=1$. We also find $\rho(0)=\langle 1/x\rangle^{2}.$

It is more complicated to find the solution of Eq. (III.11) in the case when
the energy spectrum is not degenerate. Then, due to the explicit dependence of
$\rho(z)$ on $z$, coordinate $z$ cannot be eliminated from Eq. (III.10) by a
simple rescaling. In that case, rather than rescaling $z$ we need to search
for an appropriate value of $\rho(0)$ at fixed $\epsilon_{\rm max}$, such that
$\rho(z)$ is normalized to 1.

Dynamical simulations generate numerical distribution functions, as well
Lin:2018whl . A time dependent version of this work, to be published, can
potentially do that as well. One advantage of our analytic approximation
method is scalability: dynamical simulations are limited by computational
power to relatively small galactic halos.

Finally, we note that rotation curves have simple scale invariant
representations in terms of our scaling variable and scaling function. Up to
an overall factor of dimension of velocity, $v_{0}$,

$v(r)=v_{0}\sqrt{-z\,w^{\prime}(z)}.$ (III.12)

##  IV Virial theorem in WKB

Simulations show that, in agreement with expectations, the halo is virialized
after collapse Schive:2014hza Schwabe:2016rze . Using stationary WKB wave
functions the halo does not automatically satisfy the virial theorem. However,
a time dependent generalization of the WKB approximation is expected to to
converge to a virialized state. Such a generalization is deferred to a future
publication.

The virial theorem for a self-gravitating system, without contact
interactions, is $2K+E_{g}=0$, where $K$ is the total kinetic energy and
$V_{g}$ is the total gravitational energy. Using the definitions of the
previous sections, they are

$\displaystyle K$ $\displaystyle=\frac{1}{2m_{a}^{2}}\sum_{nl}(2l+1)\int
d^{3}r\left[|\psi_{nl}^{\prime}(r)|^{2}+\frac{l(l+1)}{r^{2}}|\psi_{nl}(r)|^{2}\right],$
(IV.1) $\displaystyle E_{g}$
$\displaystyle=\frac{1}{2m_{a}}\sum_{nl}(2l+1)\int
d^{3}r|\psi_{nl}(r)|^{2}V_{g}(r),$ (IV.2)

where $V_{g}(r)$ has been defined in Eq. (I.3). Using the dimensionless wave
functions and rescaled radial parameter, $z=r/R$, we can rewrite $K$ and $V$
as

$\displaystyle K$
$\displaystyle=\frac{1}{2m_{a}^{2}R^{2}}\sum_{nl}(2l+1)M_{nl}\times\int
d^{3}z\left[\phi_{nl}^{\prime}[z]^{2}+\frac{l(l+1)}{z^{2}}|\phi_{nl}(z)|^{2}\right],$
(IV.3) $\displaystyle E_{g}$
$\displaystyle=-\frac{MG}{2R}\sum_{nl}(2l+1)M_{nl}\int
d^{3}z|\phi_{nl}(z)|^{2}w(z)$ (IV.4)

Finally, using Eq. (II.8), introducing continuous energy and angular momentum
variables as in Eq. (III.3), and expressing $G$ by size parameter $S$, we
obtain

$\displaystyle K$ $\displaystyle=2C\int
dz\int_{-w(z)}^{\epsilon_{c}}d\epsilon\int_{0}^{z^{2}[w(z)+\epsilon]}d\lambda^{2}b(\epsilon,\lambda)\left[\sqrt{w(z)+\epsilon-\frac{\lambda^{2}}{z^{2}}}+\frac{\lambda^{2}}{z^{2}\sqrt{w(z)+\epsilon-\frac{\lambda^{2}}{z^{2}}}}\right],$
(IV.5) $\displaystyle E_{g}$ $\displaystyle=-C\int
dz\int_{-w(z)}^{\epsilon_{c}}d\epsilon\int_{0}^{z^{2}[w(z)+\epsilon]}d\lambda^{2}b(\epsilon,\lambda)\frac{w(z)}{\sqrt{w(z)+\epsilon-\frac{\lambda^{2}}{z^{2}}}},$
(IV.6)

where $C=SM/(4m_{a}^{2}R^{2})$, is a constant and $\epsilon_{c}\leq 0$ is an
admissible energy cutoff parameter. Then the virial theorem applied to our
system becomes independent of all dimensional parameters and of dimensionless
size parameter, $S$. It reads as

$\displaystyle\int_{9}^{z_{\rm
vir}}dz\,z^{2}\int_{-w(z)}^{\epsilon_{c}}d\epsilon\int_{0}^{z^{2}[w(z)+\epsilon]}d\lambda^{2}\frac{b(\epsilon,\lambda)(3w(z)+4\epsilon)}{\sqrt{w(z)+\epsilon-\frac{\lambda^{2}}{z^{2}}}}=0.$
(IV.7)

In the particular case when distribution function, $b(\epsilon\lambda)\to
d_{\epsilon}$ is independent of $\lambda$, we can integrate over $\lambda$ and
the virial theorem becomes

$\displaystyle\int_{9}^{z_{\rm
vir}}dz\,z^{2}\int_{-w(z)}^{\epsilon_{c}}d\epsilon\,b(\epsilon)\sqrt{w[z]+\epsilon}(3w(z)+4\epsilon)=0.$
(IV.8)

Suppose now that an ansatz for the distribution function $b(\epsilon,\lambda)$
depends on a free parameter. A simple example is, $\epsilon_{\rm max}$,
defining a gap in the energy spectrum, $b(\epsilon,\lambda)=0$ for
$\epsilon_{\rm max}<\epsilon<0$. Another possibility is that a background
density, $\rho_{0}$, of undetermined size is subtracted from the density.
Then, using the procedure described in Sec. III, we find the numerical
gravitational scaling function, $w(z)$, as a function of $\epsilon_{\rm max}$
or $\rho_{0}$. Substituting into Eq. (IV.7) allows us to fix $\epsilon_{\rm
max}$ or, alternatively, $\rho_{0}$, which we will demonstrate in the next
section. If there is no choice for the free parameter to satisfy the condition
given by Eq. IV.7, then the distribution function $b(\epsilon,\lambda)$ does
not allow for the system to reach dynamical equilibrium and is not physically
acceptable.

##  V Examples for halos

To complete the calculation of the density (Eq. (III.7)) we have to know the
energy distribution function $b(\epsilon,\lambda)$. In this section we will
explore a variety of physical choices for $b(\epsilon,\lambda)$ to check
whether they provide acceptable density distributions. In particular, we will
pay particular attention to whether the obtained density distribution are in
agreement with the following general features of observational data and
simulations.

  * •

It is generally accepted Schive:2014dra Hui:2016ltb Schwabe:2016rze  that
after its collapse, at least asymptotically in time, the halo is virialized,
i.e. satisfies virial condition, $2K+E_{g}=0$. In simulations the system is
driven towards dynamical equilibrium. As we pointed out earlier, unlike in
simulations, $b(\epsilon,\lambda)$ must have a parameter, which can be chosen
such that the condition is satisfied.

  * •

Another important feature is that galaxies and galaxy clusters are coming in a
very wide range of sizes. As condensation $c$ controls the size of the halo,
general solutions must have solutions for a large range of $c$. However, not
all density distributions satisfy this criterion, as an example, a model which
we examined in a previous work Street:2021qxl .

  * •

Popular models for galactic halos predict that the rotation curves are
independent of the concentration. Though all observed rotation curves are
close to being flat, significant variations have been observed Zintner:2022zds
. This is especially important if one tries to describe a wide variety halos,
including galaxy clusters.

###  V.1 Thermal equilibrium

Recently, we applied the WKB expansion method to investigate a system in
thermal equilibrium Street:2021qxl . We give a simplified account of those
investigations, below. Neglecting gravitational interactions, the system is
treated as an ideal Bose-Einstein gas. The number of particles on energy
levels $E_{n}$ is $N_{n}$. Then the partition function for the for the system
is

$e^{-\beta\,g}=\prod_{n}\frac{1}{1-e^{-\beta(E_{n}-\mu)}}.$ (V.1)

The average value of $N_{n}$ is

$\langle N_{n}\rangle=\frac{1}{e^{\beta(E_{n}-\mu)}-1}$ (V.2)

The occupation number, $\langle N_{n}\rangle$, is enormous in every state,
$n$. As a result, every exponent must satisfy the inequality
$1\gg\beta\left(E_{n}-\mu\right)>0$. Consequently, we can expand them, keeping
terms to linear order, to arrive at the Rayleigh-Jeans limit of Bose-Einstein
statistics.

$\langle N_{n}\rangle=\frac{1}{\beta(E_{n}-\mu)}$ (V.3)

Using Eq. (I.13) and rescaling $\beta$ and $\mu$ as

$\displaystyle\beta$ $\displaystyle=\frac{2m_{a}R^{2}}{S}\tilde{\beta},$ (V.4)
$\displaystyle\mu=\frac{S}{2m_{a}R^{2}}\tilde{\mu}$

we substitute into Eq. (III.8), to get

$\displaystyle\rho(z)=\rho_{0}\left(\int_{-w(z)}^{0}\frac{\sqrt{w(z)+\epsilon}}{\epsilon-\tilde{\mu}}d\epsilon-\rho_{B}\right),$
(V.5)

where

$\rho_{0}=\frac{S^{3/2}}{4\pi^{2}\tilde{\beta}}$

and $\rho_{0}\rho_{B}$ is a background density term. $\rho_{B}$ is considered
to be an adjustable parameter.

Now remember that the ground state corresponds to $\epsilon=-1$, Then the
physical range of the chemical potential is $\tilde{\mu}\leq-1$. During the
collapse of the halo the system cools and the chemical potential increase
towards its critical value, $\tilde{\mu}=-1$, where Bose-Einstein condensation
starts. Integrating density $\rho(z)$ and setting $\tilde{\mu}=-1$ we obtain
Poisson equation

$\displaystyle w^{\prime\prime}(z)+\frac{2}{z}w^{\prime}(z)=-8\rho_{0}\pi$
(V.6)
$\displaystyle\left[\tan^{-1}\left(\sqrt{\frac{w(z))}{1-w(z)}}\right)\sqrt{1-w(z)}-\sqrt{w(z)}-\frac{\rho_{B}}{2}\right].$

To facilitate the numerical solution of Eq. (V.6) we rescale coordinate $z\to
x/\sqrt{8\rho_{0}\pi}$. Then integrating the Poisson equation at a series of
choices for background density parameter $\rho_{B}$, we find that at
$\rho_{B}=0.017$, the virial condition (Eq. (IV.8)) is satisfied. As at that
choice of $\rho_{B}$ $w(x_{\rm vir})=0$, where $x_{\rm vir}\simeq 4.74$. We
also calculate the concentration,

$c=\left\langle\frac{1}{x}\right\rangle x_{\rm vir}\simeq 2.78$

A similar value is obtained for $c$ if we modify the distribution function to
take into account gravitational interactions. That can be done by replacing
energy level $E_{n}$ by $H_{n}$, the contribution of that energy level to the
conserved Hamiltonian. Since observations imply $c\gtrsim 10$, the final state
of the collapse cannot be in dynamical equilibrium and thermal equilibrium at
the same time.

###  V.2 The King model

The King model 1966AJ…..71…64K  is based on the classical kinetic theory of
collision-free self-gravitating particles, assuming Maxwell velocity
distribution cut off at an escape velocity. Then the equilibrium energy
distribution takes the form

$f_{\rm King}=\begin{cases}A\left(e^{-\beta(E-E_{c})}-1\right)&\text{if }E\leq
E_{c}\\\ 0&\text{otherwise.}\end{cases}$ (V.7)

In a recent paper, Lin:2018whl , galactic halos were constructed, using the
King distribution 1966AJ…..71…64K  (and other classical distribution
functions, like the fermionic King model and the Ossipov-Merritt model
1979PAZh….5…77O ; 1985AJ…..90.1027M ), combined by wave functions of excited
states of Eq. I.1 obtained by simulations. The halos were compared to profiles
obtained by simulations of collapsing systems of ultralight bosons, with
excellent agreement. Unlike in Lin:2018whl , we do not add a contribution for
the condensate at the core of the density distribution, because our aim is to
compare the densities predicted by the King model to phenomenological models,
rather than the results of simulations. However, our result are in excellent
agreement with simulations outside of the $z\lesssim 0.1$ region where the
soliton core dominates the density distribution.

We will use our analytic WKB wave functions (Eq. (II.8)) combined with the
King distribution function to construct halos, which reduces the number of
parameters used in the construction. The fermionic King model, also used in
Lin:2018whl  give slightly better fits, as it has an extra adjustable
parameter. As we will see below, using WKB wave functions, rather than
simulations, has the advantage of scalability, which allows us to find new
features when investigating larger systems.

Adapting Eq. (V.7) to WKB wave functions we write Eq. (III.8) as

$\displaystyle\rho(z)=\frac{A\,S^{3/2}}{4\pi^{2}}\int_{-w(z)}^{\epsilon_{c}}d\epsilon\left(e^{-\beta(\epsilon-\epsilon_{c})}-1\right)\sqrt{w[z]+\epsilon},$
(V.8)

Scaling out coefficient $A\,(S/\beta)^{3/2}/4\pi^{2}$ from the density and
rescaling $z\to x$, as explained in Sec. III, we integrate over $\epsilon$ to
get the Poisson equation

$\displaystyle q^{\prime\prime}(x)$ $\displaystyle+\frac{2}{x}q^{\prime}(x)$
(V.9)
$\displaystyle=-3\left[e^{q(x)}\text{erf}\left(\sqrt{q(x)}\right)-6\sqrt{q(x)}-4q(x)^{3/2}\right],$

where

$q(x)=\beta(w(z)+\epsilon_{c}),$ (V.10)

and where

$x=z\frac{A^{1/2}S^{3/2}}{2\sqrt{6}\pi\beta^{1/4}}.$ (V.11)

The only adjustable parameter in the Poisson equation (Eq. (V.9)) is
$q(0)=\beta(1+\epsilon_{c})$, which we took be to a number of values. We found
that the concentration, $c=z_{\text{vir}}=x_{\text{vir}}\langle 1/x\rangle$,
plotted in Fig. 1, is a rapidly rising function of $q(0)$. The implication is
that the King model can describe halos of markedly different sizes.

Figure 1:  Virial radius as a function of $q(0)$

We found that in the range $7\lesssim q(0)\lesssim 10$, the King model density
distributions are quite close to NFW or Burkert distributions well in the
region where the bulk of the density is distributed, $0.1\lesssim z\lesssim
10$. As an example, we plot the WKB density distribution, $\rho(z)$,
corresponding to initial condition $q(0)=8$ in Fig. 2. We switch to the
physical, rescaled radial coordinate, $z$ and choose central density
$\rho(0)=1$ in the plot. NFW and Burkert profiles, the core density and core
radius parameters of which were fitted to the numerically calculated halo, are
also plotted. We found fitted parameters $R_{\text{NFW}}=0.437R$, and
$R_{\text{Burkert}}=0.269R$, where $R$ is the harmonic radius of the WKB
solution. Then, we were able to calculate the virial radii and concentrations
of those profiles using Eqs. (I.9) and (I.10). We found $c_{\text{NFW}}=22.3$
and $c_{\text{Burkert}}=35.7$ and physical virial radii,

$\displaystyle z_{\text{NFW}}$
$\displaystyle=c_{\text{NFW}}R_{\text{NFW}}/R=9.78,$ $\displaystyle
z_{\text{Burkert}}$
$\displaystyle=c_{\text{Burkert}}R_{\text{Burkert}}/R=9.74,$ (V.12)

which are quite close to each other. The virial radius of the WKB halo is
$z_{\text{vir}}=20.1$, as it does not have a sharp cutoff. However, in the
region $10<z<20$, $\rho(z)=\mathcal{O}\left(10^{-5}\right)$.

Figure 2: King model density distribution compared to NFW and Burkert profiles
at $q(0)=8.$

Next, we examine the predictions of the King model at a larger virial radius,
in the range where $c$ is a linearly rising function of $q(0)$, as it can be
seen in Fig. 1. Very likely that range corresponds to galaxies larger than the
Milky Way, or galaxy clusters, which are currently not accessible for
simulations. At $q(0)\gtrsim 12$ density profiles are very different from
those at $q(0)=8$, plotted in Fig. 3. In fact, they are in excellent agreement
with a slightly cored isothermal profile, with behavior $\rho(z)\propto
z^{-2}.$ We plot $\rho(z)$ at $q(0)=25$ as a function of $z$ in Fig. 3, along
with a weakly cored pseudo-isothermal profile, an NFW profile, and a Burkert
profile. All profiles are scaled to match at $z=1$, which corresponds to the
harmonic average radius of the system. Note that the bulk of contributions to
the total mass come from the range $0.01\lesssim z\lesssim 10$.

Figure 3: King model density distribution compared to isothermal, NFW and
Burkert profiles at $q(0)=25.$

Clearly, the King model profile is in excellent agreement with the isothermal
profile May:2022may Jim:2003jvt Shu:2008cis Bro:2009.bs , while the NFW and
Burkert profiles are not.

Lin et. al. Lin:2018whl  generated two halos labeled A and B, for which they
fit energy cutoff parameter, $E_{c}$ and scale parameter, $\beta_{\rm sim}$,
defining the King model. We wish to compare those with parameters we use to
construct WKB halos. However, the corresponding values of $\beta_{\rm WKB}$
and $\epsilon_{c}$ cannot be directly compared with the simulations of
Lin:2018whl , because we rescaled energy $E\to\epsilon$ and, consequently, our
$\beta$ value is also rescaled, though rescaling leaves the dimensionless
quantity $E_{c}\beta_{\rm sim}=\epsilon_{c}\beta_{\rm WKB}$, unchanged. To
check that equality, we need to use the relationship (Eq. (V.10)), which
implies of $q(0)=\beta_{\rm WKB}(1+\epsilon_{c})$. Note now, that cutoff
parameter $\epsilon_{c}$ has not been fixed in previous calculations. In other
words, $\epsilon_{c}$ is a free parameter.

Consider now that the halos generated in Lin:2018whl  are virialized, while
there is no reason why our halos would in general be. In fact, calculating
(IV.8) as

$\displaystyle 2\,K+E$ $\displaystyle\propto\int_{0}^{z_{\rm
vir}}dz\,z^{2}\int_{-w(z)}^{\epsilon_{c}}d\epsilon\sqrt{w(z)+\epsilon}$ (V.13)
$\displaystyle\times\left(e^{-\beta(\epsilon-\epsilon_{c})}-1\right)\left(\epsilon+\frac{3}{4}w(z)\right)$

we find that at every $q(0)$ we have $\langle
2\,K+E\rangle=a+b\,\beta\,\epsilon_{c}$ with varying fixed values for $a$ and
$b$. Then the vanishing of $\langle 2\,K+E\rangle$ fixes the combination
$\beta\epsilon_{c}$ at every choice of $q(0)$.

We plot the combination $\epsilon_{c}\beta$ as function of $q(0)$ along with
the values of $E_{c}\beta_{\rm sim}$ of halos $A$ and $B$, $-0.539$ and
$-0.51$, of Lin:2018whl , in Fig. 4. Those combinations are equal after
rescaling $E_{c}$ and $\beta$. The deviation of values of
$\beta\,\epsilon_{c}$ between those found in our WKB calculation and those in
simulations, are only a few percent at the relevant values of $q(0)$.

Figure 4:

Plot of King model values of $\beta\,\epsilon_{c}$ along with those of halos A
and B obtained in simulations of Lin:2018whl .

##  VI Summary

We have investigated the possibility that galactic dark matter haloes could be
made up of ultra-light bosons with wave functions that satisfy nonlinear
Schrodinger-Poisson equations. We have found eigenstates and eigenvalues of
the Schrodinger equation using the WKB approximation. The approximation method
becomes more accurate as the galaxy’s mass rises, however it has been
demonstrated that even for smaller galaxies, the leading WKB approximation may
yield accurate answers. The WKB approximation expresses the wave functions in
terms of the gravitational potential. The energy levels were determined by the
Bohr–Sommerfeld quantization condition, implied by the WKB method. To
determine how we may better explain the data we have compared two anätze for
the level occupation number distributions, the Bose-Einstein distribution and
the King model, at appropriate values of temperature and chemical potential.
The modulus square of each wave function multiplied by the occupation number,
summed over all the bound states, yielded the total density of particles. The
mass density is therefore a function of the gravitational potential $V_{g}$
since the Poisson equation connects $V_{g}$ to the particle density. This
procedure allowed us to obtain a differential equation for the gravitational
potential. Solving this equation enabled us to obtain the mass distribution of
the ULDM model of a galaxy. This technique can easily be scaled up to model
the DM halos of galaxies or galaxy clusters, well beyond the mass range
covered by simulations.

When we used the Bose Einstein density distribution in our computation of the
density profile we observed that the concentration parameter, $c$, of the
resulting galaxy was smaller in magnitude than what is observed for most
galaxies. When the King model particle distribution was utilized the
concentration parameters, on the other hand, were within the permissible range
for Milky way-like galaxies.

In future work we plan to study the dynamical collapse of boson clouds using
variational technique and study the times scales for collapse and decay in
more detail. We will also investigate how the inclusion of self-interactions
affects the the dynamics of halo.

###### Acknowledgements.

The authors are indebted to Joshua Eby for fruitful discussions. L.S. also
thanks the Department of Physics at the University of Cincinnati for financial
support in the form of the Violet M. Diller Fellowship. During part of this
research L.S. was supported by the U.S. Department of Energy (DOE), Office of
Science, Office of Workforce Development for Teachers and Scientists, Office
of Science Graduate Student Research (SCGSR) program. The SCGSR program is
administered by the Oak Ridge Institute for Science and Education (ORISE) for
the DOE. ORISE is managed by ORAU under contract number DE-SC0014664. All
opinions expressed in this paper are the authors? and do not necessarily
reflect the policies and views of DOE, ORAU, or ORISE. Research of L.C.R.W. is
partially supported by the US. Department of Energy grant DE-SC1019775.

## References

  * (1) P. J. E. Peebles, _Large-scale background temperature and mass fluctuations due to scale-invariant primeval perturbations_ , _”AtroPhys. J.”_ 263 (1982) L1.
  * (2) J. R. Bond, A. S. Szalay and M. S. Turner, _Formation of Galaxies in a Gravitino Dominated Universe_ , _Phys. Rev. Lett._ 48 (1982) 1636.
  * (3) G. R. Blumenthal, H. Pagels and J. R. Primack, _Galaxy Formation by Dissipationless Particles Heavier than Neutrinos_ , _Nature_ 299 (1982) 37.
  * (4) M. R. Baldeschi, R. Ruffini and G. B. Gelmini, _On Massive Fermions and Bosons in Galactic Halos_ , _Phys. Lett. B_ 122 (1983) 221.
  * (5) M. S. Turner, _Coherent Scalar Field Oscillations in an Expanding Universe_ , _Phys. Rev. D_ 28 (1983) 1243.
  * (6) W. H. Press, B. S. Ryden and D. N. Spergel, _Single Mechanism for Generating Large Scale Structure and Providing Dark Missing Matter_ , _Phys. Rev. Lett._ 64 (1990) 1084.
  * (7) S.-J. Sin, _Late time cosmological phase transition and galactic halo as Bose liquid_ , _Phys. Rev. D_ 50 (1994) 3650 [hep-ph/9205208].
  * (8) W. Hu, R. Barkana and A. Gruzinov, _Cold and fuzzy dark matter_ , _Phys. Rev. Lett._ 85 (2000) 1158 [astro-ph/0003365].
  * (9) F. S. Guzman and L. A. Urena-Lopez, _Evolution of the Schrodinger-Newton system for a selfgravitating scalar field_ , _Phys. Rev. D_ 69 (2004) 124033 [gr-qc/0404014].
  * (10) C. G. Boehmer and T. Harko, _Dark energy as a massive vector field_ , _Eur. Phys. J. C_ 50 (2007) 423 [gr-qc/0701029].
  * (11) D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray and A. H. G. Peter, _Cold dark matter: controversies on small scales_ , _Proc. Nat. Acad. Sci._ 112 (2015) 12249 [1306.0913].
  * (12) D. J. E. Marsh and A.-R. Pop, _Axion dark matter, solitons and the cusp–core problem_ , _Mon. Not. Roy. Astron. Soc._ 451 (2015) 2479 [1502.03456].
  * (13) L. Hui, J. P. Ostriker, S. Tremaine and E. Witten, _Ultralight scalars as cosmological dark matter_ , _Phys. Rev. D_ 95 (2017) 043541 [1610.08297].
  * (14) J.-W. Lee, _Brief History of Ultra-light Scalar Dark Matter Models_ , _EPJ Web Conf._ 168 (2018) 06005 [1704.05057].
  * (15) E. G. M. Ferreira, _Ultra-light dark matter_ , _Astron. Astrophys. Rev._ 29 (2021) 7 [2005.03254].
  * (16) H.-Y. Schive, T. Chiueh and T. Broadhurst, _Cosmic Structure as the Quantum Interference of a Coherent Dark Wave_ , _Nature Phys._ 10 (2014) 496 [1406.6586].
  * (17) H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst et al., _Understanding the Core-Halo Relation of Quantum Wave Dark Matter from 3D Simulations_ , _Phys. Rev. Lett._ 113 (2014) 261302 [1407.7762].
  * (18) B. Schwabe, J. C. Niemeyer and J. F. Engels, _Simulations of solitonic core mergers in ultralight axion dark matter cosmologies_ , _Phys. Rev. D_ 94 (2016) 043513 [1606.05151].
  * (19) J. Veltmaat and J. C. Niemeyer, _Cosmological particle-in-cell simulations with ultralight axion dark matter_ , _Phys. Rev. D_ 94 (2016) 123523 [1608.00802].
  * (20) X. Du, C. Behrens, J. C. Niemeyer and B. Schwabe, _Core-halo mass relation of ultralight axion dark matter from merger history_ , _Phys. Rev. D_ 95 (2017) 043519 [1609.09414].
  * (21) P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov et al., _Galaxy formation with BECDM – I. Turbulence and relaxation of idealized haloes_ , _Mon. Not. Roy. Astron. Soc._ 471 (2017) 4559 [1705.05845].
  * (22) D. G. Levkov, A. G. Panin and I. I. Tkachev, _Gravitational Bose-Einstein condensation in the kinetic regime_ , _Phys. Rev. Lett._ 121 (2018) 151301 [1804.05857].
  * (23) P. Mocz et al., _First star-forming structures in fuzzy cosmic filaments_ , _Phys. Rev. Lett._ 123 (2019) 141301 [1910.01653].
  * (24) M. A. Amin and P. Mocz, _Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe_ , _Phys. Rev. D_ [1902.07261].
  * (25) N. Glennon and C. Prescod-Weinstein, _Modifying PyUltraLight to model scalar dark matter with self-interactions_ , _Phys. Rev. D_ 104 (2021) 083532 [2011.09510].
  * (26) P. Mocz, A. Fialkov, M. Vogelsberger, M. Boylan-Kolchin, P. H. Chavanis, M. A. Amin, S. Bose, T. Dome, L. Hernquist and L. Lancaster, et al. Mon. Not. Roy. Astron. Soc. 521, no.2, 2608-2615 (2023) doi:10.1093/mnras/stad694 [arXiv:2301.10266 [astro-ph.CO]].
  * (27) H. Huang, H.-Y. Schive and T. Chiueh, _Cosmological Simulations of Two-Component Wave Dark Matter_ , 2212.14288.
  * (28) N. Glennon, N. Musoke and C. Prescod-Weinstein, _Simulations of multifield ultralight axionlike dark matter_ , _Phys. Rev. D_ 107 (2023) 063520.
  * (29) J. Luna Zagorac, I. Sands, N. Padmanabhan, ans R. Easther, _Schrödinger-Poisson solitons: Perturbation theory_ , _Phys. Rev. D_ 105 (2022) 103506. [1103.2054].
  * (30) S.-C. Lin, H.-Y. Schive, S.-K. Wong and T. Chiueh, _Self-consistent construction of virialized wave dark matter halos_ , _Phys. Rev. D_ 97 (2018) 103523 [1801.02320]
  * (31) F. S. Guzmán and L. A. Ureña López, _Gravitational atoms: General framework for the construction of multistate axially symmetric solutions of the Schrödinger-Poisson system_ , _Phys. Rev. D_ 101 (2020) 081302.
  * (32) P.-H. Chavanis, _Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. i. analytical results_ , _Phys. Rev.D_ [2109.01928]
  * (33) P. H. Chavanis and L. Delfini, _Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: II. Numerical results_ , _Phys. Rev. D_ 84 (2011) 043532 [1103.2054].
  * (34) J. Eby, P. Suranyi and L. C. R. Wijewardhana, _The Lifetime of Axion Stars_ , _Mod. Phys. Lett. A_ 31 (2016) 1650090 [1512.01709].
  * (35) J. F. Navarro, C. S. Frenk and S. D. M. White, _The Structure of cold dark matter halos_ , _Astrophys. J._ 462 (1996) 563 [astro-ph/9508025].
  * (36) A. Burkert, _The Structure of dark matter halos in dwarf galaxies_ , _Astrophys. J. Lett._ 447 (1995) L25 [astro-ph/9504041].
  * (37) A. Bernal, J. Barranco, D. Alic and C. Palenzuela, _Multi-state Boson Stars_ , _Phys. Rev. D_ 81 (2010) 044031 [0908.2435].
  * (38) C. Bender and S. Orszag, _Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory_ , vol. 1. 01, 1999, 10.1007/978-1-4757-3069-2.
  * (39) R. E. Langer, _On the connection formulas and the solutions of the wave equation_ , _Phys. Rev._ 51 (1937) 669.
  * (40) M. N. Sergeenko, _Semiclassical wave equation and exactness of the WKB method_ , _Phys. Rev. A_ 53 (1996) 3798.
  * (41) A. Zintner, S. Dandavate, O, Slone and M. Lisanti, _A crirical Assessment of Solutions to the Galaxy Diversity Problem_ , _JCAP_ 121 (2018) 151301 [2202.00012v2]
  * (42) L. Street, P. Suranyi and L. C. R. Wijewardhana, _Density profile of multi-state fuzzy dark matter_ , 2101.00349.
  * (43) I. R. King, _The structure of star clusters. III. Some simple dynamical models_ , _”Astron. J.”_ 71 (1966) 64.
  * (44) L. P. Osipkov, _Spherical systems of gravitating bodies with ellipsoidal velocity distribution._ , _Pisma v Astronomicheskii Zhurnal_ 5 (1979) 77.
  * (45) D. Merritt, _Spherical stellar systems with spheroidal velocity distributions_ , _”Astron. J.”_ 90 (1985) 1027.
  * (46) S. May, V. Springel, _The halo mass function and filaments in full cosmological simulations with fuzzy dark matter_ _Monthly Notices of the Royal Astronomical Society/stad2031_ [ 0302.560]
  * (47) R. Jimenez, L. Verde, T. Treu and D. Stern, _Constraints on the equation of state of dark energy and the Hubble constant from stellar ages and the CMB_ _Astrophys. J._ 593 (2003) 6220 [ 0302.560]
  * (48) C. Shu, B. Zhou, M. Bartemann, J. Comerfold, J. Huang and Y. Mellier, _Comparisons between Isothermal and NFW Mass Profiles for Strong-Lensing Galaxy Clusters_ , _Astrophys. J._ 685 (2008) 70 [ 0805.1148]
  * (49) J. Brownstein, _Modified Gravity and the Phantom of Dark Matter_ [ 0908.0040v1]


# Tuning exciton and biexciton transition energies and fine structure
splitting through hydrostatic pressure in single InGaAs quantum dots

Xuefei Wu  State Key Laboratory for Superlattices and Microstructures,
Institute of Semiconductors, Chinese Academy of Sciences, Beijing, 100083,
People’s Republic of China  Hai Wei  Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei, 230026, People’s
Republic of China  Xiuming Dou  State Key Laboratory for Superlattices and
Microstructures, Institute of Semiconductors, Chinese Academy of Sciences,
Beijing, 100083, People’s Republic of China  Kun Ding  State Key Laboratory
for Superlattices and Microstructures, Institute of Semiconductors, Chinese
Academy of Sciences, Beijing, 100083, People’s Republic of China  Ying Yu
State Key Laboratory for Superlattices and Microstructures, Institute of
Semiconductors, Chinese Academy of Sciences, Beijing, 100083, People’s
Republic of China  Haiqiao Ni  State Key Laboratory for Superlattices and
Microstructures, Institute of Semiconductors, Chinese Academy of Sciences,
Beijing, 100083, People’s Republic of China  Zhichuan Niu  State Key
Laboratory for Superlattices and Microstructures, Institute of Semiconductors,
Chinese Academy of Sciences, Beijing, 100083, People’s Republic of China  Yang
Ji  State Key Laboratory for Superlattices and Microstructures, Institute of
Semiconductors, Chinese Academy of Sciences, Beijing, 100083, People’s
Republic of China  Shushen Li  State Key Laboratory for Superlattices and
Microstructures, Institute of Semiconductors, Chinese Academy of Sciences,
Beijing, 100083, People’s Republic of China  Desheng Jiang  State Key
Laboratory for Superlattices and Microstructures, Institute of Semiconductors,
Chinese Academy of Sciences, Beijing, 100083, People’s Republic of China
Guang-can Guo  Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei, 230026, People’s Republic of China  Lixin He
<EMAIL_ADDRESS>Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei, 230026, People’s Republic of China  Baoquan
Sun<EMAIL_ADDRESS>State Key Laboratory for Superlattices and
Microstructures, Institute of Semiconductors, Chinese Academy of Sciences,
Beijing, 100083, People’s Republic of China

###### Abstract

We demonstrate that the exciton and biexciton emission energies as well as
exciton fine structure splitting (FSS) in single (In,Ga)As/GaAs quantum dots
(QDs) can be efficiently tuned using hydrostatic pressure in situ in an
optical cryostat at up to 4.4 GPa. The maximum exciton emission energy shift
was up to 380 meV, and the FSS was up to 180 $\mu$eV. We successfully produced
a biexciton antibinding-binding transition in QDs, which is the key
experimental condition that generates color- and polarization-
indistinguishable photon pairs from the cascade of biexciton emissions and
that generates entangled photons via a time-reordering scheme. We perform
atomistic pseudopotential calculations on realistic (In,Ga)As/GaAs QDs to
understand the physical mechanism underlying the hydrostatic pressure-induced
effects.

###### pacs:

78.67.Hc, 07.35.+k, 78.55.Cr, 42.50.-p

Self-assembled semiconductor quantum dots (QDs) have considerable potential
for use as fundamental building blocks in future quantum information
applications. However, so far it is impossible to use QD growth techniques for
precisely controlling QD properties, which is essential for such applications.
Therefore, externally tuning the QD properties post-growth is extremely
important. One the most prominent examples is polarization-entangled photon
pair emission through a biexciton (XX) cascade process in QDs, which requires
that the different polarized photons are energetically indistinguishable.
However, the underlying asymmetry for self-assembled (In,Ga)As/GaAs QDs leads
to splitting in degenerate bright exciton (X) states (fine structure
splitting, FSS), which is typically tens of $\mu$eV Gammon et al. (1996a);
Bayer et al. (2002); Gammon et al. (1996b); Bester et al. (2003), and much
larger than the radiative linewidth ($\sim$ 1.0 $\mu$eV); therefore, photon
entanglement is destroyed Stevenson et al. (2006); Hafenbrak et al. (2007).
Tuning techniques, such as electric Bennett et al. (2010); Gerardot et al.
(2007); Vogel et al. (2007), magnetic Hudson et al. (2007), and strain fields
Trotta et al. (2012); Ding et al. (2010); Jöns et al. (2011); Seidl et al.
(2006); Dou et al. (2008); Gong et al. (2011); Wang et al. (2012) used to
erase the FSS have been explored. Uniaxial and biaxial stresses have been used
to tune the QD structural symmetry, exciton binding energies and FSS. However,
the strain that can be generated using such techniques is limited to
approximately tens of MPa, which corresponds to a spectral shift by only
several meVs for QD peak emissions Ding et al. (2010); Jöns et al. (2011).
Herein, we report high-pressure research (up to tens of GPa) for individual
QDs using the diamond anvil cell (DAC), which has been widely used to study
metal-semiconductor transitions, electronic structures, and optical
transitions in bulk crystals and microstructures Jayaraman (1983); Ma et al.
(2004); Itskevich et al. (1998); Li et al. (1994). An exciton emission line
shift in ensemble InAs/GaAs QDs is approximately 500 meV at 8 GPa Ma et al.
(2004), which is much larger than the QD peak shift induced by a piezoelectric
actuator PMN-PT Ding et al. (2010). However, tuning the QD structural
symmetry, exciton transition states and FSS for an individual QD using the DAC
has not been reported.

In this work, we demonstrate that the exciton (X) emission energy, FSS and
biexciton (XX) binding energy can be successively tuned for extremely large
ranges using hydrostatic pressure at up to 4.4 GPa. The emission energy, FSS
and XX binding energy almost increase linearly with increased pressure. The
maximum exciton emission energy shift and FSS change can extend to 380 meV and
180 $\mu$eV, respectively, which is considerably greater than through other
techniques. By tuning the applied pressure, color-indistinguishable photons
from the biexciton and exciton emission decay through a cascade, and across
generation color coincidence for biexciton and exciton transitions are
generated. Therefore, entangled photon pairs are generated via the proposed
“time reordering” scheme Avron et al. (2008). We also perform atomistic
pseudopotential calculations on realistic (In,Ga)As/GaAs QDs under hydrostatic
stress to discern the physical mechanisms underlying the effects induced by
the hydrostatic pressure.

The investigated (In,Ga)As/GaAs QD samples with low QD density were grown
using molecular beam epitaxy (MBE) on a semi-insulating GaAs substrate with
excitonic emission energies at 1.35-1.43 eV. Figure 1(a) shows the DAC
pressure device used for tuning QD photoluminescence (PL) in situ using the
optical cryostat. To fit the QD samples into the DAC chamber [indicated in
Fig. 1(b)], the samples were mechanically thinned to a total thickness of
approximately 20 $\mu$m and then cut into pieces approximately 100$\times$100
$\mu$m2. Condensed argon was used as the pressure-transmitting medium in the
DAC, which can be used to apply hydrostatic pressure up to 9 GPa Jayaraman
(1983); Shimizu et al. (2001). The initial pressure can be adjusted at room
temperature by driving screws and can be determined in situ using the ruby R1
fluorescence line shift. To successively tune the X and XX transition energies
and fine structure splitting (FSS) by pressure at a low temperature, a novel
and easily controlled version of the DAC shown in Fig. 1(a) was developed by
combining the well-known DAC with a piezo actuator. This device can
successively generate pressures up to several GPa for the QD samples studied
at a low temperature using PL measurements, and the maximum applied pressure
depends on the actuator stroke length. The QD sample in the DAC is cooled to
20 K through a continuous-flow liquid helium cryostat and excited by a He-Ne
laser at the wavelength 632.8 nm. The excitation laser was focused to a $\sim$
2 $\mu$m spot on the sample using a microscope objective (NA: 0.35). The PL
was collected using the same objective, spectrally filtered through a 0.5 m
monochromator, and detected using a silicon-charge coupled device (CCD). A
$\lambda$/2 wave plate and linear polarizer were used to distinguish
horizontal (H) and vertical (V) linear polarization for PL components. By
carefully following the changing X and XX PL energies using the polarization
angle, we measure FSS with an $\sim$ 10 $\mu$eV accuracy by fitting the
experimental data to a sinusoidal function Ghali et al. (2012). Figure 1(c)
displays the measured pressure values and excitonic emission energies at 20 K
as a function of actuator voltage (Piezo-ceramics: PSt 150/10$\times$10/40).
We clearly show that pressure can be successively tuned in situ using an
optical cryostat from 0.5 to 4.4 GPa through a piezo actuator, and the
corresponding blue shift for the excitonic PL peak energy is $\sim$ 310 meV.

Figure 2(a) depicts the exciton emission energies as a function of the
hydrostatic pressure from 0 to 4.4 GPa for QD1-QD5. At 0 GPa, the QD exciton
emission energies are 1.401, 1.349, 1.406, 1.432 and 1.394 eV, respectively.
The exciton emission energies for the five QDs studied herein increased
linearly with the applied pressure. The blue shift for the QD1 peak energy at
4.22 GPa is approximately 330 meV, which is much larger than previously
reported shifts ( $\sim$ 10 meV) from uniaxial or biaxial stresses using
conventional methods Ding et al. (2010); Dou et al. (2008); Jöns et al.
(2011); Seidl et al. (2006); Trotta et al. (2012). The experimental data were
linearly fit, which generated pressure coefficients for QD1-QD5 of 82, 87, 93,
81 and 85 meV/GPa, respectively; such values are consistent with the reported
pressure coefficients for ensemble quantum dot Ma et al. (2004).

Figure 2(b) shows the biexciton binding energies for QD1-QD5 as a function of
hydrostatic pressure at up to 4.4 GPa. The biexciton binding energy is defined
as $E_{B}(XX)$=$E_{X}$-$E_{XX}$, where $E_{X}$ and $E_{XX}$ are the X and XX
emission energies, respectively. When $E_{B}(XX)>0$, the biexciton is in the
“binding” state, wherein the two excitons are attracted. When $E_{B}(XX)<0$,
the biexciton is in the “antibinding” state, wherein the two excitons are
repulsive. For the QDs studied herein, $E_{B}(XX)$ increases as a function of
hydrostatic pressure up to 4.4 GPa. For QD1 and QD3, the biexcitons are in an
antibinding state at zero pressure and gradually progress to the binding state
at approximately 1 and 2 GPa, respectively (i.e., $E_{B}(XX)$=0), where the
exciton and biexciton are “color-indistinguishable”.

To demonstrate the biexciton antibinding-binding transitions under pressure in
greater detail, we plotted the polarization-resolved PL spectra for the QD1 X
and XX emission lines under different hydrostatic pressures, as shown in Fig.
3(a)-(e), wherein the red and black lines correspond to the horizontal (H) and
vertical (V) polarized photons, respectively. At zero pressure, both XX
emission energies, $E(H2)$ and $E(V2)$, are higher than the X emission energy,
$E(H1)$ and $E(V1)$ [see also the scheme in Fig. 3(f)]. In addition, $E(H2)$
is slightly larger than $E(V2)$ at FSS $\sim$ 50 $\mu$eV. Under pressure, the
blue shift for the X emission energy (82 meV/GPa) is more rapid than for XX
(81 meV/GPa). Therefore, with increasing pressure, the V-polarized XX and X
emission lines first degenerate at 1.62 GPa, as shown in Fig. 3(b), and then
the H-polarized emission lines degenerate at 2.07 GPa, as shown in Fig. 3(d).
In such instances, color-indistinguishable photon pairs are generated by an
XX-X cascade emission at 1.62 GPa for V-polarized photons or at 2.07 GPa for
H-polarized photons. Therefore, it is expected that the indistinguishable two-
photon streams will be produced by adjusting a time delay between the XX and X
emissions, wherein the time delay is approximately 0.4 ns Chang et al. (2009).
Remarkably, at 1.97 GPa, across generation color coincidence for XX and X
transition energies was generated (i.e., $E(H1)$=$E(V2)$ and $E(V1)$=$E(H2)$).
This is a key condition for entangled photon generation via the proposed time
reordering scheme Avron et al. (2008). When pressure was further increased,
the separation between the XX and X emission lines again increased, as shown
in Fig. 3(e) at 3.66 GPa. Ding and coworkers demonstrated that biaxial strain
can also tune the biexciton binding energies Ding et al. (2010). However,
because their experiment generated a relatively small strain, biexciton
antibinding-binding progression was not observed.

FSS tuning by uniaxial strain has been studied experimentally Seidl et al.
(2006); Trotta et al. (2012); Kuklewicz et al. (2012) and theoreticallySingh
and Bester (2010); Gong et al. (2011); Wang et al. (2012), which has shown
that the maximum tuned FSS value is approximately 20 $\mu$eV. It is
interesting to measure the FSS change under hydrostatic pressure. Figure 2(c)
depicts the FSS for QD4 and QD5 as a function of pressure at 20 K in the range
0.5 to 4.4 GPa. The figure clearly demonstrates that increasing pressure
produces an approximately linear increase in FSS with the slope 44 and 28
$\mu$eV/GPa for QD4 and QD5, respectively, which generates a total FSS shift
as large as $\sim$ 180 and 100 $\mu$eV for QD4 and QD5, respectively. Similar
results were observed from other investigated (In,Ga)As/GaAs QDs, which
indicates that such a large shift is typical for FSS under hydrostatic
pressure.

To understand the experimental results, we calculated the electronic and
optical properties for the In1-xGaxAs/GaAs QDs under hydrostatic pressure
using an atomistic empirical pseudopotential method (EPM) Williamson et al.
(2000). The optimized QD structures are obtained by the valence force filed
method Keating (1966). We then calculate the electron/hole single-particle
energies and wave functions using the linear combination of bulk bands (LCBB)
method Wang and Zunger (1999). The exciton and biexciton energies are
calculated via the configuration interaction (CI) method Franceschetti et al.
(1999). Herein, we present results for three QDs: (i) a lens-shaped InAs/GaAs
QDs with the height $h$=1.5 nm and base diameter $b$=12 nm; (ii) a lens-shaped
In0.8Ga0.2As/GaAs QDs with $h$=1.5 nm and $b$=12, 15 nm; and (iii)
In0.8Ga0.2As/GaAs QDs with $h$=2.5 and the elliptical major (minor) axis
$a$=10 nm ($b$=7.5 nm) along the [1$\mathrm{bar}{1}$0] ([110]) crystal
direction.

The calculated exciton emission energies under pressure are shown in Fig. 4(a)
and produce blue energy shifts at approximately 76 meV/GPa, which is
consistent with the experimental values. To understand the emission energy
blue shift, we analyzed the band offsets and confinement potentials for the
QDs under pressure, which strongly depend on the strain distribution in the
dots and matrix. When hydrostatic pressure is applied, the lattice constant
for the matrix material GaAs decreases, which effectively increases the
lattice mismatch between the dot material InAs and GaAs matrix. As a result,
both the (absolute values of) isotropic and biaxial strain inside the dots
increase. The averaged isotropic strain $I$=-0.072-0.011 $P$ and biaxial
strain $\epsilon_{zz}-\epsilon_{xx}$= 0.12+0.0014 $P$, where $P$ is the
applied hydrostatic pressure in GPa. Figure 4(b) depicts the strain-modified
band offsets for the conduction band (e), heavy hole (HH), light hole (LH) and
spin-orbit (SO) bands through the dot center under P =0, 2 and 4 GPa. Whereas
the band offset change is small for holes, the band offset changes
dramatically for the conduction band. Under pressure, the electron bands move
significantly toward the higher energy. The confinement potential also
increased dramatically with increasing pressure, which is the major reason for
the observed experimental results.

Because the electron-hole Coulomb energy change is relatively small [see Fig.
4(d)], the change in exciton emission energy can be estimated using the
electron-hole single-particle gap $E_{g}$, which can be written as follows:

$E_{g}(\tensor{\epsilon})=E_{g}(0)+a_{g}I+b_{v}(\epsilon_{zz}-\epsilon_{xx})\,,$
(1)

where $\tensor{\epsilon}$ is the strain tensor inside the InAs dots,
$a_{g}$=$-$6.08 eV is the hydrostatic deformation potential for the band gap,
and $b_{v}$=$-$1.8 eV is the biaxial deformation potential for the valence
band maximum. Therefore, we estimated that the exciton PL blue shift under
hydrostatic pressure is 82 meV/GPa for pure InAs/GaAs QDs, which is consistent
with the experimental values and EPM calculations. We note that the
hydrostatic pressure is much more efficient at tuning the exciton emission
energy than uniaxial stress ($\sim$10 $\mu$eV/MPa) Jöns et al. (2011);
Kuklewicz et al. (2012); Wang et al. (unpublished).

The calculated XX binding energies $E_{B}(XX)$ are presented in Fig. 4(c). We
found that the biexciton tends toward antibinding in small QDs under zero
pressure. When the pressure increases, the $E_{B}(XX)$ for the dots calculated
increased. The binding energy tends to be saturated at very high pressure. The
XX binding energy for the In0.8Ga0.2As/GaAs QDs with $b$=12 nm and $h$=1.5 nm
is consistent with the experimental QD1.

In the calculation, we found that it is important to include many
electron/hole energy levels for the correct XX binding energies using the CI
calculations, and the XX binding energy change as a function of pressure can
be observed only using the lowest energy conduction and valence bands (i.e.,
Hartree-Fock approximation). Such observations suggest that the XX binding
energies did not change due to the correlated energies; primarily, such
changes are due to changes in the direct Coulomb integrals between the lowest
electron and hole states, as follows:

$\Delta E_{B}(XX)\approx 2\Delta J_{eh}-\Delta J_{ee}-\Delta J_{hh}\,,$ (2)

where $J_{ee}$, $J_{hh}$ and $J_{eh}$ are the direct electron-electron, hole-
hole and electron-hole Coulomb integrals, respectively. As shown in Fig. 4(d),
whereas $J_{ee}$, and $J_{eh}$ increase rapidly with pressure, $J_{hh}$ is
approximately flat. The solid purple line describes the changing exciton
binding energy calculated using Eq. (2), which is consistent with the dashed
purple line from the EPM calculations. To understand how the Coulomb integrals
change under pressure, we compared the lowest electron and hole wave functions
in Fig. 4(f) at 0.0 GPa and 4.0 GPa for the 12$\times$1.5 nm In0.8Ga0.2As/GaAs
QDs. We found that, whereas the hole wave function shape primarily does not
change, the electron becomes much more localized due to the band offset
changes shown in Fig. 4(b), which explains the Coulomb integral changes under
pressure.

Finally, we examined FSS under hydrostatic pressure. The FSS calculated as a
function of pressure is shown in Fig. 4(e), which increases dramatically with
applied pressure and is consistent with the experimental data in Fig. 2(c). It
is surprising that FSS changes under hydrostatic pressure, which does not
change the QD symmetry. However, because the electron wave functions are more
localized under pressure, an electron-hole would have a larger effective
overlap under pressure, which increases the exchange energies (e.g., the dark-
bright splitting $\Delta_{bd}$, which is also shown in Fig. 4(e)). It has been
shown that FSS can be roughly estimated as $\sim 2\eta\Delta_{bd}$ Wang et al.
(unpublished), where $\eta$ is the HH-LH mixing parameter; therefore, as
$\Delta_{bd}$ increases, FSS increases, as clearly demonstrated in Fig. 4(e).

To summarize, we experimentally and theoretically investigated the effects of
hydrostatic pressure on the exciton and biexciton transition energies as well
as FSS in single InGaAs QDs. The excitonic emission energies and FSS can be
tuned in situ by applying hydrostatic pressure in an optical cryostat for
changes over a wide energy range. The observed exciton emission energy blue
shift and FSS change were as large as $\sim$ 380 meV and $\sim$ 180 $\mu$eV,
respectively, which is greater than the values from other strain-adjusting
techniques. Tuning the QD optical properties over such a larger spectral range
yields great advantages for future QD applications, such as for generating
color-indistinguishable photon pairs from the biexciton and exciton emission
decay cascades or generating entangled photon pairs via a time-reordering
scheme. Furthermore, it expected that photon antibunching for optical
communication band QD emission can be measured using the pressure-induced blue
shift into the spectral range detected by sufficiently powerful silicon
avalanche photodiodes.

BS and LH acknowledge support from the Chinese National Fundamental Research
Program (Grant Nos. 2013CB922304, 2013CB933304, 2011CB921200, 2009CB929301),
Chinese National Natural Science Funds (Grant No. 90921015), and National
Natural Science Funds for Distinguished Young Scholars.

## References

  * Gammon et al. (1996a) D. Gammon, E. S. Snow, B. V. Shanabrook, D. S. Katzer, and D. Park, Phys. Rev. Lett. 76, 3005 (1996a).
  * Bayer et al. (2002) M. Bayer, G. Ortner, O. Stern, A. Kuther, A. A. Gorbunov, A. Forchel, P. Hawrylak, S. Fafard, K. Hinzer, T. L. Reinecke, et al., Phys. Rev. B 65, 195315 (2002).
  * Gammon et al. (1996b) D. Gammon, E. S. Snow, B. V. Shanabrook, D. S. Katzer, and D. Park, Science 5, 87 (1996b).
  * Bester et al. (2003) G. Bester, S. Nair, and A. Zunger, Phys. Rev. B 67, 161306 (2003).
  * Stevenson et al. (2006) R. M. Stevenson, R. J. Young, P. Atkinson, K. Cooper, D. A. Ritchie, and A. J. Shields, Nature 439, 179 (2006).
  * Hafenbrak et al. (2007) R. Hafenbrak, S. M. Ulrich, P. Michler, L. Wang, A. Rastelli, and O. G. Schmidt, New Journal of Physics 9, 315 (2007).
  * Bennett et al. (2010) A. J. Bennett, M. A. Pooley, R. M. Stevenson, M. B. Ward, R. B. Patel, A. Boyer de la Giroday, N. Sköld, I. Farrer, C. A. Nicoll, D. A. Ritchie, et al., Nature Phys. 6, 947 (2010).
  * Gerardot et al. (2007) B. D. Gerardot, S. Seidl, P. A. Dalgarno, R. J. Warburton, D. Granados, J. M. Garcia, K. Kowalik, O. Krebs, K. Karrai, A. Badolato, et al., Appl. Phys. Lett. 90, 041101 (2007).
  * Vogel et al. (2007) M. M. Vogel, S. M. Ulrich, R. Hafenbrak, P. Michler, L. Wang, A. Rastelli, and O. G. Schmidt, Appl. Phys. Lett. 91, 051904 (2007).
  * Hudson et al. (2007) A. J. Hudson, R. M. Stevenson, A. J. Bennett, R. J. Young, C. A. Nicoll, P. Atkinson, K. Cooper, D. A. Ritchie, and A. J. Shields, Phys. Rev. Lett. 99, 266802 (2007).
  * Trotta et al. (2012) R. Trotta, E. Zallo, C. Ortix, P. Atkinson, J. D. Plumhof, J. van den Brink, A. Rastelli, and O. G. Schmidt, Phys. Rev. Lett. 109, 147401 (2012).
  * Ding et al. (2010) F. Ding, R. Singh, J. D. Plumhof, T. Zander, V. Křápek, Y. H. Chen, M. Benyoucef, V. Zwiller, K. Dörr, G. Bester, et al., Phys. Rev. Lett. 104, 067405 (2010).
  * Jöns et al. (2011) K. D. Jöns, R. Hafenbrak, R. Singh, F. Ding, J. D. Plumhof, A. Rastelli, O. G. Schmidt, G. Bester, and P. Michler, Phys. Rev. Lett. 107, 217402 (2011).
  * Seidl et al. (2006) S. Seidl, M. Kroner, A. Högele, K. Karrai, R. J. Warburton, A. Badolato, and P. M. Petroff, Appl. Phys. Lett. 88, 203113 (2006).
  * Dou et al. (2008) X. M. Dou, B. Q. Sun, B. R. Wang, S. S. Ma, R. Zhou, S. S. Huang, H. Q. Ni, and Z. C. Niu, Chin. Phys. Lett. 25, 1120 (2008).
  * Gong et al. (2011) M. Gong, W. Zhang, G.-C. Guo, and L. He, Physical Review Letters 106, 227401 (2011).
  * Wang et al. (2012) J. Wang, M. Gong, G.-C. Guo, and L. He, Applied Physics Letters 101, 063114 (2012).
  * Jayaraman (1983) A. Jayaraman, Rev. Mod. Phys. 55, 65 (1983).
  * Ma et al. (2004) B. S. Ma, X. D. Wang, F. H. Su, Z. L. Fang, K. Ding, Z. C. Niu, and G. H. Li, J. Appl. Phys. 95, 933 (2004).
  * Itskevich et al. (1998) I. E. Itskevich, S. G. Lyapin, I. A. Troyan, P. C. Klipstein, L. Eaves, P. C. Main, and M. Henini, Phys. Rev. B 58, R4250 (1998).
  * Li et al. (1994) G. H. Li, A. R. Goñi, K. Syassen, O. Brandt, and K. Ploog, Phys. Rev. B 50, 18420 (1994).
  * Avron et al. (2008) J. E. Avron, G. Bisker, D. Gershoni, N. H. Lindner, E. A. Meirom, and R. J. Warburton, Phys. Rev. Lett. 100, 120501 (2008).
  * Shimizu et al. (2001) H. Shimizu, H. Tashiro, T. Kume, and S. Sasaki, Phys. Rev. Lett. 86, 4568 (2001).
  * Ghali et al. (2012) M. Ghali, K. Ohtani, Y. Ohno, and H. Ohno, Nat. Commun. 3, 661 (2012).
  * Chang et al. (2009) X. Y. Chang, X. M. Dou, B. Q. Sun, Y. H. Xiong, Z. C. Niu, H. Q. Ni, and D. S. Jiang, J. Appl. Phys. 106, 103716 (2009).
  * Kuklewicz et al. (2012) C. E. Kuklewicz, R. N. E. Malein, P. M. Petroff, and B. D. Gerardot, Nano Lett. 12, 3761 (2012).
  * Singh and Bester (2010) R. Singh and G. Bester, Phys. Rev. Lett. 104, 196803 (2010).
  * Williamson et al. (2000) A. J. Williamson, L.-W. Wang, and A. Zunger, Phys. Rev. B 62, 12963 (2000).
  * Keating (1966) P. N. Keating, Phys. Rev 145, 637 (1966).
  * Wang and Zunger (1999) L.-W. Wang and A. Zunger, Phys. Rev. B 59, 15806 (1999).
  * Franceschetti et al. (1999) A. Franceschetti, H. Fu, L.-W. Wang, and A. Zunger, Phys. Rev. B 60, 1819 (1999).
  * Wang et al. (unpublished) J. Wang, G.-C. Guo, and L. He (unpublished).

Figure 1: (Color online) (a) Schematic drawing of the successively applied
pressure device of the diamond anvil cell (DAC), where the DAC and piezo
actuator are assembled together via home-made copper cylinder. (b) The DAC
chamber, showing the positions of the QD sample and ruby in the DAC. (c)
Measured pressure in the DAC chamber (solid black circles) and the
corresponding QD excitonic PL peak energy (solid green circles) as a function
of actuator voltage. At zero voltage, an initial pressure of 0.5 GPa was
generated by four driven screws.  Figure 2: (Color online) (a) Exciton
emission energies as a function of pressure for QD1-QD5. (b) Biexciton binding
energies as a function of pressure for QD1-QD5, showing biexciton antibinding-
binding transitions under pressure for QD1 and QD3. (c) FSS of QD4 and QD5 as
a function of pressure.

Figure 3: (Color online) (a)-(e) Polarization-resolved PL spectra for the QD1
X and XX emission lines under different pressures, where red lines (black
lines) correspond to the horizontal (H) and vertical (V) polarized photons,
respectively. At 1.62 (b) and 2.07 (d) GPa, color-indistinguishable photon
pairs are generated by an XX-X cascade emissions for V- and H- polarized
photons, respectively. At 1.97 (c) GPa, across generation color coincidence
for XX and X transition energies is achieved. (f) Level schemes showing the
XX-X cascade at 0 GPa.  Figure 4:  (Color online) The calculated results of:
(a) The exciton emission energies of (In,Ga)As/GaAs QDs as a function of
hydrostatic pressure. (b) The band offsets of the InAs/GaAs dot under the 0.0,
2.0 and 4.0 GPa hydrostatic pressure. The black, green, blue and red lines
indicate conduction (e), heavy-hole (HH), light-hole (LH) and spin-orbit (SO)
bands, respectively. (c) The biexciton binding energies as a function of
hydrostatic pressure in (In,Ga)As/GaAs QDs. (d) The changes of direct Coulomb
integrals of the lowest electron and hole states in the 12$\times$1.5 nm
In0.8Ga0.2As/GaAs QD as a function of the applied pressure. (e) The FSS in the
10$\times$7.5$\times$2.5 nm In0.8Ga0.2As/GaAs QD as a function of the applied
pressure. (f) The wave functions of the lowest electron (e0) and hole (h0)
states under 0.0 and 4.0 GPa hydrostatic pressure in the 12$\times$1.5 nm
In0.8Ga0.2As/GaAs QD.