# Environment-Assisted Quantum Transport and Mobility Edges

Donny Dwiputra<EMAIL_ADDRESS>Theoretical Physics Laboratory,
Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung, Jl.
Ganesha 10, Bandung 40132, Indonesia  Freddy P. Zen<EMAIL_ADDRESS>Theoretical Physics Laboratory, Faculty of Mathematics and Natural Sciences,
Institut Teknologi Bandung, Jl. Ganesha 10, Bandung 40132, Indonesia
Indonesian Center for Theoretical and Mathematical Physics (ICTMP), Bandung
40132, Indonesia

###### Abstract

Environment-assisted quantum transport (ENAQT) is a unique situation where
environmental noise can, counterintuitively, enhance the transport of an open
quantum system. In this Letter, we investigate how the presence of a one-
dimensional single-particle mobility edge can generate strong ENAQT. For this
purpose, we study the energy current of a generalized Aubry-Andre-Harper (AAH)
tight binding model coupled at its edges to spin baths of differing
temperature and dephasing noise along the system. We find that the ENAQT
increases by orders of magnitude and depends on the number of localized
eigenstates and disorder strength nonmonotonically. We show that this
enhancement is the result of the cooperative interplay between population
spread and localization.

††preprint: APS/123-QED

Introduction.—Understanding and controlling the out-of-equilibrium dynamics in
open quantum systems are at the forefront of research in the field of quantum
thermodynamics. One particular phenomenon which has attracted considerable
attention, both theoretically [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] and
experimentally [14, 15, 16, 17, 18], are the environment-assisted quantum
transport (ENAQT) where moderate dephasing noise can enhance the transport
efficiency. Fueled by the quest for sustainable energy, several studies has
exploited the potential benefit from the interplay between coherence and
dissipation in quantum devices such as to enhance the power of quantum heat
engines [19, 20, 21, 22]. Nonequilibrium disordered systems with localization
properties can act as heat engine as work can be extracted from the transition
between the athermal localized states to an equilibrium extended state [23].
Systems with mobility edge on the other hand may act as a current rectifier
[24, 25] or an energy filter to obtain a strong thermoelectric response [26].
However, the link between ENAQT and mobility edge is far from clear.

ENAQT may be explained by different mechanisms: initially its origin is
understood from the destruction of Anderson localization [27]—which hinders
the transport—by dephasing in a disordered system [1, 2, 3, 4, 2, 6] such as
the photosynthetic complex [28, 29, 30, 31]. However, ENAQT persists in
ordered systems, where localization does not play role, and in fact it is
impossible only for transport in an ordered chain [7]. Thus it relies on a
beneficial competition between coherent and incoherent dynamics which results
in higher population entering the sink. In some specific cases, momentum
rejuvenation [8], line-broadening [9], and superradiance [10] may additionally
explain the high efficiency. Moreover, exposing the system to a periodic
driving may increase the efficiency even further [12, 13]. It is only revealed
recently that the uniformity of steady state population plays more universal
role in ENAQT [11]. Consider a chain of qubit sites without static disorder:
when the system is coupled only to a source and a sink (no environmental
dephasing), be it a particle pulse or a current source, transport is dominated
by the unitary dynamics and thus the corresponding wave function spreads out
with quantum interference pattern. Adding dephasing noise leads to decoherence
in site basis, which in turn reduces the fluctuation and thereby uniformize
the population. This increases the sink particle density and thus amplifying
the current. In the limit of strong dephasing, transport is suppressed under
the quantum Zeno effect [32]. For end-to-end transport, the population is
already maximally spread; dephasing will only reduce the current unless
localized eigenstates exist.

The interesting consequence of Anderson localization is the quantum phase
transition between extended (metallic) and localized (insulating) states as a
function of disorder strength. In one and two dimensions, all eigenstates are
localized for a random (uncorrelated) disorder of any strength independent of
energy while in three dimensions the localization transition occur in a
critical disorder strength forming a sharp mobility edge, i.e., a threshold
energy value that separates localized and extended eigenstates. Yet if we
replace the random disorder with a quasiperiodic potential such as the
paradigmatic Aubry-André-Harper (AAH) model [33, 34], which has been realized
in recent ultracold atom experiments with interacting and noninteracting
particles [35, 36, 37], the localization transition occurs even in one
dimension. Moreover, by introducing nonlocal hopping terms or deforming the
on-site potential, the model can manifest an exact mobility edge [38, 39, 40,
41], defined by the self-duality point.

In this Letter, we show that a thermodynamically model with single particle
mobility edge (SPME) can manifest a strong ENAQT, measured by comparing the
heat current to the one without dephasing. The population uniformization
mechanism alone is not sufficient to explain the nonmonotonic behavior of
ENAQT with the fraction of localized eigenstates; the spread of population
does not always correlate with ENAQT. Instead, both population spread and
localization must coexist in order to obtain the maximum enhancement. We also
show that the maximum ENAQT varies nonmonotonically with disorder strength. To
isolate the ENAQT from the enhancement due to interference effects, we
consider an end-to-end transport.

Model.—We study a one-dimensional generalized AAH (GAAH) model [41] of $N$
bosonic modes with nearest-neighbor hopping described by the Hamiltonian

$H_{S}=-t\sum_{n=1}^{N-1}(a_{n}^{\dagger}a_{n+1}+\text{H.c.})+\sum_{n=1}^{N}V_{n}a_{n}^{\dagger}a_{n},$
(1)

where $a_{n}(a_{n}^{\dagger})$ is the usual bosonic annihilation (creation)
operator at site $n$, $t$ is the tunneling, which sets the energy unit. We
consider zero and single-excitation maniflods given by the states $|0\rangle$
and $|n\rangle=a_{n}^{\dagger}|0\rangle$, so that the local modes
$a_{n}=|0\rangle\langle n|$. The on-site potential is

$V_{n}=\Delta+2\lambda\frac{1-\cos(2\pi\beta n+\phi)}{1+\alpha\cos(2\pi\beta
n+\phi)}.$ (2)

Here $\lambda$ sets the average disorder strength and the offset is chosen to
be $\Delta=2$ to avoid negative eigenvalues. We consider the mobility edge
parameter $\alpha\in(-1,1)$, which can be tuned to represent several limiting
cases: $\alpha=-1$ corresponds to an ordered chain, $\alpha=0$ is the original
(scaled) AAH model, and $\alpha=1$ contains singular potential. We set
$\beta=(\sqrt{5}-1)/2$ for a quasiperiodic modulation. This model has a SPME
at $E_{\mathrm{ME}}=2\text{sgn}(\lambda)(|t|-|\lambda|)/\alpha+\Delta$. The
phase $\phi$ does not vary the energy spectrum; it only shifts the position of
the localized modes [41]. From hereafter, we set $\hbar=k_{\text{B}}=1$ and
$t=1$ as the unit of energy.

The chain is coupled to two different environment: pair of thermal baths of
noninteracting spins $\\{\sigma_{\mu,l}\\}$ at different temperatures to its
edges, $l=\\{1,N\\}$; characterized by $H_{SB_{1}}$, and dephasing noise
$H_{SB_{2}}$ due to spin bath coupled to all the sites. The former generates
energy (or heat) current and the latter creates fluctuation in site energies,
which in turn destroy the coherence of the system. The total Hamiltonian is
$H=H_{S}+\sum_{i=1}^{2}(H_{B_{i}}+H_{SB_{i}})$, where
$H_{B_{1}}=\sum_{l,\mu}\varepsilon_{\mu,l}\sigma^{z}_{\mu,l}$,
$H_{B_{2}}=\sum_{m,\mu}\varepsilon_{\mu}\sigma^{z}_{\mu,m}$ with
$n=\\{1,\dots,N\\}$, and the interactions

$\displaystyle H_{SB_{1}}$ $\displaystyle=$
$\displaystyle\sum_{l,\mu}g_{\mu}^{(1)}(a_{l}+a_{l}^{\dagger})(\sigma^{+}_{\mu,l}+\sigma^{-}_{\mu,l}),$
(3) $\displaystyle H_{SB_{2}}$ $\displaystyle=$
$\displaystyle\sum_{n,\mu}g_{\mu}^{(2)}a_{n}^{\dagger}a_{n}(\sigma^{+}_{\mu,n}+\sigma^{-}_{\mu,n}),$
(4)

where $\sigma^{+}_{\mu,n}(\sigma^{-}_{\mu,n}))$ are the raising (lowering)
operators for the $\mu$-th spin coupled to site $n$. This treatment, which
explicitly considers the dynamics of the current heat baths and the
interactions, $H_{B_{1}}+H_{SB_{1}}$, guarantees a thermodynamically
consistent model [42]. In contrast to the typical models where the particle is
initiated in one site and transferred to an irreversible sink, here the
particle enters and exits the heat baths at finite temperature.

To derive the Lindblad master equation [43, 44, 45,
[Forderivationofmasterequationinsimilarsystems,
seee.g.intheAppendixof]guimaraes2016nonequilibrium,
[SeeSupplementalMaterialat[linkbypublisher]fordetailedderivationofthemasterequationandmoredetailedanalysisontheroleofpopulationspreadingtotheENAQTmechanism]supplemental]
for the density matrix $\rho$, we work in the eigenbasis
$|\eta_{k}\rangle=\sum_{n}S_{nk}|l\rangle$. In this basis, the system
Hamiltonian becomes $H_{S}=\sum_{k}\epsilon_{k}\eta^{\dagger}_{k}\eta_{k}$ and
the local mode is $a_{n}=\sum_{k}S_{nk}\eta_{k}$, where $\epsilon_{k}$ is the
$k$-th energy (arranged in increasing order) and $\eta_{k}$ and
$\eta_{k}^{\dagger}$ being its eigenmode operators. The corresponding master
equation is

$\displaystyle\frac{d\rho}{dt}$ $\displaystyle=$ $\displaystyle-i[H,\rho]$ (5)
$\displaystyle+$
$\displaystyle\sum_{l,k}|S_{lk}|^{2}J(|\epsilon_{k}|)\big{[}(1-n_{l}(\epsilon_{k}))\mathcal{D}[\eta_{k}]\rho+n_{l}(\epsilon_{k})\mathcal{D}[\eta_{k}^{\dagger}]\rho\big{]}$
$\displaystyle+$
$\displaystyle\sum_{n,k<k^{\prime}}|S_{nk}|^{2}|S_{nk^{\prime}}|^{2}\gamma(|\epsilon_{k^{\prime}}-\epsilon_{k}|)\mathcal{D}[\eta_{k}^{\dagger}\eta_{k^{\prime}}]\rho,$

where the dissipator is $\mathcal{D}[A]\rho=A\rho
A^{\dagger}-1/2\\{A^{\dagger}A,\rho\\}$ and
$n_{l}(\epsilon_{k})=(e^{\epsilon_{k}/T_{l}}+1)^{-1}$ is the spin occupation
of the bath coupled to the edge site $l$. These superoperators describe a
nonequilibrium condition in which energy is constantly pumped and absorbed at
both ends, while dephasing occurs at every eigenmodes. $J(\epsilon)$ and
$\gamma(\epsilon)$ are the spectral density of system-bath interaction
$H_{SB_{1}}$ and $H_{SB_{2}}$, respectively. Here we have assumed that $B_{2}$
is at $T=0$ so that its occupation factor $n(\epsilon)\rightarrow 0$; this
results in condition $\epsilon_{k}<\epsilon_{k^{\prime}}$ and thus the sum is
for $k<k^{\prime}$. Since we have no information about the system-bath
coupling, we assume for simplicity that $J(\epsilon_{k})=1$ and
$\gamma(|\epsilon_{k^{\prime}}-\epsilon_{k}|)=\gamma$ for all $k$ and
$k^{\prime}$. Nevertheless, Ohmic spectral density
$J(\epsilon)\propto\epsilon$ will produce similar results.

