# The Q Branch Cooling Anomaly Can Be Explained by Mergers of White Dwarfs and
Subgiant Stars

Ken J. Shen Department of Astronomy and Theoretical Astrophysics Center,
University of California, Berkeley, CA 94720, USA  Simon Blouin Department of
Physics and Astronomy, University of Victoria, Victoria, BC V8W 2Y2, Canada
Katelyn Breivik Center for Computational Astrophysics, Flatiron Institute, 162
Fifth Ave., New York, NY 10010, USA

###### Abstract

_Gaia_ ’s exquisite parallax measurements allowed for the discovery and
characterization of the Q branch in the Hertzsprung-Russell diagram, where
massive C/O white dwarfs (WDs) pause their dimming due to energy released
during crystallization. Interestingly, the fraction of old stars on the Q
branch is significantly higher than in the population of WDs that will become
Q branch stars or that were Q branch stars in the past. From this, Cheng et
al. inferred that $\sim 6\%$ of WDs passing through the Q branch experience a
much longer cooling delay than that of standard crystallizing WDs. Previous
attempts to explain this cooling anomaly have invoked mechanisms involving
super-solar initial metallicities. In this paper, we describe a novel scenario
in which a standard composition WD merges with a subgiant star. The evolution
of the resulting merger remnant leads to the creation of a large amount of
26Mg, which, along with the existing 22Ne, undergoes a distillation process
that can release enough energy to explain the Q branch cooling problem without
the need for atypical initial abundances. The anomalously high number of old
stars on the Q branch may thus be evidence that mass transfer from subgiants
to WDs leads to unstable mergers.

Binaries; Cosmochronology; Nucleosynthesis; White dwarf stars

††software: matplotlib (Hunter, 2007), MESA (Paxton et al., 2011, 2013, 2015,
2018, 2019), COSMIC (Breivik et al., 2020)

##  1 Introduction

As white dwarfs (WDs) evolve and cool, their interiors eventually undergo a
liquid to solid phase transition, and they crystallize (van Horn, 1968). The
process of crystallization releases latent heat and gravitational binding
energy due to mixing induced by the phase separation between the oxygen-rich
solid and oxygen-poor liquid (see Bauer 2023 for a recent review). These
additional luminosity sources cause the cooling evolution of a WD to slow
until the bulk of the WD has crystallized, which should result in a buildup of
WDs on the cooling sequence corresponding to the time of crystallization.

Using _Gaia_ ’s revolutionary parallax measurements (Gaia Collaboration et
al., 2016, 2018a), Gaia Collaboration et al. (2018b) and Tremblay et al.
(2019) resolved the WD crystallization overdensity in the Hertzsprung-Russell
diagram, which was nicknamed the Q branch due to its high concentration of DQ
(carbon-polluted atmosphere) WDs. However, several mysteries quickly became
apparent. For one, the hotter section of the Q branch corresponds to the
crystallization of ultramassive WDs with masses of $1.1-1.2\,M_{\odot}$.
However, WDs arising from single-star evolution are composed of C/O up to
$\simeq 1.05\,M_{\odot}$ and O/Ne above this (Siess 2007; Lauffer et al. 2018;
although see Althaus et al. 2021),111Denissenkov et al. (2013) proposed the
formation of hybrid C/O/Ne WDs near this boundary due to the quenching of the
inwardly propagating carbon-burning flame before it reaches the center of the
WD, resulting in a C/O core surrounded by an O/Ne mantle. However, Lecoanet et
al. (2016), Brooks et al. (2017), and Schwab & Garaud (2019) pointed out
several ways such a structure could not form and would not be stable even if
it did form, so we do not consider hybrid C/O/Ne WDs further in this work. and
the massive end of the Q branch coincides with the location of C/O WD
crystallization, not O/Ne WD crystallization (Camisassa et al., 2019; Bauer et
al., 2020; Camisassa et al., 2022). Furthermore, Cheng et al. (2019) realized
that the overdensity of WDs on the Q branch is higher than predicted by the
cooling delay due to the luminosity of latent heat release and the overturn of
oxygen during crystallization; another source of energy is required that
yields an extra $\sim 8$-Gyr cooling delay for $\sim 6\%$ of the WDs passing
through the Q branch.

Several recent studies have proposed solutions to some of these puzzles
(Blouin et al., 2020; Bauer et al., 2020; Caplan et al., 2020, 2021; Blouin et
al., 2021; Camisassa et al., 2021; Wu et al., 2022; Fleury et al., 2022,
2023). Of these, Blouin et al. (2021)’s work on 22Ne distillation is of
particularly interest, but this process requires a very high 22Ne abundance to
explain the multi-Gyr delay necessary to fully account for the properties of
ultramassive WDs on the Q branch. In this work, we show that accurately
accounting for compositional details during the merger of a C/O WD with the
core of a subgiant star can explain the existence of ultramassive C/O WDs that
undergo a long, multi-Gyr cooling delay during crystallization, driven by the
distillation of 22Ne and 26Mg.

##  2 Constraints on models to explain the cooling delay and a proposed
solution

The WDs that produce the overabundance on the Q branch are undergoing C/O
crystallization, because the process of O/Ne crystallization occurs at a
higher core temperature and luminosity (Camisassa et al., 2019; Bauer et al.,
2020). However, single-star evolution does not produce C/O WDs as massive as
the $\sim 1.2\,M_{\odot}$ WDs at the blue end of the Q branch. One possible
solution is the merger of two WDs that produces a $\sim 1.2\,M_{\odot}$ merger
remnant but avoids igniting carbon at any point, either explosively in a Type
Ia supernova during the merger or relatively quiescently at a later time as
the merger remnant evolves (Iben & Tutukov, 1984; Webbink, 1984; Nomoto &
Iben, 1985; Cheng et al., 2020).

During the merger, if no explosion occurs, the less massive WD is tidally
disrupted into a rapidly rotating envelope surrounding the more massive WD.
Viscous stresses in the differentially rotating material convert kinetic
energy into heat and expansion work, and the merger remnant relatively quickly
becomes a degenerate core, consisting of the initially more massive WD,
surrounded by a hot, spherical envelope made up of the less massive WD (Shen
et al., 2012; Schwab et al., 2012). This envelope then radiates and contracts,
its base temperature increases, and nuclear burning may be initiated at the
interface between the envelope and the core.

Schwab (2021) showed that any merger of two C/O WDs with a total mass $\gtrsim
1.05\,M_{\odot}$ leads to carbon ignition at the base of the contracting
envelope, which eventually destroys the carbon throughout the core; thus,
C/O+C/O WD mergers cannot be responsible for the ultramassive C/O WDs on the Q
branch. This leaves C/O + He WD mergers as a possibility, and indeed, as shown
by Wu et al. (2022), carbon-burning can be avoided in the helium-burning ashes
for C/O+He WD merger remnants with masses up to $1.2\,M_{\odot}$.

However, this progenitor scenario does not obviously explain the WDs with an
extra-long cooling delay inferred by the kinematics and quantitative
overabundance on the Q branch. The process of 22Ne distillation (Blouin et
al., 2021), in which the exclusion of 22Ne from the solid phase eventually
leads to a concentration of 22Ne in the interior of a crystallizing WD,
presents an interesting possible solution, but the implied delay is still too
short for standard compositions. The overall 22Ne abundance needs to be
several times higher than the value expected from helium-burning of CNO-
processed material with solar abundances ($X_{\rm 22Ne}\simeq 0.016$) to
explain the extra delay with 22Ne distillation. Thus, the question becomes:
can the merger of a C/O and He WD result in the production of additional 22Ne
during helium-burning?

The key to this question lies in the presence of relatively trace amounts of
other elements besides the standard helium and 14N that make up the bulk of a
He WD. In particular, during and after the merging event, it is plausible that
some amount of carbon will be dredged up from the underlying core via Kelvin-
Helmholtz instabilities and rotational mixing (see Heger et al. 2000 for a
review of such mixing processes). In addition, a small amount of hydrogen is
present on the surface of He WDs and will be incorporated into the merger
remnant’s helium envelope. As the merger remnant evolves and the envelope
contracts and becomes hotter, the hydrogen will undergo CNO burning, which
will process the dredged-up 12C primarily into 14N, since the proton capture
onto 14N is the slow step of the CNO cycle. Once all of the hydrogen is
consumed, the envelope contracts further and begins burning helium, converting
the 14N into 22Ne; see Clayton et al. (2007) and Menon et al. (2013) for
analogous considerations of the compositions of R Coronae Borealis stars
arising from lower-mass C/O + He WD merger remnants.

The resulting number abundance of 14N following the consumption of hydrogen is

$\displaystyle n_{\rm 14N,f}$ $\displaystyle=n_{\rm 14N,i}+{\rm
min}\left(\frac{n_{\rm H,i}}{2},n_{\rm 12C,i}\right),$ (1)

corresponding to a mass fraction of

$\displaystyle X_{\rm 14N,f}$ $\displaystyle=X_{\rm 14N,i}+{\rm
min}\left(7X_{\rm H,i},\frac{7X_{\rm 12C,i}}{6}\right),$ (2)

where the subscripts $i$ and $f$ refer to the states before and after
hydrogen-burning is complete, respectively. The final 14N abundance depends on
which element is the limiting factor, hence the min() function: if the initial
hydrogen abundance is too low to convert all of the 12C into 14N, the excess
12C nuclei will undergo standard $\alpha$-captures during the subsequent
helium-burning phase; conversely, if the initial carbon mass fraction is too
low to absorb all of the hydrogen, the additional protons will use the
existing 14N to go around the CNO cycle.

The solar abundance of CNO yields an initial 14N mass fraction in the helium
envelope of $X_{\rm 14N,i}\simeq 0.01$. We see that a hydrogen mass fraction
of $X_{\rm H,i}>0.001$ and a dredged-up carbon mass fraction of $X_{\rm
12C,i}>0.01$ is sufficient to double the 14N mass fraction, consequently
doubling the 22Ne in the helium-burning ash.

However, a more careful consideration of the pre-merger binary evolution
suggests that there is actually much less hydrogen available when the WDs
merge. At birth, a He WD with a mass of $0.2\,(0.3)\,M_{\odot}$ possesses a
hydrogen surface layer of $\simeq 8\,(3)\times 10^{-4}\,M_{\odot}$ (Istrate et
al., 2016). At face value, this could lead to a hydrogen mass fraction of
$0.004\,(0.001)$ in the helium-rich envelope of the merger remnant. However,
because of the low binary mass ratio and the non-degeneracy of the hydrogen-
rich surface layers, the initial phase of hydrogen-rich mass transfer in these
double WD systems is expected to be stable, leading to the transfer and
subsequent ejection of the accreted hydrogen-rich layers in classical nova-
like events. It is only after the hydrogen has been transferred and ejected
and helium-rich accretion has begun that mass transfer is expected to become
unstable during the subsequent helium novae, leading to a merger (Shen et al.,
2013; Shen, 2015). Thus, the actual hydrogen mass fraction in the merger
remnant’s envelope is likely too low to lead to a significant boost in the 14N
abundance.

Such He+C/O WD mergers will still occur and yield ultramassive C/O WDs that
crystallize along the Q branch if the total merger remnant mass is $\lesssim
1.2\,M_{\odot}$. However, they will not have the necessary additional 22Ne to
cause a long-enough cooling delay to explain the observed Q branch
overdensity.

Thus, we explore an alternative to this scenario, in the form of a C/O WD
merging with a subgiant star crossing the Hertzsprung gap. Such systems have
been considered as supersoft X-ray sources and progenitors of Type Ia
supernovae (van den Heuvel et al., 1992; Rappaport et al., 1994; Hachisu et
al., 1996; Yungelson et al., 1996; Li & van den Heuvel, 1997; Langer et al.,
2000; Han & Podsiadlowski, 2004; Toloza et al., 2023) because hydrogen-rich
mass transfer in such a system occurs on the thermal timescale of the donor,
which places it close to or within the range for thermally stable nuclear
burning on the surface of the accreting WD (Nomoto, 1982; Shen & Bildsten,
2007; Nomoto et al., 2007; Wolf et al., 2013).

Studies that find stable mass transfer and significant mass growth of the WD
typically assume that mass transfer is fully conservative (e.g., Temmink et
al. 2023) or that mass transferred above the maximum rate for stable burning
is blown away as an optically thick wind that is ejected with the specific
angular momentum of the WD (e.g., Hachisu et al. 1996; Han & Podsiadlowski
2004). However, the influence of the donor star on the formation of such a
wind has not been considered. As shown by Shen & Quataert (2022) for the
analogous case of hydrogen-rich classical novae, the gravitational influence
of the companion precludes the formation of optically thick winds for most of
the nova outburst, so that the binary’s motion becomes the primary driver of
mass loss in close binaries. This results in a destabilizing loss of orbital
angular momentum, which has strong observational support in the cataclysmic
variable mass distribution (Nelemans et al., 2016; Schreiber et al., 2016).

Furthermore, some of the aforementioned studies completely neglect the
influence of the helium novae that occur due to the buildup of helium ash from
stable hydrogen burning. Even those that do consider helium novae use helium
accumulation efficiency calculations from Kato & Hachisu (1999, 2004), which
again do not account for the influence of the companion on the formation of
the helium nova wind.

