# 8.9-hr Rotation in the Partly Burnt Runaway Stellar Remnant LP 40$-$365 (GD
492)

J. J. Hermes Department of Astronomy & Institute for Astrophysical Research,
Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA  Odelia
Putterman Department of Astronomy & Institute for Astrophysical Research,
Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA  Mark A.
Hollands Department of Physics, University of Warwick, Coventry CV4 7AL,UK
David J. Wilson McDonald Observatory, University of Texas at Austin, Austin,
TX 78712, USA  Andrew Swan Department of Physics & Astronomy, University
College London, Gower Street, London WC1E 6BT, UK  Roberto Raddi Departament
de Física, Universitat Politècnica de Catalunya, c/Esteve Terrades 5, 08860
Castelldefels, Spain  Ken J. Shen Department of Astronomy and Theoretical
Astrophysics Center, University of California, Berkeley, CA 94720, USA  Boris
T. Gänsicke Department of Physics, University of Warwick, Coventry CV4 7AL,UK

###### Abstract

We report the detection of $8.914$-hr variability in both optical and
ultraviolet light curves of LP 40$-$365 (also known as GD 492), the prototype
for a class of partly burnt runaway stars that have been ejected from a binary
due to a thermonuclear supernova event. We first detected this $1.0$%
amplitude variation in optical photometry collected by the Transiting
Exoplanet Survey Satellite. Re-analysis of observations from the Hubble Space
Telescope at the TESS period and ephemeris reveal a $5.8$% variation in the
ultraviolet of this $9800$ K stellar remnant. We propose that this $8.914$-hr
photometric variation reveals the current surface rotation rate of LP
40$-$365, and is caused by some kind of surface inhomogeneity rotating in and
out of view, though a lack of observed Zeeman splitting puts an upper limit on
the magnetic field of $<$$20$ kG. We explore ways in which the present
rotation period can constrain progenitor scenarios if angular momentum was
mostly conserved, which suggests that the survivor LP 40$-$365 was not the
donor star but was most likely the bound remnant of a mostly disrupted white
dwarf that underwent advanced burning from an underluminous (Type Iax)
supernova.

stars: individual: LP 40-365 — stars: white dwarfs — stars: rotation — stars:
kinematics and dynamics

##  1 Introduction

LP 40$-$365 (a.k.a. GD 492) is an exceptionally peculiar star: it is one of
the most metal-rich stars known, with an atmosphere dominated by oxygen and
neon that is also abundant in heavy elements from partial oxygen and silicon
burning, with no detectable hydrogen or helium (Vennes et al., 2017). Most
peculiarly, it is rapidly departing the Galaxy, with a rest-frame velocity of
$v_{\mathrm{rf}}$ $\simeq$ $852$ km s-1 that makes it unbound to the Milky
Way; unlike hypervelocity stars (Brown, 2015), LP 40$-$365 almost certainly
did not come from the Galactic center (Raddi et al., 2018a). Currently, it
appears to be an isolated star without radial-velocity variability, to a limit
of roughly $2-3$ km s-1 (Raddi et al., 2018b).

LP 40$-$365, the prototype to a small but growing class of similar objects
(Raddi et al., 2019), derives these peculiarities because it was likely
involved in a thermonuclear supernova event that ejected it from a short-
period binary system. The original discoverers suggested that stars like LP
40$-$365 are most likely the bound remnants of partially disruptive
thermonuclear supernovae, underluminous explosions often referred to as Type
Iax supernovae (Vennes et al., 2017). Type Iax supernovae are thought to arise
from an asymmetrically ignited deflagration in a Chandrasekhar-mass white
dwarf that does not fully disrupt the star and leaves behind a bound remnant
(e.g., Jordan et al., 2012; Foley et al., 2013; Fink et al., 2014; Jha, 2017).

Gaia has enabled discovery of another class of runaways, the so-called
“dynamically driven double-degenerate double-detonation” (D6) stars, which are
moving with such high velocities ($>$$1000$ km s-1) they must be the surviving
donor stars in double white dwarf systems that underwent thermonuclear
supernovae (Shen et al., 2018b). In contrast, it is hard to constrain with
certainty whether the LP 40$-$365-like stars were the donor or accretor in a
close binary, since detailed abundance analyses is complicated by highly
uncertain metal diffusion times and nucleosynthetic yields. Most aspects of
abundance analysis suggests that these objects are the bound remnants of the
actual white dwarf in which the partial deflagration occurred (Raddi et al.,
2018b). However, the LP 40$-$365-like stars have also been analyzed as the
polluted, inflated remnants of hot subdwarf donors in short-period progenitor
systems (Bauer et al., 2019).

One useful boundary condition that could help us constrain the progenitor
system for LP 40$-$365 would be the final rotation rate of the stellar
remnant. The highest-resolution spectroscopic observations of LP 40$-$365 have
only put an upper limit on the rotation velocity ($v_{\rm rot}\sin{i}<50$ km
s-1) that implies a rotation period longer than 4 hr, based on current radius
estimates (Raddi et al., 2018b).

We report here the significant detection of photometric variability in LP
40$-$365 at a period of $8.914$ hr, present in an optical light curve
collected by the Transiting Exoplanet Survey Satellite (TESS), with additional
confirmation from an ultraviolet light curve from archival observations
collected by the Hubble Space Telescope (HST) and possible detection in a
near-infrared light curve from data collected by the Wide-field Infrared
Survey Explorer (WISE).

We interpret this $8.914$-hr signal to correspond to the current surface
rotation rate of LP 40$-$365\. Section 2 and Section 3 of this Letter detail
our observations and analysis, respectively. We conclude with a discussion of
how the final rotation period we observe now in LP 40$-$365 can inform our
understanding of the progenitor conditions that led to the formation of this
partly burnt runaway stellar remnant.

##  2 Observations and Reductions

###  2.1 TESS Observations

LP 40$-$365 (TIC 198510602, $T$=15.5 mag) was observed by TESS for nearly four
months, in Sector 14 (2019 Jul 18 to 2019 Aug 14), Sector 15 (2019 Aug 15 to
2019 Sep 10), Sector 20 (2019 Dec 24 to 2020 Jan 20), and Sector 21 (2020 Jan
21 to 2020 Feb 18), all at 30-min cadence in the full-frame images (it was not
targeted by any shorter, 2-min-cadence observations). The light curves for all
four sectors were extracted using the python package eleanor (Feinstein et
al., 2019), and analyzed using the lightkurve package (Lightkurve
Collaboration et al., 2018). TESS utilizes a broad, red optical filter
centered at roughly $7865$ Å.

The second panel of Figure 1 shows the full TESS light curve. Gaps in the data
are present roughly every $13.7$ days during each regular data download per
orbit, and increased scatter due to reflected light is often present before
the data downloads when the spacecraft is closer to Earth.

Figure 1: Full light curves of LP 40$-$365 used in this analysis. Top Panel:
Five orbits of HST/STIS ultraviolet observations from 2018 May. Second Panel:
Four sectors of TESS full-frame images from 2019 July to 2020 February. Third
Panel: Periodogram of TESS observations revealing a dominant period at $8.914$
hr and its second harmonic. Bottom Panel: More than 9.5 yr of WISE infrared
photometry from 2010 April to 2019 November.

Our light curves were extracted using a fixed aperture, tuned for each sector
depending on the spacecraft roll angle. We used a 3-px aperture for Sector 14,
a 5-px aperture for Sectors 15 and 20, and a 2-px aperture for Sector 21. This
changes the amount of extra light from background sources in the large pixel
size ($>$21 arcsec pixel-1) for each sector, which makes comparing the
relative amplitude from each sector more difficult. Therefore, the relative
flux in our full light curve in Figure 1 was normalized to the median flux for
each sector. No amplitude corrections have been made for possible crowding in
the large sky area covered by these apertures, which may result in a
suppressed amplitude in the TESS bandpass caused by flux dilution. Our full
TESS light curve has 4338 points over 214.4 days, for a roughly 42% duty
cycle.

###  2.2 HST Observations

We targeted LP 40$-$365 with the Space Telescope Imagining Spectrograph (STIS)
instrument aboard HST over five orbits in 2018 May, which were first described
and analyzed in Raddi et al. (2019). The observations were proposed to analyze
the ultraviolet spectrum and better constrain the effective temperature of LP
40$-$365, as well as improve our elemental abundance analysis.

We have reanalyzed these $15{,}253$ seconds of STIS observations taken with
the G230L grating by taking advantage of their time-tagged photon collection
(Woodgate et al., 1998). We interacted with these data using the python
package stistools111https://stistools.readthedocs.io available from the Space
Telescope Science Institute. We binned the data into 3-min intervals, then
summed the spectra of each in the wavelength range with the most signal-to-
noise (S/N), $1700-3100$ Å. This effectively created an ultraviolet light
curve with 3-min exposures and an effective central wavelength of $2520$ Å. We
corrected the light curve to barycentric time using Astropy (Astropy
Collaboration et al., 2013). The full HST/STIS light curve has 82 points over
4.05 days, for a roughly 4% duty cycle, and is shown in the top panel of
Figure 1.

HST is in constant orbit around the Earth, with natural changes in the focus
of instruments like STIS due to thermal evolution (often described as
breathing, e.g., Sing et al. 2013). Since our TESS observations revealed
variability at much longer timescales than the HST orbit we have not corrected
for this effect. However, it likely has added increased scatter on timescales
shorter than $95.4$ min.

###  2.3 WISE Observations

We constructed an infrared light curve for LP 40$-$365 using data from the
Wide-field Infrared Survey Explorer (WISE) spacecraft, which is performing an
all-sky survey, visiting most sources twice per year (Wright et al., 2010). We
queried the AllWISE Multiepoch Photometry Table and NEOWISE-R Single Exposure
Source Table222Available at irsa.ipac.caltech.edu, obtaining photometry in the
W1 and W2 bands, taken during the decade 2010–2019. Here we only analyzed the
higher-S/N observations in the W1 band, which has an effective wavelength of
roughly $3.4$ µm.

The WISE scanning pattern resulted in on average 32 separate 7.7-second
exposures per visit for LP 40$-$365\. Visits are spaced about six months
apart, and last between two and eight days. Photometry is often analyzed after
combining exposures within each visit (e.g., Swan et al. 2019; Secrest &
Satyapal 2020), but we kept our observations unbinned, resulting in the
relatively large scatter seen in the bottom panel of Figure 1. Still, each
exposure used has a formal S/N $>$$4$, with an average S/N of $9$. We
corrected the light curve to barycentric time using Astropy.