The nonequilibrium steady-state (NESS) can be calculated exactly once the
Hamiltonian is diagonalized. The density matrix can be solved iteratively from
$k=N$:

$\langle\eta_{k}^{\dagger}\eta_{k}\rangle=A_{k}^{-1}\Big{(}B_{k}+\gamma\sum_{i>k}^{N}C_{ik}\langle\eta_{i}^{\dagger}\eta_{i}\rangle\Big{)},$
(6)

with the coefficients

$\displaystyle\centering A_{k}\@add@centering$ $\displaystyle=$
$\displaystyle\sum_{l}|S_{lk}|^{2}(1-2n_{l}(\epsilon_{k}))+\gamma\sum_{i<k}C_{ik},$
$\displaystyle B_{k}$ $\displaystyle=$
$\displaystyle\sum_{l}|S_{lk}|^{2}n_{l}(\epsilon_{k}),\;C_{ik}=\sum_{n}|S_{ni}|^{2}|S_{nk}|^{2},$
(7)

while the off-diagonal components are zero. Note that this state needs to be
normalized. The peculiar factor 2 in $A_{k}$ appears because the system and
the bath have different statistics. The dephasing part of Eq. (5) transfers
the population from the higher energy eigenmode $|\eta_{k^{\prime}}\rangle$ to
the lower energy $|\eta_{k}\rangle$ (since $k<k^{\prime}$). This can be
inferred from the last term in Eq. (6) in which $C_{ik}$ measures the coupling
between two eigenmodes in presence of dephasing. To determine ENAQT we
consider the NESS energy current $\mathcal{J}$ [48],

$\displaystyle\mathcal{J}$ $\displaystyle=$
$\displaystyle-\text{T}r[H\mathcal{D}_{N}\rho]$ (8) $\displaystyle=$
$\displaystyle\sum_{k}\epsilon_{k}|S_{Nk}|^{2}\left(\left(1-2n_{N}(\epsilon_{k})\right)\langle\eta_{k}^{\dagger}\eta_{k}\rangle-
n_{N}(\epsilon_{k})\right),\quad\,$

where $\mathcal{D}_{N}\cdot$ is the sum of $SB_{1}$ superoperators acting on
site $N$.

Figure 1: (a) Relative current $\mathcal{J}/\mathcal{J}_{0}$ versus the
dephasing rate $\gamma$ for different mobility edge parameters $\alpha$. Each
solid line corresponds to horizontal lines crossing the ($\alpha,\gamma$)
contour plot of unnormalized current $\mathcal{J}$ with the same colors and
symbols in (b). The parameters are $\lambda=0.4$, $\phi=\pi/3$, and the bath
temperatures $T_{h}=10^{3}$ (at $l=1$) and $T_{c}=10^{-1}$ (at $l=N$). ENAQT
occurs when $\mathcal{J}/\mathcal{J}_{0}>1$ in (a) or when a line segment in
(b) crosses increasing contour.

Results.—To capture both ENAQT and mobility edge features, we choose a
relatively short chain of $N=22$ so that dephasing does not completely
suppress the current, but still large enough to pronounce the mobility edge.
The ENAQT, defined by $\mathcal{J}/\mathcal{J}_{0}$, is measured by comparing
a current $\mathcal{J}$ relative to the one without dephasing,
$\mathcal{J}_{0}\equiv\mathcal{J}(\gamma=0)$; ENAQT is achieved when
$\mathcal{J}/\mathcal{J}_{0}>1$.

In Fig. 1(a) we report the relative current $\mathcal{J}/\mathcal{J}_{0}$ with
$\gamma$ for different $\alpha$’s corresponding to each dotted line in Fig.
1(b). We observe that between $\alpha=-1$ (uniform chain) and $\alpha=0$,
$\mathcal{J}/\mathcal{J}_{0}$ is almost always decreasing monotonically,
indicating little to no ENAQT. This behavior is expected because for an end-
to-end transport, ENAQT does not exist unless static disorder is introduced.
For negative $\alpha$, even though all eigenstates are delocalized, the weak
disorder may still reduce the spread of site population. This allows a
moderate dephasing to spread the population more evenly, producing small
ENAQT. For $\alpha>0$, the ENAQT is more pronounced through nonmonotonic
$\mathcal{J}/\mathcal{J}_{0}$. Moreover, as $\alpha$ is closer to $0.5$, the
range of $\gamma$’s for which the ENAQT appears increases significantly.

The parameter region in which ENAQT exists is illustrated more clearly in Fig.
2(a). If we consider the spread of the NESS population to be the universal
mechanism [11] underlying the ENAQT–higher spread yields higher current–then
the highest relative current should be at the highest $\alpha$, since higher
$\alpha$ gives more localized eigenstates and thus allowing the dephasing to
do more spreading (note that the spread of population in eigenbasis is same as
in the site basis as we will define the spread with basis independent
quantity). However, it turns out that the highest relative current is achieved
not at the highest $\alpha$ but at $\alpha\approx 0.6$. This indicates that
the spreading mechanism is not the whole story. Instead, the strength of ENAQT
must rely on the interplay between spreading and mobility edge.

Figure 2: (a) Contour plot of relative current $\mathcal{J}/\mathcal{J}_{0}$
as a function of $\alpha$ and $\gamma$. The dashed line corresponds to
$\mathcal{J}/\mathcal{J}_{0}=1$ and thus separates the regions with and
without ENAQT. (b) Fraction of localized eigenstates $f_{\mathrm{loc}}$ for
each $\alpha$. The shaded areas highlight the transition points at which
$f_{\mathrm{loc}}$ varies significantly. (c) Relative spread $\eta/\eta_{0}$
as a function of $\gamma$. Each line is for different value of $\alpha$ which
corresponds to horizontal lines in (a) with same color for each and increasing
thickness for increasing $\alpha$. The shaded area in (c) highlight the peaks
of $\eta/\eta_{0}$. The maximum for each $\alpha$ lines in (a) also occurs at
similar location with the maximum of $\eta/\eta_{0}$. The parameters are the
same as in Fig. 1.

To understand the role of the mobilty edge, we study the relative current
together with localization properties of the GAAH model. This is simply
revealed by the fraction of localized eigenstates $f_{\mathrm{loc}}$. A state
is localized if its energy is, for $\alpha>(<)\;0$, larger (lower) than the
SPME energy $E_{\mathrm{ME}}$. By comparing Fig. 2(a) and (b), we observe that
the contour width of ENAQT region varies quite significantly at the transition
points of $f_{\mathrm{loc}}$, that is highlighted by the shaded area around
$\alpha=0.4$ and $\alpha=0.9$. To study the interplay between spreading and
mobility edge, here we also calculate the population spread for several
$\alpha$’s, which we quantify using the normalized von Neumann entropy,
$\eta=-\sum_{k}p_{k}\log p_{k}/\log N$, where
$p_{k}=\langle\eta_{k}^{\dagger}\eta_{k}\rangle$ is given by Eq. (6). $\eta=0$
is when a single eigenvalue dominates the population and $\eta=1$ corresponds
to uniform population. This measure is has the advantage of being basis
independent. In Fig. 2(c) we compare the relative spread, $\eta/\eta_{0}$ with
$\eta_{0}\equiv\eta(\gamma=0)$, for $\alpha$ corresponding to transition
points of $f_{\mathrm{loc}}$ in Fig. 2(b). Note that the maximum ENAQT, the
peak of $\mathcal{J}/\mathcal{J}_{0}$, and the maximum $\eta/\eta_{0}$ have
different locations. The former is at $\alpha\approx 0.6$ and the latter is at
$\alpha\approx 0.3$; both occur at similar $\gamma\lesssim 10^{-2}$. Thus
higher relative spread does not guarantee higher ENAQT, since for higher
$\alpha$ the NESS is more localized, and consequently the dephasing can carry
out more spreading which results in higher ENAQT. However, for
$0.6<\alpha<0.8$ (the $f_{\mathrm{loc}}$ plateau), despite the stable
$f_{\mathrm{loc}}$, the relative spread is decreasing and thus also for the
relative current. In the plateau regime, the relative spread fully controls
the fate of ENAQT; it keeps decreasing until increases again due to the second
$f_{\mathrm{loc}}$ transition around $\alpha=0.9$. Moreover, $\eta/\eta_{0}$
also captures the $f_{\mathrm{loc}}$ transition and thus it is sensitive to
mobility edges at low $\gamma$’s (see Supplemental Material [47]). Hence, we
find that the cooperative interplay between population spread and localization
determines the fate of ENAQT.

The feature of ENAQT in this model can be explored further by varying the
phase $\phi$. This can be achieved in experiments typically by shifting the
phase of bichromatic laser field. The phase does not vary the number of
localized states, $Nf_{\mathrm{loc}}$, but it alters the spectrum $S_{nk}$.
Consequently, the position of the localized states and and the coupling
between two modes $C_{ik}$ in Eq. (Environment-Assisted Quantum Transport and
Mobility Edges) are changed, resulting a substantial change in ENAQT and
localization strength. To study the sensitivity of localization to the
dephasing we calculate the average inverse participation ratio (IPR) $\langle
I\rangle=\sum I_{k}p_{k}$ where $I_{k}=\sum_{n}|S_{nk}|^{4}$ is the IPR for
the individual eigenstates.

Figure 3:  For $\alpha=0.1,\gamma=0.006$ (red $\Diamond$; left y axis) and
$\alpha=0.9,\gamma=1$ (blue $\circ$; right y axis): (a) current $\mathcal{J}$
versus phase $\phi$, (b) maximum of coupling between two eigenmodes in
presence of dephasing $C_{ik}$, (c) spread $\eta$ as a function of $\phi$, and
(d) average inverse participation ratio $\langle I\rangle$. All plots are for
the common disorder strength $\lambda=0.4$. Shaded areas highlight the
correlations between the quantities (see text).

Figure 3(a) shows $\mathcal{J}$ as a function of $\phi$ for two cases:
$\alpha=0.1,\gamma=6\times 10^{-3}$ as the example for maximum ENAQT in the
delocalized $\alpha$ regime (red $\Diamond$; $f_{\mathrm{loc}}=0$), and
$\alpha=0.9,\gamma=1$ (blue $\circ$; $f_{\mathrm{loc}}\approx 0.5$) for the
localized regime. Both currents have peaks in different $\phi$’s and are
suppressed in overlapping range. For instance, both currents are heavily
suppressed around $\phi=\pi/2$ and $\phi=3\pi/2$. Note that here we do not
normalize $\mathcal{J}$ for ENAQT as it is already for specific $\alpha$ and
$\gamma$. For the delocalized regime, this can be explained solely by the
maximum coupling between two eigenstates, $\max(C_{ik})$. The maximum value is
always achieved for $i=k$ and thus it corresponds to
$\max(C_{kk})=\max_{k}(\sum_{n}|S_{nk}|^{4})$, i.e., the maximum IPR $I_{k}$.
We first turn our attention into the red $\Diamond$ line: in Fig. 3(b) we see
that the $\mathcal{J}$ is maximized whenever $\mathrm{max}(C_{ik})$ is small.
The current is also maximized whenever the spread $\eta$ is large (see Fig.
3(c)), confirming the previous result for low $\alpha$ regime that possesses
ENAQT. We then proceed to the blue $\circ$ line, where the $f_{loc}$ is
significantly higher than the previous one. Here the maximum coupling in Fig.
3(b) is no longer sensitive to $\phi$, since there is always at least a
localized state with IPR near unity for this case. On the other hand, the
corresponding population spread and IPR in Fig. 3(c)–3(d) simultaneously
determines the $\mathcal{J}$ peak, in contrast to the previous (delocalized)
case where $\mathcal{J}$ relies solely on the population spread. Further
analysis on the effect of $\alpha$ and $\phi$ to the spreading and
localization in extended parameter can be found in the Supplemental Material
[47].