Figure 1: Mass loss rate vs. mass of an initially $2\,M_{\odot}$ subgiant in a
binary with a $1\,M_{\odot}$ WD at an initial separation of $5\,R_{\odot}$.
Mass transfer is assumed to be fully non-conservative, with 30% of the angular
momentum in the ejected mass coming from the binary orbit.

When accounting for the influence of the donor, both of these processes – mass
transfer above the stable hydrogen-burning limit and helium novae – result in
losses of binary orbital energy and angular momentum that destabilize the
system. In Figure 1, we show the results of a binary stellar evolution
calculation with MESA 222http://mesa.sourceforge.net, version 15140 (Paxton et
al., 2011, 2013, 2015, 2018, 2019) of an initially $2\,M_{\odot}$ subgiant
with a $1\,M_{\odot}$ WD companion (modeled as a point mass) at an initial
separation of $5\,R_{\odot}$. As we are interested in the long-term evolution
of the binary, spanning many hydrogen and/or helium novae, we assume the mass
that is transferred from the subgiant is ultimately lost from the system, with
the binary orbit providing 30% of the necessary angular momentum to eject the
material (Nelemans et al., 2016). We see that when the donor mass has
decreased to $\sim 1.2\,M_{\odot}$, the mass transfer rate oscillates
dramatically and then diverges.

Figure 2: Radius of an initially $2\,M_{\odot}$ subgiant in a binary with a
$1\,M_{\odot}$ WD at an initial separation of $5\,R_{\odot}$. The initial
radius evolution, shown as a black line, corresponds to the binary calculation
shown in Figure 1. After this phase, mass is removed at a constant rate of
$10^{-6}\,M_{\odot}\,{\rm yr^{-1}}$, resulting in the red line. The blue
dotted line shows the donor’s Roche lobe radius assuming all of the angular
momentum required for ejection comes from the binary orbit.

In Figure 2, we show the donor’s radius as a function of its mass. The black
line, labeled “$\dot{M}=\dot{M}_{\rm binary}(t)$,” shows the radius during the
binary calculation described in the preceding paragraph and Figure 1. We
approximate the future evolution of the donor’s radius, shown as a red line,
by removing mass at a constant value of $10^{-6}\,M_{\odot}\,{\rm yr^{-1}}$;
while this mass loss rate is not quantitatively accurate, the corresponding
mass loss timescale is shorter than the Kelvin-Helmholtz timescale of $\sim
10^{8}$ yr and is thus qualitatively adequate for purposes of demonstration.
The point at which the helium-rich core is exposed is shown as a red circle.
The blue dotted line shows the donor’s Roche radius assuming that all of the
angular momentum required to eject the material comes from the binary orbit.

Notably, the subgiant donor’s radius does not decrease significantly as mass
is removed, and in fact increases once the non-degenerate helium-rich core is
exposed. At all times, the radius remains larger than the Roche radius. This
is in stark contrast to common envelope evolution involving a giant star: in
that case, once the degenerate core is revealed, the radius decreases by over
an order of magnitude and the donor may shrink within its Roche lobe,
resulting in a detached close binary system.

Since the helium-rich core of the subgiant is not fully degenerate and is less
condensed than the core of a giant, runaway mass transfer does not lead to the
same outcome. Instead, we expect the WD to inspiral all the way to the center
of the subgiant and become the core of a single merger remnant with a helium-
rich envelope polluted by hydrogen from the outer layers of the subgiant core
and by carbon dredged up from the WD.

The situation is analogous to the predicted mergers of low-mass WD companions
with main sequence donors due to dynamical friction induced by nova outbursts
(Nelemans et al., 2016; Schreiber et al., 2016; Metzger et al., 2021; Shen &
Quataert, 2022). However, because of the more central concentration of the
subgiant’s helium-rich core and the subgiant’s larger mass, we expect more
material to be accreted by the WD than in the case of a main sequence star +
WD merger, for which most of the main sequence star is likely ejected (Metzger
et al., 2021). We leave a quantitative calculation of the amount and
composition of the accreted material to future multi-dimensional
hydrodynamical simulations, but we conclude that it is at least plausible that
the end result of mass transfer between a $2\,M_{\odot}$ subgiant star and a
$1\,M_{\odot}$ WD is a merger between the WD and the core of the subgiant, and
the ejection of the rest of the subgiant star. Such a merger remnant would
consist of a degenerate core (the former WD) surrounded by a hot, $\sim
0.2\,M_{\odot}$ helium-rich envelope (the former subgiant core) polluted by
hydrogen from the subgiant and carbon dredged up from the WD.

##  3 The evolution of the remnant of a subgiant + white dwarf merger

In this section, we describe the evolution of such a star formed by the merger
of a $1.0\,M_{\odot}$ C/O WD and the helium-rich core of a subgiant. We assume
the core is disrupted and becomes a hot $0.2\,M_{\odot}$ envelope surrounding
the WD, in an analogous way to the expected outcome of a double WD merger
(Shen et al., 2012; Schwab et al., 2012). We again use MESA, with the
important addition of a 50-isotope nuclear network consisting of $n$, 1-2H,
3-4He, 7Li, 7,9-10Be, 8B, 12-13C, 13-15N, 14-18O, 17-19F, 18-22Ne, 21-25Na,
23-26Mg, 25-29Al, 27-30Si, 30-31P, and 31-32S. Inlists and further details are
described in the Appendix.

We first construct a $1.0\,M_{\odot}$ C/O WD, starting with a $1.0\,M_{\odot}$
helium core-burning star composed of 99% 4He and 1% 14N, by mass, as
appropriate for the stripped core of a solar-metallicity star that has
undergone previous CNO-burning. We then evolve the star, which first burns
helium in the core and then in a shell as a helium giant. Eventually, the
helium shell becomes too low in mass to support helium-burning, and the star
cools and becomes a WD.

Figure 3: Mass fractions vs. enclosed mass within a $1.0\,M_{\odot}$ C/O WD.
Panel a) shows the composition when the shell helium-burning luminosity has
decreased to $1\,L_{\odot}$. Panel b) shows profiles when the WD’s core
temperature has cooled to $10^{7}$ K, $6\times 10^{8}$ yr after the time shown
in panel a).

During the initial phases of helium-burning, the 14N captures 4He nuclei to
become 18F, which $\beta$-decays to form 18O and then subsequently captures
another alpha particle to become 22Ne. Thus, the composition of the resulting
helium-burning ash is a mixture of roughly $50/50$ C/O, by mass, plus a trace
amount of 22Ne with a mass fraction of $X_{\rm 22Ne}=0.01\times 22/14=0.016$.
However, as the core mass approaches $1.0\,M_{\odot}$, the helium-burning
shell becomes hotter, and the 22Ne undergoes an additional alpha-capture to
become 26Mg. This can be seen in panel a) of Figure 3, which shows the
composition of the newly formed WD after the helium giant phase has ended and
the helium-burning luminosity has decreased to $1\,L_{\odot}$. The 22Ne mass
fraction in the helium-burning ash begins to decrease as the core mass grows,
matching the 26Mg mass fraction after the C/O core reaches a mass of
$0.96\,M_{\odot}$.

Panel b) of Figure 3 shows the compositional profiles when the WD has cooled
to a core temperature of $10^{7}$ K, which occurs $6\times 10^{8}$ yr after
the profiles shown in panel a). Notably, because of the relatively high mean
molecular weight of the helium-burning ashes that contain 26Mg, the region
underneath the shell is unstable to thermohaline mixing (Kippenhahn et al.,
1980), which causes the outer $\sim 0.4\,M_{\odot}$ underneath the helium
layer to homogenize. However, as it will turn out, the quantitative details of
the core’s compositional profile at this time are not important because of
future mixing.

To mimic the merger of a subgiant’s core with the WD we have just constructed,
we accrete $0.2\,M_{\odot}$ of material with a composition of $X_{\rm
4He}=0.95$ and $X_{\rm 14N}=0.05$. We thus assume that the previously
discussed conversion of extra accreted hydrogen (with a mass fraction $X\geq
0.006$) and dredged-up carbon (with a mass fraction $X_{\rm 12C}\geq 0.05$)
into 14N occurs. While it would have been preferable to include this phase of
the evolution in our simulation, this proved too numerically difficult to
calculate. However, given the similar evolution calculated for lower-mass
merger remnants (Clayton et al., 2007; Menon et al., 2013), which do not
suffer from the extent of difficulties exhibited by higher-mass stars, we are
confident that such a calculation would result in the boosted 14N fraction
that we use as our initial condition.

During the accretion phase, nuclear burning is disabled. Once the desired mass
is reached, we gradually relax this constraint until nuclear reactions have
reached their proper values, at which point we begin the evolution of the
merger remnant in earnest, and the star becomes a shell helium-burning giant.
Further details of the evolution and the inlists for the calculation are
located in the Appendix.

Figure 4: Mass fractions vs. enclosed mass within a merger remnant initially
consisting of a $1.0\,M_{\odot}$ C/O core with a $0.2\,M_{\odot}$ envelope
composed of $95\%$ helium and $5\%$ 14N, by mass. Panel a) shows the
composition when the shell helium-burning luminosity has decreased to
$1\,L_{\odot}$. Panel b) shows profiles after the merger remnant has cooled to
a core temperature of $1.5\times 10^{7}$ K, $5\times 10^{8}$ yr after the time
shown in panel a).

Panel a) of Figure 4 shows the compositional profiles of the merger remnant
after the shell helium-burning phase has ended and the helium-burning
luminosity has decreased to $1\,L_{\odot}$. As expected, the 14N in the shell,
with an initial mass fraction of $X_{\rm 14N}=0.05$, has been transformed into
26Mg, with a nearly constant mass fraction of $X_{\rm 26Mg}=0.05\times
26/14=0.09$. We note that our nuclear network contains more massive elements
than Mg, so these 26Mg nuclei could in principle capture 4He, but they do not
do so under these conditions.

As before, the 26Mg in the helium-burning ashes leads to thermohaline mixing.
However, because of the larger mass of 26Mg-enriched ash and the larger
fraction of 26Mg within that ash, thermohaline mixing homogenizes the entire
core, except for a small amount of helium that diffuses to the surface. Panel
b) of Figure 4 shows the compositional profiles $5\times 10^{8}$ yr after the
time shown in panel a). By this point, which is well before crystallization
begins, the core has an essentially uniform composition. The neutron-rich
isotopes 22Ne and 26Mg have a combined mass fraction of $0.0125+0.0154=0.028$,
80% more than if the 14N present in initially solar composition material were
converted to 22Ne. This number could be increased significantly if the amounts
of accreted hydrogen and dredged-up carbon are higher than we have assumed.
Future hydrodynamic simulations will help to better quantify these initial
conditions.

##  4 A cooling delay caused by 22Ne and 26Mg distillation

Given its neutron-rich nature, 26Mg has the potential to initiate a
distillation process similar to 22Ne distillation. To investigate this
possibility, we calculate a C/O/Mg phase diagram using the semianalytic method
detailed in Medin & Cumming (2010) and Caplan et al.
(2018).333https://github.com/andrewcumming/phase_diagram_3CP Figure 5 shows a
slice of that phase diagram for a plasma coupling parameter $\Gamma_{\rm C}$
close to the crystallization temperature of a mixture with equal masses of
carbon and oxygen. The orange line represents the liquidus, the blue line
corresponds to the solidus, and the green lines indicate how the composition
changes upon freezing for different initial liquid compositions. For the low
magnesium abundances relevant here ($x_{\rm Mg}\sim 1\%$), we observe that the
solid phase is depleted in magnesium compared to the liquid: the green tie-
lines are closer to the right axis at the solidus than at the liquidus. For a
liquid with the initial $X_{\rm 26Mg}=0.0154$ of our fiducial model, the solid
is predicted to have a 26Mg mass fraction of only $\sim 0.003$.

Figure 5: C/O/Mg phase diagram. The liquidus (orange) and solidus (blue) are
connected by green lines that represent the composition difference between the
coexisting liquid and solid phases. C/O mixtures with traces of magnesium are
found close to the right axis. Note that number fractions are used in this
figure.

This 26Mg depletion is significant enough to render the solid phase buoyant,
even after considering the oxygen enrichment in the solid. We verify this
point by calculating the mass densities of the coexisting liquid and bcc solid
phases for a mixture with initially equal mass fractions of carbon and oxygen,
$X_{\rm 26Mg}=0.0154$, and $P=4\times 10^{25}\,{\rm erg}\,{\rm cm}^{-3}$,
which corresponds to the pressure at the center of a $1.2\,M_{\odot}$ C/O WD.
By employing the multi-component plasma free energy fits from Dubin (1990),
Ogata et al. (1993), and Dewitt & Slattery (2003) and accounting for the
degenerate electron gas, we determine that the mass density of the solid is
lower than that of the coexisting liquid by $0.08\%$.

