#  A counterexample to the classical Gidas Ni and Nirenberg symmetry result in
the case of Robin boundary conditions

Celentano A., Masiello A. L.*, Paoli G.

###### Abstract

Symmetry properties of solutions to elliptic quasilinear equations have been
widely studied in the context of Dirichlet boundary conditions. We provide a
counterexample to the classical symmetry result by Gidas, Ni and Nirenberg in
the case of Robin boundary conditions with positive boundary parameter
$\beta$.

MSC 2020: 35B06; 35J25; 35B35.

Keywords: Poisson problem; Robin boundary conditions; symmetry breaking.

##  1 Introduction

The task of proving symmetries of solutions to quasilinear or nonlinear PDEs
that reflect the symmetries of the domain has interested many authors. In this
context, the classical result by Gidas, Ni and Nirenberg, contained in the
celebrated paper [11, Theorem 1], is stated as follows.

######  Theorem 1.1 (Gidas-Ni-Nirenberg symmetry result).

Let $f:\operatorname{\mathbb{R}}\to\operatorname{\mathbb{R}}$ be such that
$f=f_{1}+f_{2}$, where $f_{1}$ is a locally Lipschitz function and $f_{2}$ is
non decreasing. Then, any positive solution $u\in C^{2}(\overline{B})$ to the
problem

$\begin{cases}-\Delta u=f(u)&\text{in }B,\\\ u=0&\text{on }\partial
B,\end{cases}$ (1)

where $B$ is a ball of $\operatorname{\mathbb{R}}^{n}$ of radius $R$, has to
be radial and

$\dfrac{\partial u}{\partial r}<0,\quad\text{for }\,0<r<R.$

In order to prove this result the authors make use of the method of moving
planes, first introduced by Aleksandrov in [1], and then applied by Serrin in
[19] in the context of PDEs.

In the case $n=2$, Lions in [13] gives an alternative proof of Theorem 1.1,
that allows considering weaker smoothness assumptions on $f$ and $u$, under
the additional hypothesis $f\geq 0$. The technique used in [13] relies on the
Schwartz symmetrization, the isoperimetric inequality and the Pohozaev
identity. His techniques were applied and generalized in [12] to the context
of the $n-$Laplacian, and, eventually, in [18] to the solutions to $p-$Laplace
equation for any $p$, in any dimension $n\geq 2$.

The result contained in [11] is a milestone in proving symmetry results: in
the linear case see, for instance, [4, 3, 10, 17] and, for the symmetry
results in the case of the $p-$Laplace operator, we refer to [2, 8, 9, 6]. In
all the afore-mentioned papers, it is possible to prove the radiality of
positive solution to (1) either under the regularity hypothesis on $f$ stated
in Theorem 1.1 (using the moving plane method as in [11]) or under the
assumption $f\geq 0$ (with symmetrization techniques as in [13, 18, 12]). For
a sign changing $f$, the Lipschitz continuity property cannot be relaxed to
Hölder continuity, as remarked in [11, Section 2.3]. Indeed, in this case, the
authors find a solution to (1) that is not radially symmetric, see also [7]
for the $p-$Laplacian case.

The aim of the present work is to study the behaviour of the solution to the
Robin problem

$\begin{cases}-\Delta u=f(u)&\text{in }B\\\ \displaystyle{\frac{\partial
u}{\partial\nu}+\beta u=0}&\text{on }\partial B\end{cases}$ (2)

whenever $\beta$ is a positive parameter. Up to our knowledge, in the
literature, there are few results dealing with symmetry properties of the
solutions to differential equations with Robin boundary conditions. For
instance, in [5], the authors consider the following problem

$\begin{cases}-\varepsilon^{2}\Delta u=f(u)-u&\text{in }B\\\
\displaystyle{\varepsilon\frac{\partial u}{\partial\nu}+\beta u=0}&\text{on
}\partial B\end{cases}$ (3)

where $\varepsilon>0,$
$f:\operatorname{\mathbb{R}}\to\operatorname{\mathbb{R}}$ is a continuous
function, of the form $f=f_{1}-f_{2}$ for $t\geq 0$, with $f_{1},f_{2}\geq 0$,
satisfying some structural growth conditions (for the precise details see [5,
Section 1]). Under these assumptions, it is possible to prove the existence of
a positive solution to (3) by making use of the Palais-Smale condition and the
mountain pass Lemma. The authors show that there exists a $\beta_{*}>1$, such
that, for $\beta>\beta_{*}$, the least energy solution has the maximum in the
center of the ball, while, for $\beta\leq\beta_{*}$ and $\varepsilon\to 0$,
the least energy solution has the maximum near the boundary, and,
consequently, it cannot be radially symmetric. For completeness sake, we
recall that the case $\beta=0$, i.e. the Neumann problem, and the case
$\beta=+\infty$, i.e the Dirichlet problem, have been studied, for instance,
respectively, in [14, 15] and [16]. Moreover, these results hold in the more
general case of $\Omega\subset\operatorname{\mathbb{R}}^{n}$ bounded and
smooth, as well as the one proved in [5] for the Robin boundary conditions.

We point out that in [5] the authors do not exhibit an explicit solution to
problem (3). In the present work, in dimension $n\geq 1$, we exhibit an
explicit positive solution to (2) with a right-hand side $f$ satisfying the
set of hypotheses stated in Theorem 1.1, that is not radially symmetric . In
particular, we exhibit a solution to Problem (2): see Theorem 3.1. It follows
that in the context of Robin boundary conditions this set of hypotheses on $f$
is not enough to prove the radiality.

Nevertheless, in Theorem 2.1, we show that in the one dimensional case the
symmetry result for the solution to (2) holds under the standard hypotheses of
Gidas-Ni-Nirenberg and the additional hypothesis $f\geq 0$. On the other hand,
this is not the case for $n\geq 2$, as pointed out in Corollary 3.2:

> In dimension $n\geq 2$, there exists a positive superharmonic function
> $\varphi$ that is a solution to (2) and that is not radially symmetric.

##  2 One dimensional case: the symmetry holds.

We start by analysing the one dimensional case.

######  Theorem 2.1.

Let $R>0$ and let $I=]-R,R[$ be the open ball of radius $R$. Let $u\in
C^{2}([-R,R]])$ be a solution to

$\begin{cases}-u^{\prime\prime}=f(u)&\mbox{in }I,\\\ \dfrac{\partial
u}{\partial\nu}(x)+\beta u(x)=0&\mbox{in }x=\pm R,\end{cases}$

where $\beta>0$. Let us assume that $f$ satisfies the following assumptions:

  1. (i).

$f\geq 0$ in $\operatorname{\mathbb{R}}$, $f$ is not identically zero in
$u(I)$,

  2. (ii).

$f=f_{1}+f_{2}$, where $f_{1}$ is locally Lipschitz in
$\operatorname{\mathbb{R}}$ and $f_{2}$ is non decreasing.