Figure 4:  (a) Maximum ENAQT, $\mathcal{J}/\mathcal{J}_{0}$, versus $\lambda$
with $\phi=\pi/3$. (b) Dephasing rates at which maximum ENAQT in (a) occurs.
(c) Average IPR for $\alpha_{\mathrm{opt}}$ and $\gamma_{\mathrm{opt}}$.
Highlighted area indicates the values at localization transition.

To complete the phenomenology of this system, in Fig. 4(a) we plot the maximum
ENAQT (achieved at certain $\alpha_{\mathrm{opt}}$ and
$\gamma_{\mathrm{opt}}$) with disorder strength $\lambda$. It is possible to
obtain $\mathcal{J}/\mathcal{J}_{0}\approx 400$ for $\gamma\approx 10^{-5})$
(see Fig. 4(b)) at $\lambda\approx 2.1$. However, further increase in disorder
strength reduces the maximum ENAQT as all eigenstates is maximally localized
as indicated by $\langle I\rangle_{\mathrm{opt}}\equiv\langle
I\rangle(\alpha_{\mathrm{opt}},\gamma_{\mathrm{opt}})$ near unity in Fig. 4(c)
and the fact that all eigenstates are localized for $\lambda>1$ and positive
$\alpha$. The localization transition in $\lambda=1$ to $\lambda=1.5$ is
followed by a steep increase in the system sensitivity to dephasing by orders
of magnitude, whereas giving a sharp transition for the maximum ENAQT due to
the interplay between population spread and localization. In addition, small
increments in Fig. 4(b) for $\lambda=0$ to $\lambda=1$ is followed by a series
of local maxima in $\max(\mathcal{J}/\mathcal{J}_{0})$. The fate of the
maximum ENAQT along a $\gamma_{opt}$ plateau is determined by $\alpha$ that
gives the $f_{\mathrm{loc}}$ transition as in Fig. 2(a).

Conclusions.—We have studied how the presence of dephasing noise in a model
with single particle mobility edge coupled to spin baths can manifest strong
ENAQT. While the population spread alone may explain the origin of current
enhancement, in localized regimes it cannot fully account for the nonmonotonic
dependence of the mobility edge parameter and disorder strength. Instead, both
the population spread and average IPR play a cooperative role where for highly
localized eigenstates, the latter dominates the enhancement mechanism. Strong
ENAQT thus maximizes at a specific mobility edge parameter near the
$f_{\mathrm{loc}}$ transition and it can be tuned by shifting the position of
the localized modes (by tuning the phase) and modifying the disorder strength.

It should be noted that the current is actually higher in the low $\alpha$
regime, below the mobility edge. However, if disorder or mobilty edge are
unavoidably present, or if in an engineered system the mobility edge is
desired, the intermediate dephasing regime provides better performance. We
briefly note that if particle interaction is included, the ENAQT features
should similarly persist since the many-body localization occurs in the
Hamiltonian [49, 50].

F.P.Z. thanks Ministry of Higher Education and Research of Indonesia for
Research Funding P3MI ITB 2020. The numerical results were obtained using code
written in `NumPy` [51] and `QuTiP` [52], and the figures were made using
`matplotlib` [53].

## References

  * Mohseni  _et al._ [2008] M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, Environment-assisted quantum walks in photosynthetic energy transfer, J. Chem. Phys. 129, 11B603 (2008).
  * Rebentrost  _et al._ [2009a] P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru-Guzik, Environment-assisted quantum transport, New J. Phys. 11, 033003 (2009a).
  * Rebentrost  _et al._ [2009b] P. Rebentrost, M. Mohseni, and A. Aspuru-Guzik, Role of quantum coherence and environmental fluctuations in chromophoric energy transport, J. Phys. Chem. B 113, 9942 (2009b).
  * Plenio and Huelga [2008] M. B. Plenio and S. F. Huelga, Dephasing-assisted transport: quantum networks and biomolecules, New J. Phys. 10, 113019 (2008).
  * Chin  _et al._ [2012] A. Chin, S. Huelga, and M. B. Plenio, Coherence and decoherence in biological systems: principles of noise-assisted transport and the origin of long-lived coherences, Philos. Trans. Royal Soc. A 370, 3638 (2012).
  * Chin  _et al._ [2013] A. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment–protein complexes, Nature Physics 9, 113 (2013).
  * Kassal and Aspuru-Guzik [2012] I. Kassal and A. Aspuru-Guzik, Environment-assisted quantum transport in ordered systems, New J. Phys. 14, 053041 (2012).
  * Li  _et al._ [2015a] Y. Li, F. Caruso, E. Gauger, and S. C. Benjamin, ‘momentum rejuvenation’underlies the phenomenon of noise-assisted quantum energy flow, New Journal of Physics 17, 013057 (2015a).
  * Caruso  _et al._ [2009] F. Caruso, A. W. Chin, A. Datta, S. F. Huelga, and M. B. Plenio, Highly efficient energy excitation transfer in light-harvesting complexes: The fundamental role of noise-assisted transport, J. Chem. Phys. 131, 09B612 (2009).
  * Berman  _et al._ [2015] G. P. Berman, A. I. Nesterov, G. V. López, and R. T. Sayre, Superradiance transition and nonphotochemical quenching in photosynthetic complexes, J. Phys. Chem. C 119, 22289 (2015).
  * Zerah-Harush and Dubi [2018] E. Zerah-Harush and Y. Dubi, Universal origin for environment-assisted quantum transport in exciton transfer networks, J. Chem. Phys. Lett. 9, 1689 (2018).
  * Dwiputra  _et al._ [2020] D. Dwiputra, J. S. Kosasih, A. Sulaiman, and F. P. Zen, Driving-assisted open quantum transport in qubit networks, Phys. Rev. A 101, 012113 (2020).
  * Dwiputra  _et al._ [2019] D. Dwiputra, A. Sulaiman, J. S. Kosasih, W. Hidayat, and F. P. Zen, Driving the dephasing assisted quantum transport, in  _J. Phys. Conf. Ser._ , Vol. 1245 (IOP Publishing, 2019) p. 012075.
  * Sowa  _et al._ [2017] J. K. Sowa, J. A. Mol, G. A. D. Briggs, and E. M. Gauger, Environment-assisted quantum transport through single-molecule junctions, Phys. Chem. Chem. Phys. 19, 29534 (2017).
  * Biggerstaff  _et al._ [2016] D. N. Biggerstaff, R. Heilmann, A. A. Zecevik, M. Gräfe, M. A. Broome, A. Fedrizzi, S. Nolte, A. Szameit, A. G. White, and I. Kassal, Enhancing coherent transport in a photonic network using controllable decoherence, Nat. Comm. 7, 11282 (2016).
  * Trautmann and Hauke [2018] N. Trautmann and P. Hauke, Trapped-ion quantum simulation of excitation transport: Disordered, noisy, and long-range connected quantum networks, Phys. Rev. A 97, 023606 (2018).
  * Schempp  _et al._ [2015] H. Schempp, G. Günter, S. Wüster, M. Weidemüller, and S. Whitlock, Correlated exciton transport in rydberg-dressed-atom spin chains, Phys. Rev. Lett. 115, 093002 (2015).
  * Maier  _et al._ [2019] C. Maier, T. Brydges, P. Jurcevic, N. Trautmann, C. Hempel, B. P. Lanyon, P. Hauke, R. Blatt, and C. F. Roos, Environment-assisted quantum transport in a 10-qubit network, Phys. Rev. Lett. 122, 050501 (2019).
  * Dorfman  _et al._ [2013] K. E. Dorfman, D. V. Voronine, S. Mukamel, and M. O. Scully, Photosynthetic reaction center as a quantum heat engine, Proc. Natl. Acad. Sci. 110, 2746 (2013).
  * Scully  _et al._ [2011] M. O. Scully, K. R. Chapin, K. E. Dorfman, M. B. Kim, and A. Svidzinsky, Quantum heat engine power can be increased by noise-induced coherence, Proc. Natl. Acad. Sci. 108, 15097 (2011).
  * Damanet  _et al._ [2019] F. Damanet, E. Mascarenhas, D. Pekker, and A. J. Daley, Controlling quantum transport via dissipation engineering, Phys. Rev. Lett. 123, 180402 (2019).
  * Streltsov  _et al._ [2017] A. Streltsov, G. Adesso, and M. B. Plenio, Colloquium: Quantum coherence as a resource, Rev. Mod. Phys. 89, 041003 (2017).
  * Halpern  _et al._ [2019] N. Y. Halpern, C. D. White, S. Gopalakrishnan, and G. Refael, Quantum engine based on many-body localization, Phys. Rev. B 99, 024203 (2019).
  * Balachandran  _et al._ [2019] V. Balachandran, S. R. Clark, J. Goold, and D. Poletti, Energy current rectification and mobility edges, Phys. Rev. Lett. 123, 020603 (2019).
  * Saha and Maiti [2019] M. Saha and S. K. Maiti, Particle current rectification in a quasi-periodic double-stranded ladder, J. Phys. D Appl. Phys. 52, 465304 (2019).
  * Chiaracane  _et al._ [2020] C. Chiaracane, M. T. Mitchison, A. Purkayastha, G. Haack, and J. Goold, Quasiperiodic quantum heat engines with a mobility edge, Phys. Rev. Res. 2, 013093 (2020).
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Engel  _et al._ [2007] G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mančal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems, Nature 446, 782 (2007).
  * Panitchayangkoon  _et al._ [2010] G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, Long-lived quantum coherence in photosynthetic complexes at physiological temperature, Proc. Natl. Acad. Sci. 107, 12766 (2010).
  * Duan  _et al._ [2017] H.-G. Duan, V. I. Prokhorenko, R. J. Cogdell, K. Ashraf, A. L. Stevens, M. Thorwart, and R. D. Miller, Nature does not rely on long-lived electronic quantum coherence for photosynthetic energy transfer, Proc. Natl. Acad. Sci. 114, 8493 (2017).
  * Jang and Mennucci [2018] S. J. Jang and B. Mennucci, Delocalized excitons in natural light-harvesting complexes, Rev. Mod. Phys. 90, 035003 (2018).
  * Misra and Sudarshan [1977] B. Misra and E. G. Sudarshan, The zeno’s paradox in quantum theory, J. Math. Phys. 18, 756 (1977).
  * Aubry and André [1980] S. Aubry and G. André, Analyticity breaking and anderson localization in incommensurate lattices, Ann. Israel Phys. Soc 3, 18 (1980).
  * Harper [1955] P. G. Harper, Single band motion of conduction electrons in a uniform magnetic field, Proc. Phys. Soc. London Sect. A 68, 874 (1955).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Lüschen  _et al._ [2017] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017).
  * Lüschen  _et al._ [2018] H. P. Lüschen, S. Scherg, T. Kohlert, M. Schreiber, P. Bordia, X. Li, S. D. Sarma, and I. Bloch, Single-particle mobility edge in a one-dimensional quasiperiodic optical lattice, Phys. Rev. Lett. 120, 160404 (2018).
  * Biddle  _et al._ [2011] J. Biddle, D. Priour Jr, B. Wang, and S. D. Sarma, Localization in one-dimensional lattices with non-nearest-neighbor hopping: Generalized anderson and aubry-andré models, Phys. Rev. B 83, 075105 (2011).
  * Biddle and Sarma [2010] J. Biddle and S. D. Sarma, Predicted mobility edges in one-dimensional incommensurate optical lattices: An exactly solvable model of anderson localization, Phys. Rev. Lett. 104, 070601 (2010).
  * Biddle  _et al._ [2009] J. Biddle, B. Wang, D. Priour Jr, and S. D. Sarma, Localization in one-dimensional incommensurate lattices beyond the aubry-andré model, Phys. Rev. A 80, 021603 (2009).
  * Ganeshan  _et al._ [2015] S. Ganeshan, J. Pixley, and S. D. Sarma, Nearest neighbor tight binding models with an exact mobility edge in one dimension, Phys. Rev. Lett. 114, 146601 (2015).
  * Gelbwaser-Klimovsky and Aspuru-Guzik [2017] D. Gelbwaser-Klimovsky and A. Aspuru-Guzik, On thermodynamic inconsistencies in several photosynthetic and solar cell models and how to fix them, Chem. Sci. 8, 1008 (2017).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of n-level systems, J. Math. Phys. 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Comm. Math. Phys. 48, 119 (1976).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, 2002).
  * Guimaraes  _et al._ [2016] P. H. Guimaraes, G. T. Landi, and M. J. de Oliveira, Nonequilibrium quantum chains under multisite lindblad baths, Phys. Rev. E 94, 032139 (2016).
  * [47] .
  * Wu and Segal [2008] L.-A. Wu and D. Segal, Energy flux operator, current conservation and the formal fourier’s law, J. Phys. A Math. Theor. 42, 025302 (2008).
  * Li  _et al._ [2015b] X. Li, S. Ganeshan, J. Pixley, and S. D. Sarma, Many-body localization and quantum nonergodicity in a model with a single-particle mobility edge, Phys. Rev. Lett. 115, 186601 (2015b).
  * Modak and Mukerjee [2015] R. Modak and S. Mukerjee, Many-body localization in the presence of a single-particle mobility edge, Phys. Rev. Lett. 115, 230401 (2015).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, _et al._ , Array programming with numpy, Nature 585, 357 (2020).
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, Qutip 2: A python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234 (2013).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Comput. Sci. Eng. 9, 90 (2007).