When crystallization starts at the center of the star, these 26Mg-depleted
crystals will float up and melt, displacing 26Mg-rich liquid inward in the
process. This gradual transport of 26Mg towards the center liberates
gravitational energy, just as in the case of 22Ne distillation (Isern et al.,
1991; Blouin et al., 2021). Based on the known C/O/Ne phase diagram, we also
expect the floating solids to be depleted in 22Ne. Consequently, both 22Ne and
26Mg will distill, even if the initial 22Ne mass fraction is lower than the
distillation threshold value established in Blouin et al. (2021).

The ultimate outcome of this phase separation process remains uncertain due to
the increasing abundances of 22Ne and 26Mg in the central layers, which
gradually reduce the reliability of the ternary C/O/Ne and C/O/Mg phase
diagrams. Phase diagrams including all four components would be needed, and
ideally those calculations should be performed using a simulation-based
approach (e.g., Blouin & Daligault, 2021).444Indeed, we note that significant
discrepancies between the semianalytic C/O/Ne phase diagram of Caplan et al.
(2020) and the simulations of Blouin et al. (2021) arise at high Ne mass
fractions. Nonetheless, it seems plausible that most of the star’s 22Ne and
26Mg is transported to the center. First, we already know that this is the
expected outcome of simple 22Ne distillation. Second, for 26Mg, the C/O/Mg
phase diagram of Figure 5 has a similar shape to the O/Ne/Fe phase diagram
presented in Caplan et al. (2023), where the transport of Fe to the center was
found to be possible through a combination of distillation and precipitation.

Assuming that both 22Ne and 26Mg are indeed transported to the center of the
star, we can now assess the impact of their combined distillation on WD
cooling. We can approximate this effect as the distillation of a mass fraction
of 22Ne equivalent to $0.0125+(10/12)\times 0.0154=0.025$, where the factor
$10/12$ accounts for the slightly lower neutron excess per unit charge of 26Mg
compared to 22Ne. Using the same approach as Blouin et al. (2021)555The
chemical structure of the star was adjusted to reflect the expected
stratification of the merger remnant. Note also that we assume that
distillation proceeds until $X_{\rm 22Ne}+X_{\rm 26Mg}=0.3$., we find that
distillation in this star would emit $1.8\times 10^{48}\,$erg, thereby
introducing a cooling delay of $\sim 4.2\,$Gyr.

##  5 The space density of WD + subgiant merger remnants

In this section, we calculate the local space density of subgiant + C/O WD
merger remnants in order to compare to the estimated overabundance of Q branch
WDs. We use the binary population synthesis code, COSMIC666https://cosmic-
popsynth.github.io, version 3.4.8 (Breivik et al., 2020), which implements the
single star and binary evolution prescriptions contained in Hurley et al.
(2000) and Hurley et al. (2002), respectively. We focus on mergers between C/O
WDs and subgiants such that the total mass of the WD and the subgiant’s
helium-rich core is between $1.05-1.25\,M_{\odot}$.

The population of subgiant + C/O WD mergers that produce Q branch WDs is
predominantly formed through a single formation channel in which the WD
progenitor fills its Roche lobe on the asymptotic giant branch and initiates a
common envelope evolution, which dramatically shrinks the binary, thus
allowing the companion to fill its Roche lobe as a subgiant $<1$ Gyr later.
Because of this, we simulate several populations with varying assumptions for
the common envelope ejection efficiency ranging from $\alpha=0.1-1.0$. We
initialize each population assuming solar metallicity, a $50\%$ binary
fraction, an initial stellar mass function following Kroupa (2001), a uniform
mass ratio distribution (Mazeh et al., 1992), a uniform eccentricity
distribution (Geller et al., 2019), and a log-normal orbital period
distribution (e.g., Raghavan et al., 2010).

We find between $10^{-5}-10^{-4}$ subgiant + C/O WD mergers per unit solar
mass of stars, depending on the common envelope ejection efficiency. If we
assume a constant star formation rate for the past 10 Gyr (Cukanovaite et al.,
2022) and an extra 1 Gyr of cooling to crystallization temperatures of the
ultramassive C/O WDs that arise from these WD + subgiant mergers, we can
calculate the number of C/O WDs with extra cooling delays of 4.2 Gyr currently
on the Q branch for the $4\times 10^{6}\,M_{\odot}$ of stars within 250 pc, a
distance to which _Gaia_ ’s sample of Q branch WDs is essentially complete.

We obtain between $20-200$ C/O WDs with extra cooling delays on the Q branch
at present; there are 125 such WDs if we assume $\alpha=0.3$, as suggested by
observational studies of progenitors of binaries with at least one WD (e.g.,
Zorotovic et al., 2010; Scherbak & Fuller, 2023). These numbers are comparable
to the $\sim 240$ nearby Q branch WDs with extra cooling delays estimated by
Cheng et al. (2019). Given the uncertainties in our assumptions, the rough
agreement between the number of observed Q branch WDs and our binary
population synthesis estimate suggests that mergers between subgiants and C/O
WDs could produce a significant fraction, if not all, of the Q branch WDs with
delayed cooling times.

##  6 Conclusions

In this Letter, we have proposed a novel progenitor channel to explain the
overabundance of crystallizing, ultramassive ($\sim 1.2\,M_{\odot}$) C/O WDs
on the Q branch: mergers of $\sim 1.0\,M_{\odot}$ C/O WDs with the helium-rich
cores of subgiant stars. The resulting merger remnant consists of a degenerate
C/O core surrounded by a hot, helium-rich envelope, with trace amounts of
hydrogen from the outer layers of the subgiant’s core and carbon dredged up
from the WD. These pollutants undergo CNO-burning and subsequently add to the
existing 14N mass fraction. When the envelope begins helium-burning, the 14N
is burned past 22Ne to 26Mg.

After helium-burning is completed, the merger remnant cools into an
ultramassive C/O-dominated WD. The high mean molecular weight of the 26Mg-rich
ash triggers thermohaline mixing, which homogenizes the core. Finally, when
crystallization begins, the relatively large amounts of 22Ne and 26Mg
throughout the core lead to a distillation process that yields a $\sim
4.2$-Gyr cooling delay, causing the pileup of ultramassive C/O WDs on the Q
branch. The expected number of these merger remnants is in rough agreement
with the estimated number of observed Q branch WDs with extra cooling delays.

This progenitor scenario is particularly attractive because it does not
require a super-solar initial metallicity: the large 14N abundance that
eventually leads to the creation of 26Mg, distillation, and the cooling delay
is generated during the merger and subsequent thermonuclear evolution.
However, because the abundances of hydrogen and carbon that enrich the helium
envelope following the merger are unknown, the exact amount of 26Mg produced
and thus the length of the cooling delay are uncertain. Better estimates of
these values await future multi-dimensional hydrodynamic merger simulations.

Mass transfer from a subgiant star to a C/O WD has previously been assumed to
increase the WD’s mass to the Chandrasekhar limit, resulting in a Type Ia or
Iax supernova. If instead these systems merge unstably, they will be
accompanied by much dimmer transient events and the eventual creation of
atypically long-lived, ultramassive, crystallized WDs with 22Ne- and 26Mg-rich
centers. Thus, such systems may not die in explosions but instead may linger
for billions of years as exotic probes of unstable binary mass transfer, rare
thermonuclear reactions, and chemical mixing and distillation.

We thank Evan Bauer, Sihao Cheng, and Alison Miller for helpful discussions,
the anonymous referee for useful comments, and the KITP for hosting the
program “White Dwarfs as Probes of the Evolution of Planets, Stars, the Milky
Way and the Expanding Universe” during which this research was initiated. This
work was supported by NASA through the Astrophysics Theory Program
(80NSSC20K0544) and by the NSF under Grant No. PHY-1748958. S.B. is a Banting
Postdoctoral Fellow and a CITA National Fellow, supported by the Natural
Sciences and Engineering Research Council of Canada (NSERC). The Flatiron
Institute is supported by the Simons Foundation. This research used the Savio
computational cluster resource provided by the Berkeley Research Computing
program at the University of California, Berkeley (supported by the UC
Berkeley Chancellor, Vice Chancellor of Research, and Office of the CIO). In
this Appendix, we describe the MESA inlists used in this work to perform
stellar evolution calculations. The following inlist evolves a binary
consisting of a $2\,M_{\odot}$ subgiant and a $1\,M_{\odot}$ WD (modeled as a
point mass) at an initial separation of $5\,R_{\odot}$. Any mass transferred
to the WD is assumed to be lost from the system via a combination of optically
thick winds and hydrogen- and helium-powered novae. The ejected mass is
assumed to interact with the binary components, such that 30% of the angular
momentum necessary to eject the material comes from the binary’s orbital
angular momentum. &binary_controls m1 = 2.0 ! subgiant mass in Msol m2 = 1.0 !
WD mass in Msol initial_period_in_days = -1 ! set negative to use initial
separation instead initial_separation_in_Rsuns = 5 ! assume all mass
transferred to the WD is eventually lost from the system mass_transfer_beta =
1.0 ! use other_extra_jdot routine to subtract orbital angular momentum such
that 30% of the angular momentum necessary to eject the material comes from
the orbit use_other_extra_jdot = .true. ! set other angular momentum loss
mechanisms to zero for simplicity do_jdot_gr = .false. do_jdot_mb = .false. /
! end of binary_controls namelist  In order to construct the C/O WD described
in Section 3, we begin with a $1.0\,M_{\odot}$ core helium-burning star that
then evolves into a helium shell-burning giant. Once the helium shell
decreases to a mass of $10^{-3}\,M_{\odot}$, oscillations in the outer
$10^{-4}\,M_{\odot}$ cause the timestep to decrease significantly. At this
time, we set the opacity in the convective shell to a constant value of
$0.2{\rm\,cm^{2}\,g^{-1}}$ to tame these oscillations and allow the star to
evolve past this point. This obviously changes the observational signatures of
the star, but this is an acceptable tradeoff since we are only concerned with
its compositional profile, which has essentially been set by this time. Once
the star has evolved and contracted to a radius of $0.1\,R_{\odot}$, we allow
the opacity in the outer regions to return to its normal values. At this
point, we also turn on diffusion, semiconvection, and thermohaline mixing and
stop the evolution once the core temperature reaches $10^{7}$ K. The inlist
for this calculation follows: &star_job ! change the initial composition to a
helium star; the composition is specified in &controls
relax_initial_to_xaccrete = .true. ! change the nuclear reaction network to
allow for expanded helium-burning change_net = .true. new_net_name =
’50_iso.net’ ! set the outer optical depth to be 10 x 2/3 to improve stability
set_to_this_tau_factor = 10.0 set_tau_factor = .true. / ! end of star_job
namelist &controls ! disable gold tolerances to speed up the calculation
use_gold_tolerances = .false. ! turn on MLT++ to improve stability
okay_to_reduce_gradT_excess = .true. ! set initial mass initial_mass = 1.0 !
set initial composition as appropriate for a solar composition helium star
accrete_same_as_surface = .false. accrete_given_mass_fractions = .true.
num_accretion_species = 2 accretion_species_id(1) = ’he4’
accretion_species_xa(1) = 0.99 accretion_species_id(2) = ’n14’
accretion_species_xa(2) = 0.01 ! when the helium shell decreases to 1e-3 Msol
and the timestep becomes small, limit the maximum opacity ! opacity_max = 0.2
! when the photospheric radius decreases to 0.1 Rsol, remove the opacity limit
and turn on mixing photosphere_r_lower_limit = 0.1 ! use_Ledoux_criterion =
.true. ! alpha_semiconvection = 0.1 ! semiconvection_option = ’Langer_85’ !
thermohaline_coeff = 1.0 ! do_element_diffusion = .true. / ! end of controls
namelist  To construct our $1.2\,M_{\odot}$ merger remnant, we begin with our
C/O WD and accrete $0.2\,M_{\odot}$ of material with a composition of $X_{\rm
4He}=0.95$ and $X_{\rm 14N}=0.05$. We disable nuclear burning during this
phase. After the desired mass has been accreted, we gradually turn nuclear
burning back on while keeping the core at a constant temperature of $10^{8}$
K. This is hotter than the initial core temperature of $10^{7}$ K, but such a
cold temperature leads to very small thermal diffusion timesteps due to the
steep temperature gradient. However, since the helium-burning region and
resulting ash are at an even higher temperature of $3\times 10^{8}$ K, the
core temperature does not play a significant role in the evolution. As before,
oscillations in the convective envelope set in when the helium shell has
decreased to a mass of $2\times 10^{-4}\,M_{\odot}$, and we lower the opacity
in the shell to $0.2{\rm\,cm^{2}\,g^{-1}}$. The opacity is returned to its
normal values and mixing mechanisms are turned on when the merger remnant has
contracted to a radius of $0.05\,R_{\odot}$. We use the default prescription
for thermohaline mixing, but the choice is unimportant as the core’s
composition is homogenized long before crystallization occurs. The inlist for
the remnant evolution follows: &star_job ! change the nuclear reaction network
to allow for expanded helium-burning change_net = .true. new_net_name =
’50_iso.net’ ! set the outer optical depth to be 10 x 2/3 to improve stability
set_to_this_tau_factor = 10.0 set_tau_factor = .true. / ! end of star_job
namelist &controls ! disable gold tolerances to speed up the calculation
use_gold_tolerances = .false. ! turn on MLT++ to improve stability
okay_to_reduce_gradT_excess = .true. ! accrete 14N-enriched helium shell until
the merger remnant reaches 1.2 Msol mass_change = 1d-6 max_star_mass_for_gain
= 1.2 accrete_same_as_surface = .false. accrete_given_mass_fractions = .true.
num_accretion_species = 2 accretion_species_id(1) = ’he4’
accretion_species_xa(1) = 0.95 accretion_species_id(2) = ’n14’
accretion_species_xa(2) = 0.05 ! disable nuclear burning during accretion and
relax it back to normal after reaching 1.2 Msol max_abar_for_burning = -1 !
eps_nuc_factor = 1d-15 ! dxdt_nuc_factor = 1d-15 ! keep the 1.0 Msol core at a
constant temperature with other_energy routine until nuclear burning has been
returned to normal use_other_energy = .true. ! increase varcontrol_target near
the end of helium-burning, when the helium shell has decreased to 4e-3 Msol,
to speed up the calculation ! varcontrol_target = 1d-3 ! when the helium shell
decreases to 2e-4 Msol and the timestep becomes small, limit the maximum
opacity and reset varcontrol_target back to 1d-4 ! opacity_max = 0.2 ! when
the photospheric radius decreases to 0.05 Rsol, remove the opacity limit and
turn on mixing photosphere_r_lower_limit = 0.05 ! use_Ledoux_criterion =
.true. ! alpha_semiconvection = 0.1 ! semiconvection_option = ’Langer_85’ !
thermohaline_coeff = 1.0 ! do_element_diffusion = .true. / ! end of controls
namelist