The point spread function in the W1 band has a full-width-half-maximum of 6.1
arcsec, and thus source confusion is a potential concern (Dennihy et al.,
2020). However, no excess infrared flux is observed (Raddi et al., 2018b), and
the median discrepancy between measured positions and those predicted from
Gaia astrometry is 0.3 arcsec. The full WISE/W1 light curve has 451 points
over 3480 days, for a roughly 0.001% duty cycle, and is shown in the bottom
panel of Figure 1.

##  3 Light Curve Analysis

A Lomb-Scargle periodogram of all four sectors of TESS data, shown in the
third panel of Figure 1, reveals two significant peaks: one dominant at
roughly $8.90$ hr and another peak at exactly half that period, roughly $4.45$
hr. The 0.1% significance threshold shown, at 0.39% amplitude, was
bootstrapped by keeping the time sampling identical but randomly shuffling the
flux values in the light curve, noting the maximum peak for 99.9% of all light
curves from noise alone, as described in Hermes et al. (2015).

A simultaneous nonlinear least-squares fit of all TESS data, forcing the
second frequency to be exactly twice that of the dominant signal, finds a
best-fit period of $8.91405\pm 0.00035$ hr. We have performed the same
weighted least-squares fit to the combined TESS+WISE dataset, which yields a
more precise period, which we adopt throughout: $8.914122\pm 0.000020$ hr.
From the first observed time of minimum in TESS we define the ephemeris as:

$\rm{BJD}_{\rm TDB}=2458683.7076(23)+0.37142177(83)\;E$

We fix this period and compute a linear least-squares fit for each of the
three light curves to define the amplitude of variability in each dataset. We
find an ultraviolet amplitude of $5.79\pm 0.21$% from the HST data, an optical
amplitude of $1.038\pm 0.042$% from the TESS data, and an infrared amplitude
of $1.94\pm 0.81$% from the WISE data. The uncertainties on the amplitudes are
purely formal, and do not include any systematic effects.

Figure 2: Folded light curves for LP 40$-$365, with each phased using the
ephemeris defined from the WISE+TESS observations and repeated for clarity.
Top Panel: All TESS data, which have an effective wavelength of roughly $7900$
Å, binned into 200 phase bins. Middle Panel: All HST/STIS data, unbinned,
which have an effective wavelength of roughly $2520$ Å. Bottom Panel: All
WISE/W1 data binned into 20 phase bins, which have an effective wavelength of
roughly $3.4$ µm.

We use the ephemeris to phase all light curves to the same baseline, which we
show visually in Figure 2. The best-fit amplitudes reported above are shown
visually as a faint green dashed line in each panel. The phase-folded, binned
WISE data is best fit with an amplitude larger than the full unbinned light
curve ($2.20\pm 0.87$%).

The significant detection of the second harmonic of the $8.914$-hr fundamental
period in the periodogram in Figure 1 suggests an asymmetry in the light
curve. Non-sinusoidal light curves have been observed in other white dwarfs
exhibiting photometric variability from surface spots (e.g., Hermes et al.
2017a, b; Hoard et al. 2018), and the second harmonic may be modeled in the
future to constrain the inclination to the star. The folded light curves in
Figure 2 show that the photometric variations are mostly sinusoidal; excluding
this harmonic does not change our reported ephemeris or measured amplitudes
within the stated uncertainties.

We detect $8.914$-hr variability with significance in multiple datasets; the
variability appears phase coherent over multiple years, as well as between the
optical and ultraviolet wavelengths. In 2017 we observed LP 40$-$365 for 4 hr
from the 2-m Liverpool Telescope but saw no optical variability to a limit of
roughly 0.3% amplitude (Raddi et al., 2018b); it is likely the signal we see
here was too long to detect from that light curve taken through a $V+R$-band
filter. We cannot certify that the variability observed in the WISE
observations is formally significant, although it is suggestive that the
phased data in Figure 2 shares the same photometric minimum as the TESS
observations. Only 1% of field sources vary at the wavelengths covered by W1
(Kozłowski et al., 2010; Polimera et al., 2018), so it is likely that the
periodic variation we detect in WISE is from LP 40$-$365.

We reserve a more detailed analysis of the amplitude ratios between the
ultraviolet-optical-infrared datasets until we can properly assess crowding in
the TESS observations, which could contribute to flux dilution and thus
optical amplitude suppression. It is reasonable that the ultraviolet
variations have higher amplitude than the optical variations, given that this
is a 9800 K stellar remnant (Raddi et al., 2019).

##  4 Discussion and Conclusion

We interpret the coherent $8.914$-hr photometric variability in LP 40$-$365 as
the surface rotation period of this runaway stellar remnant. Spots caused by
surface inhomogeneities on white dwarfs have been seen in a growing number of
stars observed with space-based photometry from Kepler and TESS (Hermes et
al., 2017a), especially in the presence of strong magnetic fields (Brinkworth
et al., 2013). Without knowing the topology of the magnetic field of this
white dwarf, the inferred rotation period could be off by a factor of two; our
inference of a $8.914$-hr rotation rate assumes a dipole field geometry.

We have reanalyzed the highest-resolution optical spectra collected from the
R1200R grating from the 4.2-meter William Herschel Telescope on La Palma,
described in Raddi et al. (2018b). The sharp metal lines observed in that
spectrum allowed us to put an upper limit on rotation of $v_{\rm
rot}\,\sin{i}\,<\,50$ km s-1 (we would infer a rotation velocity of roughly 22
km s-1 from an $8.9$-hr rotation period). Similarly, using the observed shape
of the calcium triplet between 8498–8662 Å (which is particularly sensitive to
small magnetic fields), we do not detect any line-profile changes caused by
Zeeman splitting from a magnetic field. This allows us to put an upper limit
on the magnetic field in LP 40$-$365 of $<$ $20$ kG. Such a small field may
not preclude formation of spots from surface inhomogeneities — for example, we
observe a spot from rotation in a pulsating white dwarf which has a global
magnetic field $<$ $10$ kG (Hermes et al., 2017b).

The best estimates of the parameters of LP 40$-$365 are detailed in Raddi et
al. (2019): $T_{\rm eff}=9800\pm 300$ K, $\log{g}=5.5\pm 0.3$ [cgs],
$L=0.20\pm 0.04$ L☉, $R=0.16\pm 0.01$ R☉, and $M=0.28^{+0.28}_{-0.14}$ M☉.
Using these atmospheric parameters, we find the gas pressure at an optical
depth of $\tau=2/3$ is $10^{4.7}$ dyne cm-2. The magnetic pressure will exceed
this gas pressure at roughly 1 kG, so our 20 kG upper limit allows for a field
that can measurably influence the atmospheric structure at the photosphere.

It is likely that the partly burnt remnant LP 40$-$365 is still highly
inflated; tracing back the kinematics yields an expected flight time of this
star away from the Galactic disk of $5.3$ Myr, which is considerably shorter
than the Kelvin-Helmholtz timescale of $32$ Myr. Considering this radius
inflation, we have explored how the current rotation period can constrain
progenitor scenarios if angular momentum in the system is mostly conserved.

We first consider if LP 40$-$365 was the white dwarf that underwent a
thermonuclear event; in this first case we assume the progenitor star had a
mass near the Chandrasekhar limit ($M_{i}$ $>$ $1.3$ $M_{\sun}$), with a
corresponding radius of $R_{i}<0.004$ R☉ (Althaus et al., 2005). In this
scenario, most of the mass from this star was lost in the thermonuclear event,
and the radius has increased both because the lower-mass white dwarf is less
degenerate but also because of eventual entropy release as the remnant evolves
(e.g., Zhang et al. 2019). Assuming angular momentum transport is efficient
and mostly conserved implies an initial rotation period of less than $2$ min.
This is very fast compared to typical isolated white dwarfs, which rotate with
periods from $0.5-2$ days (Hermes et al., 2017c), but is still much slower
than breakup velocity (which would have a spin period of order seconds).

A $2$-min initial spin period is likely faster than the expected orbital
period at explosion: LP 40$-$365 likely had an ejection velocity of roughly
$v_{\rm ej}\simeq 600$ km s-1, which suggests a progenitor orbital period of
roughly $10-30$ min for most configurations (Raddi et al., 2018a). However,
such spin-up would not be surprising for an accretor: spin-up is commonly
expected (e.g., Ritter 1985) and observed (e.g., Lopes de Oliveira et al.
2020) in cataclysmic variables with high mass-transfer rates onto accreting
white dwarfs.

We also consider if LP 40$-$365 was the donor star ejected from a once-compact
binary system; this scenario was modeled in detail by Bauer et al. (2019).
They explored two different progenitors at different orbital periods: Their
two hot subdwarf donor models (a 0.344 $M_{\sun}$ donor onto a 0.927
$M_{\sun}$ accretor and a 0.233 $M_{\sun}$ donor onto a 0.779 $M_{\sun}$
accretor) predict rapid rotation for the donor remnant $10$ Myr after
explosion; all donor remnants with final masses $>$$0.2$ $M_{\sun}$ end up
with rotation periods faster than roughly $1$ hr, much faster than $8.9$ hr.
(Shocks and interactions with high ejecta kinetic energies will strip some or
even most of the donor mass by the time it becomes a remnant, and the donor
models specified above are pre-detonation masses.)

Combining the Keplerian orbital velocity of the donor with the expectations
that it is Roche-lobe filling allowed Bauer et al. (2019) to predict the
radius of the donor at explosion in terms of its orbital velocity. The
expected range of orbital velocities for LP 40$-$365 range from $600-800$ km
s-1 (Raddi et al., 2019), which would predict a donor radius of $0.055-0.16$
$R_{\sun}$ in orbital periods ranging from roughly $5-20$ min.

The expected radius of the donor computed via the Roche-lobe radius depends on
the unknown mass ratio of the system (Eggleton, 1983) at explosion; ranges for
$M_{1}$ are unlikely to exceed $0.2-0.5$ $M_{\sun}$ for a core-He-burning
subdwarf donor, and $M_{2}$ are unlikely to exceed $0.85-1.4$ $M_{\sun}$ given
that white dwarf detonations could be sub-Chandrasekhar in a subdwarf donor
scenario (Shen et al., 2018a). All told, we cannot explain the high ejection
velocity with the comparatively slow rotation period currently observed in LP
40$-$365 if this is the donor remnant that was tidally locked in a short-
period binary; the progenitor requires significantly more radius inflation to
slow its rotation to the period we observe today.