Supplemental Material: Environment-Assisted Quantum Transport in a One-
Dimensional System with a Single-Particle Mobility Edge

Donny Dwiputra1 and Freddy P. Zen1,2

1Theoretical Physics Laboratory, Faculty of Mathematics and Natural Sciences,
Institut Teknologi Bandung, Jl. Ganesha 10, Bandung 40132, Indonesia

2Indonesian Center for Theoretical and Mathematical Physics (ICTMP), Bandung
40132, Indonesia

##  I I. Derivation of the master equation and its steady state

For the microscopic derivation of the master equation, we use the Born-Markov
and secular approximation to trace out the bath and write the Lindblad
dissipator. This is most readily done using the eigenoperators method [45]. We
first write the interaction Hamiltonians in eigenbasis $\eta_{k}$ where
$a_{n}=\sum_{k}S_{nk}\eta_{k}$ and $S$ is the $N\times N$ unitary matrix
obtained by diagonalizing $H_{S}$,

$\displaystyle H_{SB_{1}}$ $\displaystyle=$
$\displaystyle\sum_{l,\mu,k}g_{\mu}^{(1)}(S_{lk}\eta_{k}+S_{kl}^{*}\eta_{k}^{\dagger})(\sigma^{+}_{\mu,l}+\sigma^{-}_{\mu,l}),$
(S1) $\displaystyle H_{SB_{2}}$ $\displaystyle=$
$\displaystyle\sum_{n,\mu,k,k^{\prime}}g_{\mu}^{(2)}|S_{nk}|^{2}\eta_{k}^{\dagger}\eta_{k^{\prime}}(\sigma^{+}_{\mu,n}+\sigma^{-}_{\mu,n}).$
(S2)

In terms of the eigenoperators of $H_{S}$, we can write for each interaction
$H_{SB_{(1,2)}}=\sum_{\omega,k}A_{i,k}(\omega)\otimes B_{i}$ where
$i=\\{l,n\\}$ is the site index for $SB_{1}$ and $SB_{2}$ baths, respectively.
Here the decomposed operator
$A_{i,k}(\omega)\equiv\sum_{\epsilon^{\prime}-\epsilon=\omega}\Pi(\epsilon)A_{i,k}\Pi(\epsilon^{\prime})$
acts on the system, where $\Pi(\epsilon)$ is the projector onto eigenspace
belonging to the eigenvalue $\epsilon$ of $H_{S}$, and $B_{i}$ acts on the
bath. The operator $A_{i,k}(\omega)$ is chosen to satisfy

$[H,A_{i,k}(\omega)]=-\omega A_{i,k}(\omega).$ (S3)

It follows that, for the respective baths, these eigenoperator will be

$\displaystyle A_{l,k}(\omega)$ $\displaystyle=$ $\displaystyle
S_{nk}\eta_{k}\;\delta_{\omega,+\epsilon_{k}}+S_{kn}^{*}\eta_{k}^{\dagger}\;\delta_{\omega,-\epsilon_{k}},$
(S4) $\displaystyle A_{n,k}(\omega)$ $\displaystyle=$ $\displaystyle
S_{nk}S_{nk^{\prime}}\eta_{k}^{\dagger}\eta_{k^{\prime}}\;\delta_{\omega,\epsilon_{k^{\prime}}-\epsilon_{k}}.$
(S5)

while $B_{i}=\sum_{\mu}g_{\mu}^{(1,2)}(\sigma_{\mu,i}^{+}+\sigma_{\mu,i}^{-})$
for both $i=\\{l,n\\}$ [see Eq. (S1)]. Intuitively speaking, the coupling
$(\eta_{k}+\eta_{k}^{\dagger})$ and $\eta_{k}^{\dagger}\eta_{k}$ to the bath
induces transitions in the system with the allowed energy $\omega$.

For the each of above eigenoperators, the corresponding Lindblad dissipator is

$\mathcal{D}[A_{i,k}]\,\rho=\sum_{\omega}\Gamma_{i}(\omega)\left[A_{i,k}(\omega)\rho
A_{i,k}^{\dagger}(\omega)-\frac{1}{2}\\{A_{i,k}^{\dagger}(\omega)A_{i,k}(\omega),\rho\\}\right],$
(S6)

where

$\Gamma_{i}(\omega)=\int_{0}^{\infty}dt\;e^{i\omega
t}\mathrm{Tr}_{B}\left[(e^{iH_{S}t}B_{i}e^{-iH_{S}t})\,B_{i}\,\frac{e^{-H_{B}/T_{i}}}{\mathrm{Tr}(e^{-H_{B}/T_{i}})}\right]$
(S7)

is the half sided Fourier transform of bath correlation functions. Thus, for
$SB_{1}$ and $SB_{2}$ we get

$\Gamma_{i}(\omega)=\begin{cases}J_{i}(\omega)[1+n_{i}(\omega)],&\text{if
$\omega>0$}\\\ J_{i}(-\omega)n_{i}(-\omega),&\text{if $\omega<0$}\end{cases}$
(S8)

where $n_{i}(\omega)=(e^{\beta_{i}\omega}+1)^{-1}$ is the spin occupation
number for $SB_{1}$ ($SB_{2}$) bath coupled to site $l$ ($n$), and
$J_{i}(\omega)=\sum_{\mu}\pi|g_{\mu}^{(1,2)}|^{2}\delta(\omega-\epsilon_{\mu})$
is the spectral density. for $SB_{2}$, we denote the spectral density as
$\gamma_{n}(\omega)$. Thus, by writing for $\omega>0$ and $\omega<0$ in a
unified way, we can write the full dissipators as

$\displaystyle\mathcal{D}_{SB_{1}}\rho$ $\displaystyle=$
$\displaystyle\sum_{l,k}|S_{lk}|^{2}J(|\epsilon_{k}|)\Big{[}\left(1-n_{l}(\epsilon_{k})\right)(\eta_{k}\rho\eta_{k}^{\dagger}-\frac{1}{2}\\{\eta_{k}^{\dagger}\eta_{k},\rho\\})+n_{l}(\epsilon_{k})(\eta_{k}^{\dagger}\rho\eta_{k}-\frac{1}{2}\\{\eta_{k}\eta_{k}^{\dagger},\rho\\})\Big{]},$
(S9) $\displaystyle\mathcal{D}_{SB_{2}}\rho$ $\displaystyle=$
$\displaystyle\sum_{n,k,k^{\prime}}|S_{nk}|^{2}|S_{nk^{\prime}}|^{2}\gamma_{n}(|\epsilon_{k^{\prime}}-\epsilon_{k}|)\Big{[}(1-n_{n}(\epsilon_{k^{\prime}}-\epsilon_{k}))\theta(\epsilon_{k^{\prime}}-\epsilon_{k})+n_{n}(\epsilon_{k}-\epsilon_{k^{\prime}})\theta(\epsilon_{k}-\epsilon_{k^{\prime}})\Big{]}$
(S10)
$\displaystyle\times\left(\eta_{k}^{\dagger}\eta_{k^{\prime}}\,\rho\,\eta_{k^{\prime}}^{\dagger}\eta_{k}-\frac{1}{2}\left\\{\eta_{k}^{\dagger}\eta_{k^{\prime}},\rho\right\\}\right),$

where $\theta(\omega)$ is a step function. Following the main text, we assume
that the dephasing bath is at $T=0$ so that $n_{n}\rightarrow 0$. Hence, the
second term in Eq. (S10) vanishes and, since we arrange the eigenvalues
$\epsilon_{k}$ in increasing order, we can write the condition $k<k^{\prime}$
for the summation indices. To simplify even further, we take the common
$J_{l}=1$ and $\gamma_{n}=\gamma$. Hence, the $SB_{2}$ dissipator is

$\mathcal{D}_{SB_{2}}\rho=\sum_{n,k<k^{\prime}}|S_{nk}|^{2}|S_{nk^{\prime}}|^{2}\gamma\left(\eta_{k}^{\dagger}\eta_{k^{\prime}}\,\rho\,\eta_{k^{\prime}}^{\dagger}\eta_{k}-\frac{1}{2}\left\\{\eta_{k}^{\dagger}\eta_{k^{\prime}},\rho\right\\}\right).$
(S11)

To find the NESS, we consider the corresponding adjoint master equation with
$\mathcal{D}_{i,k}[A_{i,k}]^{\dagger}\,O=A_{i,k}^{\dagger}OA_{i,k}-1/2\\{A_{i,k}^{\dagger}A_{i,k},O\\}$
by replacing $O=\eta_{k}^{\dagger}\eta_{k^{\prime}}$ and solving $\langle
dO/dt\rangle=0$. The solution can be found by using only the fundamental
commutation relations which simplifies the adjoint dissipators considerably,

$\displaystyle\mathcal{D}^{\dagger}_{SB_{1}}\left(\eta_{k}^{\dagger}\eta_{k^{\prime}}\right)$
$\displaystyle=$
$\displaystyle\sum_{l}|S_{lk}|^{2}\left(n_{l}(\epsilon_{k})\delta_{kk^{\prime}}-[1-2n_{l}(\epsilon_{k})]\eta_{k}^{\dagger}\eta_{k^{\prime}}\right),$
(S12)
$\displaystyle\mathcal{D}^{\dagger}_{SB_{2}}\left(\eta_{k}^{\dagger}\eta_{k^{\prime}}\right)$
$\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{i,n}\left(|S_{nk}|^{2}+|S_{nk^{\prime}}|^{2}\right)|S_{ni}|^{2}\left(\eta_{i}^{\dagger}\eta_{i}\delta_{kk^{\prime}}\theta(\epsilon_{i}-\epsilon_{k})-\eta_{k}^{\dagger}\eta_{k^{\prime}}\right).$
(S13)

It follows that the NESS density matrix element is

$\langle\eta_{k}^{\dagger}\eta_{k^{\prime}}\rangle=\delta_{kk^{\prime}}\frac{\sum_{l}|S_{lk}|^{2}n_{l}(\epsilon_{k})+\gamma\sum_{i>k}\sum_{n}|S_{ni}|^{2}|S_{nk}|^{2}\langle\eta_{i}^{\dagger}\eta_{i}\rangle}{\sum_{l}|S_{lk}|^{2}(1-2n_{l}(\epsilon_{k}))+\gamma\sum_{i<k}\sum_{n}|S_{ni}|^{2}|S_{nk}|^{2}}.$
(S14)