## References

  * Althaus et al. (2021) Althaus, L. G., Gil-Pons, P., Córsico, A. H., et al. 2021, A&A, 646, A30, doi: 10.1051/0004-6361/202038930
  * Bauer (2023) Bauer, E. B. 2023, arXiv e-prints, arXiv:2303.10110, doi: 10.48550/arXiv.2303.10110
  * Bauer et al. (2020) Bauer, E. B., Schwab, J., Bildsten, L., & Cheng, S. 2020, ApJ, 902, 93, doi: 10.3847/1538-4357/abb5a5
  * Blouin & Daligault (2021) Blouin, S., & Daligault, J. 2021, Phys. Rev. E, 103, 043204, doi: 10.1103/PhysRevE.103.043204
  * Blouin et al. (2021) Blouin, S., Daligault, J., & Saumon, D. 2021, ApJ, 911, L5, doi: 10.3847/2041-8213/abf14b
  * Blouin et al. (2020) Blouin, S., Daligault, J., Saumon, D., Bédard, A., & Brassard, P. 2020, A&A, 640, L11, doi: 10.1051/0004-6361/202038879
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Brooks et al. (2017) Brooks, J., Schwab, J., Bildsten, L., Quataert, E., & Paxton, B. 2017, ApJ, 834, L9, doi: 10.3847/2041-8213/834/2/L9
  * Camisassa et al. (2022) Camisassa, M. E., Althaus, L. G., Koester, D., et al. 2022, MNRAS, 511, 5198, doi: 10.1093/mnras/stac348
  * Camisassa et al. (2021) Camisassa, M. E., Althaus, L. G., Torres, S., et al. 2021, A&A, 649, L7, doi: 10.1051/0004-6361/202140720
  * Camisassa et al. (2019) Camisassa, M. E., Althaus, L. G., Córsico, A. H., et al. 2019, A&A, 625, A87, doi: 10.1051/0004-6361/201833822
  * Caplan et al. (2023) Caplan, M. E., Blouin, S., & Freeman, I. F. 2023, arXiv e-prints, arXiv:2303.03409, doi: 10.48550/arXiv.2303.03409
  * Caplan et al. (2018) Caplan, M. E., Cumming, A., Berry, D. K., Horowitz, C. J., & Mckinven, R. 2018, ApJ, 860, 148, doi: 10.3847/1538-4357/aac2d2
  * Caplan et al. (2021) Caplan, M. E., Freeman, I. F., Horowitz, C. J., Cumming, A., & Bellinger, E. P. 2021, ApJ, 919, L12, doi: 10.3847/2041-8213/ac1f99
  * Caplan et al. (2020) Caplan, M. E., Horowitz, C. J., & Cumming, A. 2020, ApJ, 902, L44, doi: 10.3847/2041-8213/abbda0
  * Cheng et al. (2019) Cheng, S., Cummings, J. D., & Ménard, B. 2019, ApJ, 886, 100, doi: 10.3847/1538-4357/ab4989
  * Cheng et al. (2020) Cheng, S., Cummings, J. D., Ménard, B., & Toonen, S. 2020, ApJ, 891, 160, doi: 10.3847/1538-4357/ab733c
  * Clayton et al. (2007) Clayton, G. C., Geballe, T. R., Herwig, F., Fryer, C., & Asplund, M. 2007, ApJ, 662, 1220, doi: 10.1086/518307
  * Cukanovaite et al. (2022) Cukanovaite, E., Tremblay, P. E., Toonen, S., et al. 2022, arXiv e-prints, arXiv:2209.13919.  https://arxiv.org/abs/2209.13919
  * Denissenkov et al. (2013) Denissenkov, P. A., Herwig, F., Truran, J. W., & Paxton, B. 2013, ApJ, 772, 37, doi: 10.1088/0004-637X/772/1/37
  * Dewitt & Slattery (2003) Dewitt, H., & Slattery, W. 2003, Contributions to Plasma Physics, 43, 279, doi: 10.1002/ctpp.200310027
  * Dubin (1990) Dubin, D. H. E. 1990, Phys. Rev. A, 42, 4972, doi: 10.1103/PhysRevA.42.4972
  * Fleury et al. (2022) Fleury, L., Caiazzo, I., & Heyl, J. 2022, MNRAS, 511, 5984, doi: 10.1093/mnras/stac458
  * Fleury et al. (2023) —. 2023, MNRAS, 520, 364, doi: 10.1093/mnras/stad068
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018b, A&A, 616, A10, doi: 10.1051/0004-6361/201832843
  * Geller et al. (2019) Geller, A. M., Leigh, N. W. C., Giersz, M., Kremer, K., & Rasio, F. A. 2019, ApJ, 872, 165, doi: 10.3847/1538-4357/ab0214
  * Hachisu et al. (1996) Hachisu, I., Kato, M., & Nomoto, K. 1996, ApJ, 470, L97, doi: 10.1086/310303
  * Han & Podsiadlowski (2004) Han, Z., & Podsiadlowski, P. 2004, MNRAS, 350, 1301, doi: 10.1111/j.1365-2966.2004.07713.x
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368, doi: 10.1086/308158
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Iben & Tutukov (1984) Iben, Jr., I., & Tutukov, A. V. 1984, ApJS, 54, 335, doi: 10.1086/190932
  * Isern et al. (1991) Isern, J., Hernanz, M., Mochkovitch, R., & Garcia-Berro, E. 1991, A&A, 241, L29
  * Istrate et al. (2016) Istrate, A. G., Marchant, P., Tauris, T. M., et al. 2016, A&A, 595, A35, doi: 10.1051/0004-6361/201628874
  * Kato & Hachisu (1999) Kato, M., & Hachisu, I. 1999, ApJ, 513, L41, doi: 10.1086/311893
  * Kato & Hachisu (2004) —. 2004, ApJ, 613, L129, doi: 10.1086/425249
  * Kippenhahn et al. (1980) Kippenhahn, R., Ruschenplatt, G., & Thomas, H. C. 1980, A&A, 91, 175
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Langer et al. (2000) Langer, N., Deutschmann, A., Wellstein, S., & Höflich, P. 2000, A&A, 362, 1046, doi: 10.48550/arXiv.astro-ph/0008444
  * Lauffer et al. (2018) Lauffer, G. R., Romero, A. D., & Kepler, S. O. 2018, MNRAS, 480, 1547, doi: 10.1093/mnras/sty1925
  * Lecoanet et al. (2016) Lecoanet, D., Schwab, J., Quataert, E., et al. 2016, ApJ, 832, 71, doi: 10.3847/0004-637X/832/1/71
  * Li & van den Heuvel (1997) Li, X. D., & van den Heuvel, E. P. J. 1997, A&A, 322, L9
  * Mazeh et al. (1992) Mazeh, T., Goldberg, D., Duquennoy, A., & Mayor, M. 1992, ApJ, 401, 265, doi: 10.1086/172058
  * Medin & Cumming (2010) Medin, Z., & Cumming, A. 2010, Phys. Rev. E, 81, 036107, doi: 10.1103/PhysRevE.81.036107
  * Menon et al. (2013) Menon, A., Herwig, F., Denissenkov, P. A., et al. 2013, ApJ, 772, 59, doi: 10.1088/0004-637X/772/1/59
  * Metzger et al. (2021) Metzger, B. D., Zenati, Y., Chomiuk, L., Shen, K. J., & Strader, J. 2021, ApJ, 923, 100, doi: 10.3847/1538-4357/ac2a39
  * Nelemans et al. (2016) Nelemans, G., Siess, L., Repetto, S., Toonen, S., & Phinney, E. S. 2016, ApJ, 817, 69, doi: 10.3847/0004-637X/817/1/69
  * Nomoto (1982) Nomoto, K. 1982, ApJ, 253, 798, doi: 10.1086/159682
  * Nomoto & Iben (1985) Nomoto, K., & Iben, Jr., I. 1985, ApJ, 297, 531, doi: 10.1086/163547
  * Nomoto et al. (2007) Nomoto, K., Saio, H., Kato, M., & Hachisu, I. 2007, ApJ, 663, 1269, doi: 10.1086/518465
  * Ogata et al. (1993) Ogata, S., Iyetomi, H., Ichimaru, S., & van Horn, H. M. 1993, Phys. Rev. E, 48, 1344, doi: 10.1103/PhysRevE.48.1344
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Raghavan et al. (2010) Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1, doi: 10.1088/0067-0049/190/1/1
  * Rappaport et al. (1994) Rappaport, S., Di Stefano, R., & Smith, J. D. 1994, ApJ, 426, 692, doi: 10.1086/174106
  * Scherbak & Fuller (2023) Scherbak, P., & Fuller, J. 2023, MNRAS, 518, 3966, doi: 10.1093/mnras/stac3313
  * Schreiber et al. (2016) Schreiber, M. R., Zorotovic, M., & Wijnen, T. P. G. 2016, MNRAS, 455, L16, doi: 10.1093/mnrasl/slv144
  * Schwab (2021) Schwab, J. 2021, ApJ, 906, 53, doi: 10.3847/1538-4357/abc87e
  * Schwab & Garaud (2019) Schwab, J., & Garaud, P. 2019, ApJ, 876, 10, doi: 10.3847/1538-4357/ab113f
  * Schwab et al. (2012) Schwab, J., Shen, K. J., Quataert, E., Dan, M., & Rosswog, S. 2012, MNRAS, 427, 190, doi: 10.1111/j.1365-2966.2012.21993.x
  * Shen (2015) Shen, K. J. 2015, ApJ, 805, L6, doi: 10.1088/2041-8205/805/1/L6
  * Shen & Bildsten (2007) Shen, K. J., & Bildsten, L. 2007, ApJ, 660, 1444, doi: 10.1086/513457
  * Shen et al. (2012) Shen, K. J., Bildsten, L., Kasen, D., & Quataert, E. 2012, ApJ, 748, 35, doi: 10.1088/0004-637X/748/1/35
  * Shen et al. (2013) Shen, K. J., Guillochon, J., & Foley, R. J. 2013, ApJ, 770, L35, doi: 10.1088/2041-8205/770/2/L35
  * Shen & Quataert (2022) Shen, K. J., & Quataert, E. 2022, ApJ, 938, 31, doi: 10.3847/1538-4357/ac9136
  * Siess (2007) Siess, L. 2007, A&A, 476, 893, doi: 10.1051/0004-6361:20078132
  * Temmink et al. (2023) Temmink, K. D., Pols, O. R., Justham, S., Istrate, A. G., & Toonen, S. 2023, A&A, 669, A45, doi: 10.1051/0004-6361/202244137
  * Toloza et al. (2023) Toloza, O., Gänsicke, B. T., Guzmán-Rincón, L. M., et al. 2023, MNRAS, 523, 305, doi: 10.1093/mnras/stad1306
  * Tremblay et al. (2019) Tremblay, P.-E., Fontaine, G., Gentile Fusillo, N. P., et al. 2019, Nature, 565, 202, doi: 10.1038/s41586-018-0791-x
  * van den Heuvel et al. (1992) van den Heuvel, E. P. J., Bhattacharya, D., Nomoto, K., & Rappaport, S. A. 1992, A&A, 262, 97
  * van Horn (1968) van Horn, H. M. 1968, ApJ, 151, 227, doi: 10.1086/149432
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355, doi: 10.1086/161701
  * Wolf et al. (2013) Wolf, W. M., Bildsten, L., Brooks, J., & Paxton, B. 2013, ApJ, 777, 136, doi: 10.1088/0004-637X/777/2/136
  * Wu et al. (2022) Wu, C., Xiong, H., & Wang, X. 2022, MNRAS, 512, 2972, doi: 10.1093/mnras/stac273
  * Yungelson et al. (1996) Yungelson, L., Livio, M., Truran, J. W., Tutukov, A., & Fedorova, A. 1996, ApJ, 466, 890, doi: 10.1086/177562
  * Zorotovic et al. (2010) Zorotovic, M., Schreiber, M. R., Gänsicke, B. T., & Nebot Gómez-Morán, A. 2010, A&A, 520, A86, doi: 10.1051/0004-6361/200913658


