afqshort = AFQ, long = antiferroquadrupolar scshort = SC, long =
superconductivity qcpshort = QCP, long = quantum critical point cefshort =
CEF, long = crystal electric field ltlshort = LTL, long = Low Temperature
Laboratory hfmlshort = HFML, long = High Field Magnet Laboratory liashort =
LIA, long = lock-in amplifier ultshort = ULT, long = ultra-low temperature

# Diverse influences of hyperfine interactions on strongly correlated electron
states

Femke Bangma  High Field Magnet Laboratory (HFML-EMFL), Radboud University,
Toernooiveld 7, Nijmegen, 6525 ED, The Netherlands  Radboud University,
Institute for Molecules and Materials, Heyendaalseweg 135, 6525 AJ Nijmegen,
The Netherlands.  Lev Levitin  Department of Physics, Royal Holloway
University of London, Egham, TW20 0EX, Surrey, United Kingdom  Marijn Lucas
Department of Physics, Royal Holloway University of London, Egham, TW20 0EX,
Surrey, United Kingdom  Andrew Casey  Department of Physics, Royal Holloway
University of London, Egham, TW20 0EX, Surrey, United Kingdom  Jan Nyéki
Department of Physics, Royal Holloway University of London, Egham, TW20 0EX,
Surrey, United Kingdom  Ineke Broeders  High Field Magnet Laboratory (HFML-
EMFL), Radboud University, Toernooiveld 7, Nijmegen, 6525 ED, The Netherlands
Aaron Sutton  Department of Physics, University of Toronto, 60 St. George
Street, Toronto, M5S 1A7, Ontario, Canada  Bohdan Andraka  Department of
Physics, University of Florida, P.O. Box 118440, Gainesville, Florida
32611-8440, USA  Stephen Julian  Department of Physics, University of Toronto,
60 St. George Street, Toronto, M5S 1A7, Ontario, Canada  John Saunders
Department of Physics, Royal Holloway University of London, Egham, TW20 0EX,
Surrey, United Kingdom  Alix McCollam  High Field Magnet Laboratory (HFML-
EMFL), Radboud University, Toernooiveld 7, Nijmegen, 6525 ED, The Netherlands
<EMAIL_ADDRESS>
hyperfine interactions, ultra-low temperature, multipolar order,
superconductivity, nuclear quantum criticality

The motivation to develop materials for quantum technologies has put
exploration of novel quantum states of matter at the focus of several research
fields, with particular efforts towards understanding and controlling the
behaviour of quantum entangled and other strongly interacting electronic
states Osterloh  _et al._ (2002); Hanson  _et al._ (2007); Bar-Gill  _et al._
(2012); Kozii  _et al._ (2016); Maryenko  _et al._ (2018); Zaanen (2019);
Broholm  _et al._ (2020); Soluyanov  _et al._ (2015); Prochaska  _et al._
(2020); Paschen and Si (2020). Experimental investigation is of primary
importance, but requires measurements at ultra-low temperatures where the
quantum states of interest have long lifetimes Petta  _et al._ (2005); Bluhm
_et al._ (2010); Hanson  _et al._ (2007); Pickett and Enss (2018). Under these
conditions, low energy interactions, such as hyperfine or nuclear exchange
interactions, become relevant, and can modify electronic ground states and
their associated excitations in multiple ways that are not well understood or
characterised Eisenlohr and Vojta (2021); Schuberth  _et al._ (2016); Knapp
_et al._ (2023); Bitko  _et al._ (1996); Libersky  _et al._ (2021); Giménez-
Santamarina  _et al._ (2020). In this work, we use a recently developed
magnetic susceptibility technique, compatible with ultra-low temperatures and
high magnetic fields Bangma  _et al._ (2023); Bangma (2023), to probe the
influence of nuclear interactions on superconducting and multipole ordered
ground states in the strongly correlated electron system PrOs4Sb12. We find
that the multipole order develops a novel, entangled nuclear-electronic
character at the lowest temperatures, which significantly modifies the phase
boundary and leads to a nuclear quantum critical point H.M.Rønnow  _et al._
(2005); Eisenlohr and Vojta (2021). In the superconducting phase, we find that
hyperfine interactions suppress superconductivity in a manner that provides
evidence for superconducting pairing mediated by crystal field excitations
Fulde  _et al._ (1970); Goremychkin  _et al._ (2004); Matsumoto and Koga
(2004). Our results on PrOs4Sb12 experimentally establish a new type of non-
magnetic, nuclear quantum critical point, and give revealing insight into a
highly unusual superconducting state. They also demonstrate more generally the
feasibility of exploiting hyperfine interactions as a tuning parameter for
experimental creation and investigation of a variety of quantum states and
phenomena in correlated electron materials.

Strong correlations between electrons often lead to ordered phases in quantum
materials that we can describe in terms of the electronic variables. Nuclei
provide extra degrees of freedom but are typically regarded independently of
the electronic system because nuclear and electronic energy scales are well-
separated. There can be cases, however, where nuclear and electronic charge or
spin couple strongly, for example, through a large hyperfine interaction, and
we must then consider how nuclear interactions affect electronic ordered
phases and ground states.

To date, only a few examples have been studied where nuclear degrees of
freedom notably enhance, suppress, or induce electronic order. These examples,
including YbRh2Si2 Knapp  _et al._ (2023); Nguyen  _et al._ (2021); Schuberth
_et al._ (2016), LiHoF4 Bitko  _et al._ (1996), PrCu2 Andres (1973) and
Pr3Pd20Ge6 Steinke  _et al._ (2013); Iwakami  _et al._ (2014), cross several
material families and involve diverse orders, such as superconductivity or
magnetism. This shows the relevance of the effect in many contexts, but makes
it difficult to extract general trends in behaviour. Added to this is the fact
that the low energy scales involved usually necessitate experiments at ultra-
low temperatures, which can be technically challenging as well as limited in
dynamic range. Although temperatures below tens of millikelvin have been
achievable for many decades, experimental probes, particularly thermodynamic
probes, suitable for measurements of quantum solids in the ultra-low
temperature regime have not been widely available. This has blocked
experimental access and limited progress in understanding and characterising
the influence of nuclear interactions on electronic ordered states.

Figure 1: Phase diagram and effect of hyperfine interactions for AFQ phase of
PrOs4Sb12. (a) Energy of the crystal-electric-field (CEF) levels as a function
of magnetic field, in the presence (black lines) and absence (red lines) of
quadrupole and hyperfine interactions. The energy splitting of the CEF levels
in the presence of hyperfine interactions is too small to be resolved. (b)
Mean-field antiferroquadrupolar (AFQ) phase boundary excluding nuclear
interactions (solid line). The dashed lines are the two upper critical fields
of the superconducting phase in the absence of nuclear interactions, based on
Refs. Bauer  _et al._ (2002); Tayama  _et al._ (2003); Measson  _et al._
(2004). Gray data points are from Aoki ($\square$) Aoki  _et al._ (2002a),
Tayama ($\pentagon$) Tayama  _et al._ (2003), Sugawara ($\varhexagon$)
Sugawara  _et al._ (2005), and Rotundu ($\ocircle$) Rotundu  _et al._ (2004)
(c) Measured lower and (d) upper AFQ phase boundary at low temperatures, from
experiments at Royal Holloway ($\blacklozenge$), Toronto ($\star$) and
Nijmegen ($\triangleleft$). The black and light purple lines are the mean-
field antiferroquadrupolar phase boundaries excluding and including nuclear
interactions, respectively. At low temperatures, the experimental data (and
the light purple lines) show a dramatic departure from the mean-field phase
boundaries calculated without hyperfine interactions. The colour maps in (b)
to (d) show the occupation probability of the (hyperfine) ground state. This
level is fully and solely occupied at the lowest temperatures, whereas it
drops to $\sim\frac{1}{6}$ when the $T$ is a lot larger than the hyperfine
constant.

To consistently address the question of nuclear-electronic coupling in quantum
materials, we should ideally follow the evolution and influence of specific
nuclear-electronic interactions across ordered phases and phase boundaries in
a single material or closely related materials. This would require ultra-low
temperature exploration of a system with multiple electronic orders that can
be tuned via accessible parameters. We recently developed a high sensitivity
magnetic susceptibility technique, operating at temperatures down to $\sim 1$
mK in magnetic fields up to $\sim 7$ T Bangma  _et al._ (2023), which has
allowed us to do this in the strongly correlated electron material PrOs4Sb12.

PrOs4Sb12 is an unconventional superconductor with $T_{c}$=1.85 K and
$H_{c2}$=2.2 T Bauer  _et al._ (2002). It develops magnetic field induced afq
order between $\sim 4$ and 14 T Aoki  _et al._ (2002a), with low temperature
phase boundaries that terminate in qcp. (Fig. 1(b)) shows the phase diagram of
PrOs4Sb12 derived from experimental results (grey data points) and mean field
calculations (solid black lines) [see Methods]. Previous work indicated that
hyperfine interactions modify the afq order below 300 mK McCollam  _et al._
(2013), with the relatively high temperature scale being due to the
combination of high nuclear spin (I=5/2) and exceptionally large hyperfine
constant ($A=52$ mK) for Pr. The colour map superimposed on Fig. 1 shows the
occupation probability of the hyperfine ground state, such that the darkening
colour indicates the increasing importance of hyperfine interactions.
PrOs4Sb12 thus offers the opportunity to study how nuclear-electronic
interactions affect superconductivity, multipolar order and quantum
criticality.

We carried out magnetic susceptibility ($\chi$) measurements of PrOs4Sb12 in
magnetic field up to 6 T at temperatures down to $\sim 1$ mK, and in magnetic
field up to 30 T at temperatures down to $\sim 80$ mK (see Methods). Focussing
first on the AFQ phase, we find that the phase boundaries continuously develop
down to a few millikelvin, and are shifted to lower magnetic field below $\sim
300$ mK (Fig. 1 (c) and (d)). We ascribe this behaviour to the influence of
hyperfine interactions and the formation of a hybridised nuclear-electronic
ground state. The Pr energy levels are no longer purely electronic entities in
this region of the phase diagram, but must be described in terms of entangled
nuclear and electronic states. The associated quantum critical points
consequently have a fundamentally different character to their electronic
counterparts, and are examples of nuclear quantum critical points.

The superconducting phase is markedly suppressed at the lowest temperatures
(Fig. 3(c)). In this case, we show that the behaviour is due to hyperfine
enhanced magnetisation, and gives new information about the superconducting
pairing mechanism.

To understand the low temperature behaviour of PrOs4Sb12, we must consider the
Pr 4f2 electrons, which are subject to a tetragonal ($T_{h}$) cef from the
surrounding Sb atoms. The crystal field splits the ninefold degenerate, spin-
orbit coupled $J$=4 levels into a non-magnetic singlet ground state
$\Gamma_{1}$ ($E$ = 0 K) and an excited magnetic triplet $\Gamma_{4}^{(2)}$
($E$ = 8 K), with additional cef levels more than 130 K higher in energy. The
system can therefore be treated as a two-level singlet-triplet case Shiina and
Aoki (2004).

In the superconducting state, the singlet-triplet crystal field excitation has
a dominantly quadrupole character, and this quadrupole exciton is thought to
be integral to the superconductivity Fulde  _et al._ (1970); Goremychkin  _et
al._ (2004); Matsumoto and Koga (2004); Zwicknagl  _et al._ (2009).

On application of a magnetic field, the ground state acquires a weak magnetic
moment as the field mixes small amounts of $\Gamma_{4}^{(2)}$ into
$\Gamma_{1}$ Shiina and Aoki (2004). The $\Gamma_{4}^{(2)}$ triplet splits
into three branches, such that $\Gamma_{1}$ and the lowest $\Gamma_{4}^{(2)}$
level move closer together with increasing field (see Fig 1(a)). Once these
two lowest levels are sufficiently close in energy, the ground state becomes a
superposition of $\Gamma_{4}^{(2)}$ and $\Gamma_{1}$, and afq order develops.

The composition, occupation, and splitting of Pr cef levels are thus
fundamental to both the superconducting and afq states. It is precisely these
properties that are strongly affected by the hyperfine interaction in
PrOs4Sb12, which enables the Pr nuclear spin to profoundly influence the low
temperature phase diagram.

Figure 2: Hyperfine modified afq phase boundary. (a) Real part of the magnetic
susceptibility around the afq transition at various temperatures. The step in
$\chi^{\prime}(B)$ signals the afq transition. Superimposed on this are
quantum oscillations. The inset shows the mean field energy levels as a
function of magnetic field excluding (yellow) and including (blue) nuclear
interactions (b) Real part of the magnetic susceptibility with quantum
oscillations subtracted. Curves are vertically offset for clarity. The arrow
indicates the evolution of the AFQ phase boundary down to the lowest
temperatures. The inset shows the afq phase boundary at ultra-low
temperatures, saturating below 5 mK. Error bars in field are plotted, but are
smaller than the data points. The blue line (top $x$-axis) is the lower AFQ
phase boundary calculated by mean-field theory including nuclear interactions.