There are two assumptions that could affect our inferences. It is likely a
subdwarf donor would be close to tidally locked in these short-period binaries
(Preece et al., 2018), but breaking this assumption would affect initial spin-
period estimates. Moreover, we cannot exclude a scenario in which the donor
was highly inflated ($R$ $>$ $20$ $R_{\sun}$) after explosion and lost
significant angular momentum due to mass loss and magnetic braking. This phase
should be relatively short-lived ($<$ $10^{5}$ years, e.g. Pan et al. 2013),
but the mass loss may be significant enough to complicate the assumption of
angular momentum conservation.

That said, if angular momentum is mostly conserved, we find it less likely
that LP 40$-$365 represents the donor in a disrupted compact binary. Most
likely, it is still experiencing a dramatic radius expansion best explained by
it being the once-accreting white dwarf that itself underwent a thermonuclear
event. This radius inflation is also supported by another much hotter LP
40$-$365-like star, J1825$-$3757, which has a similar mass as LP 40$-$365 but
a radius at least a factor of three larger (Raddi et al., 2019). The event
that disrupted the binary system caused both components to sling-shot apart at
high speed. That thermonuclear event also led to partial burning in the
progenitor star to LP 40$-$365, as proposed by Vennes et al. (2017) and Raddi
et al. (2018a, b), based in part on the currently observed heavy-element
abundances.

Clarifying the role of LP 40$-$365 in the explosion is important to further
interpreting the abundances, especially to connect to expected nucleosynthetic
yields (e.g., Leung & Nomoto 2020). This connection is vital to help us better
understand the peculiar properties of the class of partially burnt stellar
remnants, the LP 40$-$365 stars (Raddi et al., 2019).

## Acknowledgements

We acknowledge helpful comments from the anonymous referee, and thank Evan
Bauer for valuable discussions. Support for this work was in part provided by
NASA TESS Cycle 2 Grant 80NSSC20K0592. M.A.H. has received funding from the
European Research Council under the European Union’s Horizon 2020 research and
innovation programme n. 677706 (WD3D). A.S. has received support from Science
and Technology Facilities Council grant ST/R000476/1. R.R. has received
funding from the postdoctoral fellowship program Beatriu de Pinós, funded by
the Secretary of Universities and Research (Government of Catalonia) and by
the Horizon 2020 program of research and innovation of the European Union
under the Maria Skłodowska-Curie grant agreement No 801370. K.J.S. was
supported in part by NASA through the Astrophysics Theory Program
(NNX17AG28G). B.T.G. was supported by the UK STFC grant ST/T000406/1 and by a
Leverhulme Research Fellowship.

This paper includes data collected by the TESS mission. Funding for the TESS
mission is provided by the NASA’s Science Mission Directorate. Based on
observations made with the NASA/ESA Hubble Space Telescope, obtained from the
Data Archive at the Space Telescope Science Institute, which is operated by
the Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-26555. These observations are associated with program #15431,
and funding in-part provided by programs #15871 and #15918. This publication
makes use of data products from the Wide-field Infrared Survey Explorer, which
is a joint project of the University of California, Los Angeles, and the Jet
Propulsion Laboratory/California Institute of Technology, funded by the
National Aeronautics and Space Administration.

## References

  * Althaus et al. (2005) Althaus, L. G., García-Berro, E., Isern, J., & Córsico, A. H. 2005, A&A, 441, 689, doi: 10.1051/0004-6361:20052996
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bauer et al. (2019) Bauer, E. B., White, C. J., & Bildsten, L. 2019, ApJ, 887, 68, doi: 10.3847/1538-4357/ab4ea4
  * Brinkworth et al. (2013) Brinkworth, C. S., Burleigh, M. R., Lawrie, K., Marsh, T. R., & Knigge, C. 2013, ApJ, 773, 47, doi: 10.1088/0004-637X/773/1/47
  * Brown (2015) Brown, W. R. 2015, ARA&A, 53, 15, doi: 10.1146/annurev-astro-082214-122230
  * Dennihy et al. (2020) Dennihy, E., Farihi, J., Fusillo, N. P. G., & Debes, J. H. 2020, ApJ, 891, 97, doi: 10.3847/1538-4357/ab7249
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368, doi: 10.1086/160960
  * Feinstein et al. (2019) Feinstein, A. D., Montet, B. T., Foreman-Mackey, D., et al. 2019, PASP, 131, 094502, doi: 10.1088/1538-3873/ab291c
  * Fink et al. (2014) Fink, M., Kromer, M., Seitenzahl, I. R., et al. 2014, MNRAS, 438, 1762, doi: 10.1093/mnras/stt2315
  * Foley et al. (2013) Foley, R. J., Challis, P. J., Chornock, R., et al. 2013, ApJ, 767, 57, doi: 10.1088/0004-637X/767/1/57
  * Hermes et al. (2017a) Hermes, J. J., Gänsicke, B. T., Gentile Fusillo, N. P., et al. 2017a, MNRAS, 468, 1946, doi: 10.1093/mnras/stx567
  * Hermes et al. (2017b) Hermes, J. J., Kawaler, S. D., Bischoff-Kim, A., et al. 2017b, ApJ, 835, 277, doi: 10.3847/1538-4357/835/2/277
  * Hermes et al. (2015) Hermes, J. J., Gänsicke, B. T., Bischoff-Kim, A., et al. 2015, MNRAS, 451, 1701, doi: 10.1093/mnras/stv1053
  * Hermes et al. (2017c) Hermes, J. J., Gänsicke, B. T., Kawaler, S. D., et al. 2017c, ApJS, 232, 23, doi: 10.3847/1538-4365/aa8bb5
  * Hoard et al. (2018) Hoard, D. W., Howell, S. B., Roettenbacher, R. M., et al. 2018, AJ, 156, 119, doi: 10.3847/1538-3881/aad238
  * Jha (2017) Jha, S. W. 2017, Type Iax Supernovae, ed. A. W. Alsabti & P. Murdin, 375, doi: 10.1007/978-3-319-21846-5_42
  * Jordan et al. (2012) Jordan, George C., I., Perets, H. B., Fisher, R. T., & van Rossum, D. R. 2012, ApJ, 761, L23, doi: 10.1088/2041-8205/761/2/L23
  * Kozłowski et al. (2010) Kozłowski, S., Kochanek, C. S., Stern, D., et al. 2010, ApJ, 716, 530, doi: 10.1088/0004-637X/716/1/530
  * Leung & Nomoto (2020) Leung, S.-C., & Nomoto, K. 2020, ApJ, 900, 54, doi: 10.3847/1538-4357/aba1e3
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python.  http://ascl.net/1812.013
  * Lopes de Oliveira et al. (2020) Lopes de Oliveira, R., Bruch, A., Rodrigues, C. V., Oliveira, A. S., & Mukai, K. 2020, ApJ, 898, L40, doi: 10.3847/2041-8213/aba618
  * Pan et al. (2013) Pan, K.-C., Ricker, P. M., & Taam, R. E. 2013, ApJ, 773, 49, doi: 10.1088/0004-637X/773/1/49
  * Polimera et al. (2018) Polimera, M., Sarajedini, V., Ashby, M. L. N., Willner, S. P., & Fazio, G. G. 2018, MNRAS, 476, 1111, doi: 10.1093/mnras/sty164
  * Preece et al. (2018) Preece, H. P., Tout, C. A., & Jeffery, C. S. 2018, MNRAS, 481, 715, doi: 10.1093/mnras/sty2091
  * Raddi et al. (2018a) Raddi, R., Hollands, M. A., Gänsicke, B. T., et al. 2018a, MNRAS, 479, L96, doi: 10.1093/mnrasl/sly103
  * Raddi et al. (2018b) Raddi, R., Hollands, M. A., Koester, D., et al. 2018b, ApJ, 858, 3, doi: 10.3847/1538-4357/aab899
  * Raddi et al. (2019) —. 2019, MNRAS, 489, 1489, doi: 10.1093/mnras/stz1618
  * Ritter (1985) Ritter, H. 1985, A&A, 148, 207
  * Secrest & Satyapal (2020) Secrest, N. J., & Satyapal, S. 2020, ApJ, 900, 56, doi: 10.3847/1538-4357/ab9309
  * Shen et al. (2018a) Shen, K. J., Kasen, D., Miles, B. J., & Townsley, D. M. 2018a, ApJ, 854, 52, doi: 10.3847/1538-4357/aaa8de
  * Shen et al. (2018b) Shen, K. J., Boubert, D., Gänsicke, B. T., et al. 2018b, ApJ, 865, 15, doi: 10.3847/1538-4357/aad55b
  * Sing et al. (2013) Sing, D. K., Lecavelier des Etangs, A., Fortney, J. J., et al. 2013, MNRAS, 436, 2956, doi: 10.1093/mnras/stt1782
  * Swan et al. (2019) Swan, A., Farihi, J., & Wilson, T. G. 2019, MNRAS, 484, L109, doi: 10.1093/mnrasl/slz014
  * Vennes et al. (2017) Vennes, S., Nemeth, P., Kawka, A., et al. 2017, Science, 357, 680, doi: 10.1126/science.aam8378
  * Woodgate et al. (1998) Woodgate, B. E., Kimble, R. A., Bowers, C. W., et al. 1998, PASP, 110, 1183, doi: 10.1086/316243
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Zhang et al. (2019) Zhang, M., Fuller, J., Schwab, J., & Foley, R. J. 2019, ApJ, 872, 29, doi: 10.3847/1538-4357/aafb34


# Deuteration removes quantum dipolar defects from KDP crystals

Bingjia Yang  Department of Chemistry, Princeton University, Princeton, NJ
08544, USA  Pinchen Xie  Program in Applied and Computational Mathematics,
Princeton University, Princeton, NJ 08544, USA  Roberto Car  Department of
Chemistry, Department of Physics, Program in Applied and Computational
Mathematics, Princeton Materials Institute, Princeton University, Princeton,
NJ 08544, USA

The structural, dielectric, and thermodynamic properties of the hydrogen-
bonded ferroelectric crystal potassium dihydrogen phosphate
($\mathbf{KH_{2}PO_{4}}$), KDP for short, differ significantly from those of
DKDP ($\mathbf{KD_{2}PO_{4}}$) [1, 2]. It is well established that deuteration
affects the interplay of hydrogen-bond switches and heavy ion displacements
that underlie the emergence of macroscopic polarization, but a detailed
microscopic model is missing [3, 4]. Here we show that all-atom path integral
molecular dynamics simulations can predict the isotope effects, revealing the
microscopic mechanism that differentiates KDP and DKDP. Proton tunneling in
the hydrogen bonds generates phosphate configurations that do not contribute
to the polarization. These dipolar defects are always present in KDP, but
disappear at low temperatures in DKDP, which behaves more classically. Quantum
disorder confers residual entropy to ferroelectric KDP, explaining its lower
spontaneous polarization and transition entropy relative to DKDP. Tunneling
should also contribute to the anomalous heat capacity observed in KDP near
absolute zero [5, 6]. The prominent role of quantum fluctuations in KDP is
related to the unusual strength of the hydrogen bonds in this system and
should be equally important in the other crystals of the KDP family, which
exhibit similar isotope effects [7].