$\mathds{H}$ $\mathds{H}$ footnotetext<EMAIL_ADDRESS>$\mathds{S}$ footnotetext<EMAIL_ADDRESS>$\mathds{N}$
footnotetext<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Unveiling Neutrino Halos with CMB Lensing

Selim C. HotinliH, William H. Miller III Department of Physics and Astronomy,
Johns Hopkins University, Baltimore, MD 21218, USA  Nashwan SabtiS, William H.
Miller III Department of Physics and Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA  Jaxon NorthN, William H. Miller III Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Marc KamionkowskiK, William H. Miller III Department of Physics and Astronomy,
Johns Hopkins University, Baltimore, MD 21218, USA

###### Abstract

The existence of a cosmic neutrino background has been inferred indirectly
from cosmological surveys through its effect on the linear-theory evolution of
primordial density perturbations, as well as from measurements of the
primordial abundances of light elements. Constraints on the masses of the
three neutrino species imply that at least two of them move non-
relativistically today. As a consequence, non-linear evolution of density
perturbations results in the formation of neutrino halos around dark-matter
halos. We study whether these neutrino halos can be detected in the
foreseeable future through measurements of weak gravitational lensing of the
cosmic microwave background, thus providing, possibly, the first beyond-
linear-theory signature of cosmic neutrinos.

##  I Introduction

A cosmic neutrino background (C$\nu$B) is an unavoidable consequence of the
hot Big-Bang model for the origin of the Universe Weinberg (1972); Dodelson
and Turner (1992). Its properties are intricately linked to the physics
governing the Universe approximately 1 second after its birth. The detection
of the C$\nu$B would therefore mark a substantial leap in unraveling the
earliest epoch of our Universe’s history.

Direct detection of these neutrinos is notoriously difficult, even by the
standards of most of neutrino physics, given that C$\nu$B neutrinos today have
energies $\lesssim 10^{-3}\,\mathrm{eV}$ and neutrino interaction strengths
scale with this energy. As early as 1962, Ref. Weinberg (1962) proposed an
inverse-beta-decay technique to directly detect neutrinos from the C$\nu$B.
Despite past and ongoing efforts in advancing this framework Cocco  _et al._
(2007), the experimental program is expected to require several more decades
to bear fruit Betti  _et al._ (2019); Alvey  _et al._ (2022a); Apponi  _et
al._ (2022). While alternative approaches to detect the C$\nu$B have been
suggested Stodolsky (1975); Langacker  _et al._ (1983); Duda  _et al._ (2001);
Eberle  _et al._ (2004); Bauer and Shergold (2021); Yoshimura  _et al._
(2015); Domcke and Spinrath (2017); Alonso  _et al._ (2019); Bauer and
Shergold (2023); Arvanitaki and Dimopoulos (2022), they often face similar
challenges in terms of sensitivity. Concurrently, an extensive laboratory
program has been established to measure fundamental properties of neutrinos,
such as their masses, with the current leading constraints set by the KATRIN
experiment, imposing an upper bound on the sum of neutrino masses of $\sum
m_{\nu}<2.4\,\mathrm{eV}$ Aker  _et al._ (2022).

From a cosmological point of view, the inferred properties of neutrinos are
contingent upon the underlying cosmological model assumed. Within the standard
$\Lambda$CDM model, neutrinos are expected to follow a Fermi-Dirac
distribution with a temperature of $1.95\,\mathrm{K}$. Within this framework,
the effects of the C$\nu$B on the linear-theory evolution of primordial
density perturbations have been seen in cosmic microwave background (CMB)
observations and galaxy surveys, see e.g. Lesgourgues and Pastor (2006) for a
review. In addition, measurements of the primordial abundances of light
elements have enabled to probe the properties of neutrinos in the very early
Universe Fields  _et al._ (2020). This has resulted in measurements of and
constraints on the relativistic and non-relativistic neutrino energy
densities, respectively, down to $O(10\%)$ precision Aghanim  _et al._ (2020);
Alvey  _et al._ (2022b). The latter measurement has yielded the strongest (but
model-dependent) constraint on the sum of neutrino masses to date, setting an
upper bound of $\sum m_{\nu}<0.12\,\mathrm{eV}$ Aghanim  _et al._ (2020).

The presence of non-zero neutrino masses has interesting implications for the
formation of cosmic structures. Measurements of solar and atmospheric
neutrinos have shown that at least two of the neutrino mass eigenstates are
heavy enough to move non-relativistically today Esteban  _et al._ (2020). As a
result, the non-linear evolution of primordial perturbations would lead to the
emergence of neutrino halos that are gravitationally bound to dark-matter
halos, but exhibit a more diffuse distribution Ringwald and Wong (2004);
Villaescusa-Navarro  _et al._ (2011); Zhang and Zhang (2018); de Salas  _et
al._ (2017); Mertsch  _et al._ (2020). Detecting these neutrino halos through
their gravitational effects would constitute a new, albeit indirect,
observational signature of the C$\nu$B. This concept has been explored in the
context of cosmic-shear measurements in Ref. Villaescusa-Navarro  _et al._
(2011), where it was shown that neutrino halos would introduce only a
(sub)percent-level effect on the shear profile, making their detection a
technical challenge.

In this paper, we study the potential of detecting neutrino halos around dark-
matter halos through their lensing influence on the CMB. On sufficiently small
angular scales, the CMB temperature/polarization can be approximated by a
gradient across the halos. Lensing of this gradient then induces a dipolar
distortion whose detailed pattern is determined by the lens mass distribution,
see e.g. (Seljak and Zaldarriaga, 2000; Lewis and King, 2006; Levy  _et al._ ,
2023; Horowitz  _et al._ , 2019; Birkinshaw and Gull, 1983; Hotinli  _et al._
, 2019, 2021). We show whether stacking CMB patches around galaxies with
future CMB experiments may enable a detection of this subtle effect, and
discuss the main challenges that need to be addressed in order to achieve this
goal.

Throughout this work, we assume a flat $\Lambda$CDM cosmology and fix the
cosmological parameters to the Planck 2018 best-fits Aghanim  _et al._ (2020):
$h=0.6727$, $\omega_{\mathrm{cdm}}=0.1202$, $\omega_{\mathrm{b}}=0.02236$,
$n_{\mathrm{s}}=0.9649$, $\ln\left(10^{10}A_{\mathrm{s}}\right)=3.045$, and
$\tau_{\mathrm{reio}}=0.0544$.

##  II Dark-matter and neutrino halos

We assume that dark-matter halos of mass $M$ follow a spherically symmetric
NFW profile, given by

$\displaystyle\rho_{\rm DM}=\frac{\rho_{\mathrm{s}}}{x(1+x)^{2}}\ ,$ (1)

where

$\displaystyle\rho_{\mathrm{s}}=\frac{M}{4\pi
r_{\mathrm{s}}^{3}\left[\ln(1+c)-c/(1+c)\right]}\ ,$ (2)

and $x=cr/r_{\mathrm{vir}}$. Here, $c$ is the concentration parameter, and
$r_{\mathrm{vir}}$ is the virial radius. We use the concentration-halo
relation from Ref. Correa  _et al._ (2015), and define the virial radius as
$r_{\mathrm{vir}}=\left(3M/(200\times
4\pi\rho_{\mathrm{crit,0}})\right)^{1/3}$. Using the Poisson equation, we can
then calculate the gravitational potential $\Phi_{\rm DM}=-4\pi
G\rho_{\mathrm{s}}r_{\mathrm{s}}^{2}{\ln(1+x)}/{x}$ and its derivative

$\displaystyle\frac{\partial\Phi_{\rm DM}}{\partial r}$ $\displaystyle=4\pi
G\rho_{\mathrm{s}}r_{\mathrm{s}}^{2}\left[\frac{\ln(1+x)}{x^{2}}-\frac{1}{x(1+x)}\right]\
,$ (3)

which we use to compute the deflection of CMB photons caused by the dark
matter.

To obtain the contribution from neutrinos to the gravitational potential, we
first solve the collisionless Boltzmann equation for the neutrino distribution
function. We do this by linearizing the Boltzmann equation, following the
approach in Refs. Singh and Ma (2003); Ringwald and Wong (2004), which after
solving gives us the neutrino density perturbation. This neutrino density
perturbation multiplied by the background density is then interpreted as the
neutrino-halo profile. In Fourier-space, the neutrino perturbation
$\hat{\delta}_{\nu}$ reads Ringwald and Wong (2004)

$\displaystyle\hat{\delta}_{\nu}(\boldsymbol{k},z,M)=\ $ $\displaystyle 4\pi
G\bar{\rho}_{\mathrm{m},0}\\!\int_{s_{\mathrm{ini}}}^{s(z)}\\!\\!\\!\\!\mathrm{d}s^{\prime}a(s^{\prime})\hat{\delta}_{\mathrm{m}}(\boldsymbol{k},s^{\prime},M)(s(z)\\!-\\!s^{\prime})$
$\displaystyle\qquad\qquad\quad\ \,\times
F\left[\frac{k(s(z)\\!-\\!s^{\prime})}{m_{\nu}}\right]\ ,$ (4)

where $s=\int\mathrm{d}t\,a^{-2}$ (and $a$ the scale factor),
$\hat{\delta}_{\mathrm{m}}$ is the matter density perturbation (which we
obtain with CLASS Blas  _et al._ (2011)), and $F$ is the Fourier transform of
the Fermi-Dirac distribution, given by

$\displaystyle
F(x)=\frac{4}{3\zeta(3)}\sum_{n=1}^{\infty}(-1)^{n+1}\frac{n}{(n^{2}+x^{2}T_{\nu}^{2})^{2}}\
,$ (5)

with $T_{\nu}=1.95(1+z)\,\mathrm{K}$. We define the real-space matter
perturbation as
$\delta_{\mathrm{m}}=\rho_{\mathrm{DM}}/\overline{\rho}_{\mathrm{m}}$, with
$\overline{\rho}_{\mathrm{m}}$ the matter background density, and use the NFW
profile as given above. The neutrino density perturbation in real space can be
obtained by Fourier transforming Eq. (II)

$\displaystyle\delta_{\nu}(r,z,M)=\frac{1}{2\pi^{2}}\int\mathrm{d}kk^{2}\frac{\sin(kr)}{kr}\hat{\delta}_{\nu}(k,z,M)\
.$ (6)

Figure 1: Left three panels: Density profiles of dark-matter (dotted) and
neutrino (solid) halos. Left-most panel shows the density profiles for three
dark-matter halo masses at redshift $z=0.5$ and for $m_{\nu}=0.3\,{\rm eV}$.
The second and third panels show the dependence of the neutrino halo profiles
on redshift and neutrino mass. Rightmost panel: The normalized lensed CMB
temperature profiles from dark matter (dotted) and neutrino (solid) halos. The
curves are normalized to unity at their maximum. CMB lensing due to neutrino
halos can be seen to peak further away from the halo center compared to the
dark-matter induced effect.

We note that while strictly speaking this approach is valid only when
$\delta_{\nu}\ll 1$, the linear approximation is generally an underestimation
of the true clustering as obtained with simulations Singh and Ma (2003);
Ringwald and Wong (2004), which makes our method here more conservative. In
what follows, we take Eq. (6) and define $\rho_{\nu}(r|z,M)\equiv f_{\nu,\rm
nl}\,\overline{\rho}_{\nu}\delta_{\nu}(r,z,M)$, where $f_{\nu,\rm nl}$ is a
fudge-factor to account for the potential mismatch between our approximation
and fits using N-1-body simulations (set to unity throughout), and
$\overline{\rho}_{\nu}$ is the background neutrino density. We then again use
the Poisson equation to calculate the neutrino halo potential and its
derivative numerically. To assess the impact of the neutrino mass scale on our
results, we consider three scenarios111Note that for
$m_{\nu}^{\mathrm{lightest}}\gtrsim 0.1\,\mathrm{eV}$ the neutrinos can be
treated as degenerate in mass Lesgourgues and Pastor (2012).: a case with
$m_{\nu}^{\mathrm{lightest}}=0.1\,\mathrm{eV}$, a higher mass scenario with
$m_{\nu}^{\mathrm{lightest}}=0.3\,\mathrm{eV}$, and the current upper limit
set by KATRIN with $m_{\nu}^{\mathrm{lightest}}=0.8\,\mathrm{eV}$ Aker  _et
al._ (2022).

We show the matter and neutrino density profiles in the first three panels of
Fig. 1. The first panel shows the profiles for three different dark-matter
halo masses $M=\\{10^{12},10^{13},10^{14}\\}\,M_{\odot}$. The neutrino halo
density can be seen to be many orders of magnitude smaller – and smoother near
the halo center – compared to the dark-matter halo. The second and third
panels show the redshift and neutrino-mass dependence of the neutrino-halo
profile. The rightmost panel shows normalized CMB lensing profiles due to dark
matter (dotted line) and neutrino (solid lines) halos, demonstrating how the
lensing signal from neutrinos peaks at larger radii compared to dark matter,
and in a way sensitive to neutrino mass.