Then, $u(x)=u(-x)$ for all $x\in[-R,R]$. Moreover,

$u^{\prime}(x)<0,\quad x\in[0,R].$

###### Proof.

We divide the proof in two steps. In the first step, we prove that the
function $u$ is strictly positive and, in the second one, we prove that we can
apply the result contained in Theorem 1.1.

Step 1. We start by proving that $u>0$ in $[-R,R]$. Since
$u^{\prime\prime}\leq 0$, $u^{\prime}$ is non increasing in $]-R,R[$, so the
minimum of $u$ on $[-R,R]$ is achieved either in $-R$ or in $R$. Let us denote
by $x_{m}$ the minimum point of $u$ in $[-R,R]$. From the Robin boundary
conditions, we have that

$-\beta u(x_{m})=\dfrac{\partial u}{\partial\nu}(x_{m})\leq 0,$

and, as a consequence, $u\geq 0$ in $[-R,R]$.

Now we want to prove that $u>0$ in $[-R,R]$. By contradiction, we assume that
$u(x_{m})=0$. If $x_{m}=-R$, the Robin boundary conditions imply

$0=\beta u(-R)=-\frac{\partial u}{\partial\nu}(-R)=u^{\prime}(-R)\geq
u^{\prime}(x),\quad\forall x\in I,$

where we have observed that $u^{\prime}$ is non increasing. This implies that
also $u$ is a non increasing function, so $-R$ should be both a minimum and a
maximum. This is not possible, since $u$ should be constant and this
contradicts the hypothesis $f\not\equiv 0$ in $u(I)$. Therefore, we have that
$x_{m}=R$ and, arguing as before, we have

$0=-\beta u(R)=\frac{\partial u}{\partial\nu}(R)=u^{\prime}(R)\leq
u^{\prime}(x),\quad\forall x\in I,$

So $u$ is a non decreasing function, the point $R$ is both a minimum and a
maximum and we get a contradiction as before.

Step 2. We prove now that $u(R)=u(-R)$. Let us assume by contradiction that
$u(R)\neq u(-R)$. Without loss of generality, we can suppose $u(R)<u(-R)$. As
a consequence of Step $1$, the function $u$ is strictly increasing in a
neighborhood of the point $-R$, so, by the continuity of $u$, there exists
$y\in I$ such that $u(y)=u(-R)$. Therefore, the following quantity is well
defined

$\lambda:=\inf\\{t\in I\,:\,u(t)=u(-R)\\}$

and $\lambda>-R$. Moreover, the continuity of $u$ also implies
$u(\lambda)=u(-R)$ and $u(x)>u(-R)$ in $(-R,\lambda)$.

We define now the function $v:=u-u(-R)$, that is a positive solution to

$\begin{cases}-\Delta v=\tilde{f}(v)&\text{in }(-R,\lambda),\\\ v=0&\text{in
}x=-R,\;x=\lambda,\end{cases}$

where $\tilde{f}(v)=f(v+u(R))$. So we can use Theorem 1.1 in the interval
$(-R,\lambda)$. We have that $u$ is symmetric with respect to the line
$x=2^{-1}(\lambda-R)$ and, as a consequence, we get

$\dfrac{du}{d\nu}(-R)=-\dfrac{du}{dx}(-R)=\dfrac{du}{dx}(\lambda)<0.$ (4)

Using (4) and the fact that $u^{\prime}$ is non increasing, we obtain

$\beta
u(R)=-\dfrac{du}{d\nu}(R)=-\dfrac{du}{dx}(R)\geq-\dfrac{du}{dx}(\lambda)=-\dfrac{du}{d\nu}(-R)=\beta
u(-R),$ (5)

and, therefore,

$u(R)\geq u(-R),$ (6)

that is a contradiction.

From Step $1$ and $2$, we have that the funtion $v=u-u(R)$ is a positive
solution to

$\begin{cases}-\Delta v=\tilde{f}(v)&\text{in }I,\\\ v=0&\text{in }x=\pm
R.\end{cases}$

So we can conclude by using Theorem 1.1.

∎

We do not know if the hypotheses $i)-ii)$ on the function $f$ are the optimal
ones to obtain the symmetry result. Nevertheless, we will show in Remark 3.1,
in the next session, that the assumption $f\geq 0$ cannot be removed.

##  3 Counterexample: symmetry breaking in dimension $n\geq 2$

In the following, we will denote by $|\cdot|$ the Euclidean norm in
$\mathbb{R}^{n}$.

######  Theorem 3.1.

Let $B_{R}\subset\mathbb{R}^{n}$, $n\geq 2$, be the ball centered at the
origin with radius $R$, let $\beta$ be a positive constant and let $x_{0}\neq
0$ in $B_{R}$. Then, there exists a the positive function $\varphi\in
C^{\infty}(\overline{B_{R}})$ that is a non radial function in $B_{R}$ (i.e.
$\varphi(x)\neq\varphi({\left|x\right|})$) and it is a solution to

$\begin{cases}-\Delta\varphi=f(\varphi)&\text{in }B_{R},\\\
\dfrac{\partial\varphi}{\partial\nu}+\beta\varphi=0&\text{on }\partial
B_{R},\end{cases}$ (7)

where

$f(t)=c_{1}t\left[c_{2}t^{\frac{1}{\beta R}}+c_{3}t^{\frac{2}{\beta
R}}\right],$ (8)

with $c_{1},c_{3}>0$, $c_{2}\in\operatorname{\mathbb{R}}$ defined as follows:

$c_{1}=2\beta R,\quad c_{2}=-2(\beta R+1)+n,\quad c_{3}=2(\beta
R+1)\alpha^{2},\quad\alpha^{2}=R^{2}-{\left|x_{0}\right|}^{2}.$

$O$$a$$x_{0}$$x$$R$$r$$\theta$ Figure 1: Construction of the function
$\varphi$

###### Proof.

We define the following quantities (see Figure 1):

  * •

$a:={\left|x_{0}\right|}$,

  * •

$\alpha^{2}:=R^{2}-a^{2}>0$.

We show the existence of a positive function
$\varphi(x)=\varphi({\left|x-x_{0}\right|})=\varphi(r)$ such that

$\dfrac{\partial\varphi}{\partial\nu}+\beta\varphi=0\qquad\text{on}\;\;\partial
B_{R}$ (9)

where $\nu$ is the unit outer normal to $\partial B_{R}$. Let us fix
$x\in\partial B_{R}$. Being

$\nabla(\varphi(x))=\varphi^{\prime}(r)\>\dfrac{x-x_{0}}{r},\qquad\nu(x)=\dfrac{x}{R},$

the Robin boundary conditions (9) becomes

$\varphi^{\prime}(r)\;\dfrac{x-x_{0}}{r}\cdot\dfrac{x}{R}+\beta\varphi(r)=0.$
(10)