Modified afq phase boundary
We first focus on the afq phase transition, and in Fig. 2(a) show the real
part of the magnetic susceptibility $\chi^{\prime}(B)$ between 3.6 T and 5.4
T, at different temperatures. $\chi^{\prime}(B)$ increases rapidly near 4 T,
signalling the lower afq phase transition, and shows de Haas-van Alphen
oscillations at higher fields. Fig. 2(b) shows a series of $\chi^{\prime}(B)$
curves with the quantum oscillations subtracted. We extracted the phase
boundary $H_{A}$ by differentiating $\chi^{\prime}(B)$ and taking the field
value corresponding to the peak of the derivative, which could be done more
accurately when the modulations from the dHvA effect were removed. $H_{A}$ vs.
temperature is plotted as red data points in Fig. 1(c). A different definition
of $H_{A}$, for example, at the foot or shoulder of the transition, shifts the
phase boundary but does not change its shape (see Methods and Extended Data
Fig1).

Remarkably, $H_{A}$ develops strongly down to 5 mK, where it saturates within
the experimental uncertainty, as shown in the inset to Fig. 2(b). Our ultra-
low temperature measurements could not be extended to the upper phase boundary
near 14 T, but Fig 1(d) shows that we observe a similar continued shift of the
upper phase boundary to at least the lowest measured temperature of 80 mK.
This behaviour is at odds with mean field models, which predict saturation of
$H_{A}$ below $\sim 400$ mK Shiina (2004); Kusunose  _et al._ (2009). A mean
field AFQ phase boundary (based on Bangma (2023); Shiina and Aoki (2004);
Shiina (2004); Kusunose  _et al._ (2009); Aoki  _et al._ (2002a); Tayama  _et
al._ (2003); Sugawara  _et al._ (2005); Rotundu  _et al._ (2004)) is shown as
a solid black line in Fig 1(b,c,d).

As discussed above, we expect nuclear degrees of freedom to be important below
a few hundred millikelvin in PrOs4Sb12. To explore the possibility that this
is responsible for the anomalous development of the AFQ phase boundary, we
incorporated hyperfine interactions into the mean field model of the AFQ phase
(see Methods).

Each cef level is a manifold of (2$I$+1)=6 nuclear levels. The hyperfine
interaction lifts their degeneracy, and causes an energy splitting
proportional to their magnetic moment. This is illustrated in the inset of
Fig. 2(a)), which shows a significantly larger splitting in the magnetic
$\Gamma_{4}^{(2)}$ state than in the barely magnetic ground state, which
consists primarily of $\Gamma_{1}$. As a result, the energy gap from
$\Gamma_{1}$ to $\Gamma_{4}^{(2)}$ diminishes as the system is cooled and
occupation of the lowest hyperfine level increases. Equivalently, we can think
of this as the hyperfine Hamiltonian mixing an increasing amount of
$\Gamma_{4}^{(2)}$ into $\Gamma_{1}$ as temperature is reduced. The AFQ order
is thus stabilised at increasingly lower magnetic field, an effect which
saturates when the lowest hyperfine level is fully occupied.

The light purple line in Fig. 1(c and d) represents the mean field AFQ phase
boundary when the nuclear interactions are included. It is in good qualitative
agreement with the experimental $H_{A}$ (red data points), as it captures the
strong shift to lower field as the temperature is reduced. The saturation of
the lower phase boundary below 5 mK is also well-described by the model (see
inset of Fig. 2(b)), indicating that the system enters its quantum ground
state below 5 mK.

The overall shape of the mean field AFQ phase boundary is not a good
quantitative match to the experimental data on the scale of a few hundred
millikelvin. This can be clearly seen in Fig.s 1(b) and (c): on the 2 K and 15
T scale of (b), the data points appear to match well to the mean-field AFQ
phase boundary, however, systematic deviations on the scale of 100 mK/0.1 T
can be clearly seen in Fig. 1(c). (See Methods for further discussion.) Even
by adjusting the parameters in the mean field model, we are not able to
accurately match the calculated phase boundary to the data at all fields and
temperatures. This is an important finding, because it shows that although a
mean field picture is extremely useful for broadly mapping and understanding
the phase diagram of PrOs4Sb12, it does not give a quantitative description at
the energy scale relevant for the hyperfine interactions and lowest
temperature behaviour.

Nuclear quantum criticality
We have shown that the AFQ order in PrOs4Sb12 is strongly temperature
dependent in the composition of the states that order, as well as in its phase
boundary. At high temperatures, the quadrupoles are formed on the (electronic)
Pr crystal field levels, but at temperatures below $\sim$ 300 mK they develop
on entangled crystal field and nuclear states. The low temperature AFQ phase
is thus an example of hybrid nuclear-electronic order, and the corresponding
quantum critical points are of a novel type, associated with transitions
between disordered and ordered nuclear-electronic phases.

At a continuous quantum phase transition, the energy of collective electronic
excitations is expected to go to zero at the quantum critical point Sachdev
(2011). In systems with significant hyperfine coupling, electronic and nuclear
modes mix, such that there are two sets of dispersive modes, one dominantly
electronic, and one at lower energy that is dominantly nuclear. At the quantum
critical point, softening of the dominantly electronic mode is incomplete and
it is the nuclear mode that becomes gapless. This behaviour leads to nuclear
quantum criticality, which cuts off electronic criticality below a certain
temperature that depends on the hyperfine coupling strength Eisenlohr and
Vojta (2021); Libersky  _et al._ (2021); McKenzie and Stamp (2018); Kovacevic
_et al._ (2016); H.M.Rønnow  _et al._ (2005)

The AFQ quantum critical points in PrOs4Sb12 are nuclear quantum critical
points, and provide a rare opportunity to explore and characterise nuclear
quantum criticality. The only previous extended investigation of nuclear
quantum criticality has been in the transverse-field Ising magnet LiHoF4
H.M.Rønnow  _et al._ (2005, 2007); Kovacevic  _et al._ (2016); McKenzie and
Stamp (2018), although the effect of nuclear spins on a number of other model
magnetic quantum phase transitions was discussed in recent theoretical work
Eisenlohr and Vojta (2021). A non-magnetic nuclear quantum critical point,
such as the AFQ nuclear quantum critical point we observe in PrOs4Sb12, has
not previously been studied either experimentally or theoretically.

In LiHoF4, the transfer of quantum criticality from the electronic to the
nuclear mode is accompanied by a marked shift of the quantum critical point
Bitko  _et al._ (1996), due to stabilisation of the magnetic phase by
interaction with the nuclear spins. The modified AFQ phase boundary we observe
in PrOs4Sb12 suggests that the general behaviour on approaching the nuclear
QCP is similar to that in LiHoF4. However, the different role of nuclear spins
in a quadrupolar system is likley to be important for a quantitative
understanding. Our results on PrOs4Sb12 strongly motivate theoretical work and
further experiments to examine nuclear quantum criticality associated with
multipole and other non-magnetic order parameters. It is of particular
interest to examine the crossover with decreasing temperature between the
regimes of electronic and nuclear quantum criticality. The behaviour near
nuclear quantum criticality in LiHoF4 has been discussed in terms of
interaction of the Ising system with a spin bath, and an interesting region of
minimum decoherence of the 2-level spin system was proposed in this crossover
region, just at the temperature where the ferromagnetic phase boundary begins
to be modfied by the hyperfine interactions H.M.Rønnow  _et al._ (2005).

Figure 3: Reduction of the superconducting upper critical field by hyperfine
interactions. (a) Imaginary and (b) real parts of the magnetic susceptibility
at various temperatures. The peak and step respectively signal the
superconducting transition. The curves are vertically offset for clarity. (c)
Comparison between the extracted $H_{c2}$(T) from (a,b) ($\blacklozenge$) and
calculated $H_{c2}$ (purple line) based on the model described in the text.
The dashed, black line shows $H_{c2}$ of the model in the absence of nuclear
interactions. The right inset shows the additional magnetization at very low
temperature caused by nuclear interactions, as calculated by mean field
theory. The left inset shows $H_{c2}$ at the lowest temperatures.

Hyperfine-induced suppression of SC
We turn now to the superconducting region of the phase diagram. Figure 3 shows
the imaginary (a) and real (b) magnetic susceptibility crossing the
superconducting transition down to 5 mK. The sharp peak in
$\chi^{\prime\prime}$ and step in $\chi^{\prime}$ signal the upper critical
field, $H_{c2}$. The slight shoulder visible in $\chi^{\prime\prime}$ at
higher field is the second superconducting transition that is typically
observed in PrOs4Sb12 Bauer  _et al._ (2002); Méasson  _et al._ (2008);
Andraka and Pocsy (2012). We focus on the temperature dependence of $H_{c2}$
below 400 mK, shown in Fig. 3(c). Interestingly, we observe a strong decrease
of $H_{c2}$ below 100 mK, until it eventually saturates below 2 mK.

The mechanism of superconductivity in PrOs4Sb12 is not fully established. The
sister compound LaOs4Sb12 is understood to be a multiband superconductor, with
conventional, phonon-mediated pairing on two bands, and $T_{c}=740$ mK Tee
_et al._ (2012). The significantly higher $T_{c}$ in PrOs4Sb12 is therefore
clearly linked to the presence of Pr $4f$ electrons, and is thought to involve
an unconventional pairing channel on one of the bands Seyfarth  _et al._
(2005). A number of scenarios have been proposed. One scenario is that a
strong interaction of the conduction electrons with the Pr $4f$ charge density
through inelastic charge scattering promotes Cooper pair formation. This
process involves exciting quadrupolar excitons (dispersive excitations of the
cef levels) Fulde  _et al._ (1970); Sugawara  _et al._ (2002); Goremychkin
_et al._ (2004); Koga  _et al._ (2006); Chang  _et al._ (2007); Namiki  _et
al._ (2007); Sato  _et al._ (2010), which could play a similar role to phonons
in conventional BCS theory. The inelastic charge scattering competes with, but
can dominate the pair-breaking 4f–conduction electron exchange interaction.
The absence of a magnetic moment in $\Gamma_{1}$ and the quadrupolar character
of the exciton in PrOs4Sb12 provides optimal conditions to realize this.

Another proposal is superconducting pairing mediated by quadrupole
fluctuations associated with the nearby afq order and quantum critical point,
similar to spin-fluctuation mediated superconductivity in the vicinity of a
magnetic quantum critical point Cox and Maple (1995); Mathur  _et al._ (1998);
Chia  _et al._ (2003); Aoki  _et al._ (2003, 2007). However, increased
quadrupolar susceptibility is reported to be associated with suppression of
superconductivity in PrOs4Sb12, rather than enhancement Tayama  _et al._
(2006), which is contrary to the scenario of quadrupole fluctuation mediated
pairing. We observed that the afq quantum critical point moves to lower
magnetic field as a result of interaction with the nuclear states, which
increases the quadrupole fluctuations at low temperature in the region of
superconductivity Eisenlohr and Vojta (2021). Calculation of the mean field
quadrupolar susceptibility at 2 T shows that it rises steeply below $\sim 100$
mK [see Methods and Extended Data Fig. 3), but we observe suppression of the
superconductivity in this temperature range. This suggests that quadrupole
fluctutations directly compete with superconductivity.

The similar temperature scale of the modified superconducting and AFQ phase
boundaries in PrOs4Sb12 points to a similar origin. We therefore look further
into the role of nuclear interactions in the superconducting state.

Independently of the superconducting or AFQ orders, the hyperfine interaction
leads to a strongly increasing magnetic moment on the Pr atoms below $\sim$200
mK, due to the mixing of $\Gamma_{4}^{(2)}$ and $\Gamma_{1}$ described above.
The change in (mean field) magnetisation $\Delta M$ per Pr atom, as a result
of this mixing, is shown as a function of temperature in the inset of Fig.
3(c)). The conduction electrons interact with the 4f magnetic moment through
the exchange interaction, which tends to break Cooper pairs. In this way, the
hyperfine interaction leads to suppression of superconductivity through
enhanced exchange scattering.

We model the effect of this exchange scattering by using the generalized
Abrikosov-Gorkov theory of Fulde and Maki Fulde and Maki (1966). This theory
describes how magnetic impurities reduce $H_{c2}$ by creating an internal
field and enhancing exchange scattering. The internal field set-up by the Pr
atoms is negligible, and the role of the magnetic impurities in our case is
played by the hyperfine-enhanced magnetic moment, i.e. $\Delta M$. We use the
mean-field $\Delta M$ to capture the amount of exchange scattering:

$\begin{split}&H_{c2}(T,A=\textrm{52\leavevmode\nobreak\ mK})=\\\
&H_{c2}(T,A=\textrm{0})-H_{c2}(T=0,A=0)\frac{\Delta
M(T)}{M_{\textrm{crit}}}\end{split}$ (1)

$H_{c2}(T,A=0)$ is the upper critical field in the absence of nuclear
interactions, for which we used a curve based on high-temperature data points
Maki  _et al._ (2007), above the temperatures relevant for hyperfine
interactions. $M_{\textrm{crit}}=0.83\leavevmode\nobreak\ \mu_{B}$/Pr was set
to ensure a match with our lowest temperature data point (see Methods for
details).

$H_{c2}(T)$ calculated from this model is shown as a purple line in Fig. 3(c).
The lower inset shows an expanded view of the lowest temperature region. The
calculated curve adequately describes the experimental data (red points): The
maximum $H_{c2}(T)$ of the model is in very close agreement with the
experimentally observed maximum at 100 mK, and the saturating $M$ of the Pr
atom below 2 mK explains the observed saturation of $H_{c2}$ below 2 mK. Thus,
we conclude that superconductivity in PrOs4Sb12 is suppressed by the
hyperfine-enhanced magnetisation of Pr as it enhances pair-breaking exchange
scattering.

In addition to enhanced quadrupole fluctuations and enhanced exchange
scattering, we have also considered the effect of the Pauli limiting field and
orbital limiting field on $H_{c2}$. We find that Pauli limiting would lead to
too small a change in $H_{c2}$ to explain the observed behaviour, and orbital
limiting would increase, rather than decrease, $H_{c2}$ (see Methods and
Extended Data Fig. 3).

Limited exchange scattering between 4$f$ and conduction electrons is a
prerequisite for superconductivity mediated by quadrupolar excitons Fulde  _et
al._ (1970), where the balance between exchange scattering and inelastic
charge scattering determines the superconducting $T_{c}$. The observed
sensitivity of superconductivity to even weakly enhanced exchange scattering
indicates the importance of this balance in PrOs4Sb12. Our ultra-low
temperature results therefore strongly suggest that Cooper pairing via
inelastic charge scattering and the formation of quadrupolar excitons
contributes significantly to superconductivity in PrOs4Sb12. This quadrupolar
exciton mechanism should be in addition to phonon-mediated pairing, comparable
to that observed in LaOs4Sb12, and can account for the increased $T_{c}$ in
the Pr compound.

Note that in the above we have focused on the first (lower field)
superconducting transition because it is the dominant feature in the
susceptibility. However, as can be seen in Fig. 3(a), the higher field
shoulder in $\chi^{\prime\prime}$, representing the second superconducting
transition, shows the same temperature dependence as the first transition, and
the two $H_{c2}$ phase boundaries have the same shape. This implies that both
superconducting transitions are affected in the same way by hyperfine
interactions. The second superconducting transition in PrOs4Sb12 is robust in
samples above a certain size Andraka and Pocsy (2012), and the double
transition has been proposed as evidence of non-unitary chiral
superconductivity in this material Kozii  _et al._ (2016). If the two
transitions delineate distinct superconducting phases Kozii  _et al._ (2016),
our results suggest that the superconducting mechanisms in the two phases are
similar.

Summary and concluding remarks
Development of a magnetic susceptibility technique compatible with ultra-low
temperature Bangma  _et al._ (2023); Bangma (2023) has allowed us to
sensitively probe the ground state of PrOs4Sb12 and examine the effect of
nuclear degrees of freedom on electronic ordered phases. We find that
hyperfine interactions have a strong effect on both the AFQ and
superconducting states, but that the consequences for the electronic order are
quite different in each case.

In the AFQ phase, the mixing of electronic and nuclear modes by the hyperfine
interaction leads to a novel, hybrid nuclear-electronic AFQ order below $\sim
300$ mK, with a strongly modified phase boundary and a new type of nuclear AFQ
quantum critical point. Our experiments have provided the first access to
nuclear quantum criticality in the ultra-low temperature regime, where we have
been able to measure the saturation of the nuclear-electronic system into the
hyperfine ground state below 5 mK. PrOs4Sb12 is thus an ideal system in which
to further explore the development of nuclear quantum criticality from
electronic quantum criticality, and may be of particular interest in the
contexts of coherent nuclear-electronic entangled states H.M.Rønnow  _et al._
(2005); Stamp and Tupitsyn (2004); Fraval  _et al._ (2004), and non-Fermi
liquid behaviour associated with AFQ quantum criticality Patri and Kim (2020);
Han  _et al._ (2022).

In the superconducting phase, we find that $H_{c2}$ is suppressed at
temperatures below $\sim 200$ mK, in a way that is well-described by a model
of hyperfine-enhanced exchange scattering. The extreme sensitivity of the the
superconducting state to exchange scattering in this material suggests a
quadrupole exciton mediated mechanism of superconducting pairing Fulde  _et
al._ (1970); Goremychkin  _et al._ (2004). This new information about the
likely pairing mechanism in PrOs4Sb12 goes some way to resolving long-standing
questions about the nature of superconductivity in this material. Further work
to understand if two-band superconductivity partly mediated by quadrupole
excitons is compatible with proposed non-unitary chiral superconductivity
Kozii  _et al._ (2016) would be valuable.

The results we present in this article demonstrate both the feasibility and
the value of thermodynamic measurements on quantum solids at ultra-low
temperatures and high magnetic fields. We observe the strong influence of
hyperfine interactions on low temperature electronic orders, and gain striking
new information about the nature of quantum criticality, the AFQ ground state
and superconductivity in PrOs4Sb12. More generally, this work provides clear
evidence that, as we search for materials to support quantum technologies,
experiments of this kind can contribute key information about the role of
nuclear degrees of freedom in phenomena such as quantum order and quantum
criticality, which shape the low temperature behaviour and stability of
potentially important quantum states Pickett and Enss (2018).

Acknowledgements
We thank Lijnis Nelemans, Michel Peters, Richard Elsom, Paul Bamford and Ian
Higgs for excellent mechanical support. The research leading to these results
has received funding from the European Union’s Horizon 2020 Research and
Innovation Programme, under Grant Agreement no 824109 (European Microkelvin
Platform). Part of this work was also supported by HFML-RU/NWO-I, a member of
the European Magnetic Field Laboratory (EMFL).

## Methods

Sample growth
Single crystals of PrOs4Sb12 were grown by a standard Sb-self-flux method
Bauer  _et al._ (2001). The sample measured at the High Field Magnet
Laboratory (HFML) and at Royal Holloway University of London (RHUL) had
dimensions 0.9 $\times$ 0.95 $\times$ 1.05 mm3, and the sample measured at the
University of Toronto had dimensions 2.5 $\times$ 2 $\times$ 2 mm3. The
samples had similar residual resistance ratios between 70 and 80.

Experiments
Magnetic susceptibility was measured using the field-modulation technique,
with the sample placed in one of two counterwound pick-up coils, mounted in
the center of a modulation coil and a larger sample magnet. The lowest-
temperature experiments were carried out at the Department of Physics, RHUL.
Measurements were performed in a nuclear demagnetization refrigerator with a
Cu nuclear stage and 9 T superconducting sample magnet. The data shown in the
main paper were acquired in a modulation field of 1.4 G at a frequency of
92.34 Hz. To investigate possible sample heating due to the modulation field,
measurements were repeated at reduced modulation amplitudes of 0.03 G, 0.14 G
and 0.55 G at the lowest temperatures (see Sample temperature below). The
voltage across the pick-up coils was amplified by a temperature-stabilized and
lead-shielded low-temperature transformer, and measured with a lock-in
amplifier. More detailed information about the experimental set-up is given
elsewhere Bangma  _et al._ (2023). Data were obtained with the sample
stabilized at a constant temperature, while ramping the sample magnet from
high to low field at $-0.05$ mT/s. To prevent hysteretic effects, the system
was always stabilized at the required temperature before entering the afq
phase by magnetisation across $H_{A}$.

Magnetic susceptibility measurements at the HFML and University of Toronto
were carried out using regular dilution refrigerators to base temperatures of
88 mK and 30 mK, respectively. The field-modulation technique was used in both
cases.

Determination of the AFQ phase boundary
We extracted the afq phase boundary $H_{A}$ by differentiating
$\chi^{\prime}(B)$ and taking the field value corresponding to the peak of the
derivative, which could be done more accurately when the modulations from the
de Haas–van Alphen effect were removed. A different definition of $H_{A}$
shifts the phase boundary but does not qualitatively change its shape. For
example, if we take the foot or the shoulder of the transition, defined to be
at the low- or high-field half-heights of Gaussian fits to
$\text{d}\chi^{\prime}/\text{d}B$, we obtain the phase boundary shown in
Extended Data Fig. 1.

Sample temperature
Noise thermometry Casey  _et al._ (2014) and melting-curve thermometry
Greywall (1985) were used during the lowest temperature measurements (RHUL).
Both thermometers were located outside of the high field region, at a distance
of $\sim$ 70 cm from the experiment. To obtain an estimate of the thermal
gradients in the system, an extensive thermal model was developed and tested,
details are given in ref. Bangma  _et al._ (2023). The main source of heat
close to the sample is eddy current heating due to the modulation field.
Repeated measurements of the afq phase boundary at ultra-low temperatures show
no difference when performed at 0.6 G or 1.4 G, and thus give no indication of
heating by eddy currents in the vicinity of the afq phase. In addition, the
saturation of the afq phase boundary below 5 mK is in good agreement with
mean-field models, which predict saturation below 7 mK. We therefore conclude
that there is no significant eddy current heating of the sample at the
modulation fields used.

Sample temperature in the superconducting phase and determination of ${\bf
H_{c2}}$
As discussed above, no heating of the sample was observed in the afq phase.
When the sample is in the superconducting phase, however, it is subject to
additional heating due to the motion of flux vortices. The vortex motion is
driven by both the modulation field and by sweeping the DC field, and is a
function of the penetration field ($B_{p}$) Bean (1964). $B_{p}$ scales
linearly with the critical current density $j_{c}(B)$. We therefore expect
this heating to be insignificant close to $H_{c2}$, where $j_{c}$ is small,
and to grow as we enter further into the superconducting phase. Here we show
that this heating is negligible very close to the normal state.

The superconducting phase transition is recognized in the AC susceptibility as
a step in $\chi^{\prime}$ and a peak in $\chi^{\prime\prime}$, see Extended
Data Fig. 2(a), (b). We define $H_{c2}$ as some fraction of the step down in
$\chi^{\prime}(B)$. Choosing different fractions leads to small apparent
shifts in $H_{c2}$ (see Extended Data Fig.2). Smaller apparent values of
$H_{c2}$ (larger fractions) are further into the superconducting state and
imply larger values of $j_{c}$ and associated vortex heat.

Based on mean-field calculations, we expect to fully occupy the quantum ground
state of PrOs4Sb12 at 2 T in the low-mK temperature range (see discussion and
Fig. 3(c) in the main paper). In this regime, saturation of $H_{c2}$ should
occur. At higher temperatures, changing occupations of the hyperfine levels
will result in temperature dependence of $H_{c2}$ (see discussion in main
paper). In Extended Data Fig. 2(c)-(f) we compare the temperature dependence
of $H_{c2}$ between 0.5 mK and 5 mK obtained at $b_{mod}$ = 1.4 G. For
$H_{c2}$ defined at fractions of 0.01 and 0.03 of the step down in
$\chi^{\prime}(B)$ (Extended Data Fig. 2(c) - (d)), $H_{c2}$ indeed shows
saturation below 3 mK. This behaviour is also observed in measurements
obtained at lower modulation fields of 0.03 G and 0.14 G (Extended data Fig.
2(g)-(n)), and indicates that vortex motion does not heat the sample.
Definitions of $H_{c2}$ with significantly higher fractions, e.g. 0.5 and 0.9
in Extended Data Fig. 2(e) - (f), which probe the sample further into the
superconducting state, show a continuously evolving $H_{c2}$ down to the
lowest temperatures, indicating heating We therefore defined the upper
critical field $H_{c2}$, as shown in Fig. 3, as 1% of the step down in
$\chi^{\prime}(B)$, where we know the sample is not subject to large heating,
see Extended Data Fig. 2(b).

Understanding the superconducting upper critical field boundary in PrOs4Sb12
To understand the anomalous behaviour of the upper critical field boundary in
PrOs4Sb12, we use a model derived from the generalized Abrikosov-Gorkov theory
of Fulde and Maki Fulde and Maki (1966), which we modify to make applicable to
PrOs4Sb12. The original model describes how magnetic impurities reduce
$H_{c2}$ via (1) spin-flip scattering of conduction electrons by magnetic
impurities and (2) spin polarization of conduction electrons by the exchange
field generated by magnetic impurities:

$\begin{split}&H_{c2}(x,T)=\\\
&H_{c2}(0,T)-H_{c2}(0,0)\Big{(}\frac{x}{x_{crit}}+\frac{P}{P_{crit}}\Big{)}\end{split}$
(1)

in which $x$ is the concentration of magnetic impurities, $H_{c2}(x,T)$ is the
upper critical field at temperature $T$, and $P$ is the Pauli polarization
term. Fulde and Maki (1966)
In the case of PrOs4Sb12, we assume that the nuclear spins act as magnetic
impurities. As limited exchange scattering is considered to be key for
superconductivity in PrOs4Sb12, we limit ourselves to pair-breaking mechanism
(1) above. In adapting the Fulde-Maki model, we convert the concentration of
magnetic impurities $x$ to an increase in magnetization $\Delta M$, caused by
the hyperfine interaction of strength $A$. Our expression for the upper
critical field is then:

$\begin{split}&H_{c2}(A=52\leavevmode\nobreak\ \textrm{mK},T)=\\\
&H_{c2}(0,T)-H_{c2}(0,0)\Big{(}\frac{\Delta M(T)}{M_{crit}}\Big{)}\end{split}$
(2)

in which $\Delta M(T)=M(A=52\leavevmode\nobreak\ \textrm{mK},T)-M(A=0,T)$ is
the difference in mean-field calculated magnetizations including ($A=52$ mK)
and excluding ($A=0$) the hyperfine interaction. $\Delta M$ is expressed in
$\mu_{B}/\textrm{Pr atom}$ by comparison to measured magnetization Tayama  _et
al._ (2003). $H_{c2}(A=0,T)$ is taken from Maki et al. Maki  _et al._ (2007),
based on measurements of PrOs4Sb12 at high temperature, which are largely
independent of effects related to Pr nuclei. The $H_{c2}$ values in Maki et
al. were larger than the values we measured. This could be due to differences
in sample quality or because of a different definition of $H_{c2}$. We
therefore rescaled $H_{c2}(A=0,T)$ from Maki  _et al._ (2007) to match our
high temperature data points above 300 mK. This results in a value of
$H_{c2}(A=0,T=0)=2.245$ T. $M_{\textrm{crit}}=0.68\leavevmode\nobreak\
\mu_{B}$/Pr was set to ensure a match with our lowest temperature data point
at 1 mK.

Enhanced exchange scattering is not the only mechanism that could suppress
superconductivity. We also considered how hyperfine effects on Pauli limiting
or orbital limiting could modify $H_{c2}$.

In the case of Pauli limiting, we convert $\Delta M$, determined above, to an
internal field. When this is added to the external applied field, we find that
the Pauli-limiting field changes by $\sim$ 1 mT. This change is orders of
magnitude too small to account for the observed suppression of $H_{c2}$.

We calculate the effect on the orbital-limiting field by considering the
temperature dependence of the effective mass $m_{e}$ of the quasiparticles:
$H_{\text{orb}}\propto(1/\xi^{2})\propto m_{e}^{2}$. A suggested mass-
enhancement mechanism of quasiparticles in PrOs4Sb12 is through the aspherical
coulomb scattering, captured in terms of the quadrupolar susceptibility
Zwicknagl  _et al._ (2009):

$\frac{m_{e}}{m_{0}}=1+\alpha\chi_{Q}$ (3)

where $m_{0}$ is the free electron mass. $\alpha$ contains various pre-factors
and is arbitrarily set to 1. We can thus model the qualitative effect of mass-
enhancement on the orbital-limiting field using:

$\begin{split}&H_{\text{orb}}(A=52\ \textrm{mK},T)=\\\
&H_{\text{orb}}(0,T)\left(\frac{1+\alpha\chi_{Q}(A=52\
\textrm{mK})}{1+\alpha\chi_{Q}(A=0\ \textrm{mK})}\right)\end{split}$ (4)

The quadrupolar susceptibility $\chi_{Q}$ increases as temperature is
decreased, see inset of Extended Data Fig. 3. This would enhance
$H_{\text{orb}}$, and thus cannot account for the observed suppression of
superconductivity.

Enhanced spin susceptibilities in the vicinity of magnetic phase transitions
can stabilize superconductivity Monthoux  _et al._ (2007). We examine whether
the enhanced quadrupolar susceptibility $\chi_{Q}$ can similarly affect
superconductivity in PrOs4Sb12. As $Q_{zx}$ is the ordering quadrupole moment
for $B\parallel$[100] (which is the direction of the magnetic field in our
experiments), we consider $\chi_{Q_{zx}}$, given by Zwicknagl  _et al._
(2009):

$\begin{split}&\chi_{Q_{zx}}=\\\
&\sum_{m=0}^{m=6}\sum_{n=6}^{n=24}\frac{2\lvert\exp{(-E_{m}/T)}\langle\lambda_{n}\rvert\hat{O}_{zx}\lvert\lambda_{m}\rangle\rvert^{2}}{E_{n}-E_{m}}\\\
&\times\sum_{m=0}^{m=6}\frac{1}{\exp{(-E_{m}/T)}}\end{split}$ (5)

in which $E_{i}$ is the mean-field energy of the hyperfine levels, and the
sums are over all six hyperfine ground states $m$ and all 18 hyperfine excited
states $n$. Following Monthoux et al., this susceptibility could set the size
of the superconducting pairing interaction strength Monthoux  _et al._ (2007).
We assume $\chi_{Q}$ to be proportional to $H_{c2}$, and simulate the effect
of hyperfine interactions on $H_{c2}$ through this quadrupole susceptibility
using:

$\begin{split}&H_{\chi_{Q}}(A=52\textrm{mK},T)=\\\
&H_{c2}(0,T)\left(\frac{\chi_{Q}(A=52\ \textrm{mK})}{\chi_{Q}(A=0\
\textrm{mK})}\right)\end{split}$ (6)

As hyperfine interactions cause $\chi_{Q_{zx}}$ to increase with decreasing
temperature (see inset of Extended Data Fig. 3), we expect $H_{c2}$ to
increase with decreasing temperature. This is in conflict with experimental
results, and thus cannot explain the observed low-temperature suppression of
$H_{c2}$.

We summarize the effects of hyperfine-modified Pauli and orbital limiting and
quadrupolar susceptibility on $H_{c2}$ in Extended Data Fig. 3.

Mean-field theory
To determine the mean-field behaviour of PrOs4Sb12, we have closely followed
the work by Shiina Shiina (2004), Shiina and Aoki Shiina and Aoki (2004),
Kohgi et al. Kohgi  _et al._ (2003) and Kusunose et al. Kusunose  _et al._
(2009). In agreement with their work, dipole-dipole interactions are not
included in our model because the (crystal field) excitons in PrOs4Sb12 are
quadrupolar in character Kuwahara  _et al._ (2005); Kaneko  _et al._ (2007)
and including significant spin-spin interactions does not give an accurate
description of the behaviour of the system Kaneko  _et al._ (2007). Similar to
McCollam et al. McCollam  _et al._ (2013), we have added nuclear interactions
to the mean-field Hamiltonian and extended the basis to include the Pr nuclear
spin.

Our mean-field Hamiltonian consists of the crystal electric field (CEF)
Hamiltonian, mean-field quadrupole-quadrupole interaction, electronic Zeeman,
and nuclear interactions. The terms in the nuclear Hamiltonians are given
below:

$\begin{split}\hat{\mathcal{H}}_{\textrm{HF}}&=A\hat{\vec{I}}\cdot\hat{\vec{J}}=A\left(\hat{I}_{z}\hat{J}_{z}+\frac{1}{2}(\hat{I}_{+}\hat{J}_{-}+\hat{I}_{-}\hat{J}_{+})\right)\\\
\hat{\mathcal{H}}_{\textrm{NQI}}&=P\hat{I}_{z}^{2}-\frac{P}{3}\hat{I}^{2}\\\
\hat{\mathcal{H}}_{\textrm{Z,
nuclear}}&=\frac{g_{N}\mu_{N}}{k_{B}}\hat{\vec{I}}\cdot\vec{B}\end{split}$ (7)

in which $\hat{\mathcal{H}}_{\textrm{HF}}$ is the hyperfine interaction with
hyperfine constant $A=52$ mK Kondo (1961); Bleaney (1963);
$\hat{\mathcal{H}}_{\textrm{NQI}}$ is the nuclear quadrupole interaction, $P$
= 7 mK Aoki  _et al._ (2002b); and $\hat{\mathcal{H}}_{\textrm{Z, nuclear}}$
is the nuclear Zeeman interaction. With the Pr nuclear g-factor $g_{N}$ =
1.7102 Macfarlane  _et al._ (1982), the prefactor amounts to only
$(g_{N}\mu_{N})/k_{B}$=0.6 mK/T. Thus, nuclear interactions in magnetic field
are dominated by the hyperfine interaction. For an overview of the basis
vectors, Hamiltonians, and specific values of parameters used in our model, we
refer to Bangma (2023). In all calculations and measurements discussed in this
paper, the magnetic field was oriented along the $x$-axis, $B\parallel[100]$.
Experiments with $B\parallel[110]$ also show the relevance of hyperfine
effects McCollam  _et al._ (2013). Indeed, calculations with magnetic field
along [110] or [001] yield qualitatively similar results to those described
here, but the precise low-temperature shape of the afq phase boundary changes
slightly for different field directions. This is because the composition of
the energy levels depends on magnetic field orientation, due to absence of
four-fold symmetry on the Pr atom in PrOs4Sb12.

The mean-field quadrupole-quadrupole interaction constant was set to 240 mK,
to give a match between the lower afq phase boundaries from our experimental
data and the mean-field model (including hyperfine and nuclear terms) at 100
mK (see Fig. 1(c)). However, a match at 100 mK leads to disagreement between
the experimental and mean-field phase boundaries at other temperatures. Most
significantly, our mean-field model overestimates the maximum critical
temperature of the afq phase by a factor of 1.9 K/1.2 K. A similar discrepancy
is reported by other authors for PrOs4Sb12 Kohgi  _et al._ (2003); Shiina
(2004); Kusunose  _et al._ (2009), and is also observed in mean-field models
of other materials Tallon  _et al._ (2011); Dollberg  _et al._ (2022). We were
unable to find a value of the quadrupole-quadrupole interaction strength that
resulted in a good match of experimental and mean-field phase boundaries at
all temperatures: a weaker (stronger) quadrupole-quadrupole interaction
reduces (increases) the maximum critical temperature, but results in a strong
mismatch of the low-(high-)temperature phase boundaries. This mismatch is
likely to be because the high temperature AFQ phase boundary is strongly
suppressed by fluctuations, which are not taken into account in the mean-field
theory.

## Extended data

Extended Data Fig. 1:  Extraction of lower afq phase boundary. (a)
$\chi^{\prime}$ at 1 mK from which quantum oscillations have been subtracted.
(b) Field-derivative of $\chi^{\prime}$ from (a) with a fit to a Gaussian
function (red dashed line). The diamond indicates the center of the Gaussian
fit, and the downward and upward triangles represent the half-height points of
the same Gaussian at low and high field, respectively, and are also shown in
(a). The triangles are close to the foot and shoulder of the transition. (c)
Temperature dependence of the afq phase boundary based on the center
(diamonds) or low/high field half-heights (downwards/upwards triangles) of a
Gaussian fit. The different definitions shift the phase boundary, but do not
qualitatively change its shape. Extended Data Fig. 2:  Superconducting
transition at various modulation fields for different definitions of
$\bf{H_{c2}}$. Real (a) and imaginary (b) magnetic susceptibility at 70 mK.
Data are normalized to the amplitude of the modulation field. A smaller
modulation field moves the step in $\chi\prime$ and peak in $\chi\prime\prime$
to higher magnetic field (e.g. Ge  _et al._ (2014); Müller (1989)). The
resulting different upper critical fields are indicated by red diamonds and
are defined at 1% of the step down in $\chi^{\prime}$ (grey, dashed lines for
1.4 G data). (c) - (n): $H_{c2}(T)$ for different modulation fields of 1.4 G
((c) - (f)); 0.14 G ((g) - (j)); and 0.03 G ((k) - (n)), defined at fractions
of 1%, 3%, 50% and 90% of the step down in $\chi^{\prime}$. See text.
Extended Data Fig. 3: Expected changes in the upper critical field from
various mechanisms due to hyperfine interactions. $H_{\textrm{orb}}$ and
$H_{{\chi}_{Q}}$ are estimates, with shaded areas added to give an idea of the
uncertainty in the estimate. The inset shows the quadrupolar susceptibility vs
temperature, expression 5 Zwicknagl  _et al._ (2009).