It is apparent that only the highest energy population,
$\langle\eta_{N}^{\dagger}\eta_{N}\rangle$ is independent while lower energy
populations should be solved via iterative backward calculation. Once the NESS
is found, we can proceed to calculate the energy current with a similar
calculation.

##  II II. Population spread and localization properties of the steady state

In Fig. 2 of the main text we have studied the interplay between localization
and population spread to ENAQT. Here we extend the discussion using by
considering the average IPR $\langle I\rangle$ and the relative spread
$\eta/\eta_{0}$ over all $\alpha$’s and $\gamma$’s. The transitions of
dephasing-independent $f_{\mathrm{loc}}$ in Fig. 2(b) of the main text, at
$\alpha\approx 0.4$ (indicated by horizontal dotted lines) and less
prominently at $\alpha\approx 0.9$, is nicely depicted here with the average
IPR in the low dephasing regime until the quantum Zeno effect dominates the
transport for $\gamma\gg 10^{-2}$.

Figure S1:  Contour plots of (a) relative current
$\mathcal{J}/\mathcal{J}_{0}$, (b) average IPR $\langle I\rangle$, (c)
relative spread $\eta/\eta_{0}$, and (d) current $\mathcal{J}$ in logarithmic
scale as a function of $\alpha$ and $\gamma$. Vertical lines indicate the
value of $\gamma$ where the relative spread is maximized and horizontal lines
indicate $\alpha$ at which the average IPR (and also $f_{\mathrm{loc}}$)
undergo a transition. The parameters are the same as in Fig. 1 and 2 of the
main text ($N=22,\lambda=0.4$, $T_{h}=10^{3}$, $T_{l}=10^{-1}$, and
$\phi=\pi/3$).

The value of $\gamma$ when the relative spread $\eta/\eta_{0}$ is maximized is
indicated by the vertical dotted lines. By comparing Figs. S1(a) and S1(c), we
can see that the maximum of $\eta/\eta_{0}$ for each $\alpha$ is in general
coincides with the maximum ENAQT. To see this correlation more clearly, we can
infer from the fact that in Fig. S1(d) the vertical dotted line crosses the
inflection points of current $\mathcal{J}$ which indicate the maximum
enhancement for each constant $\alpha$ slice. We can simultaneously observe
that $\eta/\eta_{0}$ is sensitive to the mobility edge as it captures the
signature of $\langle I\rangle$ (and $f_{\mathrm{loc}})$ transition at
$\alpha\approx 0.4$ for the low $\gamma$ regime.

Figure S2:  Contour plots of current $\mathcal{J}$ in logarithmic scale with
$\lambda=0.4$ and $\phi=\pi/3$ for (a) $\alpha=0.1$ (delocalized regime;
$f_{\mathrm{loc}}=0$) and (b) $\alpha=0.9$ (localized regime;
$f_{\mathrm{loc}}\approx 0.5$). The vertical lines in (a) and (b) respectively
correspond to current for the red ($\gamma=0.06$) and blue ($\gamma=1$) lines
plotted in Fig. 3(a) of the main text.

Figures S2(a) and S2(b) show the current $\mathcal{J}$ for delocalized and
localized regimes, respectively. For clarity, we plot the current in
logarithmic scale—in contrast to linear scale as in Fig. 3(a) of the main
text. For the delocalized regime, current is maximized at low $\gamma$’s while
for the localized regime it is at high $\gamma$’s. This contrasting behavior
shows that the quantum Zeno effect may result in different consequences with
the presence of mobility edges.

##  III III. Comparison with the other measure of population spread

This section clarify the reason of why we choose normalized entropy as the
measure for population spread. In Fig. S3 we compare our measure $\eta$ with
the population spread $\Delta_{n}$ from Ref. [11], which is defined as

$\Delta_{n}=1-\left(\frac{1}{N}\sum_{i}n_{i}-n_{\mathrm{ext}}\right)^{2}=1-\left(\frac{1}{N}-\sum_{k}|S_{Nk}|^{2}p_{k}\right)^{2},$
(S15)

where $n_{\mathrm{ext}}$ is the site basis population at the extraction site
$N$. In energy basis it is $n_{\mathrm{ext}}=\sum_{k}|S_{Nk}|^{2}p_{k}$ where
$p_{k}$ is $k$-th energy population. Although $\Delta_{n}$ can mimic the
$\mathcal{J}$ (see Fig. S1(d)), its relative spread $\Delta_{n}/\Delta_{n0}$
($\Delta_{n0}\equiv\Delta_{n}(\gamma=0)$) fails to contain the maximum and
$f_{\mathrm{loc}}$ transition features of the relative current in Fig. S1(a)
particularly for $\alpha>0.4$—which $\eta/\eta_{0}$ does so relatively well.
Both features are crucial as they characterize the role of mobility edges in
sytems with ENAQT. In addition, the value $\eta$ gives a broader variance
while still bounded in $[0,1]$.

Figure S3:  Contour plots of (a) $\eta$, (b) $\Delta_{n}$, and (c) its
relative spread $\Delta_{n}/\Delta_{n0}$. The parameters are the same as in
Fig. S1

.


# Exotic Compact Objects with Two Dark Matter Fluids

Marie Cassing Institut für Theoretische Physik, Max-von-Laue-Straße 1, 60438
Frankfurt am Main, Germany  Alexander Brisebois Present Affiliation:
Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario
N2L 3G1, Canada  Division of Engineering Science, University of Toronto,
Toronto, Ontario M5S 2E4, Canada  Muhammad Azeem Department of Physics and
Astronomy, York University, Toronto, Ontario M3J 1P3, Canada  Jürgen
Schaffner-Bielich Institut für Theoretische Physik, Max-von-Laue-Straße 1,
60438 Frankfurt am Main, Germany

###### Abstract

The generic properties of compact objects made of two different fluids of dark
matter are studied in a scale invariant approach. We investigate compact
objects with a core-shell structure, where the two fluids are separated, and
with mixed dark matter components, where both dark matter fluids are immersed
within each other. The constellations considered are combinations of
incompressible fluids, free and interacting Fermi gases, and equations of
state with a vacuum term, i.e. self-bound dark matter. We find novel features
in the mass-radius relations for combined dark matter compact objects which
distinguishes them from compact objects with a single dark matter fluid and
compact stars made of ordinary baryonic matter, as white dwarfs, neutron stars
and quark stars. The maximum compactness of certain combined dark matter stars
can reach values up to the causality limit for compact stars but not beyond
that limit if causality of the dark matter fluids is ensured.

##  1 Introduction

The phenomena of dark matter are so far unexplained features from cosmological
observations and astrophysical data from the small and large scale structure
of the universe. Possible candidates of dark matter are weakly interacting
massive particles WIMPs, which are fermions, or axions and axion-like
particles. However, the standard cold dark matter (CDM) model of collisionless
free particles is in tension with the small scale structure observations of
galaxies as e.g. the core-cusp problem and the satellite problem.
Selfinteracting dark matter can provide a possible explanation of these
observations, see (Tulin & Yu, 2018).

If dark matter consists of selfinteracting particles they can form spheres of
a fluid which are bound by gravity, i.e. compact stars of dark matter. Boson
stars made of bosonic dark matter particles have been discussed for a long
time, first by a noninteracting scalar field (Wheeler, 1955; Kaup, 1968;
Ruffini & Bonazzola, 1969) and then by introducing selfinteractions (Colpi et
al., 1986) or a solitonic potential (Lee & Pang, 1987; Lee, 1987; Friedberg et
al., 1987a, b). The case for fermions, the properties of compact stars made of
selfinteracting fermionic dark matter compact stars, have been investigated in
addition, see e.g. (Narain et al., 2006; Kouvaris & Nielsen, 2015; Gresham &
Zurek, 2019). Stable compact stars which are bound without gravity can be
generated by a specific form of the equation of state where the pressure
vanishes at a nonvanishing energy density. These are so called self-bound
stars which have been discussed for absolutely stable strange quark matter
(Witten, 1984; Alcock et al., 1986; Haensel et al., 1986), for a review see
(Weber, 2005). As shown in (Chang et al., 2019) compact stars made of
interacting dark matter particles can be produced by cosmological structure
formation on various different scales depending on the interaction strength
and masses of the dark matter particles.

Dark matter can also be present in more sophisticated forms, as in the form of
multi-component dark matter which arises naturally in models beyond the
standard model (Zurek, 2009). There could be different kinds of particles, for
example in the form of composite dark matter (Khlopov & Kouvaris, 2008a, b).
Recently, the exploration of the properties of compact stars made of two
different fluids of dark matter has been put forward by studying dark white
dwarfs (Ryan & Radice, 2022) and exotic cores (Zöllner & Kämpfer, 2022).

The measurements of gravitational waves from binary mergers of compact objects
has the potential to observe and study exotic compact objects made of dark
matter (Cardoso et al., 2016; Maselli et al., 2017a, a, b; Mark et al., 2017;
Wystub et al., 2021). The LIGO-Virgo collaboration has published a catalog of
gravitational wave sources from binary mergers from their third run, GWTC-3
(Abbott et al., 2021). The catalog exhibits several masses of compact objects
which are lying in the mass gap between the mass of the lightest black hole,
set to $5M_{\odot}$, and the most massive neutron stars with a maximum mass of
about $2M_{\odot}$ from pulsar data (Cromartie et al., 2019; Fonseca et al.,
2021). A unique identification of a neutron star merger has been so far only
be possible for GW170817 by the observation of a jet and the optical afterglow
in the form of a kilonova (Abbott et al., 2020). Other neutron star-neutron
star or neutron star-black hole merger candidates are classified by their
masses. It is not clear if these gravitational wave events emerge from compact
objects involving ordinary neutron stars. The gravitational wave event
GW190814, for example, involves a compact object with a mass of $2.6M_{\odot}$
(Abbott et al., 2020), well above the most massive neutron star known at
present.

In this work we study the properties of spheres of dark matter for several
cases of two different kinds of fluids of dark matter. To keep our discussion
more general, we use scaling relations so that our results can be adopted for
arbitrary scales, be it masses, vacuum energies or interaction strengths. In
Section 2 the Tolman-Oppenheimer-Volkoff equation for gravitating spheres of
fluids is introduced in scale invariant form and its scaling relations are
discussed. The generic equations of state for an interacting gas of fermions,
self-bound matter and the analytic solution for an incompressible fluid are
summarized then. In Section 3 the results of our investigations are presented
for the following specific cases which can be seen in table 1: i) Two
incompressible fluids, ii) incompressible fluid core and an interacting Fermi
gas shell, iii) interacting Fermi gas core and a free Fermi gas shell, iv)
self-bound matter core and an interacting Fermi gas shell, and v) two
homogeneously mixed interacting Fermi fluids. The results of our study are
summarized in Section 4.

Model | Section | Description
---|---|---
1 | 3.1 | Two incompressible fluids
2 | 3.2 | Incompressible fluid core and an interacting Fermi gas shell
3 | 3.3 | Interacting Fermi gas core and a free Fermi gas shell
4 | 3.4 | Self-bound matter core and an interacting Fermi gas shell
5 | 3.5 | Star of two homogeneously mixed interacting Fermi fluids
Table 1: The different models of exotic compact objects consisting of two
fluids of dark matter studied in this work

##  2 Theoretical background

For a perfect fluid, the equations for hydrostatic equilibrium in general
relativity are denoted as the Tolman–Oppenheimer–Volkoff (TOV) equations and
can be derived by solving the Einstein-equations for a spherically symmetric
and static metric. The latter are given by the following differential
equations:

$\displaystyle\frac{dP}{dr}$
$\displaystyle=\frac{-Gm}{r^{2}}\rho\left(1+\frac{P}{\rho}\right)\left(1+\frac{4\pi
r^{3}P}{m}\right)\left(1-\frac{2Gm}{r}\right)^{-1}$ (1)
$\displaystyle\frac{dm}{dr}$ $\displaystyle=4\pi r^{2}\rho,$ (2)