Denoting, now, by $\theta$ the angle between the vectors $x_{0}$ and $x$, we
have the following relation

$\cos(\theta)=\dfrac{R^{2}+a^{2}-r^{2}}{2aR}.$ (11)

So, from (10) and (11), recalling that $\alpha^{2}=R^{2}-a^{2}$, we have

$\frac{\varphi^{\prime}(r)}{rR}\left(R^{2}-Ra\cos(\theta)\right)+\beta\varphi(r)=\frac{\varphi^{\prime}(r)}{2rR}\left(r^{2}+\alpha^{2}\right)+\beta\varphi(r)=0,$

and, therefore,

$\dfrac{\varphi^{\prime}(r)}{\varphi(r)}=-\left(2\beta
R\right)\dfrac{r}{r^{2}+\alpha^{2}}.$ (12)

Integrating (12), we get

$\varphi(r)=\dfrac{c}{(r^{2}+\alpha^{2})^{\beta R}}.$ (13)

If we choose $c=1$, we have

$\displaystyle\varphi^{\prime}(r)$ $\displaystyle=-\dfrac{2\beta
Rr}{(r^{2}+\alpha^{2})^{\beta R+1}}$ $\displaystyle\varphi^{\prime\prime}(r)$
$\displaystyle=-(-\beta R-1)\dfrac{4\beta Rr^{2}}{(r^{2}+\alpha^{2})^{\beta
R+2}}-\dfrac{2\beta R}{(r^{2}+\alpha^{2})^{\beta R+1}}.$

and, consequently,

$\displaystyle-\Delta\varphi$
$\displaystyle=-\varphi^{\prime\prime}(r)-\dfrac{n-1}{r}\varphi^{\prime}(r)=$
$\displaystyle=\dfrac{2\beta R}{(r^{2}+\alpha^{2})^{\beta R+1}}\left[n-2(\beta
R+1)+(\beta R+1)\dfrac{2\alpha^{2}}{r^{2}+\alpha^{2}}\right]$
$\displaystyle=2\beta R\varphi(r)^{\frac{1}{\beta R}+1}\left[n-2(\beta
R+1)+2(\beta R+1)\alpha^{2}\varphi(r)^{\frac{1}{\beta R}}\right]$
$\displaystyle=f(\varphi(r)),$

where $f$ is the function defined in (8). So, we have proved the desired
claim, since we have found a non radial function of the form
$\varphi(x)=\varphi({\left|x-x_{0}\right|})=\varphi(r)$, defined in (13), that
satisfies (7). ∎

As a consequence of Theorem 3.1, we obtain the following Corollary.

######  Corollary 3.2.

Let $n\geq 2$. There exists a positive superharmonic function $\varphi$ that
is a solution to (2) and that is not radially symmetric.

###### Proof.

In the case $n=2$, the right-hand side of (8) becomes

$f(t)=4\beta Rt\left(-\beta Rt^{\frac{1}{\beta R}}+\alpha^{2}(1+\beta
R)t^{\frac{2}{\beta R}}\right).$

We notice that

$f(t)\geq 0,\qquad{\rm if}\;t\geq\left(\dfrac{\beta R}{\alpha^{2}(1+\beta
R)}\right)^{\beta R},$ (14)

so the function $f\circ\varphi$ is positive if $\varphi$ satisfies (14) for
all $x\in B_{R}$, and this follows by imposing the following geometric
constraint

$\beta\leq\dfrac{R-a}{R(R+a)}.$

If $n\geq 3$, we can choose the constant $c_{2}\geq 0$, by imposing the
condition

$\beta\leq\dfrac{n-2}{2R}$ (15)

and, under these assumptions, we have that $f(t)\geq 0$ for $t\geq 0$.

Therefore, we can see that, by imposing the geometrical constraints (14) and
(15) respectively for $n=2$ and $n\geq 3$, the function $\varphi$ defined in
(13) is an example of positive superharmonic function, that is non radial and
that satisfies (7). ∎

We conclude with a remark for the one dimensional case.

######  Remark 3.1.

The function $\varphi$ defined in Theorem 3.1, in the case $n=1$, satisfies
the problem

$\begin{cases}-\varphi^{\prime\prime}=f(\varphi)&\text{in }(-R,+R),\\\
\dfrac{\partial\varphi}{\partial\nu}+\beta\varphi=0&\text{ in}\;x=\pm
R.\end{cases}$ (16)

We note that $\varphi\in C^{\infty}\left([-R,R]\right)$ and $f$ is a locally
Lipschitz function, but, $f$ does not satisfy the hypothesis $i)$, that is the
positiveness. Indeed, by straightforward computations, we obtain that
$f\circ\varphi$ is a sign changing function in $\varphi([-R,+R])$ for every
$\beta>0$ and $R>a>0$.

## Statements and Declarations

Funding: All the authors are partially supported by Gruppo Nazionale per
l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of
Istituto Nazionale di Alta Matematica (INdAM). Gloria Paoli is supported by
the Alexander von Humboldt Foundation through an Alexander von Humboldt
research fellowship.

Conflict of interst: There is no conflict of interest to disclose