## References

  * Osterloh  _et al._ (2002) A. Osterloh, L. Amico, G. Falci, and R. Fazio, Nature 416, 608 (2002).
  * Hanson  _et al._ (2007) R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen, Reviews of Modern Physics 79, 1217 (2007).
  * Bar-Gill  _et al._ (2012) N. Bar-Gill, L. Pham, C. Belthangady, D. L. Sage, P. Cappellaro, J. Maze, M. Lukin, A. Yacoby, and R. Walsworth, Nature Communications 3 (2012), 10.1038/ncomms1856.
  * Kozii  _et al._ (2016) V. Kozii, J. Venderbos, and L. Fu, Science Advances 2 (2016), 10.1126/sciadv.1601835.
  * Maryenko  _et al._ (2018) D. Maryenko, A. McCollam, J. Falson, Y. Kozuka, J. Bruin, U. Zeitler, and M. Kawasaki, Nature Communications 9 (2018), 10.1038/s41467-018-06834-6.
  * Zaanen (2019) J. Zaanen, SciPost Physics 6 (2019), 10.21468/scipostphys.6.5.061.
  * Broholm  _et al._ (2020) C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Science 367, 263 (2020).
  * Soluyanov  _et al._ (2015) A. Soluyanov, D. Gresch, Z. Wang, Q. Wu, M. Troyer, X. Dai, and A. Bernevig, Nature 527, 495 (2015).
  * Prochaska  _et al._ (2020) L. Prochaska, X. Li, D. C. MacFarland, A. M. Andrews, M. Bonta, E. F. Bianco, S. Yazdi, W. Schrenk, H. Detz, A. Limbeck, Q. Si, E. Ringe, G. Strasser, J. Kono, and S. Paschen, Science 367, 285 (2020).
  * Paschen and Si (2020) S. Paschen and Q. Si, Nature Reviews Physics 3, 9 (2020).
  * Petta  _et al._ (2005) J. Petta, A. Johnson, J. Taylor, E. Laird, A. Yacoby, M. Lukin, C. Marcus, M. Hanson, and A. Gossard, Science 309, 2180 (2005).
  * Bluhm  _et al._ (2010) H. Bluhm, S. Foletti, I. Neder, M. Rudner, D. Mahalu, V. Umansky, and A. Yacoby, Nature Physics 7, 109 (2010).
  * Pickett and Enss (2018) G. Pickett and C. Enss, Nature Reviews Materials 3 (2018), 10.1038/natrevmats.2018.12.
  * Eisenlohr and Vojta (2021) H. Eisenlohr and M. Vojta, Physical Review B 103 (2021), 10.1103/physrevb.103.064405.
  * Schuberth  _et al._ (2016) E. Schuberth, M. Tippmann, L. Steinke, S. Lausberg, A. Steppke, M. Brando, C. Krellner, C. Geibel, R. Yu, Q. Si, and F. Steglich, Science 351, 485 (2016).
  * Knapp  _et al._ (2023) J. Knapp, L. Levitin, J. Nyéki, A. Ho, B. Cowan, J. Saunders, M. Brando, C. Geibel, K. Kliemt, and C. Krellner, Physical Review Letters 130 (2023), 10.1103/PhysRevLett.130.126802.
  * Bitko  _et al._ (1996) D. Bitko, T. F. Rosenbaum, and G. Aeppli, Physical Review Letters 77, 940 (1996).
  * Libersky  _et al._ (2021) M. Libersky, R. McKenzie, D. Silevitch, P. Stamp, and T. Rosenbaum, Physical Review Letters 127 (2021), 10.1103/physrevlett.127.207202.
  * Giménez-Santamarina  _et al._ (2020) S. Giménez-Santamarina, S. Cardona-Serra, J. M. Clemente-Juan, A. Gaita-Ariño, and E. Coronado, Chemical Science 11, 10718 (2020).
  * Bangma  _et al._ (2023) F. Bangma, L. Levitin, M. Lucas, A. Casey, J. Nyéki, J. Saunders, S. Julian, and A. McCollam, To be published x (2023).
  * Bangma (2023) F. Bangma,  _Novel ordered phases and quantum criticality of $4f$-electron systems in high magnetic fields and ultra-low temperatures_, Ph.D. thesis, Radboud University (2023).
  * H.M.Rønnow  _et al._ (2005) H.M.Rønnow, R. Parthasarathy, J. Jensen, G. Aeppli, T. Rosenbaum, and D. McMorrow, Science 308, 389 (2005).
  * Fulde  _et al._ (1970) P. Fulde, L. L. Hirst, and A. Luther, Zeitschrift für Physik A Hadrons and nuclei 230, 155 (1970).
  * Goremychkin  _et al._ (2004) E. A. Goremychkin, R. Osborn, E. D. Bauer, M. B. Maple, N. A. Frederick, W. M. Yuhasz, F. M. Woodward, and J. W. Lynn, Physical Review Letters 93 (2004), 10.1103/physrevlett.93.157003.
  * Matsumoto and Koga (2004) M. Matsumoto and M. Koga, Journal of the Physical Society of Japan 73, 1135 (2004).
  * Bauer  _et al._ (2002) E. D. Bauer, N. A. Frederick, P.-C. Ho, V. S. Zapf, and M. B. Maple, Physical Review B 65 (2002), 10.1103/physrevb.65.100506.
  * Tayama  _et al._ (2003) T. Tayama, T. Sakakibara, H. Sugawara, Y. Aoki, and H. Sato, Journal of the Physical Society of Japan 72, 1516 (2003).
  * Measson  _et al._ (2004) M.-A. Measson, D. Braithwaite, J. Flouquet, G. Seyfarth, J. P. Brison, E. Lhotel, C. Paulsen, H. Sugawara, and H. Sato, Physical Review B 70 (2004), 10.1103/physrevb.70.064516.
  * Aoki  _et al._ (2002a) Y. Aoki, T. Namiki, S. Ohsaki, S. R. Saha, H. Sugawara, and H. Sato, Journal of the Physical Society of Japan 71, 2098 (2002a).
  * Sugawara  _et al._ (2005) H. Sugawara, M. Kobayashi, S. Osaki, S. R. Saha, T. Namiki, Y. Aoki, and H. Sato, Physical Review B 72 (2005), 10.1103/physrevb.72.014519.
  * Rotundu  _et al._ (2004) C. R. Rotundu, H. Tsujii, Y. Takano, B. Andraka, H. Sugawara, Y. Aoki, and H. Sato, Physical Review Letters 92 (2004), 10.1103/physrevlett.92.037203.
  * Nguyen  _et al._ (2021) D. H. Nguyen, A. Sidorenko, M. Taupin, G. Knebel, G. Lapertot, E. Schuberth, and S. Paschen, Nature Communications 12 (2021), 10.1038/s41467-021-24670-z.
  * Herrmannsdörfer and Pobell (2005) T. Herrmannsdörfer and F. Pobell, in  _Frontiers in Superconducting Materials_ (Springer-Verlag, 2005) pp. 71–107.
  * Andres (1973) K. Andres, Physical Review B 7, 4295 (1973).
  * Steinke  _et al._ (2013) L. Steinke, K. Mitsumoto, C. F. Miclea, F. Weickert, A. Dönni, M. Akatsu, Y. Nemoto, T. Goto, H. Kitazawa, P. Thalmeier, and M. Brando, Physical Review Letters 111 (2013), 10.1103/physrevlett.111.077202.
  * Iwakami  _et al._ (2014) O. Iwakami, Y. Namisashi, S. Abe, K. Matsumoto, G. Ano, M. Akatsu, K. Mitsumoto, Y. Nemoto, N. Takeda, T. Goto, and H. Kitazawa, Physical Review B 90 (2014), 10.1103/physrevb.90.100402.
  * McCollam  _et al._ (2013) A. McCollam, B. Andraka, and S. R. Julian, Physical Review B 88 (2013), 10.1103/physrevb.88.075102.
  * Shiina and Aoki (2004) R. Shiina and Y. Aoki, Journal of the Physical Society of Japan 73, 541 (2004).
  * Zwicknagl  _et al._ (2009) G. Zwicknagl, P. Thalmeier, and P. Fulde, Physical Review B 79 (2009), 10.1103/physrevb.79.115132.
  * Shiina (2004) R. Shiina, Journal of the Physical Society of Japan 73, 2257 (2004).
  * Kusunose  _et al._ (2009) H. Kusunose, M. Matsumoto, and M. Koga, Journal of the Physical Society of Japan 78, 094713 (2009).
  * Sachdev (2011) S. Sachdev,  _Quantum Phase Transitions (2nd edition)_ (Cambridge University Press, Cambridge, 2011).
  * McKenzie and Stamp (2018) R. McKenzie and P. Stamp, Physical Review B 97 (2018), 10.1103/PhysRevB.97.214430.
  * Kovacevic  _et al._ (2016) I. Kovacevic, P. Babkevich, M. Jeong, J. Piatek, G. Boero, and H.M.Rønnow, Physical Review B 94 (2016), 10.1103/PhysRevB.94.214433.
  * H.M.Rønnow  _et al._ (2007) H.M.Rønnow, J. Jensen, R. Parthasarathy, G. Aeppli, T. Rosenbaum, D. McMorrow, and C. Kraemer, Physical Review B 75 (2007), 10.1103/PhysRevB.75.054426.
  * Méasson  _et al._ (2008) M.-A. Méasson, D. Braithwaite, G. Lapertot, J.-P. Brison, J. Flouquet, P. Bordet, H. Sugawara, and P. C. Canfield, Physical Review B 77, 134517 (2008).
  * Andraka and Pocsy (2012) B. Andraka and K. Pocsy, Journal of Applied Physics 111, 07E115 (2012).
  * Tee  _et al._ (2012) X. Y. Tee, H. G. Luo, T. Xiang, D. Vandervelde, M. B. Salamon, H. Sugawara, H. Sato, C. Panagopoulos, and E. E. M. Chia, Physical Review B 86, 064518 (2012).
  * Seyfarth  _et al._ (2005) G. Seyfarth, J. P. Brison, M.-A. Méasson, J. Flouquet, K. Izawa, Y. Matsuda, H. Sugawara, and H. Sato, Physical Review Letters 95, 107004 (2005).
  * Sugawara  _et al._ (2002) H. Sugawara, S. Osaki, S. R. Saha, Y. Aoki, H. Sato, Y. Inada, H. Shishido, R. Settai, Y. Ōnuki, H. Harima, and K. Oikawa, Physical Review B 66 (2002), 10.1103/physrevb.66.220504.
  * Koga  _et al._ (2006) M. Koga, M. Matsumoto, and H. Shiba, Journal of the Physical Society of Japan 75, 014709 (2006).
  * Chang  _et al._ (2007) J. Chang, I. Eremin, P. Thalmeier, and P. Fulde, Physical Review B 76 (2007), 10.1103/physrevb.76.220510.
  * Namiki  _et al._ (2007) T. Namiki, Y. Aoki, H. Sato, C. Sekine, I. Shirotani, T. D. Matsuda, Y. Haga, and T. Yagi, Journal of the Physical Society of Japan 76, 093704 (2007).
  * Sato  _et al._ (2010) H. Sato, K. Tanaka, D. Kikuchi, T. Namiki, R. Higashinaka, Y. Aoki, M. Kohgi, H. Sugawara, K. Kuwahara, and K. Iwasa, Physica C: Superconductivity and its Applications 470, S525 (2010).
  * Cox and Maple (1995) D. L. Cox and M. B. Maple, Physics Today 48, 32 (1995).
  * Mathur  _et al._ (1998) N. D. Mathur, F. M. Grosche, S. R. Julian, I. R. Walker, D. M. Freye, R. K. W. Haselwimmer, and G. G. Lonzarich, Nature 394, 39 (1998).
  * Chia  _et al._ (2003) E. E. M. Chia, M. B. Salamon, H. Sugawara, and H. Sato, Physical Review Letters 91 (2003), 10.1103/physrevlett.91.247003.
  * Aoki  _et al._ (2003) Y. Aoki, A. Tsuchiya, T. Kanayama, S. R. Saha, H. Sugawara, H. Sato, W. Higemoto, A. Koda, K. Ohishi, K. Nishiyama, and R. Kadono, Physical Review Letters 91 (2003), 10.1103/physrevlett.91.067003.
  * Aoki  _et al._ (2007) Y. Aoki, T. Tayama, T. Sakakibara, K. Kuwahara, K. Iwasa, M. Kohgi, W. Higemoto, D. E. MacLaughlin, H. Sugawara, and H. Sato, Journal of the Physical Society of Japan 76, 051006 (2007).
  * Tayama  _et al._ (2006) T. Tayama, T. Sakakibara, H. Sugawara, and H. Sato, Journal of the Physical Society of Japan 75, 043707 (2006).
  * Fulde and Maki (1966) P. Fulde and K. Maki, Physical Review 141, 275 (1966).
  * Maki  _et al._ (2007) K. Maki, M.-A. Measson, D. Parker, W. Li, and H. Won, Journal of Magnetism and Magnetic Materials 310, 631 (2007).
  * Stamp and Tupitsyn (2004) P. Stamp and I. Tupitsyn, Physical Review B 69 (2004), 10.1103/PhysRevB.69.014401.
  * Fraval  _et al._ (2004) E. Fraval, M. Sellars, and J. Longdell, Physical Review Letters 92 (2004), 10.1103/PhysRevLett.92.077601.
  * Patri and Kim (2020) A. S. Patri and Y. B. Kim, Physical Review X 10 (2020), 10.1103/PhysRevX.10.041021.
  * Han  _et al._ (2022) S. Han, D. J. Schultz, and Y. B. Kim, Physical Review B 106 (2022), 10.1103/PhysRevB.106.155155.
  * Bauer  _et al._ (2001) E. D. Bauer, A. Slebarski, E. J. Freeman, C. Sirvent, and M. B. Maple, Journal of Physics: Condensed Matter 13, 4495 (2001).
  * Casey  _et al._ (2014) A. Casey, F. Arnold, L. V. Levitin, C. P. Lusher, J. Nyéki, J. Saunders, A. Shibahara, H. van der Vliet, B. Yager, D. Drung, T. Schurig, G. Batey, M. N. Cuthbert, and A. J. Matthews, Journal of Low Temperature Physics (2014), 10.1007/s10909-014-1147-z.
  * Greywall (1985) D. S. Greywall, Physical Review B 31, 2675 (1985).
  * Bean (1964) C. P. Bean, Reviews of Modern Physics 36, 31 (1964).
  * Monthoux  _et al._ (2007) P. Monthoux, D. Pines, and G. G. Lonzarich, Nature 450, 1177 (2007).
  * Kohgi  _et al._ (2003) M. Kohgi, K. Iwasa, M. Nakajima, N. Metoki, S. Araki, N. Bernhoeft, J.-M. Mignot, A. Gukasov, H. Sato, Y. Aoki, and H. Sugawara, Journal of the Physical Society of Japan 72, 1002 (2003).
  * Kuwahara  _et al._ (2005) K. Kuwahara, K. Iwasa, M. Kohgi, K. Kaneko, N. Metoki, S. Raymond, M.-A. Méasson, J. Flouquet, H. Sugawara, Y. Aoki, and H. Sato, Physical Review Letters 95 (2005), 10.1103/physrevlett.95.107003.
  * Kaneko  _et al._ (2007) K. Kaneko, N. Metoki, R. Shiina, T. D. Matsuda, M. Kohgi, K. Kuwahara, and N. Bernhoeft, Physical Review B 75 (2007), 10.1103/physrevb.75.094408.
  * Kondo (1961) J. Kondo, Journal of the Physical Society of Japan 16, 1690 (1961).
  * Bleaney (1963) B. Bleaney, Journal of Applied Physics 34, 1024 (1963).
  * Aoki  _et al._ (2002b) Y. Aoki, T. Namiki, T. D. Matsuda, K. Abe, H. Sugawara, and H. Sato, Physical Review B 65 (2002b), 10.1103/physrevb.65.064446.
  * Macfarlane  _et al._ (1982) R. M. Macfarlane, D. P. Burum, and R. M. Shelby, Physical Review Letters 49, 636 (1982).
  * Tallon  _et al._ (2011) J. L. Tallon, J. G. Storey, and J. W. Loram, Physical Review B 83 (2011), 10.1103/physrevb.83.092502.
  * Dollberg  _et al._ (2022) T. Dollberg, J. C. Andresen, and M. Schechter, Physical Review B 105 (2022), 10.1103/physrevb.105.l180413.
  * Ge  _et al._ (2014) J. Ge, J. Gutierrez, J. Li, J. Yuan, H.-B. Wang, K. Yamaura, E. Takayama-Muromachi, and V. V. Moshchalkov, Applied Physics Letters 104, 112603 (2014).
  * Müller (1989) K.-H. Müller, Physica C: Superconductivity 159, 717 (1989).