##  III CMB lensing

CMB photons from the last-scattering surface get deflected by potential
gradients along our line-of-sight direction $\hat{\boldsymbol{n}}$ (see Ref.
(Lewis and Challinor, 2006) for a review). The lensed CMB temperature and
polarization satisfy
${X}^{l}(\hat{\boldsymbol{n}})=X^{u}(\hat{\boldsymbol{n}}+\boldsymbol{\alpha}(\hat{\boldsymbol{n}}))$,
where $X\in\\{T,Q,U\\}$ represents the dimensionless CMB temperature $T$ or
the two Stokes parameters $Q$ or $U$ (normalized by the mean CMB temperature).
Throughout this work, we will use $l$ ($u$) to denote lensed (unlensed) CMB
fluctuations. At leading order in perturbations, the lensed CMB can be
expanded as

$\begin{split}X^{u}(\hat{\boldsymbol{n}}+\boldsymbol{\alpha}(\hat{\boldsymbol{n}}))\simeq
X^{u}(\hat{\boldsymbol{n}})+\boldsymbol{\alpha}(\hat{\boldsymbol{n}})\cdot\boldsymbol{\nabla}X^{u}(\hat{\boldsymbol{n}})\
,\\\ \end{split}$ (7)

and the lensing deflection $\boldsymbol{\alpha}$ can be assumed to be a pure
gradient
$\boldsymbol{\alpha}(\hat{\boldsymbol{n}})=\boldsymbol{\nabla}\phi(\hat{\boldsymbol{n}})$,
such that

$\displaystyle\boldsymbol{\alpha}(\hat{\boldsymbol{n}})=-2\,\boldsymbol{\nabla}\\!\\!\int_{0}^{\chi_{\star}}\mathop{{\rm
d}\chi}\frac{D_{\rm ds}}{D_{\rm s}}\Phi(\chi\hat{\boldsymbol{n}})\ .$ (8)

Here, $\Phi(\chi\hat{\boldsymbol{n}})$ is the gravitational potential, $\chi$
is the comoving distance, $\chi_{\star}$ is the comoving distance to the
recombination surface, $D_{\rm ds}\equiv\chi_{\star}-\chi$, and $D_{\rm
s}\equiv\chi_{\star}$. We approximate the gravitational potential near a dark-
matter halo to be spherically symmetric around the halo center,
$\boldsymbol{\nabla}\Phi(r)=\hat{\mathbf{r}}\,\partial\Phi(r)/\partial r$, and
define the unit vector $\hat{\boldsymbol{r}}$ with origin at the halo center.

The CMB modulations due to lensing sourced by a single halo can be written in
the flat-sky approximation as

$\displaystyle\begin{split}X^{l}(\boldsymbol{r}\\!_{\perp})\\!=\\!-\frac{4}{c^{2}}[\boldsymbol{\nabla}X^{u}(\boldsymbol{r}\\!_{\perp})\\!\cdot\\!\boldsymbol{r}_{\perp}]\frac{D_{\rm
ds}}{D_{\rm s}}\\!\int_{r_{\perp}}^{\infty}\\!\\!\\!\\!\frac{{\rm
d}r}{\sqrt{{r}^{2}-{r}^{2}_{\perp}}}\frac{\partial\Phi}{\partial r}\
,\end{split}$ (9)

where $\boldsymbol{r}\\!_{\perp}$ is the transverse radial distance away from
the halo center, and $\boldsymbol{\nabla}X^{u}(\boldsymbol{r}\\!_{\perp})$ is
the unlensed CMB fluctuation gradient. Here we assume both the CMB temperature
gradient and the parameter $D_{\rm ds}$ are constant in the integration
region. Before describing our stacking analysis, we first detail the galaxy
survey configuration we consider.

##  IV Galaxy survey

For the galaxy survey, we consider the Vera Rubin Observatory (LSST) deep
sample following Ref. (Ferraro  _et al._ , 2022). We assume that the galaxy
catalog can be separated into halo mass and redshift bins, and calculate the
lensing signal on an evenly-spaced grid of 10 halo mass bins within
$\log_{10}M\in[12,15]\,M_{\odot}$ by 5 redshift bins within $z\in[0.1,3]$. The
total galaxy number density at each redshift roughly satisfies $n(z)=n_{\rm
gal}({z}/{z_{0}})^{2}\exp(-z/z_{0})/(2z_{0})$, where $z_{0}=0.3$ and we set
the galaxy number to roughly twice the number anticipated for the LSST gold
sample: $n_{\rm gal}\simeq 80\,{\rm arcmin}^{-2}$ (Abell  _et al._ , 2009;
Ferraro  _et al._ , 2022). We compute the anticipated number count of halos in
each mass-redshift bin using the standard Sheth–Tormen halo mass function
(Sheth and Tormen, 1999), normalized to the expected total galaxy number
$n(z)$ at each redshift in a volume corresponding to a sky coverage
appropriate for the LSST survey.

The individual halo masses are expected to be measured imperfectly, with
around $40\%$ error in $\ln M$, from combinations of lensing and Sunyaev-
Zeldovich (SZ) measurements. In addition, redshift determinations are subject
to photo-$z$ errors (Palmese  _et al._ , 2020; Murata  _et al._ , 2018;
Ballardini  _et al._ , 2019); for example, the anticipated LSST photo-$z$
error is $\sigma_{z}(z)=0.03(1+z)$ (Abell  _et al._ , 2009). As we show next,
using information on halo mass and redshift would lead to a more optimal
estimation of the lensing profile, and can potentially break degeneracies
between parameters that affect CMB lensing in similar ways. Here, we assume
that the effect of the measurement uncertainties of halo mass and redshift can
be accounted for by taking galaxy catalog mass-redshift bin sizes to be
sufficiently larger than these errors222Future experiments such as MegaMapper
(Schlegel  _et al._ , 2019), which is expected to have precise redshift
information for all halos observable with LSST, could allow for a much finer
binning in redshift, significantly increasing the prospects of probing
neutrino halos at higher precision..

Figure 2: Detection signal-to-noise ratio of neutrino halos for next-
generation CMB experiments and a galaxy catalog representative of the LSST
deep sample. The dashed lines are obtained by assuming frequency-dependent
foregrounds can be removed from maps, leaving the lensed black-body CMB signal
and black-body foregrounds only, while the solid ones correspond to the
detection SNR assuming standard ILC-cleaning as described in the text. In the
left (right) panel, we assume the beam of the residual CMB noise is
1.4-arcminute (15-arcsecond), appropriate to upcoming (future) CMB
experiments. The horizontal axis corresponds to the ILC-cleaned white-noise
RMS of CMB experiments, with gray shaded regions depicting anticipated noise
levels matching CMB-S4/Advanced SO (left panel) and CMB-HD (right panel). The
blue, green, and yellow lines in both panels correspond to neutrino masses
$m_{\nu}=\\{0.1,0.3,0.8\\}\,{\rm eV}$, respectively. Here we set $f_{\rm
sky}=0.4$ and $f_{\nu,\rm nl}=1$.

##  V Stacking analysis

To extract the lensing signal, we stack patches of the CMB around dark-matter
halos oriented along the unlensed CMB gradients. We apply a matched filter to
each patch that returns the unbiased minimum-variance lensing signal to
maximize the prospects of detection.

We begin by writing the observed real-space intensity map around a dark-matter
halo as composed of the lensing signal,
$\boldsymbol{\alpha}(\boldsymbol{r})\equiv\hat{\boldsymbol{r}}_{\perp}\alpha(\boldsymbol{r})$,
and all other contributions

$\tilde{X}^{l}({\boldsymbol{r}})=-\boldsymbol{\nabla}X^{u}({\boldsymbol{r}})\cdot\hat{\boldsymbol{r}}\\!_{\perp}\,{\alpha}(\boldsymbol{r})+\tilde{X}({\boldsymbol{r}})\
,$ (10)

where $\alpha(\boldsymbol{r})=\alpha_{\nu}(\boldsymbol{r})+\alpha_{\rm
DM}(\boldsymbol{r})+\alpha_{\rm b}(\boldsymbol{r})$ is the total lensing
deflection profile including neutrinos, dark matter, and baryons,
respectively. We assume that the observed lensed CMB maps can be filtered to
get unbiased and minimum-variance estimates for the norm of the unlensed CMB
gradient $\widehat{|\boldsymbol{\nabla}X^{u}|}$ for each halo. The
corresponding estimator noise
$N^{X}\equiv\langle(\widehat{|\boldsymbol{\nabla}X^{u}|}-{|\boldsymbol{\nabla}X^{u}|})^{2}\rangle$
from neutrino-lensing can be shown to satisfy (Hotinli  _et al._ , 2021)

$\displaystyle N^{X}_{\nu}(M,z)=\frac{1}{2}\int\frac{{\rm
d}^{2}\boldsymbol{\ell}}{(2\pi)^{2}}\frac{|{\alpha}_{\nu}(\boldsymbol{\ell})|^{2}}{\tilde{C}_{\ell}^{XX}}\
,$ (11)

where $\alpha_{\nu}(\boldsymbol{\ell})$ is the neutrino-halo lensing
deflection profile in Fourier space, and $\tilde{C}_{\ell}^{XX}$ is the ILC-
cleaned CMB spectrum, consisting of the lensed primary CMB signal,
foregrounds, and CMB experiment noise333Here we assumed that the neutrino-halo
lensing can be isolated from the CMB patches. In practice, lensing due to
dark-matter halos and baryons will bias the measurement of the unlensed CMB
gradients, roughly proportional to
${\sim}|\alpha_{\nu}(\boldsymbol{\ell})\alpha_{\rm DM}(\boldsymbol{\ell})|$
and ${\sim}|\alpha_{\nu}(\boldsymbol{\ell})\alpha_{\rm
b}(\boldsymbol{\ell})|$, affecting prospects to unambiguously detect the
neutrino-halo lensing. We take a fist step towards investigating this bias in
what follows..

We calculate the lensed CMB black-body temperature and polarization power-
spectra using CAMB Lewis  _et al._ (2000). The Stokes parameters spectra
$C_{\ell}^{QQ}$ and $C_{\ell}^{UU}$ are derived from $C_{\ell}^{EE}$ and
$C_{\ell}^{BB}$. We account for the foregrounds following the prescription in
the appendix of Ref. (Hotinli  _et al._ , 2021). In particular, we add the
following contributions to the temperature power spectrum: the frequency-
dependent cosmic infrared background, the thermal Sunyaev-Zeldovich (tSZ)
effect, radio sources, and the reionization and late-time kinetic Sunyaev-
Zeldovich (kSZ) effects (which are black-body). We assume that the
polarization spectra are foreground-free. The CMB noise is modeled at each
frequency channel and includes the pink and white noise components as
$N_{\ell}=\Delta_{X}^{2}\exp\left[\ell(\ell+1){\theta_{\rm
FWHM}^{2}}/{(8\log(2))}\right][1+(\ell_{\rm knee}/\ell)^{\alpha_{\rm knee}}]$,
with $\theta_{\rm FWHM}$ the beam full-width at half-maximum (FWHM),
$\Delta_{X}$ the noise root mean square (RMS), and $\\{\alpha_{\rm
knee},\ell_{\rm knee}\\}=\\{-3,100\\}$ parameters that model the effect of the
Earth’s atmosphere for upcoming ground-based experiments. For our forecast, we
consider next-generation CMB experiments, ranging from the upcoming Simons
Observatory (SO) Abitbol  _et al._ (2019) to the future CMB-S4 Abazajian  _et
al._ (2016, 2019), Advanced SO Abitbol  _et al._ (2019) and CMB-HD Sehgal  _et
al._ (2019); Aiola  _et al._ (2022), whose survey parameters we take from
table 2 of Ref. (Hotinli  _et al._ , 2021). We combine the anticipated CMB
observations at different frequency channels to calculate the minimum-variance
standard internal linear combination444Ongoing work on mitigating the effects
of foregrounds for lensing measurements, as well as delensing, and development
of more optimal lensing-potential reconstruction techniques suggests the
fidelity of lensing reconstruction in the future may improve significantly
compared to the ILC-cleaned CMB spectrum $\tilde{C}_{\ell}^{XX}$ which we
consider here (see e.g. (Namikawa  _et al._ , 2013; Millea  _et al._ , 2019;
Horowitz  _et al._ , 2019; Hadzhiyska  _et al._ , 2019; Darwish  _et al._ ,
2023; Sailer  _et al._ , 2020, 2023; Parker  _et al._ , 2022; Foreman  _et
al._ , 2023; Li  _et al._ , 2022; Han and Sehgal, 2022; Hotinli  _et al._ ,
2022)). (ILC) cleaned (Tegmark  _et al._ , 2003) CMB spectrum
$\tilde{C}_{\ell}^{XX}$.