KDP is the prototypical member of a class of isomorphous hydrogen-bonded
molecular crystals with formula unit $\mathrm{MH_{2}XO_{4}}$ (M=K, Rb, Cs,
$\mathrm{NH_{4}}$; X=P, As) [7]. The $\mathrm{XO_{4}}$ groups are held
together by strong hydrogen bonds (H-bonds) (Fig.1(a)). The H-bonds are
correlated with the polar distortion of the $\mathrm{XO_{4}}$ groups giving
rise to ferroelectricity when M=K, Rb, Cs, and to anti-ferroelectricity when
M=$\mathrm{NH_{4}}$. In these materials, $T_{c}$, the paraelectric to
ferroelectric (antiferroelectric) phase transition (PT) temperature, decreases
with increasing pressure tending to zero at some critical pressure [7]. When
hydrogen is substituted with its heavier isotopes, changes of $T_{c}$ of the
order of $100$K are observed [7]. At ambient pressure, KDP and DKDP have
$T_{c}$ equal to 123K and 229K, respectively. The isotope effects modulate
piezoelectric, electro-optical, and nonlinear optical properties of KDP/DKDP
[8, 9], leading to different applications such as electro-optic modulators,
frequency converters, and piezoelectric devices [10, 11]. Isotope effects on
thermodynamic properties are a manifestation of nuclear quantum effects (NQE)
in the atomic dynamics.

Empirical observations indicate that NQE make strong H-bonds stronger [12], an
effect that should favor more symmetric bond configurations in KDP than in
DKDP. Structural changes across the PT can be related to two non-exclusive
mechanisms: displacive and order-disorder [13]. In a displacive transition,
spontaneous polarization emerges from an almost uniform structural distortion
combining irreversible off-centering of the H/D ions with a polar distortion
of the $\text{PO}_{4}$ groups induced by a displacement of the P ions as in
Fig.1(b). In the paraelectric phase, the H-bonds are uniformly symmetrized and
the $\text{PO}_{4}$ groups are non-polar. The displacive mechanism underlies
the tunneling model [14] and its variants [15, 3], which are defined by
effective Hamiltonians that treat the protons as two-level systems (TLS) and
ignore or substantially reduce the degrees of freedom of the heavy ions (K, P,
O). These models attribute the PT to the softening of a polar phonon mode
combining collective proton tunneling with a polar distortion of the
$\mathrm{PO_{4}}$ groups. The displacive mechanism was challenged by
spectroscopic experiments [3, 16, 17], which found that, in the paraelectric
phase, the O-H…O bridges are asymmetric but disordered, so that the
$\mathrm{PO_{4}}$ groups are polarized in random directions with distortions
not substantially weaker than in the ferroelectric phase. This scenario
attributes spontaneous polarization to the establishment of long-range order,
consistent with neutron diffraction experiments that show a lack of
displacive-type structural distortions across the PT [18]. Indeed, the spatial
distribution of the H/D ion relative to the neighboring O ions shows off-
centered peaks in both ferroelectric and paraelectric phases, as illustrated
in Fig.1(c). It is of great interest to understand precisely how NQE affect
the disordering of the H-bond network. The insight gained on KDP should help
understand the PT in the entire KDP family of crystals [7]. It may also
elucidate the isotope shifts larger than $100$K observed in hydrogen-bonded
$\mathrm{H_{2}O_{4}C_{4}}$ and some $\mathrm{PbHPO_{4}}$-type crystals [19,
20].

Figure 1: (a) Schematic representation [21] of the Fdd2 equilibrium structure
of KDP. The z-axis is parallel to the $c$ crystallographic axis. (b) Sketch of
a ferroelectric KDP unit. K (shared), H (shared), O, and P ions are
represented by purple, light blue, red, and yellow spheres. The grey arrow
shows the direction of the molecular dipole moment $\mathbf{d}$. (c) Left:
schematics of an O-H…O bridge. Colored shadows suggest the quantum
delocalization of the particles. The geometric center of the oxygen pair is
the origin of the coordinates. $\boldsymbol{\delta}_{xy}$ is the projection on
the $xy$ plane of the position vector of the proton. Right: contour plot of
the probability density function $n(\boldsymbol{\delta}_{xy})$, below and
above $T_{c}$, inferred from neutron diffraction experiments [2, 18]. Darker
contour lines indicate higher probability density. $\Delta_{xy}(T)$, the
distance of the maximum of $n(\boldsymbol{\delta}_{xy})$ from the origin,
quantifies proton off-centering. (d) 16 types of KDP units are divided into
four classes with decreasing energetic stability from left to right. Units in
the same class have similar energetic stability. K ions and H-bonded H ions
are not shown.

Key aspects of the order-disorder mechanism were elucidated by Slater [22] and
Takagi [23, 24], who pointed out that the local H-bond structure could be
associated with the 16 phosphate configurations depicted in Fig. 1(d), which
groups the phosphates into four classes, i.e., Ferro, Slater, Takagi, and
Nonpolar units, in order of decreasing energy stability from left to right.
The more stable Ferro and Slater units satisfy the “ice rules”, because each
phosphate has two donor and two acceptor H-bonds, in analogy with
$\mathrm{H_{2}O}$ ice. The less stable Takagi and Nonpolar units break the
“ice rules” as they have unequal numbers of donor and acceptor bonds.
Attributing nominal charges to the ions and taking $\mathrm{K^{+}}$ into
account, Ferro and Slater units ($\mathrm{H_{2}PO_{4}^{-}}$) satisfy local
charge neutrality while Takagi and Nonpolar units violate it. Within
electronic density functional theory (DFT), a molecular dipole moment
$\mathbf{d}$ can be assigned to each phosphate group sharing H and K ions with
its neighbors, as depicted in Fig. 1(b) [25]. The local dipole $\mathbf{d}$ is
defined in terms of the positions and charges of the ions (nuclei plus core
electrons) and the centroids of the valence electrons (see Methods). The
shared ions (K, H) carry half charge, the non-shared ions (P, O) and
electronic centroids carry full charge. The product of polarization
$\boldsymbol{\mathcal{P}}$ and volume $V$ for a periodic supercell is equal to
the sum of the dipoles $\mathbf{d}$ contained in that supercell modulo an
immaterial constant that does not affect the changes of
$\boldsymbol{\mathcal{P}}$ [26, 25]. In view of its connection with the
macroscopic polarization, the average $\mathbf{d}$ is a proper local order
parameter for the PT. We find that Nonpolar units have nearly vanishing
$\mathbf{d}$, Takagi and Slater units have non-zero $\mathbf{d}$ along
directions approximately orthogonal to $z$, while Ferro units have non-zero
$\mathbf{d}$ along $z$. For $T<T_{c}$, broken symmetry Ferro units dominate
and the other units act as dipolar defects destabilizing the ordered phase.
For $T>T_{c}$, extensive excitation of dipolar defects occurs.

Figure 2:  (a) The longitudinal spontaneous polarization $|\mathcal{P}_{z}|$
as a function of $T$. The experimental data are for pure KDP and
KD1.96H0.04PO4 [1]. (b) Predicted and experimental [18] lattice constants $a$
and $b$ as functions of $T$. (c) The joint PDFs $n(d_{x},0,d_{z})$ (upper) and
$n(d_{x},d_{y},0)$ (lower) for KDP and DKDP at the temperatures indicated in
panel (a).

Early models [4] explained semi-quantitatively the PT by reducing the degrees
of freedom of each phosphate unit to the 16 discrete states of Fig. 1(d),
ignoring the fluctuations that connect distorted configurations. At low
temperatures, these fluctuations should be negligible in classical but not in
quantum ferroelectrics. Thermal and quantal fluctuations, including the
effects of tunneling and anharmonicity [27, 28], can be simulated with
controllable accuracy with all-atom path-integral molecular dynamics (PIMD).
In this approach, $N$ distinguishable quantum ions are mapped onto $N$
classical ring polymers made of an equal number, $L$, of beads, subject to
intra- and inter-polymer interactions. The intra-polymer interaction is
harmonic and describes the quantum delocalization of a particle, while the
inter-polymer interaction, given by the potential energy surface of the system
rescaled by $L$, describes the physical interactions of the ions. For large
enough $L$, the equilibrium values of position-dependent observables estimated
by sampling the ring polymer distribution reproduce the corresponding quantum
values (see Methods). Modeling the polarization evolution requires knowledge
of the polarization surface, which gives the dependence of the polarization on
the ion coordinates. This is possible with ab initio PIMD [29], an approach in
which the potential energy and polarization surfaces are computed on the fly
within Kohn-Sham DFT [30]. So far only a small periodic system of four KDP
units could be studied in this way [31], due to the computational cost of ab
initio PIMD. In these simulations, the accuracy of the DFT approximation is
crucial. Ref. [28] found that semi-local functionals tend to overestimate the
quantum delocalization of the hydrogens, destabilizing or even completely
eliminating the ferroelectric phase, a difficulty likely associated to the
self-interaction or delocalization error [32, 33] of DFT approximations. Due
to the exponential dependence of tunneling on the barrier height, a small
underestimation of the latter can make the system a quantum paraelectric by
removing the broken symmetry O-H…O configurations of the ferroelectric phase.

In this work, we overcome the limitations of ab initio PIMD and establish DFT-
based models for KDP/DKDP that agree well with experiments. Specifically, we
use neural network representations of the potential energy [34, 35] and
polarization surfaces [36] (see Methods), to boost numerical efficiency and
scalability. The neural networks are trained on DFT data generated with the
Strongly Constrained and Appropriately Normed (SCAN) functional approximation
[37]. Most importantly, we correct the excessive quantum delocalization of the
H/D ions by stiffening the springs of the associated ring polymers, while
leaving the other polymers unaffected. In practice, we multiply the spring
constant $k=\frac{mL}{\beta^{2}\hbar^{2}}$, with $m$ the ion mass and $\beta$
the inverse temperature, by a stiffening factor $\mu$ whose value is fixed by
requiring that the H/D off-centering in the $xy$ plane, i.e., $\Delta_{xy}$ in
Fig. 1(c), match the value estimated by neutron diffraction experiments [18,
38, 39, 40]. This procedure leads to $\mu=2.5$ and $\mu=16$, respectively, for
H and D ions (see Methods). This empirical fix leaves the DFT potential energy
surface unaffected, but lowers the zero point energy of H/D to compensate for
the underestimation of the double well barrier in the H-bonds, enhancing the
weight of the asymmetric O-H…O configurations. It turns out that static
position-dependent observables are predicted with remarkable accuracy with
this procedure. However, momentum-dependent observables, like the momentum
distribution measured in deep inelastic neutron scattering experiments [41,
42] are compromised and should not be calculated with this approach.