# Solutions of the Einstein equations for a black hole surrounded by a
galactic halo

R. A. Konoplya Research Centre for Theoretical Physics and Astrophysics,
Institute of Physics, Silesian University in Opava,
Bezručovo náměstí 13, CZ-74601 Opava, Czech Republic  A. Zhidenko Centro de
Matemática, Computação e Cognição (CMCC), Universidade Federal do ABC (UFABC),
Rua Abolição, CEP: 09210-180, Santo André, SP, Brazil

###### Abstract

Various profiles of matter distribution in galactic halos (such as the
Navarro-Frenk-White, Burkert, Hernquist, Moore, Taylor-Silk, and others) are
considered here as the source term for the Einstein equations. We solve these
equations and find exact solutions that represent the metric of a central
black hole immersed in a galactic halo. Even though in the general case the
solution is numerical, very accurate general analytical metrics, which include
all the particular models, are found in the astrophysically relevant regime,
when the mass of the galaxy is much smaller than the characteristic scale in
the halo.

Astrophysical black holes(98) — General relativity(641) — Galaxy dark matter
halos(1880)

##  1 Introduction

Almost every large galaxy has a supermassive black hole in its center Kormendy
& Ho (2013). Galactic matter is usually modeled by an anisotropic fluid with
some density distribution, which implies an almost spherical halo dominated by
dark matter Benson (2010). Depending on the size, mass, and form of a galaxy
one or another distribution is preferable. The generic density distribution of
a galactic halo has the following form (see, e. g., Taylor & Silk (2003)):

$\rho(r)=2^{(\gamma-\alpha)/k}\rho_{a}(r/a)^{-\alpha}(1+r^{k}/a^{k})^{-(\gamma-\alpha)/k},$
(1)

which interpolates between the slope $\alpha$ near the galactic center and the
slope $\gamma$ at large distance $r\gg a$. Here $a$ is the characteristic
scale of the galactic halo. For a dwarf galaxy, composed of about a thousand
up to several billion stars, the Burkert model Burkert (1995); Salucci &
Burkert (2000) ($\alpha=1$, $\gamma=3$, $k=2$) is suitable, while for galaxies
with the largest content of dark matter, the Navarro-Frenk-White model Navarro
et al. (1995, 1997) ($\alpha=1$, $\gamma=3$, $k=1$) is mostly used. The
Hernquist profile Hernquist (1990) ($\alpha=1$, $\gamma=4$, $k=1$) is applied
for modeling the Sérsic profiles observed in bulges and elliptical galaxies.
When supposing the cold dark-matter halos that form within cosmological N-body
simulations, the Moore model Moore et al. (1998) ($\alpha=7/5$, $\gamma=14/5$,
$k=7/5$) is considered. Within super-symmetric models the lightest neutralino
is an excellent candidate to form the universe’s cold dark matter, and they
can be observed indirectly owing to annihilation in regions of high dark-
matter density, such as centers of galactic halos Feng (2010). When studying
signals from such annihilation events, the Taylor-Silk model Taylor & Silk
(2003) ($\alpha=3/2$, $\gamma=3$, $k=3/2$) is suggested. Dynamical constraints
on such dark-matter models of galaxies were studied in Lacroix (2018), while
the impact of relativistic corrections on the detectability of dark-matter
spikes with gravitational waves were considered in Speeney et al. (2022). A
general relativistic description of a black hole surrounded by a central
region of a galaxy was given in Sadeghian et al. (2013).

The natural question in this context is whether we can ascribe a general
relativistic metric to such galactic distribution of matter that includes the
spacetime of a central black hole. One way is to consider an isolated black-
hole spacetime that is matched to some distribution of matter via the mass
function Xu et al. (2018); Zhang et al. (2021, 2022); Liu et al. (2021);
Jusufi et al. (2020); Hou et al. (2018); Konoplya (2019). In contrast to such
a cut-and-paste approach, a straightforward solution has been recently
suggested in Cardoso et al. (2021) where the problem of general relativistic
description of a central black hole immersed in the (Hernquist distribution)
galactic halo was considered self-consistently, i. e., via a solution of the
corresponding Einstein equations with the energy-momentum tensor representing
the galactic matter. Quasinormal modes, scattering, and optical phenomena for
this solution have been studied in Konoplya (2021); Stuchlík & Vrba (2021),
while an exact solution for a different equation of state of the galactic
matter was proposed in Jusufi (2022) in a similar fashion.

In the present paper we propose a general approach of this kind and find exact
solutions of the Einstein equations with the energy-momentum tensor
corresponding to various distributions of the galactic medium. Even though the
analytical solutions can be obtained only in particular cases, we show that a
very good analytical approximation can be obtained in the general case by
expanding the accurate solution in terms of the small parameter $M/a$, where
$M$ is the mass of a galaxy.

Thus, for a spherically symmetric line element,

$ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{1-2m(r)/r}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}),$ (2)

where $m(r)<r/2$ is the mass function and $f(r)>0$ is the redshift function,
we find an analytical approximate form of the metric. In particular, we will
show that the approximate metric takes the following compact form for the
Navarro-Frenk-White model Navarro et al. (1995, 1997),

$\displaystyle
f(r)=\left(1-\dfrac{r_{0}}{r}\right)\left(1-\eta\left(\dfrac{a}{r}\ln\dfrac{a}{r+a}+\mu\right)\right),$
(3) $\displaystyle
1-\dfrac{2m(r)}{r}=\left(1-\dfrac{r_{0}}{r}\right)\left(1-\eta\dfrac{a}{r+a}-\eta\dfrac{a}{r}\ln\dfrac{a}{r+a}\right),$
$\displaystyle\mu=\dfrac{a}{a+s},\qquad\eta=\dfrac{2Ms(a+s)}{a(s-r_{0})(s+(a+s)\ln(a/(a+s)))},$

and for the Burkert model Burkert (1995); Salucci & Burkert (2000),

$\displaystyle
f(r)=\left(\\!1\\!-\\!\dfrac{r_{0}}{r}\\!\right)\left(\\!1\\!-\\!\eta\\!\left(\\!\dfrac{a}{2r}\ln\dfrac{r^{2}+a^{2}}{a^{2}}+\mu-\arctan\dfrac{r}{a}\\!\right)\\!\right)\\!,$
$\displaystyle
1-\dfrac{2m(r)}{r}=\left(1-\dfrac{r_{0}}{r}\right)\left(1-\eta\dfrac{a}{2r}\ln\dfrac{r^{2}+a^{2}}{a^{2}}\right),$
(4)
$\displaystyle\mu=\arctan\dfrac{s}{a},\qquad\eta=\dfrac{4Ms}{a(s-r_{0})\ln(1+s^{2}/a^{2})}.$

Here $s$ is the radius of the galactic halo and $r_{0}$ is the radius of the
event horizon.

While for the model of Taylor & Silk (2003) the approximate analytic
expression has a rather cumbersome form, for a Taylor-Silk-like model
($\alpha=3/2$, $\gamma=3$, $k=1$), which provides the same slopes for the
density distribution in the central and far regions (though with a slightly
different interpolation between them), the corresponding metric functions are

$\displaystyle
f(r)=\left(\\!1\\!-\\!\dfrac{r_{0}}{r}\\!\right)\left(\\!1\\!-\\!\eta\\!\left(\\!\dfrac{a}{r}\operatorname{arcsinh}\sqrt{\dfrac{r}{a}}+\mu-\sqrt{\dfrac{r+a}{r}}\\!\right)\\!\right)\\!,$
$\displaystyle
1-\dfrac{2m(r)}{r}=\left(\\!1\\!-\\!\dfrac{r_{0}}{r}\\!\right)\left(\\!1\\!-\\!\eta\dfrac{a}{r}\\!\left(\\!\operatorname{arcsinh}\sqrt{\dfrac{r}{a}}\\!-\\!\sqrt{\dfrac{r}{r+a}}\\!\right)\\!\right)\\!,$
$\displaystyle\mu=\sqrt{\dfrac{s}{a+s}},\quad\eta=\dfrac{2Ms}{a(s-r_{0})(\operatorname{arcsinh}\sqrt{s/a}-\mu)}.$
(5)

We will show that for the general case given by Eq. (1), the analytic
approximation can be found in the form of the hypergeometric function.

##  2 Black hole surrounded by the galactic halo

We assume that Eq. (2) is the solution to the Einstein equations with the
stress-energy tensor corresponding to the anisotropic matter with the density
$\rho(r)$ and only the tangential pressure $P(r)$,

$T_{0}^{0}=-\rho(r),\qquad T_{2}^{2}=T_{3}^{3}=P(r).$ (6)

The Einstein equations imply

$m^{\prime}(r)=4\pi
r^{2}\rho(r),\qquad\dfrac{f^{\prime}(r)}{f(r)}=\dfrac{2m(r)}{r^{2}-2rm(r)},$
(7)