where $P$ and $\rho$ denote the pressure and energy density, respectively, and
$G$ is the gravitational constants. Note that throughout the paper we are
using natural units, i.e. $\hbar=c=1$. For a given equation of state, which
provides a relation between $P$ and $\rho$, we can solve the above set of
equations employing appropriate initial and boundary conditions. We denote the
total radius of a star by $R$, which is found using the condition that the
pressure vanishes at the surface of the star ($P(R)=0$). The mass $M(r=0)$
must be zero at $r=0$ and $M(R)$ gives the total mass of the star at $r=R$. In
the model cases investigated here these quantities are found numerically by
solving the differential equations (1) and (2) for a star containing different
types of matter and therefore different equations of state.

###  2.1  Scaling the TOV equation

The TOV equations contain a dimensional quantity $G$, which can expressed in
terms of the Planck mass ($G=M_{p}^{-2})$. To transform the TOV equation to
scale independent variables we use the rescaling as introduced in (Narain et
al., 2006), see also (Schaffner-Bielich, 2020) (Throughout the paper we will
denote dimensional quantities with a prime). The rescaling factor
$\epsilon_{0}$ can be calculated for a given equation of state. The
dimensionless form of the equations facilitates the computational solutions
since one does not have to care about units,

$P=\epsilon_{0}P^{\prime},\qquad\epsilon=\epsilon_{0}\epsilon^{\prime}\,.$ (3)

Pressure and energy-density have the same units and thus can be rescaled by
the constant $\epsilon_{0}$ (for a Fermi gas of mass $m_{f}$ we use
$\epsilon_{0}=m_{f}^{4}$). Then the radius and the mass also have to be
rescaled:

$r=ar^{\prime},\qquad M=bM^{\prime}$ (4)

Plugging the new quantities into the TOV-equation we obtain a dimensionless
TOV-equation and relations for $a$ and $b$ which we can use below to obtain
these quantities in physical units,

$a=(G\epsilon_{0})^{-1/2},\qquad b=(G^{3}\epsilon_{0})^{-1/2}\,.$ (5)

We will now discuss the different types of equations of state (EOS) used as
input for the TOV-equation.

###  2.2 Constant Density Solution

An analytical solution for the above set of TOV-equations for the interior of
a star $r\leq R$ can be obtained for an incompressible fluid, i.e. for
constant density $\rho(r)=\rho_{0}$. In this case the solution for Eq. (2),
for $r\leq R$, is:

$M(r)=\frac{4\pi}{3}r^{3}\rho_{0}.$ (6)

With (6) and the boundary condition $P(r=R)=0$ the solution of the TOV-
equation is given by the Schwarzschild interior solution (Schwarzschild, 1916)
(see also (Misner et al., 1973)) :

$P(r)=\rho_{0}\left[\frac{R\sqrt{R-2GM}-\sqrt{R^{3}-2GMr^{2}}}{\sqrt{R^{3}-2GMr^{2}}-3R\sqrt{R-2GM}}\right].$
(7)

The central pressure (at $r=0$) is related to the mass $M$ by

$P(r=0)=\rho_{0}\left[\frac{\sqrt{1-\frac{2GM}{R}}-1}{1-3\sqrt{1-\frac{2GM}{R}}}\right].$
(8)

We define the compactness of a star by the mass to radius ratio,

$\displaystyle C=\frac{GM}{R}.$ (9)

Note that the central pressure, eq. (8), diverges for
$R\rightarrow\frac{9GM}{4}$, which is known as the Buchdahl-limit (Buchdahl,
1959) giving the maximum compactness $C_{max}=\frac{4}{9}$. Objects with a
compactness larger than the Buchdahl limit will be black holes. Another
constraint on the maximal mass and compactness is the causality limit
$C_{s}\approx 3/8$ (Schaffner-Bielich, 2020). It arises from the fact that the
speed of sound in the star has to be smaller then the speed of light $c$ to
fulfill causality and thus $c_{s}^{2}\leq 1$, i.e. by setting $P(r=0)=\rho$.
The speed of sound squared $c^{2}_{s}$ in a medium is defined by the
derivative of the pressure with respect to the energy density $\rho$ at
constant entropy $S$,

$c_{s}^{2}=\left(\frac{\partial P}{\partial\rho}\right)_{S}.$ (10)

Considering a polytropic equation of state the polytropic index $\Gamma=2$
will lead to the stiffest possible EoS and the limit $P=\rho$ with
$c_{s}^{2}=1$. The speed of sound for a fermionic type of matter – which is
used in four of our models – is $c_{s}=1/\sqrt{3}$ at high densities, thus
obeying causality.

###  2.3 Star with step density profile

In most cases studied in this paper, we consider a fluid density profile for a
star containing two fluids. We investigate stars with a core of fluid density
$\rho_{0}$ and dimensionless core radius $r_{c}$. Then the density profile
reads as follows:

$\displaystyle\rho(r)=\begin{cases}\rho_{0},&0\leq r\leq r_{c}\\\
\rho_{1},&r_{c}\leq r\leq R,\\\ 0&\text{else}\end{cases}.$ (11)

Here $\rho_{0}$ and $\rho_{1}$ can be a function of $P$ depending on the
equation of state for each fluid. For such a profile there is no longer an
analytical solution for the TOV equations (1) and (2) and the latter have to
be solved numerically (e.g. by the Runge-Kutta method).

###  2.4 Fermi gas with interactions

We use the same equation of state as in Ref. (Narain et al., 2006) for a Fermi
gas of mass $m_{f}$ with interactions at zero temperature. This EoS is given
as an explicit expression of $P$ and $\rho$:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{g}{2\pi^{2}}\int_{0}^{p_{F}}dpEp^{2}+\left(\frac{1}{3}\pi^{2}\right)^{2}y^{2}p_{F}^{6}=m_{f}^{4}\rho^{\prime}$
(12) $\displaystyle P$ $\displaystyle=$
$\displaystyle\frac{g}{2\pi^{2}}\int_{0}^{p_{F}}dp\frac{p^{4}}{3E}+\left(\frac{1}{3}\pi^{2}\right)^{2}y^{2}p_{F}^{6}=m_{f}^{4}P^{\prime}$
(13)

with $y=m_{f}/m_{I}$ being the interaction strength, where $m_{I}$ is the
scale of the interaction, and $p$ the momentum of the particles. The spin-
factor for $s=1/2$ fermions is $g=2$.

For comparison, the mass of a neutron is $m_{n}=0.938\,GeV/c^{2}$ and we use
the neutron mass or fractions of it to describe the different types of matter.
Using $m_{1}=m_{n}$ and some type of dark matter with $m_{2}=0.5m_{1}$ would
lead to a configuration similar to a white dwarf. For the rescaling of the
interacting Fermi gas we calculate the rescaling factor using equations (4)
and (5) to $\epsilon_{0}=m_{f}^{4}$ for dimensionless units.

###  2.5 Conformal equation of state with vacuum energy

Models with an non-vanishing vacuum energy are used to describe self-bound
stars. A well known example of such a vacuum model is the MIT-bag model which
is often used for the description of cold quark stars, see (Weber, 2005) for a
review. The MIT-bag model is a phenomenological model for mimicking quark
confinement and was introduced to describe hadron properties in terms of a bag
of quarks with a bag constant $B$. The vacuum energy density in this model
equals to 4 times the bag constant.

$\rho_{vac}=4B,\quad\text{ with }B^{1/4}\approx 145\,\text{MeV}.$ (14)

It allows self-bound stars to be stabilized by itself as without gravity. The
equation of state for the MIT bag model is

$P=\frac{1}{3}(\rho-\rho_{vac})\,.$ (15)

Such an equation of state is used in the fourth model case of our studies for
the description of a core of self-bound matter. It is related to the constant
speed of sound parameterization used for describing quark matter in compact
stars (Alford et al., 2013) with a speed of sound of $c_{s}^{2}=1/3$ and can
be motivated from perturbative QCD calculations for quark matter at high
density (Fraga et al., 2001). We will use eq. (15) generically as an EOS for
self-bound dark matter below. We will set $\rho_{vac}$ as the scaling variable
$\epsilon_{0}$, so that the results can be rescaled accordingly without the
need to specify the vacuum energy scale explicitly. The form of the EOS given
in eq. (15) emerges also from conformal field theories and is also known as
the conformal EOS.

##  3 Combined dark matter compact objects

###  3.1 Incompressible fluid stars

Figure 1: The dimensionless pressure $P^{\prime}$ of a star with a 2-step
density profile and a fixed central pressure $P^{\prime}_{c}=1$ displayed as a
function of the radial coordinate $r^{\prime}$ for two cases: i) for varying
core radius $r_{c}$ with fixed ratio $\rho_{1}/\rho_{0}$ (left) and ii) for
varying density ratios $\rho_{1}/\rho_{0}$ with fixed $r_{c}$ (right).

Figure 2: The mass-radius-relation for different central pressures $0\leq
P_{c}\leq 1$ for a star with a 2-step density profile containing two
incompressible fluids with a density-ratio of $\rho_{1}/\rho_{0}=0.5$ and
varying core-radii (left), and with a core-radius of $r_{c}=0.2$ and varying
density-ratios $\rho_{1}/\rho_{0}$ (right). The dashed and solid black lines
represent the causality and Buchdahl limits, respectively.

In this case we examine a star consisting of two layers of incompressible
fluids with densities $\rho_{0}$ and $\rho_{1}$. The density changes at a
chosen critical radius denoted by $r_{c}$. Two separate cases are explored in
this section (Fig. 1): i) the ratio between the constant densities is fixed to
$\rho_{1}/\rho_{0}=0.5$ while the core-radius $r_{c}$ is varied; ii) the
critical radius is fixed to $r_{c}=0.2$ while the ratio of the densities
$\rho_{1}/\rho_{0}$ is varied. To demonstrate the shape of a single pressure
profile throughout the star, the central pressure is fixed at $P_{c}=1$ and
the profiles for the two different cases are presented in Fig. 1.

The pressure profiles for case i) with a fixed density-ratio and varying
$r_{c}$ is shown on the left of Fig. 1 and are continuous throughout the star.
The slope of the pressure deviates from that of a single fluid star at the
chosen values of $r_{c}$ forming a kink. In case ii) with fixed $r_{c}$ and
varying density-ratio $\rho_{1}/\rho_{0}$ (r.h.s. of the figure) the slope of
the pressure similarly deviates at $r=r_{c}$, while the magnitude of the slope
decreases with dropping ratios $\rho_{1}/\rho_{0}$.

In Fig. 2 we display the relation between the total mass $M^{\prime}$ and the
total radius $R^{\prime}$ of the star for the two cases i) and ii), but for
different values of the central pressure $P^{\prime}_{c}$. For some
combinations of $r_{c}$ and $\rho_{1}/\rho_{0}$ we reach the causality limit,
which results in very compact configurations.

From the mass-radius relation of the $r_{c}$-varying case i) (left of Fig. 2),
lower values of $r_{c}$ result in a higher compactness for configurations with
high central pressures $P_{c}$ (some reaching the causality limit). In case
ii) (right of Fig. 2) with varying ratio $\rho_{1}/\rho_{0}$, a higher
$\rho_{1}/\rho_{0}$ will result in a star with a higher compactness, which
implies that configurations with larger shell density are more likely to reach
the causality limit. One also notices that configurations with a small density
ratio reach higher maximal masses and larger radii for high central pressure.

###  3.2 Incompressible fluid and Fermi gas star