Using the above PIMD scheme we estimate the phase transition temperature
$T_{c}$ from the temperature evolution of the spontaneous polarization
$|\mathcal{P}_{z}|$, reported with the corresponding experimental data in Fig.
2(a). The estimated $T_{c}$ is $\approx 120$K for KDP, and $\approx 230$K for
DKDP, with an error relative to experiments of less than 10K in both cases.
The predicted saturated polarizations, for $T<T_{c}-20$K, closely match their
experimental counterparts. Notably, the experimental observation that in KDP
$|\mathcal{P}_{z}|$ is $\approx 25\%$ smaller than in DKDP, is reproduced by
the simulations, further supporting the validity of the model. Due to the
finite size of the simulations, which include 512 KDP/DKDP units, the observed
onset of spontaneous polarization is smoother than in experiments. A finite
size analysis can be found in the Supplementary Information (SI), where we
also report the calculated temperature dependence of the electric
susceptibility, which shows approximate Curie-Weiss behavior near $T_{c}$ for
both KDP and DKDP.

Experimentally, the PT in KDP/DKDP is first-order [43], structurally
associated with a shear distortion of the lattice in the ferroelectric phase.
The corresponding tetragonal to orthorhombic transformation is illustrated by
the changes of the $a$ and $b$ lattice constants in Fig. 2(b). The agreement
of theory and experiment is excellent. The models reproduce the measured
thermal expansion for $T>T_{c}$. For $T<T_{c}$, they slightly underestimate
the orthorhombic distortion, while predicting a less sharp transition than in
experiments due to finite-size effects. The shift of the $a$ and $b$ curves
upon deuteration is well described. In the SI, we also analyze the isotope
effects on the lattice constant $c$, the unit cell volume $\Omega$, and the
H-bond length $d_{\mathrm{O}...\mathrm{H}}$. The simulation reproduces the
Ubbelohde effect [44, 45], i.e., the elongation of the H-bonds following
deuteration, a fine effect in the order of the hundredths of an angstrom.

Figure 3:  (a) Schematic representation of the rules for the classification of
each KDP/DKDP replica. Each P-O-H (yellow-red-blue) structure can be donor
(left) or acceptor (right) type. The grey sphere indicates the location of the
electronic centroid. (b) The joint probability distribution of $\alpha$ and
$r_{\mathrm{O}\rightarrow\mathrm{H}}$ for KDP and DKDP, respectively, in the
ferroelectric ($T=T_{c}-10$K) and the paraelectric phase (insets,
$T=T_{c}+20$K). The vertical grey dashed line shows $r^{*}=1.25\text{\AA}$.
(c) Percent population of the different classes of KDP/DKDP units vs
temperature. The vertical grey dashed line shows $T_{c}$. (d) RDF of the
Takagi defects in KDP at 50K. The distance $r$ between the defects is the
spatial separation of the respective P ions.

The simulation gives access to microscopic properties that it would be
difficult, if not impossible, to observe in experiments, such as the
probability density function (PDF) of the local dipole moments
$\mathbf{d}=(d_{x},d_{y},d_{z})$. The distributions $n(d_{x},d_{z})\equiv
n(d_{x},0,d_{z})$ and $n(d_{x},d_{y})\equiv n(d_{x},d_{y},0)$ at two different
temperatures are displayed in Fig. 2(c) for DKDP and KDP, respectively. In
DKDP, breaking of the $z$-reflection symmetry of $n(d_{x},d_{z})$ for
$T<T_{c}$ is evident in the upper $\mathrm{c_{1}}$ panel. This symmetry is
restored for $T>T_{c}$ (upper $\mathrm{c_{2}}$ panel). The magnitude of
$d_{z}$ in the peak of the distribution in panel $\mathrm{c_{1}}$ closely
matches the dipole of the Ferro units (Fig. 1(d)). The peak magnitude of
$d_{z}$ does not change appreciably across the transition ($\mathrm{c_{2}}$
panel), as expected for a PT dominated by order-disorder. In the lower
$\mathrm{c_{1}}$ panel, $n(d_{x},d_{y})$ is uniformly close to zero,
reflecting the $z$ alignment of the dipoles in the ferroelectric phase. In the
paraelectric phase (lower $\mathrm{c_{2}}$ panel), $n(d_{x},d_{y})$ shows
eight off-center spots with octagonal symmetry, associated to the Slater units
of Fig. 1(d). Each one of the four Slater units in that figure contributes to
two off-center spots due to the glide reflection symmetry of the crystal. The
central peak in the lower $\mathrm{c_{2}}$ panel originates from the diffusive
part of the two dominant peaks of the upper $\mathrm{c_{2}}$ panel. The dipole
distributions of KDP are reported in panels $\mathrm{c_{3}}$ and
$\mathrm{c_{4}}$. Spontaneous symmetry breaking is again evident in
$n(d_{x},d_{z})$, with a dipole magnitude that does not change much across the
transition, consistent with order-disorder character. Relative to DKDP, the
KDP distributions have a more strongly diffusive nature, signaling the
presence of significantly more disorder above and below $T_{c}$. This is
particularly evident in the paraelectric $n(d_{x},d_{y})$ distribution, which,
in KDP, shows roughly circular symmetry (lower $\mathrm{c_{4}}$ panel). As we
will see, this is due to the prominent role played by Takagi units generated
by quantum fluctuations. While in DKDP Slater units dominate the paraelectric
phase near $T_{c}$, in KDP Takagi and Slater units are equally important,
blurring the octagonal pattern of the lower $\mathrm{c_{2}}$ panel.

In PIMD simulations each configuration of the beads at equal imaginary time
defines a distinct replica of the system (see Methods). Each KDP/DKDP unit in
a replica can be assigned to one of the four classes displayed in Fig. 1(d),
depending on the character of its four OH bonds, each one of which can have
donor or acceptor H-bond character. The length and character of each OH bond
is linearly correlated, to excellent approximation, with the length and
character of the OP bond and the position of the electronic centroid about the
oxygen (see Methods), as illustrated in Fig. 3(a). It is because of this
correlation that the dipole moment of a phosphate unit depends on the class to
which it belongs. In view of this correlation, the two types of P-O-H(D)
structure of Fig. 3(a) can be distinguished by the length of
$r_{\mathrm{O}\rightarrow\mathrm{H}}$. In the following, we choose the value
$r^{*}=1.25\text{\AA}$ for the divider between donor and acceptor type
structures for both KDP and DKDP below and above $T_{c}$. This choice is
supported by the calculated joint PDF of $\alpha$ and
$r_{\mathrm{O}\rightarrow\mathrm{H}}$, i.e.,
$n(\alpha,r_{\mathrm{O}\rightarrow\mathrm{H}})$, reported in Fig. 3(b) for
DKDP and KDP at $T=T_{c}-10$K, and at $T=T_{c}+20$K in the insets. In all
cases, $n(\alpha,r_{\mathrm{O}\rightarrow\mathrm{H}})$ is strongly bimodal
with the peaks corresponding to typical donor and acceptor type P-O-H(D)
structures, respectively. Moreover,
$n(\alpha,r_{\mathrm{O}\rightarrow\mathrm{H}})$ changes little across the PT.
In DKDP, the divider $r^{*}$, indicated by the grey dashed line in the figure,
separates well donor and acceptor configurations. In KDP, quantum fluctuations
make donor and acceptor distributions overlap, but symmetric H-bond
configurations, in which H lies approximately in the middle of the O-O bond,
have a substantially smaller statistical weight than asymmetric
configurations, even at the lowest temperature of the simulations ($T=50$K),
indicating that our choice of $r^{*}$ is meaningful.

The temperature variation of the population of phosphate units, classified
using $r^{*}$, is plotted in Fig. 3(c). For $T<T_{c}$, broken symmetry Ferro
units dominate both DKDP and KDP. Long-range order is lost in the paraelectric
phase where the Ferro units, while still present, are not dominant and the
majority of the units are Slater and/or Takagi. In this phase, the frequent
interconversions of Ferro, Slater and Takagi units, restore the
$z\rightarrow-z$ symmetry of the Ferro units. The population of Nonpolar units
is negligible at all temperatures. The most important message of Fig. 3(c) is
that Slater and Takagi populations are markedly different in DKDP and KDP. In
DKDP, when $T<150$K almost all units are Ferro. Dipolar defects, mostly
associated with Slater units, appear only for $T>150$K, causing a gradual
reduction of the spontaneous polarization. By contrast, in KDP, a substantial
fraction of Takagi, rather than the more stable Slater units, persists down to
$T=50K$, whereby Ferro units only account for $\approx 75\%$ of the
phosphates. On the time scale of the PIMD simulations, the Takagi units are
dynamic and fluctuate with H-bond switches every several picoseconds. By
contrast, H-bond switches require tens of picoseconds at low temperatures in
DKDP. Given that the molecular dipole of an ideal Ferro unit is only
marginally affected by deuteration, the roughly $25\%$ smaller polarization of
KDP relative to DKDP can only originate from the dipolar defects present in
KDP, and is mostly absent in DKDP at low temperatures. Dipolar disorder also
explains why in experiments the excess transition entropy is $\approx 30\%$
smaller in KDP than in DKDP [46].

The Takagi defects are spatially correlated as shown by their radial
distribution function (RDF) in KDP at $T=50$K in Fig. 3(d). Due to averaging
over the replicas in imaginary time, this RDF is already well converged in a
single PIMD snapshot. The first peak corresponds to pairs of adjacent defects
created mostly by switches of the H-bonds between neighboring Ferro units. The
weak temperature dependence of the Takagi fraction in Fig. 3(c) suggests that
the effect should originate from the quantum tunneling of the protons observed
in deep inelastic neutron scattering [41]. Interestingly, the same experiments
do not provide tunneling evidence in DKDP [42]. The prominent second peak of
the RDF indicates that an additional defect is likely to reside near a defect
pair, while the non-monotonic decay of the RDF with the Euclidean distance
reflects a monotonic decay with the graph distance defined as the minimal
number of O-H…O bridges connecting pairs of Takagi defects. The defects are
also dynamically correlated as suggested by analyses of their joint PDF
reported in the SI. Recombination and migration to a neighboring site are two
sources of dynamical correlation. The Takagi defects are nominally charged and
correlation acts to reduce local charge fluctuations.