and the tangential pressure has the form,

$P(r)=\frac{r\rho(r)}{2}\frac{m(r)}{r^{2}-2rm(r)}.$ (8)

Thus, once the density distribution is specified, all the functions can be
determined by solving Eqs. (7) with the following conditions:

$m(0)=0,\qquad f(\infty)=1.$ (9)

We will consider the density distribution (1), where the constant
$\rho_{a}\equiv\rho(a)$ fixes the total mass of the galaxy

$M=m(s)=4\pi\intop_{0}^{s}\rho(r)r^{2}dr,$ (10)

and $s$ is the radius of the halo, such that $s>a\gg M$. We notice that for
$\gamma>3$ the galaxy size can be taken infinite, since, as $s\to\infty$, the
improper integral (10) converges. When $s$ is finite, in order to have a
finite total mass, we suppose that the space is empty outside the galactic
halo, i.e. $\rho(r>s)=0$.

Figure 1: Comparison of the numerical solution of equation (7) (blue) and the
analytic expression (15) (red) for $a=4M_{t}=40r_{0}$, $s=50a$ ($\alpha=1$,
$\gamma=3$, $k=1$).

When there is a black hole in the center of the galaxy, the density
distribution is modified near the event horizon, located at $r_{0}\ll M$, in
such a way that the galactic distribution of matter is reproduced in the far
zone. In the general case we consider

$\rho(r)\to\widetilde{\rho}(r)=b(r)\rho(r),$ (11)

where the prefactor $b(r)$ approaches unity for $r\gg r_{0}$. For this purpose
we can define the function $b(r)$ through the following expansion:

$b(r)=1+C_{1}\frac{r_{0}}{r}+C_{2}\frac{r_{0}^{2}}{r^{2}}+C_{3}\frac{r_{0}^{3}}{r^{3}}+\ldots.$
(12)

In particular, the choice

$b(r)=(1-r_{0}/r)^{n+1},$ (13)

sets to zero the density and its first $n$ derivatives at the event horizon.
By solving Eqs. (7) with the following conditions (cf. 9),

$m(r_{0})=r_{0}/2,\qquad f(\infty)=1,$ (14)

we obtain numerically the accurate metric functions describing the galactic
halo with the central black hole of radius $r_{0}$. For the Hernquist-type
density distribution ($\alpha=1$, $\gamma=4$, $k=1$) and $n=0$ in (13), the
resulting metric has been obtained in an analytic form in Cardoso et al.
(2021).

In order to simplify analysis in the general case, we introduce the new
functions, $A(z)$ and $B(z)$, which are finite at the horizon,

$\begin{array}[]{rcl}f(r)&=&\left(1-\dfrac{r_{0}}{r}\right)A\left(\dfrac{r}{a}\right),\\\
1-\dfrac{2m(r)}{r}&=&\left(1-\dfrac{r_{0}}{r}\right)B\left(\dfrac{r}{a}\right).\end{array}$
(15)

The functions $A(z)$ and $B(z)$ are dimensionless and must depend on the
following small dimensionless parameters: $\frac{M}{a}$, $\frac{r_{0}}{a}$ and
$\frac{r_{0}}{s}$. For our purposes, we can safely ignore the dependence on
the two latter parameters, since the black hole size is negligible comparing
to the size of the galaxy. Therefore, we have

$\begin{array}[]{rcl}A(z)&=&1-\dfrac{2M}{a}\widetilde{A}(z)+{\cal
O}\left(\dfrac{M}{a}\right)^{2},\\\
B(z)&=&1-\dfrac{2M}{a}\widetilde{B}(z)+{\cal
O}\left(\dfrac{M}{a}\right)^{2}.\end{array}$ (16)

We notice that in the region near the black hole, i. e., for $r\simeq r_{0}\ll
a$, the redshift function gains the factor

$A\left(\dfrac{r}{a}\right)\approx A(0)=1-\dfrac{2M}{a}\widetilde{A}(0),$

which corresponds to the dominant redshift correction to the frequencies due
to the galactic halo. Further, we shall calculate the value of
$\widetilde{A}(0)$ explicitly, and also note that $\widetilde{B}(0)=0$.

By solving Eqs. (7), taking the dominant order in $1/a$ and neglecting the
terms of order $r_{0}/a$, we find that $\widetilde{B}(z)$ does not depend on
the particular choice of the prefactor $b(r)$ in (11),

$\displaystyle M\widetilde{B}(z)$ $\displaystyle=$
$\displaystyle\frac{4\pi\rho_{a}a^{3}}{3-\alpha}2^{(\gamma-\alpha)/k}\times$
$\displaystyle\times$
$\displaystyle\operatorname{{}_{2}F_{1}}\left(\dfrac{3-\alpha}{k},\dfrac{\gamma-\alpha}{k},\dfrac{3-\alpha}{k}+1;-z^{k}\right)z^{2-\alpha},$

where

$\operatorname{{}_{2}F_{1}}\left(a,b,c;x\right)\equiv\sum_{n}\frac{\Gamma(a+n)\Gamma(b+n)\Gamma(c)}{\Gamma(a)\Gamma(b)\Gamma(c+n)}\cdot\frac{x^{n}}{n!}$

is the hypergeometric function. Substituting (15) into (7) and neglecting the
black-hole size as compared to the galactic scales, we obtain

$A^{\prime}(z)=A(z)(B^{-1}(z)-1)/z,$ (18)

which, after substitution of the expansion (16), leads in the dominant order
to the following relation:

$\widetilde{A}^{\prime}(z)=-\widetilde{B}(z)/z.$ (19)

Using (2) in (19) one can find explicitly $A(z)$ in terms of the generalized
hypergeometric functions,

$\displaystyle
M\widetilde{A}(z)=\mu-\frac{4\pi\rho_{a}a^{3}}{3-\alpha}2^{(\gamma-\alpha)/k}\times$
(20)
$\displaystyle\times\int\operatorname{{}_{2}F_{1}}\left(\dfrac{3-\alpha}{k},\dfrac{\gamma-\alpha}{k},\dfrac{3-\alpha}{k}+1;-z^{k}\right)z^{1-\alpha}dz.$

The constant of integration $\mu$ is fixed in order to match the Schwarzschild
geometry for $r>s$,

$g_{tt}=-g_{rr}^{-1}=-1+2M_{t}/r,\quad r>s,$ (21)

where $M_{t}=M+r_{0}/2$ is the total asymptotic mass. Thus, the values of
$\rho_{a}$ and $\mu$ are determined in terms of the asymptotic mass $M_{t}$
and the cutoff parameter $s$ by matching the Schwarzschild metric (21) at
$r=s$,

$\displaystyle\widetilde{A}\left(\dfrac{s}{a}\right)$ $\displaystyle=$
$\displaystyle\dfrac{a}{2M}\left(1-\dfrac{s-2M_{t}}{s-r_{0}}\right)\approx\dfrac{a}{s},$
(22) $\displaystyle\widetilde{B}\left(\dfrac{s}{a}\right)$ $\displaystyle=$
$\displaystyle\dfrac{a}{2M}\left(1-\dfrac{s-2M_{t}}{s-r_{0}}\right)\approx\dfrac{a}{s}.$
(23)

Even for quite large values of $M/a$ the resulting analytic expressions
approximate very well the accurate metric functions, which can be found only
numerically (see Fig. 1). For the particular cases of the Navarro-Frenk-White,
Burkert, and Taylor-Silk-like models, the above hypergeometric functions take
a relatively simple form leading to equations (3-5) for the metric functions.
The analytic approximation for the metric functions is available in the
Wolfram Mathematica® ancillary file111The ancillary file is available from
https://arxiv.org/src/2202.02205/anc..

We would like to note that within our approach the cutoff occurs not in an
arbitrary place, but at the radius of the galactic halo. Outside this radius
the (conditionally) empty space is described by the Schwarzschild metric
produced by the total mass of the halo. Therefore, the results depend on the
cutoff parameter exactly in the same way, as they depend on the size of the
galaxy. If we fix the total mass and change the size of the galaxy (and
consequently the value of $s$), we change the density of the halo, and the
observables are changed correspondingly.

##  3 Circular photon orbit and ISCO

| $M=5r_{0}$ | $M=10r_{0}$ | $M=50r_{0}$
---|---|---|---
$a/r_{0}$ | accurate | approximation | accurate | approximation | accurate | approximation
$1000$ | $0.4940-0.1840i$ | $0.4940-0.1840i$ | $0.4916-0.1831i$ | $0.4916-0.1831i$ | $0.4719-0.1758i$ | $0.4711-0.1755i$
$500$ | $0.4916-0.1831i$ | $0.4916-0.1831i$ | $0.4866-0.1813i$ | $0.4865-0.1812i$ | $0.4479-0.1668i$ | $0.4443-0.1654i$
$100$ | $0.4723-0.1758i$ | $0.4715-0.1754i$ | $0.4485-0.1668i$ | $0.4451-0.1653i$ | $0.2741-0.1013i$ | $0.0657-0.0164i$
Table 1: Fundamental ($n=0$, $\ell=1$) quasinormal mode of the electromagnetic field (in units of $r_{0}$) calculated for the accurate metric in Konoplya (2021) compared to time-domain profile values found using the analytic approximation for the Hernquist model ($\alpha=1$, $\gamma=4$, $k=1$). The QNMs are computed via the WKB method Schutz & Will (1985); Iyer & Will (1987); Konoplya (2003); Matyjasek & Opala (2017); Konoplya et al. (2019).  $k=1$ | $M=5r_{0}$ | $M=10r_{0}$ | $M=50r_{0}$
---|---|---|---
$a/r_{0}$ | accurate | approximation | accurate | approximation | accurate | approximation
$1000$ | $0.4940-0.1840i$ | $0.4940-0.1840i$ | $0.4935-0.1838i$ | $0.4935-0.1838i$ | $0.4815-0.1793i$ | $0.4812-0.1792i$
$500$ | $0.4935-0.1838i$ | $0.4935-0.1838i$ | $0.4905-0.1827i$ | $0.4905-0.1827i$ | $0.4666-0.1738i$ | $0.4653-0.1733i$
$100$ | $0.4816-0.1793i$ | $0.4813-0.1792i$ | $0.4668-0.1738i$ | $0.4655-0.1732i$ | $0.3543-0.1316i$ | $0.3118-0.1152i$
$k=2$ | $M=5r_{0}$ | $M=10r_{0}$ | $M=50r_{0}$
$a/r_{0}$ | accurate | approximation | accurate | approximation | accurate | approximation
$1000$ | $0.4950-0.1844i$ | $0.4950-0.1844i$ | $0.4934-0.1838i$ | $0.4934-0.1838i$ | $0.4808-0.1791i$ | $0.4805-0.1790i$
$500$ | $0.4934-0.1838i$ | $0.4934-0.1838i$ | $0.4902-0.1826i$ | $0.4902-0.1826i$ | $0.4652-0.1733i$ | $0.4638-0.1728i$
$100$ | $0.4808-0.1791i$ | $0.4805-0.1789i$ | $0.4652-0.1732i$ | $0.4640-0.1727i$ | $0.3460-0.1286i$ | $0.3001-0.1111i$
Table 2: Fundamental ($n=0$, $\ell=1$) quasinormal mode of the electromagnetic
field (in units of $r_{0}$) calculated via the WKB approach for the accurate
metric numeric metric ($b(r)=1-r_{0}/r$) compared to the values found using
the analytic approximation ($\alpha=1$, $\gamma=3$, $s=10a$) for the Navarro-
Frenk-White ($k=1$) and Burkert ($k=2$) models.

The shadow radius $R_{sh}$ depends only on the redshift function $f(r)$,
corresponding to the minimum

$R_{sh}=\min_{r>r_{0}}\frac{r}{\sqrt{f(r)}}=\frac{r_{ph}}{\sqrt{f(r_{ph})}},$
(24)

where $r_{ph}$ is the radius of the circular photon orbit. Substituting (15)
into (24) we find that

$r_{ph}/r_{0}=3/2+{\cal O}\left(r_{0}/a\right),$ (25)

where we neglect the radius of the black hole as compared to the
characteristic scale of the galaxy $a$. Therefore, we obtain the following
expression for the shadow radius:

$\displaystyle R_{sh}$ $\displaystyle=$
$\displaystyle\frac{3r_{0}}{2\sqrt{f(3r_{0}/2)}}+{\cal
O}\left(\frac{r_{0}}{a}\right)$ $\displaystyle=$
$\displaystyle\frac{3\sqrt{3}r_{0}}{2}\left(1+\frac{M}{a}\widetilde{A}(0)+{\cal
O}\left(\frac{r_{0}}{a}\right)+{\cal O}\left(\frac{M}{a}\right)^{2}\right).$