## References

  * [1] A. D. Aleksandrov.  Uniqueness theorems for surfaces in the large. V.  Amer. Math. Soc. Transl. (2), 21:412–416, 1962.
  * [2] M. Badiale and E. Nabana.  A note on radiality of solutions of $p$-Laplacian equation.  Appl. Anal., 52(1-4):35–43, 1994.
  * [3] H. Berestycki and L. Nirenberg.  Monotonicity, symmetry and antisymmetry of solutions of semilinear elliptic equations.  J. Geom. Phys., 5(2):237–275, 1988.
  * [4] H. Berestycki and L. Nirenberg.  On the method of moving planes and the sliding method.  Bol. Soc. Brasil. Mat. (N.S.), 22(1):1–37, 1991.
  * [5] H. Berestycki and J. Wei.  On singular perturbation problems with Robin boundary condition.  Ann. Sc. Norm. Super. Pisa Cl. Sci. (5), 2(1):199–230, 2003.
  * [6] F. Brock.  Radial symmetry for nonnegative solutions of semilinear elliptic equations involving the $p$-Laplacian.  In Progress in partial differential equations, Vol. 1 (Pont-à-Mousson, 1997), volume 383 of Pitman Res. Notes Math. Ser., pages 46–57. Longman, Harlow, 1998.
  * [7] F. Brock.  Continuous rearrangement and symmetry of solutions of elliptic problems.  Proc. Indian Acad. Sci. Math. Sci., 110(2):157–204, 2000.
  * [8] L. Damascelli and F. Pacella.  Monotonicity and symmetry results for $p$-Laplace equations and applications.  Adv. Differential Equations, 5(7-9):1179–1200, 2000.
  * [9] L. Damascelli and B. Sciunzi.  Regularity, monotonicity and symmetry of positive solutions of $m$-Laplace equations.  J. Differential Equations, 206(2):483–515, 2004.
  * [10] L. E. Fraenkel.  An introduction to maximum principles and symmetry in elliptic problems, volume 128 of Cambridge Tracts in Mathematics.  Cambridge University Press, Cambridge, 2000.
  * [11] B. Gidas, W.-M. Ni, and L. Nirenberg.  Symmetry and related properties via the maximum principle.  Comm. Math. Phys., 68(3):209–243, 1979.
  * [12] S. Kesavan and F. Pacella.  Symmetry of positive solutions of a quasilinear elliptic equation via isoperimetric inequalities.  Appl. Anal., 54(1-2):27–37, 1994.
  * [13] P.-L. Lions.  Two geometrical properties of solutions of semilinear problems.  Applicable Anal., 12(4):267–272, 1981.
  * [14] W.-M. Ni and I. Takagi.  On the shape of least-energy solutions to a semilinear Neumann problem.  Comm. Pure Appl. Math., 44(7):819–851, 1991.
  * [15] W.-M. Ni and I. Takagi.  Locating the peaks of least-energy solutions to a semilinear Neumann problem.  Duke Math. J., 70(2):247–281, 1993.
  * [16] W.-M. Ni and J. Wei.  On the location and profile of spike-layer solutions to singularly perturbed semilinear Dirichlet problems.  Comm. Pure Appl. Math., 48(7):731–768, 1995.
  * [17] E. Rosset.  An approximate Gidas-Ni-Nirenberg theorem.  Math. Methods Appl. Sci., 17(13):1045–1052, 1994.
  * [18] J. Serra.  Radial symmetry of solutions to diffusion equations with discontinuous nonlinearities.  J. Differential Equations, 254(4):1893–1902, 2013.
  * [19] J. Serrin.  A symmetry problem in potential theory.  Arch. Rational Mech. Anal., 43:304–318, 1971.

E-mail address, A. Celentano<EMAIL_ADDRESS>
E-mail address, A.L. Masiello* (corresponding author):
<EMAIL_ADDRESS>
Dipartimento di Matematica e Applicazioni “R. Caccioppoli”, Università degli
studi di Napoli Federico II, Via Cintia, Complesso Universitario Monte S.
Angelo, 80126 Napoli, Italy.

E-mail address, G. Paoli<EMAIL_ADDRESS>
Department Mathematik, Chair in Dynamics, Control, Machine Learning and
Numerics (Alexander von Humboldt-Professorship), Cauerstr. 11, 91058 Erlangen,
Germany.


# Broad-Range Directional Detection of Light Dark Matter in Cryogenic Ice

Nora Taufertshöfer  Molecular Foundry Division, Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA  Materials Sciences Division,
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
Institute for Theoretical Physics, University of Würzburg, Am Hubland, D-97074
Würzburg, Germany  Maurice Garcia-Sciveres  Physics Division, Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA  Sinéad M.
Griffin<EMAIL_ADDRESS>Molecular Foundry Division, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA  Materials Sciences
Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720,
USA

###### Abstract

We propose hexagonal ice (H2O) as a new target for light dark matter (DM)
direct detection. Ice, a polar material, is suitable for single phonon
detection through DM scattering for which we consider light dark photon and
light scalar mediator models. We report a rate sensitivity down to a DM mass
of $\sim$ keV, constituting a broader mass range than other promising
candidates. We find better sensitivity for near-term experimental thresholds
from the presence of high-frequency phonons. These advantages, and ice’s
availability, make it highly promising for single-phonon detection.

New proposals for the direct detection of dark matter (DM) have steadily
pushed into the sub-GeV ‘light’ DM rangeEssig  _et al._ (2012). In particular,
for DM masses in the kev-GeV range, freeze-in DM scattering is well motivated
through both a light scalar mediator and a dark photon mediatorHall  _et al._
(2010); Bernal  _et al._ (2017). The smaller kinetic energy for this lower
mass DM motivates the exploration of low-energy excitations in condensed
matter systems as possible direct detection channelsKahn and Lin (2022).
Recent proposals for DM detection in crystal targets include quasiparticle
excitations ranging from phonon and electron scattering in conventional
semiconductorsEssig  _et al._ (2016); Knapen  _et al._ (2018) to more exotic
quantum materials such as Dirac semimetalsHochberg  _et al._ (2018), and
axions in topological insulators and multiferroicsMarsh  _et al._ (2019);
Røising  _et al._ (2021).

Of these proposals for DM direct detection via quasiparticle scattering,
single phonon production presents a promising detection scheme for several
reasonsKnapen  _et al._ (2018); Griffin  _et al._ (2018). Typical optical
phonons in crystals range from 10-100 meV, making them kinematically matched
with sub-MeV DM. Dark photon mediators can mix with standard model photons
resulting in a coupling to optical phonon modes in polar semiconductors.
Furthermore, common polar semiconductors (e.g. GaAs) have sizeable bandgaps
resulting in small screening. The intrinsic anisotropy of polar semiconductors
allows for a directional response giving signal discrimination based on
daily/annual modulation of the DM windDrukier  _et al._ (1986). Importantly,
high-quality detector-grade crystals of a range of polar semiconductors are
readily available owing to their use in microelectronics and quantum computing
(e.g. GaAsKnapen  _et al._ (2018), SiCGriffin  _et al._ (2021)), and phonon
sensing technologies are steadily improving towards the low sensitivities
required for detectionPyle  _et al._ (2015).

With the next-generation of low-mass DM experiments in planning, it remains to
identify the best possible target material for each detection scheme. Previous
work addressed this at the DM-target interaction level by performing a
comparison of DM-phonon reach, isolating materials-specific factors that
optimize for low-mass sensitivity and greatest cross sectionGriffin  _et al._
(2020a). Beyond this, other materials factors for target selection should also
be considered that take into account current synthesis/fabrication
capabilities and sensing constraints for near-term experiments. While
experimental thresholds for phonon sensing continue to improve towards the
single phonon excitation regime ($\sim$10-100 meV), near-term experiments with
larger thresholds could access higher energy responses. For instance, beyond
the single phonon excitation regime, multiphonon responses can be accessed for
DM scattering above the single phonon spectra though at a suppressed
rateCampbell-Deem  _et al._ (2020). However, an alternative approach to this
is to search for single phonon and other quasiparticle responses that lie
within the current sensing thresholds (>100 meV).

In this Letter, we propose cryogenic water ice (H2O) as a particularly
appealing target material for direct detection of low-mass DM from single-
phonon excitations. We find that the mixed nature of bonding in H2O results in
a broad range of phonon frequencies suitable for both lower-mass detection,
and for higher frequency single phonon excitations that are within the range
of current sensing thresholds. Moreover, with the ubiquity and long history of
fabrication and characterization of solid H2O, it can be readily fabricated as
large single crystals comprising earth-abundant (inexpensive) elements.