The estimated signal-to-noise ratio (SNR) from CMB field $X$ for a halo of
mass $M$ at redshift $z$ is then
$\mathrm{SNR}=\sqrt{{\langle|\boldsymbol{\nabla}X^{u}|^{2}\rangle}/{N^{X}_{\nu}(M,z)}}$,
where we set the fiducial values of the dimensionless CMB gradient RMS
$\\{\sqrt{\langle|\boldsymbol{\nabla}T^{u}|^{2}\rangle},\sqrt{\langle|\boldsymbol{\nabla}Q^{u}|^{2}\rangle},\sqrt{\langle|\boldsymbol{\nabla}U^{u}|^{2}\rangle}\\}$
equal to555Or, equivalently,
$\\{\sqrt{\langle|\boldsymbol{\nabla}T^{u}|^{2}\rangle},\sqrt{\langle|\boldsymbol{\nabla}Q^{u}|^{2}\rangle},\sqrt{\langle|\boldsymbol{\nabla}U^{u}|^{2}\rangle}\\}$
= $\\{13,1.3,1.3\\}$ in ${\mu{\rm K}\,{\rm arcmin}^{-1}}$ units.
$\\{4.7,0.47,0.47\\}\\!\times\\!10^{-6}$. We first calculate the single-halo
detection SNRs from the three components $X\in\\{T,Q,U\\}$ for each halo mass
and redshift bin in our grid. We then weight the SNR values with the number
count of halos at each mass-redshift bin expected from LSST. Finally, we sum
these in inverse quadrature to obtain the total detection SNR from the
catalog.

##  VI Results

The SNR is shown in Fig. 2 as a function of the CMB temperature white noise
amplitude after ILC-cleaning. We also consider forecasts assuming all
frequency-dependent foregrounds can be removed from the CMB (dashed lines),
leaving only the black-body CMB signal and foregrounds. We find that at the
present cosmological mass limit within $\Lambda$CDM the SNR is too low for the
neutrino-halo-lensing effect to be detectable in the near future, unless
$f_{\nu,\mathrm{nl}}\gtrsim 10$. Nevertheless, we see that larger neutrino
masses lead to a significant increase in the detection prospects, suggesting
measurements of CMB lensing could potentially be used for putting upper bounds
on neutrino masses. Notably, future cosmological surveys may have the
sensitivity to place such an upper bound to high significance if the
frequency-dependent foregrounds can be cleaned efficiently from the CMB.

##  VII Discussion and Conclusion

In this work, we considered a new way to detect cosmic neutrinos beyond their
linear-theory effects. In particular, we studied the statistical power of
upcoming CMB and large-scale structure surveys to detect neutrino halos from
oriented stacking of lensed CMB temperature and polarization patches around
galaxy clusters. We showed that this could be feasible if the neutrino masses
are sufficiently large.

So far, we have omitted uncertainties in the modeling of the dark-matter halo
profiles and the effects of baryonic matter on the lensing profiles, both of
which would influence the SNR of an unambiguous neutrino-halo detection. While
the lensing signal caused by neutrinos have a different shape compared to dark
matter (see rightmost panel in Fig. 1) and baryons (see e.g. (Peirani  _et
al._ , 2017; Sorini  _et al._ , 2022; Chung  _et al._ , 2020)), we anticipate
that joint modeling of these effects may make it somewhat more difficult to
isolate the impact of neutrinos. As a preliminary analysis, we assess the
reduction in the detection SNR of neutrino halos when the dark-matter halo
mass $M$ and concentration parameter $c$ are modeled as free parameters. We
perform an ensemble-statistics analysis where we calculate the information
matrix $\mathcal{F}_{\alpha\beta}^{(z)}=\sum_{\ell}\\!\frac{f_{\rm
sky}({2\ell+1})}{2}[\partial_{\alpha}C_{\ell}^{TT,\nu(z)}\mathcal{C}^{-1}\partial_{\beta}C_{\ell}^{TT,\nu(z)}]$.
Here, $\tilde{C}_{\ell}^{TT,\nu(z)}$ is the neutrino-induced lensing profile
at redshift $z$, $\mathcal{C}=\tilde{C}_{\ell}^{TT,\nu(z)}+N_{\nu}^{T}$ is the
corresponding lensing estimator covariance, and indices $\alpha,\beta$ for
partial derivatives vary over parameters $\\{m_{\nu},M,c\\}$. For a given halo
of mass $M\sim 10^{13}\,M_{\odot}$ and $m_{\nu}=0.2\,{\rm eV}$, we calculate
$\mathcal{F}_{\alpha\beta}^{(z)}$ at each of the five redshift bins we
consider within $z\in[0.1,3.0]$ . We then sum these matrices after inverse-
weighting them by the number-count of halos anticipated to be in the LSST
catalog in these mass-redshift bins. We find that marginalizing over dark-
matter halo parameters leads to a modest reduction in the neutrino-halo
detection SNR by ${\sim}20\%$ compared to parameterizing the lensing signal
with $m_{\nu}$ alone, i.e.,
$\sqrt{(\mathcal{F}^{-1})_{m_{\nu}m_{\nu}}/(\mathcal{F}_{m_{\nu}m_{\nu}})^{-1}}\sim
0.8$. In a follow-up work, we will investigate the impact of modeling dark-
matter and baryonic effects jointly with neutrinos in further detail.

Another assumption we made was that the unlensed CMB temperature and
polarization gradients can be known perfectly when aligning the patches. This
is a reasonable assumption as the unlensed CMB gradients are coherent on small
scales due to Silk damping (over ${\sim}90$ percent of the RMS of
$\langle|\boldsymbol{\nabla}X^{u}|^{2}\rangle$ comes from scales $\ell\lesssim
1000$ (Hui  _et al._ , 2007)) and the observed lensed CMB fields can be
delensed very effectively for upcoming CMB surveys that cover large scales
(see e.g. (Hotinli  _et al._ , 2022; Green  _et al._ , 2017)). Nevertheless,
in practice, our knowledge of the unlensed CMB gradients will not be perfect
due to residual lensing and other interactions of CMB photons with the
intervening structure along their trajectories. We leave a more detailed
analysis of this effect on the neutrino halo detection prospects to upcoming
work.

Finally, in computing the neutrino halo profile, we solved the linearized
Boltzmann equation, which underestimates the true clustering magnitude of
massive neutrinos. To enhance the accuracy of our results, a more advanced
method, such as one based on N-1-body simulations, would capture the non-
linear behavior and result in a clustering boost by a factor of a few. We will
implement such an improved methodology in future work to allow us to refine
our findings and obtain more precise results.

The relic neutrino background is an integral part of our cosmological model,
yet its direct detection has remained elusive. While laboratory experiments
are currently being designed to directly detect these species, cosmology has
already provided us with probes to indirectly infer their existence and
properties. Neutrinos, being massive, cluster together to create neutrino
halos, a mechanism that extends beyond their linear-theory effect in
cosmology. In this study, we investigated the prospects of detecting neutrino
halos through a joint analysis of CMB lensing and galaxy surveys. Such a
detection would open up new windows to the origin of our Universe.

Acknowledgements. — We thank Niayesh Afshordi, Neal Dalal, Miguel Escudero,
Matthew Johnson, Neelima Sehgal, and Kendrick Smith for useful discussions.
SCH and NS were supported by a Horizon Fellowship from Johns Hopkins
University. MK was supported by NSF Grant No. 2112699 and the Simons
Foundation.