Because of quantum disorder, the ferroelectric state of KDP is fundamentally
different from that of DKDP. Linear extrapolation in Fig. 3(c) suggests that a
substantial fraction ($\approx 15\%$) of fluctuating Takagi defects should be
present near absolute zero, providing a possible realization of the TLS that
explains the anomalous heat capacity of insulating glasses near $T=0$K [47,
48]. This speculation is supported by the experimental findings of an
anomalous heat capacity varying approximately linearly with $T$ near absolute
zero in KDP [49, 5, 6]. The TLS in KDP should be an intrinsic equilibrium
property rather than the outcome of quenched glassy disorder, but it would be
difficult to test this hypothesis as it is hard to control the various sources
of extrinsic disorder present in the experimental samples [49].

The huge effect of quantum fluctuations in KDP requires a fine balance between
the tunneling barrier and the zero point energy. This balance, which is broken
by deuteration with the ensuing large isotope effects, is realized in the KDP
family of crystals due to the presence of unusually strong H-bonds. H-bonds
can be made stronger by pressure. In water ice, where the H-bonds at ambient
conditions are relatively strong but significantly weaker than in KDP, isotope
effects are relatively small, but they become large at pressures in the
hundreds of kbar range [50]. In KDP, experiments show that the temperature of
the ferroelectric transition is reduced with increasing pressure, until, at
pressure in excess of 17kbar the PT disappears and the system becomes a
quantum paraelectric [51, 52]. We understand this phenomenology as the result
of an increasing Takagi population following increasing H-bond symmetrization
with pressure. We speculate that the phase point at which the PT disappears
($T=0$K, $P=17$kbar) may be associated to quantum criticality [53, 54], an
hypothesis consistent with the large dielectric susceptibility measured in the
vicinity of this point [52]. In future studies, it would be of interest to
look for possible non-classical temperature dependence of the inverse
dielectric susceptibility at very low temperature [54].

## Methods

### Density functional theory and atomistic models

We adopt SCAN, the strongly constrained and appropriately normed functional
approximation for exchange and correlation of Ref. [37]. It has been reported
that SCAN outperforms PBE and LDA in modeling the structural properties of KDP
[55]. We perform electronic structure calculations with Quantum Espresso [56],
using norm-conserving pseudopotentials [57] that include in the valence the
$2s$ and $2p$ states of O, the $3s$ and $3p$ states of P, and the $3s$, $3p$,
and $4s$ states of K. We find a stable classical ferroelectric structure of
space group Fdd2 at $T=0\text{K}$ (Fig.1(a)), with lattice constants
$a=10.62\text{\AA}$, $b=10.67\text{\AA}$, and $c=6.86\text{\AA}$. The length
of the O-H covalent bond is $1.042\text{\AA}$, that of the O…H H-bond is
$1.438\text{\AA}$, and the H off-centering is equal to
$\Delta_{xy}^{\text{Fdd2}}=0.198\text{\AA}$. The above results agree well with
Ref. [55] with minor differences due to the different pseudopotentials and
numerical methods adopted in the two calculations.

We use DeePMD-Kit [35] to train a DP model [34] for the potential energy
surface within SCAN-DFT. The typical error of the model relative to DFT is
less than $1\text{meV/atom}$. Training details and error distributions are
reported in the SI. In addition, we train a Deep Wannier (DW) model [36] for
the polarization surface in terms of the Wannier decomposition of the valence
electronic structure [26]. Maximally localized Wannier functions are computed
with Wannier90 [58] for a subset of the atomic configurations in the DP
training set (see SI). Each Wannier distribution accommodates two spin-
degenerate electrons. In all the atomic configurations four Wannier centers
are uniquely associated to each O and four others to each K atom, forming
approximate atom-centered tetrahedral structures. The Wannier centroids (WCs),
i.e., the geometric centers of the four Wannier centers associated to O or K
atoms provide all the necessary information to compute the polarization [36].
The DW model gives the environmental dependence of the WCs. While the oxygen
WCs depend crucially on the chemical environment, as illustrated in Fig. 3(a),
the potassium WCs, associated with the $3s$ and $3p$ semi-core states, are
essentially independent of the environment and their contribution to
polarization changes is negligible. We do not include them in the DW model,
i.e., to calculate the polarization we treat the semi-core states as frozen
core electrons of the $\mathrm{K^{+}}$ ion. The training of the DW model is
done with DeePMD-Kit. The details are reported in the SI.

Following the theory of polarization [26, 25], we assign to each elementary
KDP/DKDP unit ( Fig. 1(b)) a molecular dipole $\mathbf{d}$ given by:

$\begin{split}\mathbf{d}=\sum_{i=1}^{2}q_{\mathrm{K}}(\mathbf{r}_{\mathrm{K}}^{(i)}-\mathbf{r}_{\mathrm{P}})+\sum_{i=1}^{4}q_{\mathrm{O}}(\mathbf{r}_{\mathrm{O}}^{(i)}-\mathbf{r}_{\mathrm{P}})\\\
+\sum_{i=1}^{4}q_{\mathrm{H}}(\mathbf{r}_{\mathrm{H}}^{(i)}-\mathbf{r}_{\mathrm{P}})+\sum_{i=1}^{4}q_{\mathrm{WC}}(\mathbf{r}_{\mathrm{WC}}^{(i)}-\mathbf{r}_{\mathrm{P}}).\end{split}$
(1)

Here, $\mathbf{r}_{\mathrm{K}}^{(i)}$, $\mathbf{r}_{\mathrm{O}}^{(i)}$,
$\mathbf{r}_{\mathrm{H}}^{(i)}$, and $\mathbf{r}_{\mathrm{P}}$ are the
position vectors of the K, O, H/D, and P ions belonging to the elementary
unit, and $q_{\mathrm{K}}=0.5e$, $q_{\mathrm{O}}=6e$, and
$q_{\mathrm{H}}=0.5e$ are the charges carried by the K, O, and H/D ions,
respectively. Since the K and H/D ions are shared with a neighboring unit,
they carry half charge. The P ion, carrying a charge of $5e$, is taken as the
local reference and does not contribute to the dipole. Finally,
$\mathbf{r}_{\mathrm{WC}}^{(i)}$ are the position vectors of the WCs
associated with the O ions and $q_{\mathrm{WC}}=-8e$ are their charges. When
all the ionic and electronic charges are taken into account, each elementary
unit is electrically neutral. The positions of the ions and of the WCs
corresponding to each atomic configuration
$\mathbf{R}=\\{\mathbf{r}_{1},\cdots,\mathbf{r}_{N}\\}$ visited by PIMD in a
periodic simulation box of volume $V$ and $N$ physical particles are provided
by the DP and DW models. The global dipole of the configuration $\mathbf{R}$
is $\mathbf{D}(\mathbf{R})=\sum_{j}\mathbf{d}_{j}$, where the sum extends to
all the local molecular dipoles $\mathbf{d}_{j}$ contained in the simulation
box. The polarization surface $\boldsymbol{\mathcal{P}}(\mathbf{R})$, i.e.,
the polarization of the configuration $\mathbf{R}$, is given, modulo an
immaterial constant, by
$\boldsymbol{\mathcal{P}}(\mathbf{R})=\mathbf{D}(\mathbf{R})/V$ [26].

### Path integral molecular dynamics

Within Feynman’s path integral formulation of quantum statistical mechanics
[59], the canonical partition function of $N$ distinguishable quantum
particles can be approximated with a Trotter factorization in imaginary time:

$Z_{L}\approx\int\mathcal{D}\mathbf{R}e^{-\beta\mathcal{H}^{\text{eff}}_{L}},$
(2)

where

$\mathcal{D}\mathbf{R}\equiv\prod_{i=1}^{N}\left(\frac{m_{i}L}{2\pi\beta\hbar^{2}}\right)^{\frac{3L}{2}}\text{d}\mathbf{r}_{i}^{(1)}\cdots\text{d}\mathbf{r}_{i}^{(L)}.$
(3)

The effective Hamiltonian $\mathcal{H}^{\text{eff}}_{L}$ is

$\mathcal{H}^{\text{eff}}_{L}=\sum^{L}_{s=1}\Bigg{(}\sum_{i=1}^{N}\frac{k_{i}}{2}\left(\mathbf{r}^{(s)}_{i}-\mathbf{r}^{(s+1)}_{i}\right)^{2}+\frac{1}{L}U({\mathbf{R}}^{(s)})\Bigg{)},$
(4)

where $k_{i}=\frac{m_{i}L}{\beta^{2}\hbar^{2}}$ and
$\mathbf{R}^{(s)}=\\{\mathbf{r}_{1}^{(s)},\cdots,\mathbf{r}_{N}^{(s)}\\}$. The
index $s$ labels imaginary times, the positive integer $L$ is the total number
of imaginary time slices, and the condition $L+1=1$ applies in the argument of
the sum over the beads so that $\mathbf{r}^{(L+1)}_{i}=\mathbf{r}^{(1)}_{i}$.
Thus, a discretized Feynman path in imaginary time,
$\\{\mathbf{r}_{i}^{(1)},\cdots,\mathbf{r}_{i}^{(L)}\\}$, is equivalent to a
ring polymer. $U$ is the potential energy surface provided by the DP model.
$Z_{L}$ recovers the exact quantum partition function in the limit
$L\rightarrow\infty$. For finite $L$, Eq. (2) can be seen as the
configurational partition function of a classical system of $N$ ring polymers
each made of $L$ beads and described by the effective Hamiltonian
$\mathcal{H}^{\text{eff}}_{L}$, in which the constant $k_{i}$ defines the
spring that links neighboring beads in imaginary time, and
$\frac{1}{L}U({\mathbf{R}}^{(s)})$ is the potential of interaction between
different ring polymers at equal imaginary time. The ensemble of all the beads
at equal imaginary time constitutes a replica of the system. PIMD [60] is a
molecular dynamics scheme for sampling the configurations of the ring polymers
with Boltzmann weights proportional to
$e^{-\beta\mathcal{H}^{\text{eff}}_{L}}$. Then, a canonical average of a
position-dependent observable of the quantum system is estimated from a
classical canonical average of the ring-polymer system. For example, the
average macroscopic polarization is calculated as follows:

$\begin{split}\boldsymbol{\mathcal{P}}=\frac{1}{L}\left\langle\sum_{s=1}^{L}\boldsymbol{\mathcal{P}}(\mathbf{R}^{(s)})\right\rangle_{Z_{L}}=\\\
\frac{\int\mathcal{D}\mathbf{R}\sum_{s=1}^{L}\boldsymbol{\mathcal{P}}(\mathbf{R}^{(s)})e^{-\beta\mathcal{H}^{\text{eff}}_{L}}}{L\int\mathcal{D}\mathbf{R}e^{-\beta\mathcal{H}^{\text{eff}}_{L}}},\end{split}$
(5)

where $\boldsymbol{\mathcal{P}}(\mathbf{R}^{(s)})$, the polarization of the
$s$-th replica, is provided by the DW model.

### PIMD simulations

We use I-PI [61], LAMMPS [62], and DeePMD-Kit to perform NPT-PIMD simulations
on a supercell with 512 KDP/DKDP units in the temperature range
$[50,300]\text{K}$. The temperature is controlled with the PILE thermostat
[63] and the pressure is maintained with an anisotropic barostat [64]. The
spring constants $k_{i}$ (Eq. (4)) associated with K, P, and O ions take their
physical values, but the $k_{i}$ associated with H is multiplied by a
stiffening factor $\mu=2.5$, while the one corresponding to D is multiplied by
$\mu=16$. These values are chosen to fit the experimental neutron data on the
off-center displacement of H/D, as explained in the next section.

To control convergence with $L$, we use the estimator $K_{L}$ of the
thermodynamic quantum kinetic energy per particle suggested in Ref. [65]:

$K_{L}=\frac{3}{2\beta}+\frac{1}{2LN}\Bigg{\langle}\sum_{s=1}^{L}\sum_{i=1}^{N}\left({\mathbf{r}_{i}^{(s)}}-\overline{{\mathbf{r}_{i}}}\right)\cdot\frac{\partial
U}{\partial{\mathbf{r}_{i}^{(s)}}}\Bigg{\rangle}_{Z_{L}}.$ (6)

Here, $\overline{{\mathbf{r}_{i}}}$ is the centroid of all the beads of
particle $i$. In the classical limit, $L=1$ and $K_{L}=3k_{B}T/2$, the
classical kinetic energy per particle. By monitoring the convergence of
$K_{L}$ with $L$ (see SI) we estimate $|K_{L}-K_{\infty}|$ at different
temperatures. In the production runs we adopt $L$ values for which
$|K_{L}-K_{\infty}|$ is smaller than $1\text{meV/atom}$, which is also the
approximate error of the DP representation. A study of the effect of finite
$N$ is impractical with PIMD simulations as $L$ can be as large as $128$. For
that reason, we only report in the SI a study of the effect of $N$ on the
corresponding classical system ($L=1$).

Figure 4: (Upper) The relation between $T_{c}$ and H/D off-centering. The
black stars mark the experimental measurements of $T_{c}$ and
$\tilde{\Delta}_{xy}^{\mathrm{EXP}}$ for KDP and DKDP at different external
pressures [18, 38, 39, 40]. Blue spheres mark the numerical $T_{c}$ and
$\Delta_{xy}(T\rightarrow T_{c}^{+})$ obtained with PIMD with different
$m^{*}$ at atmospheric pressure. The error bar of $T_{c}$ is $\pm 10$K,
reflecting the $T$ spacing adopted in the calculations of the order parameter
near $T_{c}$. The error bar of $\Delta_{xy}(T\rightarrow T_{c}^{+})$ is $\pm
0.005\text{\AA}$, reflecting the space resolution of the histogram
representing $n(\boldsymbol{\delta}_{xy})$ from which we extract
$\Delta_{xy}(T\rightarrow T_{c}^{+})$. Blue shadows are a guide for the eye,
suggesting an almost linear $T_{c}$-$\Delta_{xy}(T\rightarrow T_{c}^{+})$
relation. (Lower) The heatmap of $n(\boldsymbol{\delta}_{xy})$ computed at
$T>T_{c}$. The left figure is for $m^{*}=2.5$Da and $T=150$K. The right figure
is for $m^{*}=32$Da and $T=250$K.

### The spring stiffening factor $\mu$ for H and D ions

In KDP, SCAN-DFT predicts a ferroelectrically ordered classical ground state
[66], but this structure is destabilized by NQE. Absence of ferroelectricity
was also found in ab initio PIMD simulations of KDP with the PBE-DFT
functional [31]. This occurs because common functionals tend to underestimate
the double-well potential energy barrier in the H-bonds, wiping out the
classical ferroelectric ground state [28].

In the presence of NQE, a correct description of the ferroelectric phase
requires a delicate balance of tunneling barrier and zero-point energy of the
protons, which are both strongly correlated with the geometric distortion of
the phosphate groups, as shown in numerous experimental [18, 38, 39, 40, 19]
and theoretical [67, 27, 68] studies.

Experiments find that, in KDP and DKDP, $T_{c}$ at different pressure
conditions is to very good approximation linearly correlated with
$\tilde{\Delta}_{xy}^{\mathrm{EXP}}$, the off-center displacement of the H/D
ions in the paraelectric phase just above the PT [18, 38, 39, 40]. This
correlation, shown in Fig. 4, strongly supports the notion that H/D off-
centering is the key geometric parameter controlling $T_{c}$, and hints at a
possible way of improving the experimental agreement of DFT-based PIMD
simulations by fixing the predicted H/D off-centering. In KDP, SCAN-based PIMD
always yields unimodal proton distributions with $\Delta_{xy}(T)=0$, but the
excessive quantum delocalization of H/D, responsible for the effect, can be
reduced as desired by multiplying the H/D spring constant
$k=\frac{mL}{\beta^{2}\hbar^{2}}$ by a stiffening factor $\mu$ greater than 1,
a procedure equivalent to assigning an effective mass $m^{*}=\mu m$ to H/D.
The displacements $\Delta_{xy}(T\rightarrow T_{c}^{+})$ and the corresponding
$T_{c}$ calculated with PIMD for three different values of $m^{*}$ at ambient
pressure (AP) are reported in Fig. 4, showing a linear behavior of
$\Delta_{xy}$ with $T_{c}$ that closely mimics the experimental one. For
$m^{*}=2.5\text{Da}$, i.e. $\mu=2.5$, and for $m^{*}=32\text{Da}$, i.e.,
$\mu=16$, $\Delta_{xy}(T\rightarrow T_{c}^{+})$ and $T_{c}$ match the
experimental displacements and transition temperatures of KDP and DKDP,
respectively. This is the choice of $\mu$ that we adopt in our PIMD
simulations of KDP and DKDP.

The displacement $\Delta_{xy}(T\rightarrow T_{c}^{+})$ is estimated from
thermal averages at $T\approx T_{c}+20$K to reduce the large impact of the
fluctuations near $T_{c}$. The probability distribution functions
$n(\boldsymbol{\delta}_{xy})$ for H and D are also reported in Fig. 4, showing
a bimodal character in qualitative agreement with neutron diffraction
experiments [2]. Relative to KDP, DKDP has a more localized
$n(\boldsymbol{\delta}_{xy})$ and a larger $\Delta_{xy}$, consistent with the
more classical behavior of the deuterated H-bonds. In fact, in the classical
limit, when all the atoms in PIMD have infinite mass, $T_{c}$ is approximately
equal to $250$K, slightly above the predicted transition temperature of DKDP,
consistent with a predominantly classical character of the transition in that
system.

## Acknowledgements

We thank Yifan Li, Chunyi Zhang, Xifan Wu, Chenyuan Li, Chenxing Luo, and
Linfeng Zhang for fruitful discussions. All authors were supported by the
Computational Chemical Sciences Center: Chemistry in Solution and at
Interfaces (CSI) funded by DOE Award DE-SC0019394. The authors are pleased to
acknowledge that the work reported in this paper was performed largely using
the Princeton Research Computing resources at Princeton University. This
research also used resources of the National Energy Research Scientific
Computing Center (NERSC) operated under Contract No. DE-AC02-05CH11231 using
NERSC award ERCAP0021510.