This paper is laid out as follows; we summarize the calculation of spin-
independent DM scattering for single-phonons and the features of ice H2O
(specifically its XIh polymorph) of relevance to DM-phonon interactions. We
then present our calculations of the reach for a light dark photon mediator
and a light scalar mediator with comparisons to best-performing targets for
these DM models, and their calculated directional dependence. Finally we
discuss other material limiting considerations for target materials and
prospects for its use for near-term and next-generation DM detection.

First, we briefly summarize the formalism for calculating spin-independent DM
scattering processes including a light dark photon mediator and a light scalar
mediatorTrickle  _et al._ (2020). For a detector with target density
$\rho_{\text{T}}$ the total rate per target mass is given by

$\displaystyle R=\frac{1}{\rho_{\text{T}}}\frac{\rho_{\chi}}{m_{\chi}}\int
d^{3}vf_{\chi}(\mathbf{v})\Gamma(\mathbf{v})$ (1)

where $\rho_{\chi}$ is the local DM energy density, $m_{\chi}$ the DM mass,
and $f_{\chi}(\mathbf{v})$ is the DM velocity distribution in the lab frame
which is modelled with a boosted Maxwell-Boltzmann distribution 111Due to
Earth’s rotation the orientation of the incoming DM particles with regard to
the target crystal changes over a day, which leads to a directional modulation
of the rate. We choose the z-axis in the crystal frame to be parallel to
Earth’s velocity at daytime $t=0$Lewin and Smith (1996).. $\Gamma(\mathbf{v})$
is the scattering rate per non-relativistic DM particle which is defined by an
integral over the momentum transfer
$\mathbf{q}=\mathbf{p}-\mathbf{p}^{\prime}$ with initial DM momentum
$\mathbf{p}=m_{\chi}\mathbf{v}$:

$\displaystyle\Gamma(\mathbf{v})=\frac{\pi\bar{\sigma}}{\mu^{2}}\int\frac{d^{3}q}{(2\pi)^{3}}(q_{0}/q)^{4}S(\mathbf{q},\omega_{\mathbf{q}})$
(2)

Here, $\bar{\sigma}$ is a model-dependent reference cross section and we
define $\bar{\sigma}:=\frac{\mu^{2}}{\pi}\overline{|\mathcal{M}(q_{0})|^{2}}$
with $\mathcal{M}$ the target-independent $2\rightarrow 2$ scattering matrix
element (see Supplemental Material). $\mu$ is the reduced mass of an electron
or nucleon and DM particle for light-dark-photon- or light-scalar-mediated
scattering, respectively. $q_{0}$ is a reference momentum and
$S(\mathbf{q},\omega_{\mathbf{q}})$ the target-dependent dynamic structure
factor.

For single phonon excitations we recall the dynamical structure factor:

$\displaystyle\begin{split}&S(\mathbf{q},\omega_{\mathbf{q}})=\\\
&\frac{\pi}{\Omega}\sum_{\nu}\frac{1}{\omega_{\nu,\mathbf{k}}}\Big{|}\sum_{j}\frac{e^{-W_{j}(\mathbf{q})}}{\sqrt{m_{j}}}\,e^{i\mathbf{G}\cdot\mathbf{x}^{0}_{j}}\,(\mathbf{Y}_{j}\cdot\bm{\epsilon}^{\ast}_{\nu,j,\mathbf{k}})\Big{|}^{2}\delta(\omega_{\mathbf{q}}-\omega_{\nu,\mathbf{k}})\end{split}$
(3)

where $\omega_{\mathbf{q}}=\mathbf{q}\cdot\mathbf{v}-\frac{q^{2}}{2m_{\chi}}$
explicitly depends on the energy deposition. The sums run over all phonon
branches $\nu$ and over all ions $j$ in the primitive cell. The ionic masses
are $m_{j}$ with equilibrium positions $\mathbf{x}^{0}_{j}$ and $\Omega$ is
the primitive cell volume. $\mathbf{\epsilon}_{\nu,j,\mathbf{k}}$ defines the
phonon polarization vectors. $\mathbf{G}$ is the reciprocal lattice vector
that satisfies $\mathbf{G}=\mathbf{q}-\mathbf{k}$ with $\mathbf{k}$ within the
first Brillouin zone. DM couplings come in with the model-specific
$\mathbf{Y}_{j}$ terms defined below. In the continuum limit for $\mathbf{k}$,
the Debye-Waller factor is given by:

$\displaystyle
W_{j}(\mathbf{q})=\frac{\Omega}{4m_{j}}\sum_{\nu}\int_{1\text{BZ}}\frac{d^{3}k}{(2\pi)^{3}}\frac{|\mathbf{q}\cdot\bm{\epsilon}_{\nu,j,\mathbf{k}}|^{2}}{\omega_{\nu,\mathbf{k}}}.$
(4)

The materials-dependent quantities, namely the phonon frequencies
$\omega_{\nu,\mathbf{k}}$ for branch $\nu$ and momentum $\mathbf{k}$, and the
phonon polarization vectors $\mathbf{\epsilon}_{\nu,j,\mathbf{k}}$, can be
calculated using first-principles methods such as Density Functional
TheoryGeorge  _et al._ (2020).

Finally, we consider two well-motivated models of light DM, namely a
kinematically mixed dark photon and a light scalar mediatorKnapen  _et al._
(2018). In the former case, the standard model photon kinematically mixes with
the dark photon owing to their same quantum numbers, resulting in
‘millicharged’ DM under $U(1)$ of the standard model. In the $q\rightarrow 0$
limit valid for light DM scatteringTrickle  _et al._ (2020), $\mathbf{Y}_{j}$
is given by

$\displaystyle\mathbf{Y}_{j,\text{dark
photon}}=-\frac{q^{2}}{\mathbf{q}\cdot\varepsilon_{\infty}\cdot\mathbf{q}}(\mathbf{q}\cdot\mathbf{Z}^{\ast}_{j})$
(5)

where $\varepsilon_{\infty}$ is the high-frequency dielectric tensor and
$\mathbf{Z}^{\ast}_{j}$ is the Born effective charge tensor of the $j$th atom
in the unit cell, both of which can be calculated using first-principles
methods.

We also consider DM that only couples to nucleons via a light scalar mediator
with identical coupling for proton and neutron. The now scalar-valued
$\mathbf{Y}_{j}$ is given by

$\displaystyle\mathbf{Y}_{j,\text{hadrophilic scalar}}=A_{j}F_{N_{j}}(q)$ (6)

with $A_{j}$ the atomic mass number and $F_{N_{j}}$ an isotropic nuclear form
factor of the $j$th atom, respectively. The derivation of (6) is valid for the
$q\rightarrow 0$ limit, which is fulfilled for the single phonon excitation
regime. Further, for this type of scattering $F_{N_{j}}$ can be set to one.