## References

  * Weinberg (1972) Steven Weinberg,  _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_ (John Wiley and Sons, New York, 1972).
  * Dodelson and Turner (1992) Scott Dodelson and Michael S. Turner, “Nonequilibrium neutrino statistical mechanics in the expanding universe,” Phys. Rev. D 46, 3372–3387 (1992).
  * Weinberg (1962) Steven Weinberg, “Universal Neutrino Degeneracy,” Phys. Rev. 128, 1457–1473 (1962).
  * Cocco  _et al._ (2007) Alfredo G. Cocco, Gianpiero Mangano, and Marcello Messina, “Probing low energy neutrino backgrounds with neutrino capture on beta decaying nuclei,” JCAP 06, 015 (2007), arXiv:hep-ph/0703075 .
  * Betti  _et al._ (2019) M. G. Betti _et al._ (PTOLEMY), “Neutrino physics with the PTOLEMY project: active neutrino properties and the light sterile case,” JCAP 07, 047 (2019), arXiv:1902.05508 [astro-ph.CO] .
  * Alvey  _et al._ (2022a) James Alvey, Miguel Escudero, Nashwan Sabti, and Thomas Schwetz, “Cosmic neutrino background detection in large-neutrino-mass cosmologies,” Phys. Rev. D 105, 063501 (2022a), arXiv:2111.14870 [hep-ph] .
  * Apponi  _et al._ (2022) A. Apponi _et al._ (PTOLEMY), “Heisenberg’s uncertainty principle in the PTOLEMY project: A theory update,” Phys. Rev. D 106, 053002 (2022), arXiv:2203.11228 [hep-ph] .
  * Stodolsky (1975) Leo Stodolsky, “Speculations on Detection of the Neutrino Sea,” Phys. Rev. Lett. 34, 110 (1975), [Erratum: Phys.Rev.Lett. 34, 508 (1975)].
  * Langacker  _et al._ (1983) Paul Langacker, Jacques P. Leveille, and Jon Sheiman, “On the Detection of Cosmological Neutrinos by Coherent Scattering,” Phys. Rev. D 27, 1228 (1983).
  * Duda  _et al._ (2001) Gintaras Duda, Graciela Gelmini, and Shmuel Nussinov, “Expected signals in relic neutrino detectors,” Phys. Rev. D 64, 122001 (2001), arXiv:hep-ph/0107027 .
  * Eberle  _et al._ (2004) Birgit Eberle, Andreas Ringwald, Liguo Song, and Thomas J. Weiler, “Relic neutrino absorption spectroscopy,” Phys. Rev. D 70, 023007 (2004), arXiv:hep-ph/0401203 .
  * Bauer and Shergold (2021) Martin Bauer and Jack D. Shergold, “Relic neutrinos at accelerator experiments,” Phys. Rev. D 104, 083039 (2021), arXiv:2104.12784 [hep-ph] .
  * Yoshimura  _et al._ (2015) M. Yoshimura, N. Sasao, and M. Tanaka, “Experimental method of detecting relic neutrino by atomic de-excitation,” Phys. Rev. D 91, 063516 (2015), arXiv:1409.3648 [hep-ph] .
  * Domcke and Spinrath (2017) Valerie Domcke and Martin Spinrath, “Detection prospects for the Cosmic Neutrino Background using laser interferometers,” JCAP 06, 055 (2017), arXiv:1703.08629 [astro-ph.CO] .
  * Alonso  _et al._ (2019) Rodrigo Alonso, Diego Blas, and Peter Wolf, “Exploring the ultra-light to sub-MeV dark matter window with atomic clocks and co-magnetometers,” JHEP 07, 069 (2019), arXiv:1810.00889 [hep-ph] .
  * Bauer and Shergold (2023) Martin Bauer and Jack D. Shergold, “Limits on the cosmic neutrino background,” JCAP 01, 003 (2023), arXiv:2207.12413 [hep-ph] .
  * Arvanitaki and Dimopoulos (2022) Asimina Arvanitaki and Savas Dimopoulos, “The Cosmic Neutrino Background on the Surface of the Earth,” (2022), arXiv:2212.00036 [hep-ph] .
  * Aker  _et al._ (2022) M. Aker _et al._ (KATRIN), “Direct neutrino-mass measurement with sub-electronvolt sensitivity,” Nature Phys. 18, 160–166 (2022), arXiv:2105.08533 [hep-ex] .
  * Lesgourgues and Pastor (2006) Julien Lesgourgues and Sergio Pastor, “Massive neutrinos and cosmology,” Phys. Rept. 429, 307–379 (2006), arXiv:astro-ph/0603494 .
  * Fields  _et al._ (2020) Brian D. Fields, Keith A. Olive, Tsung-Han Yeh, and Charles Young, “Big-Bang Nucleosynthesis after Planck,” JCAP 03, 010 (2020), [Erratum: JCAP 11, E02 (2020)], arXiv:1912.01132 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Alvey  _et al._ (2022b) James Alvey, Miguel Escudero, and Nashwan Sabti, “What can CMB observations tell us about the neutrino distribution function?” JCAP 02, 037 (2022b), arXiv:2111.12726 [astro-ph.CO] .
  * Esteban  _et al._ (2020) Ivan Esteban, M. C. Gonzalez-Garcia, Michele Maltoni, Thomas Schwetz, and Albert Zhou, “The fate of hints: updated global analysis of three-flavor neutrino oscillations,” JHEP 09, 178 (2020), arXiv:2007.14792 [hep-ph] .
  * Ringwald and Wong (2004) Andreas Ringwald and Yvonne Y. Y. Wong, “Gravitational clustering of relic neutrinos and implications for their detection,” JCAP 12, 005 (2004), arXiv:hep-ph/0408241 .
  * Villaescusa-Navarro  _et al._ (2011) Francisco Villaescusa-Navarro, Jordi Miralda-Escudé, Carlos Peña Garay, and Vicent Quilis, “Neutrino Halos in Clusters of Galaxies and their Weak Lensing Signature,” JCAP 06, 027 (2011), arXiv:1104.4770 [astro-ph.CO] .
  * Zhang and Zhang (2018) Jue Zhang and Xin Zhang, “Gravitational clustering of cosmic relic neutrinos in the Milky Way,” Nature Commun. 9, 1833 (2018), arXiv:1712.01153 [astro-ph.CO] .
  * de Salas  _et al._ (2017) P. F. de Salas, S. Gariazzo, J. Lesgourgues, and S. Pastor, “Calculation of the local density of relic neutrinos,” JCAP 09, 034 (2017), arXiv:1706.09850 [astro-ph.CO] .
  * Mertsch  _et al._ (2020) P. Mertsch, G. Parimbelli, P. F. de Salas, S. Gariazzo, J. Lesgourgues, and S. Pastor, “Neutrino clustering in the Milky Way and beyond,” JCAP 01, 015 (2020), arXiv:1910.13388 [astro-ph.CO] .
  * Seljak and Zaldarriaga (2000) Uros Seljak and Matias Zaldarriaga, “Lensing induced cluster signatures in cosmic microwave background,” Astrophys. J. 538, 57–64 (2000), arXiv:astro-ph/9907254 .
  * Lewis and King (2006) Antony Lewis and Lindsay King, “Cluster masses from cmb and galaxy weak lensing,” Phys. Rev. D 73, 063006 (2006), arXiv:astro-ph/0512104 .
  * Levy  _et al._ (2023) Kevin Levy, Srinivasan Raghunathan, and Kaustuv Basu, “A Foreground-Immune CMB-Cluster Lensing Estimator,” (2023), arXiv:2305.06326 [astro-ph.CO] .
  * Horowitz  _et al._ (2019) Benjamin Horowitz, Simone Ferraro, and Blake D. Sherwin, “Reconstructing Small Scale Lenses from the Cosmic Microwave Background Temperature Fluctuations,” Mon. Not. Roy. Astron. Soc. 485, 3919–3929 (2019), arXiv:1710.10236 [astro-ph.CO] .
  * Birkinshaw and Gull (1983) M. Birkinshaw and S. F. Gull, “A test for transverse motions of galaxies and clusters,” Nature 302, 315–317 (1983).
  * Hotinli  _et al._ (2019) Selim C. Hotinli, Joel Meyers, Neal Dalal, Andrew H. Jaffe, Matthew C. Johnson, James B. Mertens, Moritz Münchmeyer, Kendrick M. Smith, and Alexander van Engelen, “Transverse Velocities with the Moving Lens Effect,” Phys. Rev. Lett. 123, 061301 (2019), arXiv:1812.03167 [astro-ph.CO] .
  * Hotinli  _et al._ (2021) Selim C. Hotinli, Matthew C. Johnson, and Joel Meyers, “Optimal filters for the moving lens effect,” Phys. Rev. D 103, 043536 (2021), arXiv:2006.03060 [astro-ph.CO] .
  * Correa  _et al._ (2015) Camila A. Correa, J. Stuart B. Wyithe, Joop Schaye, and Alan R. Duffy, “The accretion history of dark matter haloes – III. A physical model for the concentration–mass relation,” Mon. Not. Roy. Astron. Soc. 452, 1217–1232 (2015), arXiv:1502.00391 [astro-ph.CO] .
  * Singh and Ma (2003) Shwetabh Singh and Chung-Pei Ma, “Neutrino clustering in cold dark matter halos : Implications for ultrahigh-energy cosmic rays,” Phys. Rev. D 67, 023506 (2003), arXiv:astro-ph/0208419 .
  * Blas  _et al._ (2011) Diego Blas, Julien Lesgourgues, and Thomas Tram, “The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes,” JCAP 07, 034 (2011), arXiv:1104.2933 [astro-ph.CO] .
  * Lesgourgues and Pastor (2012) Julien Lesgourgues and Sergio Pastor, “Neutrino mass from Cosmology,” Adv. High Energy Phys. 2012, 608515 (2012), arXiv:1212.6154 [hep-ph] .
  * Lewis and Challinor (2006) Antony Lewis and Anthony Challinor, “Weak gravitational lensing of the CMB,” Physics Reports 429, 1–65 (2006), arXiv:astro-ph/0601594 [astro-ph] .
  * Ferraro  _et al._ (2022) Simone Ferraro, Emmanuel Schaan, and Elena Pierpaoli, “Is the Rees-Sciama effect detectable by the next generation of cosmological experiments?” (2022), arXiv:2205.10332 [astro-ph.CO] .
  * Abell  _et al._ (2009) Paul A. Abell _et al._ (LSST Science, LSST Project), “LSST Science Book, Version 2.0,” (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Sheth and Tormen (1999) Ravi K. Sheth and Giuseppe Tormen, “Large scale bias and the peak background split,” Mon. Not. Roy. Astron. Soc. 308, 119 (1999), arXiv:astro-ph/9901122 .
  * Palmese  _et al._ (2020) A. Palmese _et al._ (DES), “Stellar Mass as a Galaxy Cluster Mass Proxy: Application to the Dark Energy Survey redMaPPer Clusters,” Mon. Not. Roy. Astron. Soc. 493, 4591–4606 (2020), arXiv:1903.08813 [astro-ph.CO] .
  * Murata  _et al._ (2018) Ryoma Murata, Takahiro Nishimichi, Masahiro Takada, Hironao Miyatake, Masato Shirasaki, Surhud More, Ryuichi Takahashi, and Ken Osato, “Constraints on the mass-richness relation from the abundance and weak lensing of SDSS clusters,” Astrophys. J. 854, 120 (2018), arXiv:1707.01907 [astro-ph.CO] .
  * Ballardini  _et al._ (2019) Mario Ballardini, William Luke Matthewson, and Roy Maartens, “Constraining primordial non-Gaussianity using two galaxy surveys and CMB lensing,” Mon. Not. Roy. Astron. Soc. 489, 1950–1956 (2019), arXiv:1906.04730 [astro-ph.CO] .
  * Schlegel  _et al._ (2019) David J. Schlegel _et al._ , “Astro2020 APC White Paper: The MegaMapper: a z $>$ 2 Spectroscopic Instrument for the Study of Inflation and Dark Energy,” (2019), arXiv:1907.11171 [astro-ph.IM] .
  * Lewis  _et al._ (2000) Antony Lewis, Anthony Challinor, and Anthony Lasenby, “Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models,” Astrophys. J. 538, 473–476 (2000), arXiv:astro-ph/9911177 [astro-ph] .
  * Hotinli  _et al._ (2021) Selim C. Hotinli, Kendrick M. Smith, Mathew S. Madhavacheril, and Marc Kamionkowski, “Cosmology with the moving lens effect,” Phys. Rev. D 104, 083529 (2021), arXiv:2108.02207 [astro-ph.CO] .
  * Abitbol  _et al._ (2019) Maximilian H. Abitbol _et al._ (Simons Observatory), “The Simons Observatory: Astro2020 Decadal Project Whitepaper,” Bull. Am. Astron. Soc. 51, 147 (2019), arXiv:1907.08284 [astro-ph.IM] .
  * Abazajian  _et al._ (2016) Kevork N. Abazajian _et al._ (CMB-S4), “CMB-S4 Science Book, First Edition,” (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Abazajian  _et al._ (2019) Kevork Abazajian _et al._ , “CMB-S4 Science Case, Reference Design, and Project Plan,” (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Sehgal  _et al._ (2019) Neelima Sehgal _et al._ , “CMB-HD: An Ultra-Deep, High-Resolution Millimeter-Wave Survey Over Half the Sky,” (2019), arXiv:1906.10134 [astro-ph.CO] .
  * Aiola  _et al._ (2022) Simone Aiola _et al._ (CMB-HD), “Snowmass2021 CMB-HD White Paper,” (2022), arXiv:2203.05728 [astro-ph.CO] .
  * Namikawa  _et al._ (2013) Toshiya Namikawa, Duncan Hanson, and Ryuichi Takahashi, “Bias-Hardened CMB Lensing,” Mon. Not. Roy. Astron. Soc. 431, 609–620 (2013), arXiv:1209.0091 [astro-ph.CO] .
  * Millea  _et al._ (2019) Marius Millea, Ethan Anderes, and Benjamin D. Wandelt, “Bayesian delensing of CMB temperature and polarization,” Phys. Rev. D 100, 023509 (2019), arXiv:1708.06753 [astro-ph.CO] .
  * Hadzhiyska  _et al._ (2019) Boryana Hadzhiyska, Blake D. Sherwin, Mathew Madhavacheril, and Simone Ferraro, “Improving Small-Scale CMB Lensing Reconstruction,” Phys. Rev. D 100, 023547 (2019), arXiv:1905.04217 [astro-ph.CO] .
  * Darwish  _et al._ (2023) Omar Darwish, Blake D. Sherwin, Noah Sailer, Emmanuel Schaan, and Simone Ferraro, “Optimizing foreground mitigation for CMB lensing with combined multifrequency and geometric methods,” Phys. Rev. D 107, 043519 (2023), arXiv:2111.00462 [astro-ph.CO] .
  * Sailer  _et al._ (2020) Noah Sailer, Emmanuel Schaan, and Simone Ferraro, “Lower bias, lower noise CMB lensing with foreground-hardened estimators,” Phys. Rev. D 102, 063517 (2020), arXiv:2007.04325 [astro-ph.CO] .
  * Sailer  _et al._ (2023) Noah Sailer, Simone Ferraro, and Emmanuel Schaan, “Foreground-immune CMB lensing reconstruction with polarization,” Phys. Rev. D 107, 023504 (2023), arXiv:2211.03786 [astro-ph.CO] .
  * Parker  _et al._ (2022) Liam Parker, Dongwon Han, Pablo Lemos Portela, and Shirley Ho, “Recovering Galaxy Cluster Convergence from Lensed CMB with Generative Adversarial Networks,” arXiv e-prints , arXiv:2211.08990 (2022), arXiv:2211.08990 [astro-ph.CO] .
  * Foreman  _et al._ (2023) Simon Foreman, Selim C. Hotinli, Mathew S. Madhavacheril, Alexander van Engelen, and Christina D. Kreisch, “Subtracting the kinetic Sunyaev-Zeldovich effect from the cosmic microwave background with surveys of large-scale structure,” Phys. Rev. D 107, 083502 (2023), arXiv:2209.03973 [astro-ph.CO] .
  * Li  _et al._ (2022) Peikai Li, Ipek Ilayda Onur, Scott Dodelson, and Shreyas Chaudhari, “High-Resolution CMB Lensing Reconstruction with Deep Learning,” arXiv e-prints , arXiv:2205.07368 (2022), arXiv:2205.07368 [astro-ph.CO] .
  * Han and Sehgal (2022) Dongwon Han and Neelima Sehgal, “Mitigating foreground bias to the CMB lensing power spectrum for a CMB-HD survey,” Phys. Rev. D 105, 083516 (2022), arXiv:2112.02109 [astro-ph.CO] .
  * Hotinli  _et al._ (2022) Selim C. Hotinli, Joel Meyers, Cynthia Trendafilova, Daniel Green, and Alexander van Engelen, “The benefits of CMB delensing,” JCAP 04, 020 (2022), arXiv:2111.15036 [astro-ph.CO] .
  * Tegmark  _et al._ (2003) Max Tegmark, Angelica de Oliveira-Costa, and Andrew Hamilton, “A high resolution foreground cleaned CMB map from WMAP,” Phys. Rev. D 68, 123523 (2003), arXiv:astro-ph/0302496 .
  * Peirani  _et al._ (2017) S. Peirani, Y. Dubois, M. Volonteri, J. Devriendt, K. Bundy, J. Silk, C. Pichon, S. Kaviraj, R. Gavazzi, and M. Habouzit, “Density profile of dark matter haloes and galaxies in the horizon–agn simulation: the impact of AGN feedback,” Mon. Not. Roy. Astron. Soc. 472, 2153–2169 (2017), arXiv:1611.09922 [astro-ph.GA] .
  * Sorini  _et al._ (2022) Daniele Sorini, Romeel Dave, Weiguang Cui, and Sarah Appleby, “How baryons affect haloes and large-scale structure: a unified picture from the Simba simulation,” Mon. Not. Roy. Astron. Soc. 516, 883–906 (2022), arXiv:2111.13708 [astro-ph.GA] .
  * Chung  _et al._ (2020) Eegene Chung, Simon Foreman, and Alexander van Engelen, “Baryonic effects on CMB lensing and neutrino mass constraints,” Phys. Rev. D 101, 063534 (2020), [Erratum: Phys.Rev.D 102, 109903 (2020)], arXiv:1910.09565 [astro-ph.CO] .
  * Hui  _et al._ (2007) Lam Hui, Enrique Gaztañaga, and Marilena Loverde, “Anisotropic magnification distortion of the 3D galaxy correlation. I. Real space,” Phys. Rev. D 76, 103502 (2007), arXiv:0706.1071 [astro-ph] .
  * Green  _et al._ (2017) Daniel Green, Joel Meyers, and Alexander van Engelen, “CMB Delensing Beyond the B Modes,” JCAP 12, 005 (2017), arXiv:1609.08143 [astro-ph.CO] .