Figure 3: The dimensionless maximal mass $M^{\prime}$ as a function of the
radius $R^{\prime}$ for different initial values of the central pressure
$P^{\prime}_{c}$ (left) and the dimensionless central pressure
$P^{\prime}_{c}$ as a function of the radius of the star $R^{\prime}$ (right)
for a star combined out of an incompressible fluid up to different core radii
$0.1\leq r_{c}\leq 0.3$ and a Fermi gas with interactions ($y=10$) as a shell.
Again the dashed and solid lines represent the causality and Buchdahl limit.

In this model we consider an incompressible fluid core and a shell with an
interacting Fermi gas. In our example the mass of the fermion particles is
$m_{f}$, which might be some kind of dark matter with interaction strength
$y=10$. In the core the energy density of the incompressible fluid is set to
be equal to the maximal possible value of the central pressure
$P_{c,max}=\epsilon_{0}=1$. The TOV equations then are solved for different
values of the central pressure, and for small $P_{c}$ we would see a first
order phase transition in the energy density. An analytical solution of the
TOV equations can be found for the core of the star. The mass as a function of
the radius then can be calculated as follows.

$M=\frac{4\pi}{3}\epsilon_{0}r^{3}\,.$ (16)

The compactness at the core should not be larger than the Buchdahl limit
$M_{c}/r_{c}\leq 4/9$ and inserting eq. (16) we get a limit on the possible
core radii in dimensionless units $r^{\prime}_{c}\leq\sqrt{1/(3\pi)}\approx
0.33$. Another constraint on the core radius comes from the causality limit
$M_{c}/r_{c}\leq 3/8$. In dimensionless units the core radius should be
smaller than $r^{\prime}_{c}\leq\sqrt{9/(32\pi)}\approx 0.3$. An analytical
solution for the pressure in the core of the star can also be obtained in the
following form:

$P_{1}(r)=\epsilon\left[\frac{\epsilon_{0}+P_{c}-(3P_{c}+\epsilon_{0})\sqrt{1-\frac{8\pi}{3}\epsilon_{0}r^{2}}}{(\epsilon_{0}+3P_{c})\sqrt{1-\frac{8\pi}{3}\epsilon_{0}r^{2}}-3(\epsilon_{0}+P_{c})}\right]\,.$
(17)

Now for a specific value of the central pressure $P_{c}$ and a selected core
radius $r_{c}$ the pressure at the core radius should be continuous and we
would start to solve for the interacting Fermi gas starting from the pressure
$P_{1}(r_{c})$. In Fig. 3 we see the typical shape of the mass-radius-curves.
We obtain a swirl and an unstable branch for large central pressures due to
the behaviour of the interacting Fermi gas. To the right of the maximal mass
$M_{TOV}$ we drop to the stable branch with smaller central pressures. It can
be seen that a larger core radius and therefore a larger amount of the
incompressible fluid leads to a smaller maximal mass $M_{TOV}$. It can be
explained such that interactions in the Fermi gas allow for more compact
configurations and increasing the core radius decreases the amount of the
surrounding Fermi gas. For low pressures the mass has a turning point and
curves to the left in the direction of small radii $R$. This branch appears
due to the presence of the incompressible fluid core which results in the
Fermi solution for larger central pressures. For larger core radii this
transition happens smoothly and for smaller core-radii we find a sharper
turning point. On the right hand side of Fig. 3 we see that for some small
values of the central pressure $P_{c}$ the total radius $R$ increases very
fast, such that we obtain an almost horizontal line for the central pressure.
Therefore for small central pressures we first obtain small total radii and
then large total radii of the star. At the point where the Fermi gas dominates
over the incompressible fluid and the mass-curve has a turning point the total
radius of the star $R$ decreases for increasing $P_{c}$. It then turns up and
back to a bit larger value of $R$ staying almost constant afterwards for
increasing $P_{c}$.

###  3.3 Two Fermi gases star

Figure 4: The dimensionless maximal mass $M^{\prime}$ as a function of the
radius $R^{\prime}$ for different initial central pressure $P^{\prime}_{c}$
(left) and the dimensionless central pressure $P^{\prime}_{c}$ as a function
of the radius of the star $R^{\prime}$ (right) for a star combined out of a
Fermi gas $m_{1}$ with interactions ($y=10$) as a core up to $r^{\prime}_{c}$
followed by a free Fermi gas of dark matter particles $m_{2}=10^{-6}m_{1}$.
The solid black line represents the Buchdahl limit and the dashed one the
causality limit.

In this model we consider a star with an interacting Fermi gas core and a free
Fermi gas shell. For this combination we use an interaction strength of $y=10$
and a particle-mass-ratio of $m_{2}/m_{1}=10^{-6}$ where the mass of the first
matter type $m_{1}$ is our scaling unit, i.e. the nucleon-mass, and the mass
of the particles in the shell $m_{2}$ is very small and could be considered as
axinos or some kind of warm dark matter. In Fig. 4 we present the
dimensionless maximal masses $M^{\prime}$ and total radii $R^{\prime}$ for
different central pressures $P_{c}$. The mass-radius-relation is only
displayed for the stable solutions up to the maximal mass of the star.

We notice that for the variation of the particle-mass-ratio $m_{2}/m_{1}$ at a
fixed interaction-strength $y$ the mass-radius-relation or behaviour of the
central pressure almost doesn’t change. The interaction parameter $y$
dominates and determines the behaviour as can be seen from the theoretical
expression (13).

From the mass-radius-relation we find that for small matching radii $r_{c}$
the star almost behaves as a free Fermi gas and the slope of the mass radius
relation roughly corresponds to the one for an EOS with a polytropic index of
$\Gamma=5/3$, i.e. to the one of a fermion star without interactions. When
increasing the core-radius $r_{c}$ the slope of the mass changes and we are
able to obtain more compact configurations. For large values of $r_{c}$ the
interacting Fermi gas is the dominant component in the star and we find a
change in the slope of the mass-radius curve. Then the slope comes close to
that of a relativistic Fermi gas with $\Gamma=4/3$.

###  3.4 Self-bound matter and Fermi gas star

Figure 5: The dimensionless mass $M^{\prime}$ (left) and pressure
$P^{\prime}$(right) as a function of the radius $R^{\prime}$ for different
initial values of the central pressure $P_{c}$ for a star combined out of a
self-bound matter core and a Fermi gas with interactions ($y=10$) as a shell
for different core-radii $r^{\prime}_{c}$. Dashed and solid black lines are
the causality and Buchdahl limit.

In this example we investigate a core described by the equation of state of
self-bound matter, eq. (15). The shell contains an interacting Fermi gas with
particle mass $m_{2}=10^{-6}m_{1}$ and interaction strength $y=10$ which could
represent some type of dark matter. The core can be seen as dark matter coming
from e.g. some dark QCD theory. Then some dark quarks are the components of
the dark self-bound stars core. In Fig. 5 the dimensionless mass and pressure
can be seen as a function of the total radius of the star for different core
radii $r_{c}$. We present only the stable solutions of this configuration. The
total mass $M^{\prime}$ slowly increases for increasing central pressure
$P^{\prime}_{c}$ whereas the total radius $R^{\prime}$ increases fast similar
to the behaviour seen in model case two. In this branch the star is dominated
by the self-bound matter component then smoothly turning into the branch
dominated by the interacting Fermi gas. There the total radius slowly
decreases again for high pressures while the mass keeps increasing up to
$M_{TOV}$. Afterwards the mass radius relation would produce a swirl due to
the fermionic components of the shell resulting in an unstable branch which is
not shown in Fig. 5. An unusual behaviour is found in the pressure for small
core-radii as for the example of $r_{c}=1$. There the pressure increases very
slow but the total radius of the star increases a lot and turns back in the
direction of smaller radii for increasing central pressures. A similar
behaviour has been seen for strange dwarfs, see (Glendenning et al., 1995a, b)
which consist of a core of absolutely stable strange quark matter and a halo
of white dwarf matter. For a discussion of the stability analysis of strange
dwarfs we refer to (Alford et al., 2017; Di Clemente et al., 2022) noting that
it does not directly apply to the case studied here with a constant core
radius.

###  3.5 Mixed Fermi star

Figure 6: The maximal mass (left) and central pressure (right) as a function
of the total radius $R^{\prime}$ for the mixed fluid containing an interacting
Fermi gas with the particle-mass of $m_{2}=10^{-6}m_{1}$ and an interacting
Fermi gas with particle mass $m_{1}$ and different interaction strengths
$y=0,\,10,\,100,\,1000$ . The solid lines represent the stable branches while
the dashed lines correspond to the unstable ones. The solid black line is the
Buchdahl limit and the dashed black one is the causality limit.

Another case studied is a star with two homogeneously mixed components. The
number densities of the two components are taken to be equal $n_{1}=n_{2}$.
Such a case can be realized for the scenario of dark fermions with opposite
dark charges, as studied for dark white dwarfs in (Ryan & Radice, 2022). In
our case these two components are interacting Fermi gases, where one has a
particle mass $m_{1}$ and a small particle mass $m_{2}=10^{-6}m_{1}$. In Fig.
6 the mass-radius-relation and the pressure are presented in dimensionless
units on a double-logarithmic scale. The mass radius relation reproduces the
expected behaviour of a single interacting Fermi gas star. The stable branches
with increasing mass up to $M_{TOV}$ can be well seen and the total mass of
the star as well as $M_{TOV}$ increases for increasing interaction parameter
$y$. In this model, we kept the unstable solutions in the Fig. 6 (dashed
lines) to show the typical behaviour of the Fermi gases and the similarities
of the solutions for high central pressures. The mass radius relation - making
a swirl at high central pressures - corresponds to the wiggles seen in the
pressure on the right. Although the particle masses of the two components
greatly differ, the behaviour is similar to that of a single interacting Fermi
gas. Also, changing the ratio of the particle masses doesn’t affect the
behaviour since it is dominantly determined by the interaction strength.

###  3.6 Compactness

Figure 7: The compactness at $M_{TOV}$ the maximal mass for the five different
models with different core radii $r_{c}$ and/or mass ratios.

Fig. 7 summarizes our results. There the dimensionless compactness $C$ is
presented as a function of the total radius of the star $R^{\prime}$ for all
five model cases studied in this work. In several models we have investigated
stars of two different types of matter combined at different values of the
core radius $r_{c}$. For each star with different core-radius we present here
the compactness at the maximal total mass $M_{TOV}$ which corresponds to the
most compact configuration for different $r_{c}$.

It can be seen that the stars combined of two incompressible fluids are the
most compact ones. As we have already observed before for very small core
radii $r_{c}$ at a fixed density ratio $\rho_{1}/\rho_{0}$, and large density
ratios $\rho_{1}/\rho_{0}$ at fixed core radius the compactness almost reaches
the causality limit of $C=3/8$. It is a novel feature and finding that such a
simple two-step-density profile configuration can reach such high values of
the compactness.

The second most compact combination is that of model case 3 where the star has
an interacting Fermi gas core and a free Fermi gas shell consisting of very
light particles. This shows that even with a light component we can obtain
compact configurations and we find that the behaviour is dominated more by the
interaction strength than by the particle mass.

In the list of most compact configurations follows model case 4 which
corresponds to a star containing a self-bound matter core and an interacting
Fermi gas shell. The configuration with the largest amount of self-bound
matter has the highest compactness. This model is followed by model 2 which
corresponds to a a star with an incompressible fluid core and an interacting
Fermi gas shell. The configuration with the smallest core radius is the most
compact one.

The least compact configurations are obtained in model 5 in which the two
interacting Fermi gases with a large particle mass ratio are mixed together.
The mixed fluid is very insensitive to the particle mass but changes with the
interaction strength of the Fermi gas as can be seen from the results in
Figure 6. Not all explored cases of the mixed star are shown in Fig. 6. The
missing cases with a very high interaction strength of $y=50,\,100$ have very
large total radii of the star. We observe that the radius increases greatly
with interaction strength, but the compactness stays almost the same and
increases only slightly.

##  4 Summary