Fig. 1: (a) Primitive unit cell of hexagonal ice XIh and its first Brillouin
zone with high-symmetry points labelled. (b) Calculated phonon band structure
for hexagonal ice XIh. The cartoons on the right hand side show the movement
of ions for three selected phonon modes representative of the types of optical
phonon modes in ice: translational modes (bottom), libration/bending (middle),
and stretch (top).

The pervasiveness of liquid and solid H2O has generated extensive studies of
the phase diagram of ice. So far, twenty polymorphs of crystalline H2O have
been identified under various conditions of temperature and pressure, with the
most recent, ice XIX reported in 2021Hansen (2021). Common to all of these is
the local atomic-scale arrangement of hydrogen atoms in bonded units that
fulfil the ‘Bernal-Fowler ice rules’: each oxygen forms a covalent bond with
two hydrogens and a weaker van-der-Waals bond with two other hydrogensBernal
and Fowler (1933). Since each oxygen is also shared between two hydrogen
atoms, this results in an average oxygen bonding environment comprising ‘two-
in’ (strong) and ‘two-out’ (weak) bonds with hydrogen. Such an atomic
arrangement is geometrically frustrated on a tetrahedral latticeGriffin and
Spaldin (2017), resulting in disordered relative arrangements of the strong
and weak bonding networks known as proton disorder. Common ice (Ih) forms such
a hydrogen-disordered network of H-O-H bonds at ambient pressure, giving rise
to a residual entropy as first predicted by PaulingPauling (1935).

However, spontaneous ordering of the hydrogen networks can occur on cooling,
or by the introduction of a dopant to overcome the kinetic barrier to
ordering. For example KOH doping of Ih results in the formation of a fully
ordered phase of hexagonal ice, XIh below $\sim$ 72 KLeadbetter  _et al._
(1985); Tajima  _et al._ (1982). In fact, the hydrogen ordering in the XIh
structure gives it a net dipole moment, making it potentially
ferroelectricParkkinen  _et al._ (2014). Since XIh is the stable form of ice
at ambient pressures, and can be synthesized in its ordered form, we focus on
this polymorph for the remaining discussion, noting, however, that other
structures of ice should have comparable DM reach for single phonon-based
detection owing to their similar bonding networks. We also considered heavy
ice D2O XIh but we did not find a significant difference to H2O (results
reported in the Supplemental Material).

We now present our DFT-calculated structural and phonon properties of ice XIh.
The full calculation details can be found in the Supplemental Material. With
its combination of covalent bonding, hydrogen bonding, and van-der-Waals
bonding, the challenge of treating ice H2O with DFT functionals has been
explored extensivelySantra  _et al._ (2013). We benchmark our choice of
exchange correlation functional for accurate structural and vibrational
properties. Consistent with previous workRaza  _et al._ (2011), we find that
the van-der-Waals corrected nonlocal functional (optPBE-vdw) of Klimeš et al.
performs bestKlimeš  _et al._ (2009, 2011), resulting in the lattice constants
$a=4.470$ Å and $c=7.212$ Å. This corresponds to 0.6$\%$ and 1.5$\%$ deviation
from experimental measurements at T = 2 KFortes (2018).

In Fig. 1(b) we plot the calculated phonon dispersion for ice XIh. All
polymorphs of ice exhibit a large range of phonon frequencies as a result of
the varied nature of bonding in H2O crystalsCherubini  _et al._ (2021); Yuan
_et al._ (2017, 2022). The low-energy optical phonon range is made up of
‘translational’ modes where individual H2O molecules behave as atom-like
clusters, and vibrate with respect to one another. Such low-frequency phonon
modes are common in molecular crystals such as H2O where molecular units are
weakly bonded to each other to form the crystal networkBrown-Altvater  _et
al._ (2016). The next highest energy set of phonon modes corresponds to
libration and bending of H-O bonds, with the highest-frequency range caused by
the stretching of the O-H bond. As the lightest element in the periodic table,
hydrogen sets an upper limit on the range of these high-frequency stretch
modes that can be found in a material.

Fig. 2: Projected reach for light-dark-photon-mediated DM scattering via
single phonon excitations with 1kg-year exposure and no background. Solid,
long-dashed and dashed lines refer to a 1 meV, 20 meV and 100 meV detector
threshold. For comparison, additional reach data Griffin  _et al._ (2020a) is
shown. Stellar constraints from red giants (RG) and the horizontal branch (HB)
are taken from Chang  _et al._ (2021), the freeze-in as referenced in Griffin
_et al._ (2020a).

We present our calculated reach for light-dark-photon-mediated scattering of
single phonons in ice XIh in Fig. 2, including comparisons to top-performing
targets previously studiedGriffin  _et al._ (2020a). All reach calculations
were done using the code PhonoDarkTrickle  _et al._ (2022); Trickle (2020).
The reach of ice XIh extends further into both the low-mass DM range (well
into the constrained regions) and has better reach in the high-mass DM range
than any other previously suggested single target material. To understand the
exceptionally broadband sensitivity of ice XIh, we evaluate its performance
with respect to our previously suggested quality factors for dark photon
mediators in the high- and low-frequency phonon rangesGriffin  _et al._
(2020a). The lowest DM mass accessible is determined by the energy of the
lowest-frequency optical mode where
$m_{\chi}\sim\frac{1}{3}\omega_{0}^{min}\times 10^{6}$, suggesting that
materials with low-lying optical phonons such as CsI ($\omega_{0}^{min}\sim 7$
meV) have the best low-mass reach. In our case, the van-der-Waals bonded
molecular units also have very low lying optical modes ($\omega_{0}^{min}\sim
6$ meV), making ice XIh optimal for the low-mass DM range. We expect all
similar molecular crystals with such low-frequency translational phonon modes
to have competitive reach for low-mass DM.

For the high-mass range, a quality factor, $Q$, is defined asGriffin  _et al._
(2020a):

$Q\equiv\frac{Z^{\ast
2}}{A_{1}A_{2}\varepsilon_{\infty}^{2}}\left(\frac{\text{meV}}{\omega_{LO}}\right)$
(7)

where $Z^{\ast}$ are the Born effective charges, $A_{1}$, $A_{2}$ are the
atomic masses, $\varepsilon_{\infty}$ is the high-frequency dielectric
constant, and $\omega_{LO}$ is the longitudinal optical phonon frequency. A
high $Q$ corresponds to better reach in the high-mass DM regime. We calculate
each of these values for ice XIh using DFT, with the full set of results given
in the Supplemental Material. The best performing materials previously
proposed include LiF with $Q=270\times 10^{-7}$ and SiO2 with $Q=200\times
10^{-7}$. For ice XIh, we find $Q=533\times 10^{-7}$ for its highest-frequency
stretch modes (>375 meV), steadily increasing for lower-frequency
$\omega_{LO}$ clusters. The substantially enhanced $Q$ of ice XIh in the high-
mass range is due to the low masses of hydrogen and oxygen. This optimized
reach is dominated by hydrogen’s extremely low mass, despite these small
masses also resulting in higher frequency phonon modes. However, this quality
factor does not take the detection threshold into account which is especially
important for near-term experiments.