## References

  * [1] Samara, G. A.  The effects of deuteration on the static ferroelectric properties of kh2p04 (kdp).  _Ferroelectrics_ 5, 25–37 (1973).
  * [2] Nelmes, R.  Structural studies of kdp and the kdp-type transition by neutron and x-ray diffraction: 1970–1985.  _Ferroelectrics_ 71, 87–123 (1987).
  * [3] Tokunaga, M. & Matsubara, T.  Review on tunneling model for kh2po4.  _Ferroelectrics_ 72, 175–191 (1987).
  * [4] Schmidt, V. H.  Review of order-disorder models for kdp-family crystals.  _Ferroelectrics_ 72, 157–173 (1987).
  * [5] Lawless, W. & Lawless, T.  Specific heat studies of small kdp crystals at low temperatures.  _Ferroelectrics_ 45, 149–155 (1982).
  * [6] Foote, M. & Anderson, A.  Low-temperature specific heats of batio3 and kh2po4.  _Ferroelectrics_ 62, 11–15 (1985).
  * [7] Lines, M. E. & Glass, A. M.  _Principles and applications of ferroelectrics and related materials_ (Oxford university press, 2001).
  * [8] Mason, W.  The elastic, piezoelectric, and dielectric constants of potassium dihydrogen phosphate and ammonium dihydrogen phosphate.  _Phys. Rev._ 69, 173 (1946).
  * [9] Wang, D. _et al._ Characteristics of nonlinear optical absorption and refraction for kdp and dkdp crystals.  _Opt. Mater. Express_ 7, 533–541 (2017).
  * [10] Baisden, P. _et al._ Large optics for the national ignition facility.  _Fusion Sci. Technol._ 69, 295–351 (2016).
  * [11] Tressler, J. F., Alkoy, S. & Newnham, R. E.  Piezoelectric sensors and sensor materials.  _J. Electroceram._ 2, 257–272 (1998).
  * [12] Li, X.-Z., Walker, B. & Michaelides, A.  Quantum nature of the hydrogen bond.  _Proc. Natl. Acad. Sci._ 108, 6369–6373 (2011).
  * [13] Rabe, K. M., Ahn, C. H. & Triscone, J.-M.  _Physics of ferroelectrics: a modern perspective_ , vol. 105 (Springer Science & Business Media, 2007).
  * [14] Blinc, R.  On the isotopic effects in the ferroelectric behaviour of crystals with short hydrogen bonds.  _J. Phys. Chem. Sol._ 13, 204–211 (1960).
  * [15] De Gennes, P.  Collective motions of hydrogen bonds.  _Solid State Commun._ 1, 132–137 (1963).
  * [16] Tominaga, Y., Kasahara, M., Urabe, H. & Tatsuzaki, I.  Raman spectral of low-lying longitudinal mode of kh2po4 and kd2po4 in ferroelectric phase.  _Solid State Commun._ 47, 835–837 (1983).
  * [17] Tokunaga, M. & Tatsuzaki, I.  Light scattering spectra of polarization fluctuations and models of the phase transition in kdp type ferroelectrics.  _Phase Transitions_ 4, 97–155 (1984).
  * [18] Nelmes, R., Tun, Z. & Kuhs, W.  A compilation of accurate structural parameters for kdp and dkdp, and a users’ guide to their crystal structures.  _Ferroelectrics_ 71, 125–141 (1987).
  * [19] McMahon, M. _et al._ Geometric effects of deuteration on hydrogen-ordering phase transitions.  _Nature_ 348, 317–319 (1990).
  * [20] Řezina, B., Smutnỳ, F. & Fousek, J.  Ferroelectricity in pbhaso4 and pbdaso4.  _Czech. J. Phys. B_ 25, 1411–1412 (1975).
  * [21] Wang, X., Tian, T. & Aschauer, U.  Beautiful Atoms: A python library for the manipulation and visualization of atomistic structure using blender (2022).  URL https://github.com/beautiful-atoms/beautiful-atoms.
  * [22] Slater, J. C.  Theory of the transition in kh2po4.  _J. Chem. Phys._ 9, 16–33 (1941).
  * [23] Takagi, Y.  Theory of the transition in kh2po4.(i).  _Phys. Soc. Jpn._ 3, 271–272 (1948).
  * [24] Takagi, Y.  Theory of the transition in kh2po4,(ii).  _Phys. Soc. Jpn._ 3, 273–274 (1948).
  * [25] Vanderbilt, D.  _Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators_ (Cambridge University Press, 2018).
  * [26] Resta, R. & Vanderbilt, D.  Theory of polarization: a modern approach.  In _Physics of Ferroelectrics_ , 31–68 (Springer, 2007).
  * [27] Koval, S., Kohanoff, J., Migoni, R. & Tosatti, E.  Ferroelectricity and isotope effects in hydrogen-bonded kdp crystals.  _Phys. Rev. Lett._ 89, 187602 (2002).
  * [28] Menchón, R. _et al._ Ab initio study of the structure, isotope effects, and vibrational properties in kdp crystals.  _Phys. Rev. B_ 98, 104108 (2018).
  * [29] Marx, D. & Parrinello, M.  Ab initio path integral molecular dynamics: Basic ideas.  _J. Chem. Phys._ 104, 4077–4082 (1996).
  * [30] Kohn, W. & Sham, L. J.  Self-consistent equations including exchange and correlation effects.  _Phys. Rev._ 140, A1133 (1965).
  * [31] Srinivasan, V. & Sebastiani, D.  The isotope-effect in the phase transition of kh2po4: new insights from ab initio path-integral simulations.  _J. Phys. Chem. C_ 115, 12631–12635 (2011).
  * [32] Perdew, J. P. & Zunger, A.  Self-interaction correction to density-functional approximations for many-electron systems.  _Phys. Rev. B_ 23, 5048 (1981).
  * [33] Mori-Sánchez, P., Cohen, A. J. & Yang, W.  Localization and delocalization errors in density functional theory and implications for band-gap prediction.  _Phys. Rev. Lett._ 100, 146401 (2008).
  * [34] Zhang, L., Han, J., Wang, H., Car, R. & Weinan, E.  Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics.  _Phys. Rev. Lett._ 120, 143001 (2018).
  * [35] Wang, H., Zhang, L., Han, J. & Weinan, E.  Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics.  _Comput. Phys. Commun._ 228, 178–184 (2018).
  * [36] Zhang, L. _et al._ Deep neural network for the dielectric response of insulators.  _Phys. Rev. B_ 102, 041121 (2020).
  * [37] Sun, J., Ruzsinszky, A. & Perdew, J. P.  Strongly constrained and appropriately normed semilocal density functional.  _Phys. Rev. Lett._ 115, 036402 (2015).
  * [38] Tibballs, J., Nelmes, R. & McIntyre, G.  The crystal structure of tetragonal kh2po4 and kd2po4 as a function of temperature and pressure.  _J. Phys. C: Solid State Phys._ 15, 37 (1982).
  * [39] Tibballs, J. & Nelmes, R.  The pt dependence of the crystal structure of kdp and dkdp above tc.  _J. Phys. C: Solid State Phys._ 15, L849 (1982).
  * [40] Nelmes, R.  On the structural evidence for a direct proton tunneling effect in the kh2po4-type transition.  _J. Phys. C: Solid State Phys._ 21, L881 (1988).
  * [41] Reiter, G., Mayers, J. & Platzman, P.  Direct observation of tunneling in kdp using neutron compton scattering.  _Phys. Rev. Lett._ 89, 135505 (2002).
  * [42] Reiter, G., Shukla, A., Platzman, P. & Mayers, J.  Deuteron momentum distribution in kd2po4.  _New J. Phys._ 10, 013016 (2008).
  * [43] Reese, W.  Studies of phase transitions in order-disorder ferroelectrics. iii. the phase transition in k h 2 p o 4 and a comparison with k d 2 p o 4.  _Phys. Rev._ 181, 905 (1969).
  * [44] Robertson, J. M. & Ubbelohde, A. R. J. P.  Structure and thermal properties associated with some hydrogen bonds in crystals i. the isotope effect.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 170, 222–240 (1939).
  * [45] Ubbelohde, A. & Woodward, I.  Isotope effect in potassium dihydrogen phosphate.  _Nature_ 144, 632–632 (1939).
  * [46] Strukov, B., Amin, M. & Kopchik, V.  Comparative investigation of the specific heat of kh2po4 (kdp) and kd2po4 (dkdp) single crystals.  _Phys. Status Solidi B_ 27, 741–749 (1968).
  * [47] Anderson, P. W., Halperin, B. I. & Varma, C. M.  Anomalous low-temperature thermal properties of glasses and spin glasses.  _Philos. Mag._ 25, 1–9 (1972).
  * [48] Phillips, W. A.  Two-level states in glasses.  _Rep. Prog. Phys._ 50, 1657 (1987).
  * [49] Lawless, W. N.  ${T}^{\frac{3}{2}}$ contribution to the specific heat of ferroelectrics at low temperatures.  _Phys. Rev. Lett._ 36, 478–479 (1976).  URL https://link.aps.org/doi/10.1103/PhysRevLett.36.478.
  * [50] Pruzan, P.  Pressure effects on the hydrogen bond in ice up to 80 gpa.  _J. Mol. Struct._ 322, 279–286 (1994).
  * [51] Samara, G.  Pressure dependence of the static and dynamic properties of kh2po4 and related ferroelectric and antiferroelectric crystals.  _Ferroelectrics_ 71, 161–182 (1987).
  * [52] Endo, S., Deguchi, K. & Tokunaga, M.  Quantum paraelectricity in ${\mathrm{kd}}_{2}{\mathrm{po}}_{4}$ and ${\mathrm{kh}}_{2}{\mathrm{po}}_{4}$ under high pressure.  _Phys. Rev. Lett._ 88, 035503 (2002).
  * [53] Sachdev, S.  Quantum phase transitions.  _Phys. World_ 12, 33 (1999).
  * [54] Rowley, S. _et al._ Ferroelectric quantum criticality.  _Nat. Phys._ 10, 367–372 (2014).
  * [55] Zhang, Y., Sun, J., Perdew, J. P. & Wu, X.  Comparative first-principles studies of prototypical ferroelectric materials by lda, gga, and scan meta-gga.  _Phys. Rev. B_ 96, 035143 (2017).
  * [56] Giannozzi, P. _et al._ Quantum espresso: a modular and open-source software project for quantum simulations of materials.  _J. Condens. Matter Phys._ 21, 395502 (2009).
  * [57] Hamann, D.  Optimized norm-conserving vanderbilt pseudopotentials.  _Phys. Rev. B_ 88, 085117 (2013).
  * [58] Pizzi, G. _et al._ Wannier90 as a community code: new features and applications.  _J. Condens. Matter Phys._ 32, 165902 (2020).
  * [59] Feynman, R. P., Hibbs, A. R. & Styer, D. F.  _Quantum mechanics and path integrals_ (Courier Corporation, 2010).
  * [60] Tuckerman, M.  _Statistical mechanics: theory and molecular simulation_ (Oxford university press, 2010).
  * [61] Kapil, V. _et al._ i-pi 2.0: A universal force engine for advanced molecular simulations.  _Comput. Phys. Commun._ 236, 214–223 (2019).
  * [62] Thompson, A. P. _et al._ Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  _Comput. Phys. Commun._ 271, 108171 (2022).
  * [63] Ceriotti, M., Parrinello, M., Markland, T. E. & Manolopoulos, D. E.  Efficient stochastic thermostatting of path integral molecular dynamics.  _J. Chem. Phys._ 133, 124104 (2010).
  * [64] Martyna, G. J., Hughes, A. & Tuckerman, M. E.  Molecular dynamics algorithms for path integrals at constant pressure.  _J. Chem. Phys._ 110, 3275–3290 (1999).
  * [65] Ceriotti, M. & Manolopoulos, D. E.  Efficient first-principles calculation of the quantum kinetic energy and momentum distribution of nuclei.  _Phys. Rev. Lett._ 109, 100604 (2012).
  * [66] Zhong, W., Vanderbilt, D. & Rabe, K. M.  First-principles theory of ferroelectric phase transitions for perovskites: The case of ${\mathrm{batio}}_{3}$.  _Phys. Rev. B_ 52, 6301–6312 (1995).
  * [67] Ichikawa, M., Amasaki, D., Gustafsson, T. & Olovsson, I.  Structural parameters determining the transition temperature of tetragonal kh 2 po 4-type crystals.  _Phys. Rev. B_ 64, 100101 (2001).
  * [68] Koval, S., Kohanoff, J., Lasave, J., Colizzi, G. & Migoni, R.  First-principles study of ferroelectricity and isotope effects in h-bonded kh2po4 crystals.  _Phys. Rev. B_ 71, 184102 (2005).