We have investigated 5 different model cases for the composition of exotic
stars and obtained novel results for the mass-radius relation of compact stars
made of fluids of dark matter. In the model cases we have combined two
different types of matter with different particle mass ratios at different
core-radii. We focused on scenarios with large particle mass ratios where the
lighter component of the star can be thought as some sort of another dark
matter particle. Such a light component, for example axinos, and the core-
shell-structure of the star have an impact on the total mass and the
compactness of the star.

We find that the mass-radius relation for exotic compact objects with two
different fluids of dark matter can exhibit unusual and sophisticated
patterns. For comparison with our results specific for each model studied in
this work listed below we note that the characteristic mass-radius relation
for white dwarfs and neutron stars follows usually simple typical relations.
For white dwarfs and low mass neutron stars the mass simply decreases with the
radius cubed as dictated from the properties of a free Fermi gas equation of
state. Single fluid selfbound stars, for an equation of state with a vacuum
term, show as a simple characteristic an increasing mass with radius cubed.

First we have considered a model with two incompressible fluids with different
density-ratios and different core radii which correspond to a two-step density
profile. There we have obtained exotic mass-radius relations and find
extremely compact configurations almost reaching the causality limit. We have
investigated the case in which we add an interacting Fermi gas shell on the
core of an incompressible fluid. The Fermi gas has particles with mass $m_{f}$
and is strongly interacting. We have found interesting shapes of the mass-
radius-relation with very small total radii and masses for small central
pressures. After the maximal total radius, the branch dominated by the
incompressible fluid then turns to a branch dominated by the interacting Fermi
gas which shows the typical behaviour of a Fermi gas. We have found that
increasing the core radius, i.e. the amount of incompressible fluid, the
maximal mass of the star $M_{TOV}$ decreases.

This happens due to the larger amount of the non-interacting incompressible
fluid leading to less compact configurations compared to a simple interacting
Fermi gas star.

The third model case we have investigated contains two Fermi gases. The core
has components with a particle mass of $m_{1}$ and strong interactions. The
shell is constructed of non-interacting particles with a small mass
$m_{2}=10^{-6}m_{1}$. Such a combination of the different components leads to
typical behaviours of the compact star configuration related to the Fermi gas
in the relativistic and non-relativistic limit. It is new that we obtain
configurations with a high compactness for some values of the core radius
$r_{c}$. The total radius of the star is decreasing for an increasing value of
the central pressure.

In model case four we investigated a star with a self-bound matter core,
described by a nonvanishing vacuum energy density and an interacting Fermi gas
shell with strength $y=10$. For small central pressures the star is dominated
by the self-bound component and the total radius increases very fast while the
central pressure and total mass increase slowly. This behaviour is similar to
model case two in which for small central pressures the star is dominated by
the component of the core and the central pressure looks almost horizontally
for small values. A similar behaviour has already been seen in strange dwarfs.
At some larger total radius of the star the total mass and pressure turn to
smaller radii and increase in value. There the star is dominated by the Fermi
gas.

The last model case explored consists of two interacting Fermi gases with
different interaction strength mixed together. One component has a particle
mass of $m_{1}$ and the second component of the star in contrast has a small
particle mass $m_{2}=10^{-6}m_{1}$. We have found that adding a component with
such a small particle mass still reproduces the typical known behaviour of a
Fermi star with an interacting Fermi gas. The total mass of the star and
$M_{TOV}$ increases with increasing interaction strength.

All the cases considered here lead to quite compact configurations. The cases
containing two incompressible fluids give the most compact configurations,
even getting very close to the causality limit. In conclusion, we can say that
combining two different types of matter where either one or both of them could
be dark matter gives interesting configurations in mass and radius when
compared to known single fluid dark matter models. The additional dark matter
component has an impact on the total mass and the total radius of the star,
therefore making more compact configurations possible. The mass-radius
relation can be substantially changed in some cases studied leading to unique
features distinctly different to the mass-radius relation for ordinary matter,
such as neutron stars and white dwarfs. We point out that this is an
exploratory work for the study of two-component dark matter stars which are
encouraging. Certainly, the whole plethora of combinations of two dark matter
fluids in exotic compact objects has just been touched upon. Further
investigations of exotic compact objects with different types of dark matter
will be considered in future work.

This work emerged as a follow-up from the German-Canadian student research
collaboration in theoretical astroparticle physics ’EXPeriential Learning
Opportunity through Research and Exchange’ EXPLORE of the York University, the
University of Toronto, and the Goethe University of Frankfurt in the summer
term of 2021. We gratefully thank Nassim Borzognia, Saeed Rastgoo, Laura
Sagunski, Sean Tulin for initiating and setting up this research collaboration
and all participants of EXPLORE for discussions and a fruitful collaboration.
We thank Sean Tulin for a careful reading of the manuscript. We acknowledge
support by the Academic Innovation Fund of the York University, Toronto,
Canada, by the QSL (Quality Assurance in Teaching) of Goethe University
Frankfurt, by the Research Cluster ’ELEMENTS: Exploring the Universe from
Microscopic to Macroscopic Scales’ of the State of Hesse, and by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR
211 ’Strong-interaction matter under extreme conditions’ – project number
315477589-TRR-211.

## References

  * Abbott et al. (2020) Abbott, R., et al. 2020, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2021) —. 2021.  https://arxiv.org/abs/2111.03606
  * Alcock et al. (1986) Alcock, C., Farhi, E., & Olinto, A. 1986, Astrophys. J., 310, 261, doi: 10.1086/164679
  * Alford et al. (2013) Alford, M. G., Han, S., & Prakash, M. 2013, Phys. Rev. D, 88, 083013, doi: 10.1103/PhysRevD.88.083013
  * Alford et al. (2017) Alford, M. G., Harris, S. P., & Sachdeva, P. S. 2017, Astrophys. J., 847, 109, doi: 10.3847/1538-4357/aa8509
  * Buchdahl (1959) Buchdahl, H. A. 1959, Phys. Rev., 116, 1027, doi: 10.1103/PhysRev.116.1027
  * Cardoso et al. (2016) Cardoso, V., Hopper, S., Macedo, C. F. B., Palenzuela, C., & Pani, P. 2016, Phys. Rev. D, 94, 084031, doi: 10.1103/PhysRevD.94.084031
  * Chang et al. (2019) Chang, J. H., Egana-Ugrinovic, D., Essig, R., & Kouvaris, C. 2019, JCAP, 03, 036, doi: 10.1088/1475-7516/2019/03/036
  * Colpi et al. (1986) Colpi, M., Shapiro, S. L., & Wasserman, I. 1986, Phys. Rev. Lett., 57, 2485, doi: 10.1103/PhysRevLett.57.2485
  * Cromartie et al. (2019) Cromartie, H. T., et al. 2019, Nature Astron., 4, 72, doi: 10.1038/s41550-019-0880-2
  * Di Clemente et al. (2022) Di Clemente, F., Drago, A., Char, P., & Pagliara, G. 2022.  https://arxiv.org/abs/2207.08704
  * Fonseca et al. (2021) Fonseca, E., et al. 2021, Astrophys. J. Lett., 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Fraga et al. (2001) Fraga, E. S., Pisarski, R. D., & Schaffner-Bielich, J. 2001, Phys. Rev. D, 63, 121702, doi: 10.1103/PhysRevD.63.121702
  * Friedberg et al. (1987a) Friedberg, R., Lee, T. D., & Pang, Y. 1987a, Phys. Rev. D, 35, 3640, doi: 10.1103/PhysRevD.35.3640
  * Friedberg et al. (1987b) —. 1987b, Phys. Rev. D, 35, 3658, doi: 10.1103/PhysRevD.35.3658
  * Glendenning et al. (1995a) Glendenning, N. K., Kettner, C., & Weber, F. 1995a, Phys. Rev. Lett., 74, 3519, doi: 10.1103/PhysRevLett.74.3519
  * Glendenning et al. (1995b) —. 1995b, Astrophys. J., 450, 253, doi: 10.1086/176136
  * Gresham & Zurek (2019) Gresham, M. I., & Zurek, K. M. 2019, Phys. Rev. D, 99, 083008, doi: 10.1103/PhysRevD.99.083008
  * Haensel et al. (1986) Haensel, P., Zdunik, J. L., & Schaeffer, R. 1986, Astron. Astrophys., 160, 121
  * Kaup (1968) Kaup, D. J. 1968, Phys. Rev., 172, 1331, doi: 10.1103/PhysRev.172.1331
  * Khlopov & Kouvaris (2008a) Khlopov, M. Y., & Kouvaris, C. 2008a, Phys. Rev. D, 77, 065002, doi: 10.1103/PhysRevD.77.065002
  * Khlopov & Kouvaris (2008b) —. 2008b, Phys. Rev. D, 78, 065040, doi: 10.1103/PhysRevD.78.065040
  * Kouvaris & Nielsen (2015) Kouvaris, C., & Nielsen, N. G. 2015, Phys. Rev. D, 92, 063526, doi: 10.1103/PhysRevD.92.063526
  * Lee (1987) Lee, T. D. 1987, Phys. Rev. D, 35, 3637, doi: 10.1103/PhysRevD.35.3637
  * Lee & Pang (1987) Lee, T. D., & Pang, Y. 1987, Phys. Rev. D, 35, 3678, doi: 10.1103/PhysRevD.35.3678
  * Mark et al. (2017) Mark, Z., Zimmerman, A., Du, S. M., & Chen, Y. 2017, Phys. Rev. D, 96, 084002, doi: 10.1103/PhysRevD.96.084002
  * Maselli et al. (2017a) Maselli, A., Pnigouras, P., Nielsen, N. G., Kouvaris, C., & Kokkotas, K. D. 2017a, Phys. Rev. D, 96, 023005, doi: 10.1103/PhysRevD.96.023005
  * Maselli et al. (2017b) Maselli, A., Völkel, S. H., & Kokkotas, K. D. 2017b, Phys. Rev. D, 96, 064045, doi: 10.1103/PhysRevD.96.064045
  * Misner et al. (1973) Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W. H. Freeman)
  * Narain et al. (2006) Narain, G., Schaffner-Bielich, J., & Mishustin, I. N. 2006, Phys. Rev. D, 74, 063003, doi: 10.1103/PhysRevD.74.063003
  * Ruffini & Bonazzola (1969) Ruffini, R., & Bonazzola, S. 1969, Phys. Rev., 187, 1767, doi: 10.1103/PhysRev.187.1767
  * Ryan & Radice (2022) Ryan, M., & Radice, D. 2022.  https://arxiv.org/abs/2201.05626
  * Schaffner-Bielich (2020) Schaffner-Bielich, J. 2020, Compact star physics (Cambridge: Cambridge University Press), doi: 10.1017/9781316848357
  * Schwarzschild (1916) Schwarzschild, K. 1916, in Sitzungsber. K. Preuss. Akad. Wiss., 424–434
  * Tulin & Yu (2018) Tulin, S., & Yu, H.-B. 2018, Phys. Rept., 730, 1, doi: 10.1016/j.physrep.2017.11.004
  * Weber (2005) Weber, F. 2005, Prog. Part. Nucl. Phys., 54, 193, doi: 10.1016/j.ppnp.2004.07.001
  * Wheeler (1955) Wheeler, J. A. 1955, Phys. Rev., 97, 511, doi: 10.1103/PhysRev.97.511
  * Witten (1984) Witten, E. 1984, Phys. Rev. D, 30, 272, doi: 10.1103/PhysRevD.30.272
  * Wystub et al. (2021) Wystub, S., Dengler, Y., Christian, J.-E., & Schaffner-Bielich, J. 2021.  https://arxiv.org/abs/2110.12972
  * Zöllner & Kämpfer (2022) Zöllner, R., & Kämpfer, B. 2022.  https://arxiv.org/abs/2203.17228
  * Zurek (2009) Zurek, K. M. 2009, Phys. Rev. D, 79, 115002, doi: 10.1103/PhysRevD.79.115002