We conclude that for optimal high-mass DM single-phonon-based detection with
dark photon mediators, chemistries that include light elements such as H and
He, and especially H-H bonds, will give the best reach in the high-mass range
owing to their high-frequency single phonon modes. In the low-mass range, we
find that both compounds containing heavy elements and molecular crystals that
possess low-frequency optical phonon modes will give the best sensitivity.
With ice combining both of these properties (molecular units and H-H bonds),
it displays better sensitivity for dark photon mediators than previously
suggested targets across a broad range of DM masses with CsI being the only
exception for a very small low-mass range (Fig. 2).

The reach for light-hadrophilic-scalar-mediated scattering of single phonons
in ice XIh is plotted in Fig. 3 for several detector thresholds, and compared
to previous top candidates from Ref. Griffin  _et al._ (2020a). For the light
scalar mediator, we also find that ice XIh has broad-band sensitivity for both
low- and high-mass DM. Looking at the reach curves, we find a competition
between which of diamond or ice is favorable for different mass ranges and
thresholds. For all thresholds considered, we find that ice has sensitivity to
lower masses than diamond, which become comparable by a threshold of 100 meV.
In the low mass range, the lowest DM mass sensitivity is generally governed by
$\omega_{min}$/$c_{s}$ where $\omega_{min}$ is the detector threshold and
$c_{s}$ is the material’s speed of sound. Therefore, materials with higher
$c_{s}$ such as diamond will have the best sensitivity to lower-mass DM given
low detector thresholds. We calculated the directional averaged speed of sound
of ice XIh to be 4376 ms-1 (details given in the Supplemental Material), which
is a factor of three smaller than in diamond. However, for materials with
smaller $c_{s}$ the acoustic phonons can be kinematically inaccessible,
especially as the threshold increases to higher energies. When this happens,
the lowest DM sensitivity is then set by the lowest optical phonon available.
For ice XIh we find that the minimum reachable DM mass from kinematic
considerations is $\sim$ 29 keV for a 1 meV threshold, while the lowest
optical phonon is 6 meV. This manifests as a change in the slope as seen in
Fig. 3 where optical phonons, rather than acoustic phonons, dominate the low
mass response.

Fig. 3: Projected reach for light-scalar-mediated DM scattering via single
phonon excitations with 1kg-year exposure and no background. Solid, long-
dashed, dashed and dotted lines refer to a 1 meV, 20 meV, 100 meV and 400 meV
detector threshold. Additional reach data is taken from Griffin  _et al._
(2020a), light-colored lines show the expected curve progression for scale
ranges not shown in the reference data.

In the high-mass, high-threshold range (detector threshold of 400 meV), we
only include ice as it is the only material with single phonons above
$\sim$160 meV. We note that the definition of $\bar{\sigma}_{n}\propto
m_{\chi}^{-4}$ leads to the generally better reach at high massesGriffin  _et
al._ (2020b).

Fig. 4: Daily modulation for H2O XIh ice with (a) a light dark photon and (b)
a light scalar mediator for various $m_{\chi}$.

We finally calculate the daily modulation rate assuming that the z-axis of the
target crystal is aligned with the DM wind at $t=0$. Each rate is normalised
to the average rate $\bar{R}$ over one day. For a dark photon mediator we find
a strong modulation for a mass range of $m_{\chi}=10^{4}-10^{5}$ eV (Fig. 4
(a)). We also find a significant daily modulation for the scalar mediator case
(Fig. 4 (b)), especially in the low 10 keV range. Finally, we extract the
phonon lifetimes $\tau$ from measured phonon linewidths $\Gamma_{\text{ph}}$
of ice XICherubini  _et al._ (2021); Lory  _et al._ (2017) from
$\tau=(\pi\Gamma_{\text{ph}})^{-1}$ giving a lifetime range of 0.5 - 2.1 ps
(see Supplemental Material), consistent with other systems with high-frequency
optical phononsKang  _et al._ (2008).

Our results suggest that ice H2O is a competitive candidate target for single-
phonon based light DM detection both in the low-mass and high-mass range (here
we considered ice XIh, but our conclusions are generally applicable to other
ice polymorphs). Molecular units provide ultra-low-mass optical phonons (down
to 6 meV) giving excellent reach in the low-mass DM range for dark photon
mediators. We expect similar reach to low masses for polar semiconductors
containing heavy ions (as e.g. CsI) or molecular crystals where the relative
effective masses of the ions/molecules result in low-frequency optical
phonons. For near-term experiments, we find that the broad range of single
phonons in ice up to 400 meV gives it enhanced reach for higher experimental
thresholds compared to alternative proposals such as multiphonon excitations.
However, the lower density of ice due to its van-der-Waals bonding results in
a reduced reach compared to other solid-state targets; this can be addressed
by the growth of large single crystals of ice which is already well
exploredKhusnatdinov and Petrenko (1996). Finally, while ice H2O is cheap,
earth-abundant, and has been extensively characterized, it remains an
engineering challenge as to the best way to incorporate a material that is a
liquid at room temperature into a (cryogenic) detector architecture.

We thank Dan Carney, Tanner Trickle, Kathryn Zurek, Katherine Inzani and
Thomas Harrelson for helpful discussions. This work was supported by the US
Department of Energy under contract DEAC02-05CH11231 and Quantum Information
Science Enabled Discovery (QuantISED) for High Energy Physics grant KA2401032.
Computational resources were provided by the National Energy Research
Scientific Computing Center and the Molecular Foundry, DOE Office of Science
User Facilities supported by the Office of Science, U.S. Department of Energy
under Contract No. DEAC02-05CH11231. The work performed at the Molecular
Foundry was supported by the Office of Science, Office of Basic Energy
Sciences, of the U.S. Department of Energy under the same contract.