Taking into account that Eq. (2) implies $\widetilde{B}(0)=0$ ($\alpha<2$), we
obtain the Lyapunov exponent

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\sqrt{\left(1-\frac{2m(r_{ph})}{r_{ph}}\right)\frac{2f(r_{ph})-f^{\prime\prime}(r_{ph})}{2r_{ph}^{2}}}$
$\displaystyle=$
$\displaystyle\frac{2}{3\sqrt{3}r_{0}}\left(1-\frac{M}{a}\widetilde{A}(0)+{\cal
O}\left(\frac{r_{0}}{a}\right)+{\cal O}\left(\frac{M}{a}\right)^{2}\right).$

The radius of the innermost stable circular orbit $r_{ISCO}$ satisfies the
relation

$\frac{3f^{\prime}(r_{ISCO})}{r_{ISCO}}-\frac{2f^{\prime}(r_{ISCO})^{2}}{f(r_{ISCO})}+f^{\prime\prime}(r_{ISCO})=0.$
(28)

Substituting (15) into (28) we obtain

$r_{ISCO}/r_{0}=3+{\cal O}\left(r_{0}/a\right)^{2-\alpha}.$ (29)

Neglecting the black-hole size as compared to the characteristic scale of the
galaxy we find the corresponding frequency at ISCO,

$\displaystyle\Omega_{ISCO}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{f^{\prime}(r_{ISCO})}{2r_{ISCO}}}$ $\displaystyle=$
$\displaystyle\frac{1}{3\sqrt{6}r_{0}}\left(1-\frac{M}{a}\widetilde{A}(0)+{\cal
O}\left(\frac{r_{0}}{a}\right)+{\cal O}\left(\frac{M}{a}\right)^{2}\right).$

Figure 2: The redshift factor $\widetilde{A}(0)$ as a function of the galactic
size $s$ for $\gamma=3$: $\alpha=1$, $k=1$ (red) and $k=2$ (magenta),
$\alpha=3/2$, $k=1$ (black) and $k=3/2$ (blue).

It is well known that the high-frequency (eikonal) quasinormal modes of test
fields and, at least in a great number of cases, of gravitational
perturbations are fully determined by the circular frequency and Lyapunov
exponent of a null ray orbiting around the black hole Cardoso et al. (2009);
Konoplya & Stuchlík (2017). Thus, the quasinormal frequencies in the eikonal
regime ($\ell\rightarrow\infty$) and the ISCO frequency gain the same redshift
due to the galactic halo with the factor $\widetilde{A}(0)$, which depends on
$s$,

$\widetilde{A}(0)\approx\dfrac{(3-\alpha)a}{(2-\alpha)s}\dfrac{\operatorname{{}_{2}F_{1}}\left(\frac{2-\alpha}{k},\frac{\gamma-\alpha}{k},\frac{2-\alpha}{k}+1;-\left(\frac{s}{a}\right)^{k}\right)}{\operatorname{{}_{2}F_{1}}\left(\frac{3-\alpha}{k},\frac{\gamma-\alpha}{k},\frac{3-\alpha}{k}+1;-\left(\frac{s}{a}\right)^{k}\right)}.$
(31)

For $\gamma=3$, $\widetilde{A}(0)$ goes to zero as $s$ grows (see Fig. 2),
because the constant halo mass leads to the vanishing density in this limit.
For $\gamma=4$ in the limit $s\to\infty$ Eq. (31) reads

$\lim_{s\to\infty}\widetilde{A}(0)=\dfrac{(3-\alpha)\Gamma\left(\frac{2}{k}\right)\Gamma\left(\frac{2-\alpha}{k}+1\right)}{(2-\alpha)\Gamma\left(\frac{1}{k}\right)\Gamma\left(\frac{3-\alpha}{k}+1\right)}.$
(32)

##  4 Accuracy of the approximation

First of all, we will compare the analytic solution for the particular case of
$\alpha=1$, $\gamma=4$, $k=1$ (cf. (6) of Cardoso et al. (2021)) and our
approximation (16), yielding

$m(r)=\frac{r_{0}}{2}+\frac{Mr^{2}}{(r+a)^{2}}\left(1-\frac{r_{0}}{r}\right),$
(33)

which corresponds to the following density distribution (cf. (10) of Cardoso
et al. (2021))

$4\pi\rho(r)=\frac{m^{\prime}(r)}{r^{2}}=\frac{M}{r(r+a)^{3}}\left(2a+r_{0}-\frac{ar_{0}}{r}\right).$
(34)

Notice that $\rho(r_{0})\neq 0$, because we have neglected some terms
proportional to the black-hole radius. The corresponding approximation for the
redshift function takes the simple form

$f(r)=\left(1-\frac{r_{0}}{r}\right)\left(1-\frac{2M}{r+a}\right),$ (35)

which coincides with (7) of Cardoso et al. (2021) within the considered
approximation. The redshift factor (32) for $\alpha=1$, $\gamma=4$, $k=1$ is
unity, so that for the eikonal quasinormal modes and the ISCO frequency we
have (cf. (12) and (16) of Cardoso et al. (2021))

$\Omega=\Omega_{0}(1-M/a+{\cal O}\left(r_{0}/a\right)+{\cal
O}\left(M/a\right)^{2}).$

Comparison of the metric coefficients of the accurate numerical or analytical
solution and the approximate one is not meaningful, because the metric
coefficients are not observable gauge invariant characteristics. Instead we
will compare the dominant proper oscillation frequencies, called _quasinormal
modes_ (QNMs) Konoplya & Zhidenko (2011); Kokkotas & Schmidt (1999), which are
sensitive to the near-horizon behavior. From the table 1 we see that the
approximation provides good estimations for the quasinormal modes for the
electromagnetic perturbations, even for large black holes ($M=5r_{0}$) as long
as $M/a$ is small. A similar behavior we observe for the other models,
examples of which are shown on Tables 2 for the Navarro-Frenk-White and
Burkert profiles.

##  5 Conclusions

When constructing the metric of a supermassive black hole immersed in the
galactic halo, a cut-and-paste approach is usually used, which simply matches
the Schwarzschild solution with the weak field regime matter distribution via
the mass function. On the contrary to this approach, here we developed the
fully general relativistic approach and found self-consistent solutions to the
Einstein equations describing a black hole immersed in some general
distribution of matter (1) which includes various profiles used for modeling
the galactic halo. In the astrophysically motivated range of parameters the
general analytical expression for the metric functions has been obtained in
the form of the hypergeometric functions and the excellent accuracy of this
expression is confirmed via analysis of electromagnetic quasinormal modes,
frequencies at ISCO and the radius of the black hole shadow. Even though the
influence of the galactic environment is relatively small for the radiation
processes around central black holes, they might be potentially observable in
future, for example, when detecting quasinormal modes, due to many cycles of
rotation of a binary system before the merger in the galactic medium Cardoso
et al. (2021) or in optical phenomena owing to the dark-matter spikes in the
central region Nampalliwar et al. (2021). The current and expected in the near
future sensitivity of the gravitational wave detectors is certainly not
sufficient to detect the influence of the galactic environment.

A. Z. was supported by Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq). R. K. would like to acknowledge support of the grant
19-03950S of Czech Science Foundation (GAČR).

## References

  * Benson (2010) Benson, A. J. 2010, Phys. Rept., 495, 33, doi: 10.1016/j.physrep.2010.06.001
  * Burkert (1995) Burkert, A. 1995, Astrophys. J. Lett., 447, L25, doi: 10.1086/309560
  * Cardoso et al. (2021) Cardoso, V., Destounis, K., Duque, F., Macedo, R. P., & Maselli, A. 2021.  https://arxiv.org/abs/2109.00005
  * Cardoso et al. (2009) Cardoso, V., Miranda, A. S., Berti, E., Witek, H., & Zanchin, V. T. 2009, Phys. Rev. D, 79, 064016, doi: 10.1103/PhysRevD.79.064016
  * Feng (2010) Feng, J. L. 2010, Ann. Rev. Astron. Astrophys., 48, 495, doi: 10.1146/annurev-astro-082708-101659
  * Hernquist (1990) Hernquist, L. 1990, Astrophys. J., 356, 359, doi: 10.1086/168845
  * Hou et al. (2018) Hou, X., Xu, Z., Zhou, M., & Wang, J. 2018, JCAP, 07, 015, doi: 10.1088/1475-7516/2018/07/015
  * Iyer & Will (1987) Iyer, S., & Will, C. M. 1987, Phys. Rev. D, 35, 3621, doi: 10.1103/PhysRevD.35.3621
  * Jusufi (2022) Jusufi, K. 2022.  https://arxiv.org/abs/2202.00010
  * Jusufi et al. (2020) Jusufi, K., Jamil, M., & Zhu, T. 2020, Eur. Phys. J. C, 80, 354, doi: 10.1140/epjc/s10052-020-7899-5
  * Kokkotas & Schmidt (1999) Kokkotas, K. D., & Schmidt, B. G. 1999, Living Rev. Rel., 2, 2, doi: 10.12942/lrr-1999-2
  * Konoplya (2003) Konoplya, R. A. 2003, Phys. Rev. D, 68, 024018, doi: 10.1103/PhysRevD.68.024018
  * Konoplya (2019) —. 2019, Phys. Lett. B, 795, 1, doi: 10.1016/j.physletb.2019.05.043
  * Konoplya (2021) —. 2021, Phys. Lett. B, 823, 136734, doi: 10.1016/j.physletb.2021.136734
  * Konoplya & Stuchlík (2017) Konoplya, R. A., & Stuchlík, Z. 2017, Phys. Lett. B, 771, 597, doi: 10.1016/j.physletb.2017.06.015
  * Konoplya & Zhidenko (2011) Konoplya, R. A., & Zhidenko, A. 2011, Rev. Mod. Phys., 83, 793, doi: 10.1103/RevModPhys.83.793
  * Konoplya et al. (2019) Konoplya, R. A., Zhidenko, A., & Zinhailo, A. F. 2019, Class. Quant. Grav., 36, 155002, doi: 10.1088/1361-6382/ab2e25
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, Ann. Rev. Astron. Astrophys., 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Lacroix (2018) Lacroix, T. 2018, Astron. Astrophys., 619, A46, doi: 10.1051/0004-6361/201832652
  * Liu et al. (2021) Liu, D., Yang, Y., Wu, S., et al. 2021, Phys. Rev. D, 104, 104042, doi: 10.1103/PhysRevD.104.104042
  * Matyjasek & Opala (2017) Matyjasek, J., & Opala, M. 2017, Phys. Rev. D, 96, 024011, doi: 10.1103/PhysRevD.96.024011
  * Moore et al. (1998) Moore, B., Governato, F., Quinn, T. R., Stadel, J., & Lake, G. 1998, Astrophys. J. Lett., 499, L5, doi: 10.1086/311333
  * Nampalliwar et al. (2021) Nampalliwar, S., Kumar, S., Jusufi, K., et al. 2021, Astrophys. J., 916, 116, doi: 10.3847/1538-4357/ac05cc
  * Navarro et al. (1995) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1995, Mon. Not. Roy. Astron. Soc., 275, 720, doi: 10.1093/mnras/275.3.720
  * Navarro et al. (1997) —. 1997, Astrophys. J., 490, 493, doi: 10.1086/304888
  * Sadeghian et al. (2013) Sadeghian, L., Ferrer, F., & Will, C. M. 2013, Phys. Rev. D, 88, 063522, doi: 10.1103/PhysRevD.88.063522
  * Salucci & Burkert (2000) Salucci, P., & Burkert, A. 2000, Astrophys. J. Lett., 537, L9, doi: 10.1086/312747
  * Schutz & Will (1985) Schutz, B. F., & Will, C. M. 1985, Astrophys. J. Lett., 291, L33, doi: 10.1086/184453
  * Speeney et al. (2022) Speeney, N., Antonelli, A., Baibhav, V., & Berti, E. 2022.  https://arxiv.org/abs/2204.12508
  * Stuchlík & Vrba (2021) Stuchlík, Z., & Vrba, J. 2021, JCAP, 11, 059, doi: 10.1088/1475-7516/2021/11/059
  * Taylor & Silk (2003) Taylor, J. E., & Silk, J. 2003, Mon. Not. Roy. Astron. Soc., 339, 505, doi: 10.1046/j.1365-8711.2003.06201.x
  * Xu et al. (2018) Xu, Z., Hou, X., Gong, X., & Wang, J. 2018, JCAP, 09, 038, doi: 10.1088/1475-7516/2018/09/038
  * Zhang et al. (2022) Zhang, C., Zhu, T., Fang, X., & Wang, A. 2022.  https://arxiv.org/abs/2201.11352
  * Zhang et al. (2021) Zhang, C., Zhu, T., & Wang, A. 2021, Phys. Rev. D, 104, 124082, doi: 10.1103/PhysRevD.104.124082