## References

  * Essig  _et al._ (2012) R. Essig, J. Mardon, and T. Volansky, Phys. Rev. D 85, 076007 (2012).
  * Hall  _et al._ (2010) L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, Journal of High Energy Physics 2010, 1 (2010).
  * Bernal  _et al._ (2017) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, International Journal of Modern Physics A 32, 1730023 (2017).
  * Kahn and Lin (2022) Y. Kahn and T. Lin, Reports on Progress in Physics (2022).
  * Essig  _et al._ (2016) R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto, T. Volansky, and T.-T. Yu, Journal of High Energy Physics 2016, 1 (2016).
  * Knapen  _et al._ (2018) S. Knapen, T. Lin, M. Pyle, and K. M. Zurek, Physics Letters B 785, 386 (2018).
  * Hochberg  _et al._ (2018) Y. Hochberg, Y. Kahn, M. Lisanti, K. M. Zurek, A. G. Grushin, R. Ilan, S. M. Griffin, Z.-F. Liu, S. F. Weber, and J. B. Neaton, Physical Review D 97, 015004 (2018).
  * Marsh  _et al._ (2019) D. J. E. Marsh, K. C. Fong, E. W. Lentz, L. Šmejkal, and M. N. Ali, Phys. Rev. Lett. 123, 121601 (2019).
  * Røising  _et al._ (2021) H. S. Røising, B. Fraser, S. M. Griffin, S. Bandyopadhyay, A. Mahabir, S.-W. Cheong, and A. V. Balatsky, Physical Review Research 3, 033236 (2021).
  * Griffin  _et al._ (2018) S. Griffin, S. Knapen, T. Lin, and K. M. Zurek, Phys. Rev. D 98, 115034 (2018).
  * Drukier  _et al._ (1986) A. K. Drukier, K. Freese, and D. N. Spergel, Physical Review D 33, 3495 (1986).
  * Griffin  _et al._ (2021) S. M. Griffin, Y. Hochberg, K. Inzani, N. Kurinsky, T. Lin, and T. C. Yu, Physical Review D 103, 075002 (2021).
  * Pyle  _et al._ (2015) M. Pyle, E. Figueroa-Feliciano, and B. Sadoulet, arXiv preprint arXiv:1503.01200 (2015).
  * Griffin  _et al._ (2020a) S. M. Griffin, K. Inzani, T. Trickle, Z. Zhang, and K. M. Zurek, Phys. Rev. D 101, 055004 (2020a).
  * Campbell-Deem  _et al._ (2020) B. Campbell-Deem, P. Cox, S. Knapen, T. Lin, and T. Melia, Physical Review D 101, 036006 (2020).
  * Trickle  _et al._ (2020) T. Trickle, Z. Zhang, K. M. Zurek, K. Inzani, and S. M. Griffin, Journal of High Energy Physics 2020 (2020), 10.1007/jhep03(2020)036.
  * Lewin and Smith (1996) J. Lewin and P. Smith, Astroparticle Physics 6, 87 (1996).
  * George  _et al._ (2020) J. George, G. Hautier, A. P. Bartók, G. Csányi, and V. L. Deringer, The Journal of Chemical Physics 153, 044104 (2020).
  * Hansen (2021) T. C. Hansen, Nature Communications 12, 3161 (2021).
  * Bernal and Fowler (1933) J. D. Bernal and R. H. Fowler, The Journal of Chemical Physics 1, 515 (1933), https://doi.org/10.1063/1.1749327 .
  * Griffin and Spaldin (2017) S. M. Griffin and N. A. Spaldin, Journal of Physics: Condensed Matter 29, 343001 (2017).
  * Pauling (1935) L. Pauling, Journal of the American Chemical Society 57, 2680 (1935).
  * Leadbetter  _et al._ (1985) A. J. Leadbetter, R. C. Ward, J. W. Clark, P. A. Tucker, T. Matsuo, and H. Suga, The Journal of Chemical Physics 82, 424 (1985), https://doi.org/10.1063/1.448763 .
  * Tajima  _et al._ (1982) Y. Tajima, T. Matsuo, and H. Suga, Nature 299, 810–812 (1982).
  * Parkkinen  _et al._ (2014) P. Parkkinen, S. Riikonen, and L. Halonen, J. Phys. Chem. C 118, 26264 (2014).
  * Santra  _et al._ (2013) B. Santra, J. Klimeš, A. Tkatchenko, D. Alfè, B. Slater, A. Michaelides, R. Car, and M. Scheffler, The Journal of Chemical Physics 139, 154702 (2013).
  * Raza  _et al._ (2011) Z. Raza, D. Alfè, C. G. Salzmann, J. Klimeš, A. Michaelides, and B. Slater, Phys. Chem. Chem. Phys. 13, 19788 (2011).
  * Klimeš  _et al._ (2009) J. Klimeš, D. R. Bowler, and A. Michaelides, Journal of Physics: Condensed Matter 22, 022201 (2009).
  * Klimeš  _et al._ (2011) J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).
  * Fortes (2018) A. D. Fortes, Acta Crystallographica Section B 74, 196 (2018).
  * Cherubini  _et al._ (2021) M. Cherubini, L. Monacelli, and F. Mauri, The Journal of Chemical Physics 155, 184502 (2021), https://doi.org/10.1063/5.0062689 .
  * Yuan  _et al._ (2017) Z.-Y. Yuan, P. Zhang, S.-k. Yao, Y.-B. Lu, H.-Z. Yang, H.-W. Luo, and Z.-J. Zhao, RSC Adv. 7, 36801 (2017).
  * Yuan  _et al._ (2022) X.-Q. Yuan, X.-H. Yu, X.-L. Zhu, X.-C. Wang, X.-Y. Liu, J.-W. Cao, X.-L. Qin, and P. Zhang, ACS Omega 7, 2970 (2022).
  * Brown-Altvater  _et al._ (2016) F. Brown-Altvater, T. Rangel, and J. B. Neaton, Phys. Rev. B 93, 195206 (2016).
  * Chang  _et al._ (2021) J. H. Chang, R. Essig, and A. Reinert, Journal of High Energy Physics 2021, 141 (2021).
  * Trickle  _et al._ (2022) T. Trickle, Z. Zhang, and K. M. Zurek, Phys. Rev. D 105 (2022), 10.1103/PhysRevD.105.015001.
  * Trickle (2020) T. Trickle, “Phonodark,” https://github.com/tanner-trickle/PhonoDark (2020).
  * Griffin  _et al._ (2020b) S. M. Griffin, K. Inzani, T. Trickle, Z. Zhang, and K. M. Zurek, Phys. Rev. D 101, 055004 (2020b).
  * Lory  _et al._ (2017) P.-F. Lory, S. Pailhès, V. M. Giordano, H. Euchner, H. D. Nguyen, R. Ramlau, H. Borrmann, M. Schmidt, M. Baitinger, M. Ikeda, P. Tomeš, M. Mihalkovič, C. Allio, M. R. Johnson, H. Schober, Y. Sidis, F. Bourdarot, L. P. Regnault, J. Ollivier, S. Paschen, Y. Grin, and M. de Boissieu, Nature Communications 8, 491 (2017).
  * Kang  _et al._ (2008) K. Kang, D. Cahill, T. Ozel, and M. Shim, in  _APS March Meeting Abstracts_ (2008) pp. S29–007.
  * Khusnatdinov and Petrenko (1996) N. Khusnatdinov and V. Petrenko, Journal of Crystal Growth 163, 420 (1996).