# Quadrilateral Particle Arrangement within Shocks in a Two-Dimensional Dusty
Plasma

Anton Kananovich  Department of Physics and Astronomy, Appalachian State
University, Boone, North Carolina 28608, USA  Department of Physics and
Astronomy, University of Iowa, Iowa City, Iowa 52242, USA  J. Goree
Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa
52242, USA

###### Abstract

The microscopic structure within a two-dimensional shock was studied using
data from a dusty plasma experiment. A single layer of charged microparticles,
levitated in a glow-discharge plasma, was perturbed by an electrically
floating wire that was moved at a steady supersonic speed to excite a
compressional shock. A rearrangement of particles was observed, from a
hexagonal lattice in the preshock into a quadrilateral microstructure within
the shock. This quadrilateral structure would not be stable in a monolayer of
identical repulsive particles, under equilibrium conditions. Glaser-Clark
polygon analysis of the microstructure helped in identifying quadrilaterals.
Voronoi analysis was used to characterize the defect fraction behind the
shock, as an indication of shock-induced melting.

##  I Introduction

Shock waves (shocks), which can be defined as disturbances propagating at
supersonic speed, exist in all kinds of media, including gases, solids, and
plasmas. One type of plasma, which can also sustain shocks, is two-dimensional
(2D) dusty plasma. This type of plasma consists of a layer of microparticles
levitating in the sheath above the lower electrode in the gas discharge [1,
2]. The ground state for such a 2D layer is a triangular lattice with six-fold
symmetry.

Shocks in 2D dusty plasma have recently seen an increase in interest,
including experimental and molecular dynamics (MD) simulations studies [3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 3, 4, 5, 6, 7, 10,
11, 12, 13, 15, 16, 17, 18, 19, 20]. Among recently reported MD simulations,
we can mention studies of shock-front speed and hydrodynamic quantities [7],
hydrodynamic quantities and velocity distribution functions along with
analytical expression between the thermal and the drift velocities [10], the
evolution of blast waves [15], and a head-on collision of shocks [16].

Due to its unique properties, 2D dusty plasma is an excellent medium for
experimentally studying the microscopic structure of shocks. The sound speed
in that medium is the order of centimeters per second so that video microscopy
diagnostics allow the resolving and tracking of individual microparticles
[21]. These particle-level measurements allow experimental study of phenomena
such as the finite thickness of the shock front [14], and they also allow
direct comparisons to MD simulations. The results may also be applicable to
other 2D physical systems that similarly exhibit an equilibrium state with
hexagonal structure [20].

The microscopic properties of shocks in 2D dusty plasma have been studied
using the method of Voronoi diagrams, with an input of particle-position data
obtained from MD simulations. In Ref. [12], the authors studied melting, which
was reported to occur when the moving exciter had a speed exceeding a
threshold, as determined by Voronoi diagrams that were examined qualitatively
and were also used to quantify defect fraction. A leap of the defect ratio,
which can indicate the occurrence of a phase transition from the solid to
liquid states, was observed. In Ref. [13], a postshock region of 2D crystal,
with both liquid and solid coexisting regions, was observed when the exciter
speed was sufficiently high. For this phenomenon, again Voronoi diagrams were
used to examine disorder and quantify defect fraction.

In this paper, we report an experimental investigation of a point that has not
been much studied yet: the microscopic structure at the shock front and
immediately behind it. We find that the six-fold triangular geometry of a
ground-state crystal does not persist, under the sudden compression of a
shock, but instead the particles tend to rearrange into a quadrilateral
microscopic structure. The quadrilateral microstructure is unusual, because
for a 2D crystal with isotropic repulsion the triangular (also called
hexagonal) lattice is the only close-packed lattice in two dimensions [22],
under the conditions of steady equilibrium. A quadrilateral lattice in 2D
dusty plasma under steady conditions can exist only in special cases, such as
binary mixtures [23] or when the crystal buckles [24, 25]. In this paper, we
observe quadrilateral microscopic structure under the conditions of shock
compression that are not a steady equilibrium.

To study this microscopic structure, we analyze our experimental particle-
position data using the methods of Voronoi diagrams [26] and the polygon
construction method of Glaser and Clark [27, 28, 29]. The polygon construction
method is particularly effective in revealing the quadrilateral microscopic
structure that we report, for the nonequilibrium conditions of shock
compression in a 2D substance.

##  II Method

###  II.1 Experiment

In this paper, we report a further analysis of experimental runs that we
previously reported Ref. [9]. In that experiment, as an exciter was moved at a
steady speed, there was a localized compression of the 2D dust layer in front
of the exciter. This compression was strong and supersonic, indicating that it
was a shock. The main result in that previous report was the relationship of
two speeds, for the shock itself and for the exciter. Here we analyze the
particle-position data for a different purpose, to assess the microscopic
structure, particularly within the shock and immediately behind it. We briefly
review that experiment next; further details are provided in Ref. [9].

The experiment was performed in our Kuda-Topf chamber [9, 30, 31]. To sustain
a gas discharge in argon, 13.56 MHz radio-frequency (rf) power with a peak-to-
peak voltage of 184 V was applied to the lower electrode through a coupling
capacitor, while the outer wall of the chamber was grounded. As a result, a dc
bias of -85.5 V developed on the lower electrode. A 2D layer of monodisperse
(8.69 µm in diameter) melamine-formaldehyde microparticles was levitated in
the plasma sheath above the lower electrode. The particle charge was
$-1.4\times 10^{4}$ elementary charges. In the absence of the perturbation of
a shock, the lattice constant was $b=0.80$ mm, and correspondingly the 2D
Wigner-Seitz radius was $a=0.42$ mm, while the nominal 2D dust plasma
frequency was $\omega_{pd}=49$ s-1.

Compressional shocks were generated by continuously moving an exciter at a
controlled velocity. The exciter was a nickel wire, shown in the Supplemental
Material 111The Supplemental Material for this preprint will be published,
along with the article itself, in a future journal publication.. This exciter
wire was aligned perpendicular to its direction of motion but parallel to the
2D microparticle cloud. Runs were performed at five different speeds of the
exciter at various speeds, ranging from 50.8 to 101.6 mm/s all above the sound
speed of 19 mm/s. In dimensionless units, this range was from 2.5 to 5.0
$a\omega_{pd}$, significantly exceeding the exciter speeds in the simulation
of Sun et al. [15].

For this paper, we analyzed the microscopic structure for all of these runs;
we report snapshots of the data for three runs in the main text, and more
detailed video data for all the runs in the Supplemental Material [32].

In the experiment, the 2D microparticle cloud was imaged at 870 frames/s using
video microscopy [33, 34, 35]. Afterwards, positions of individual
microparticles were measured in each video frame, using the moment method,
which has a sub-pixel precision [33, 36]. We restricted our analysis to
particles within a $25.4\times 10.2$ mm region of interest (ROI), as described
in Supplemental Material of Ref. [9]. The microparticle positions within the
ROI, in each video frame, served as the input for the analysis that we
describe next.

###  II.2 Polylgon analysis

Figure 1: Microscopic structure of the 2D layer of microparticles. Particle
positions are shown in (a) and (b). In the region of interest (ROI) shown
here, the exciter wire moved from left to right. The left column is for an
early time when the shock had not yet disturbed the ROI, while the right
column at 248 ms later had a shock at the location marked with a black line.
Voronoi diagrams in (c) and (d), and polygons in (e) and (f), reveal regions
of order and disorder differently. In Voronoi diagrams, a topological defect
corresponds to a non-6-fold coordination, while in a polygon construction map,
a geometric defect is recognized as a polygon with more than four sides. The
polygon construction map sensitively reveals a non-hexagonal arrangement, near
the shock, where there is an abundance of quadrilateral polygons. These data
are from a run, corresponding to Fig. 2 in Ref. [9], with speeds of 101.6 ms/s
for the exciter and 117.5 mm/s for the shock. In the postshock for this run,
the crystal has clearly melted, as indicated by the considerable disorder that
is typical of a liquid. The snapshots shown here are specific frames from
videos r_A.mp4, v_A.mp4, and g_A.mp4 in the Supplemental Material [32].

Two methods of structural analysis were used for this paper: Voronoi analysis
[26] and the Glaser-Clark polygon analysis method [28, 29]. Both methods have
been used in the literature to characterize defects, for the melting of a 2D
crystalline lattice. Both of these methods start with particle position
measurements as the input for a Delaunay triangulation calculation, which
identifies bonds as line segments connecting nearby particles. In an ideal 2D
crystal under equilibrium conditions at zero temperature, these bonds will be
arranged as equilateral triangles, i.e., with six-fold symmetry. However, in
actual conditions, they will not be perfectly equilateral.

The Voronoi method has been used to study 2D melting with experimental data
since at least the 1987 colloid experiments of Ref. [37]. The Voronoi
algorithm is to form a polygon having its edges being the perpendicular
bisectors of the bonds. In an equilibrium state of a 2D crystal, these Voronoi
polygons will be hexagons, which is to say the local geometry has six-fold
coordination. Five- and seven-sided polygons, and any other polygon that is
not six-fold, are identified as a disclination, which is a defect. This kind
of defect analysis has an essentially binary outcome: the vertex at the center
of a polygon is either six-fold or it is not. Sufficient displacements of
particles from the ground-state arrangement can result in a pair of polygons
changing from hexagons to a five-seven pair of Voronoi cells. Sometimes only a
tiny displacement of a single particle can result in such a change, from non-
defect to defect. This binary classification, defect vs non-defect, is a
disadvantage of the Voronoi method.

The Glaser-Clark polygon analysis method, which was developed in 1990 [28],
has a distinctive difference compared to the Voronoi approach because it does
not have a simple binary outcome of defects vs non-defects. Instead, Glaser-
Clark polygons have various orders: triangular, quadrilateral, pentagonal, or
even higher order polygons. A triangular polygon is non-defective. A
quadrilateral is the least severe defect, followed by a pentagon and so on for
more severe defects [38, 27]. This gradation of defects is an advantage of the
Glaser-Clark method. The algorithm for producing these polygons from the
Delaunay triangulation is to remove bonds if their angles deviate sufficiently
from an equilateral shape. The criterion we will use, as in Refs [27, 28, 29],
is to remove bonds opposite an unusually large angle of 75°. We used the
source code for Glaser-Clark polygon analysis provided in the Supplementary
Material of Ref. [27]. Despite its advantages, this method has not been used
as much as Voronoi analysis in the literature. Glaser and Clark originally
demonstrated it using simulation data for a 2D lattice under steady conditions
[28, 29]; it was then used with experimental data to characterize the sudden
melting of a 2D dusty plasma crystal [27, 30], the melting process in 2D
Yukawa systems [38], and structure and dynamics in screened 2D dipole systems
[39].

##  III Results

In Figs. 1-3 we present 2D maps of the microstructure. The particle positions
shown in the upper rows of these figures were the input to a Delaunay
triangulation calculation, for preparing the Voronoi polygons in the middle
rows and the Glaser-Clark polygons in the lower rows. These three figures are
for different runs, each with a different speed for the exciter wire, which
moved in the positive $x$ direction. In Fig. 1 we also provide an example of
the conditions of the microparticle cloud before it was disturbed by the
movement of the exciter. In the Supplemental Material [32] we provide videos
of the particle positions, Voronoi cells and Glaser-Clark polygons; these
videos are for six experimental runs, including the three runs shown in the
snapshots of Figs. 1-3.

Figure 2: Microscopic structure, with data from a run with a slower speed of
88.9 mm/s for the exciter and 110.5 mm/s for the shock. In this run, there was
again an abundance of quadrilaterals near the shock. In the postshock, the
rearrangement relaxed more nearly toward a hexagonal lattice, meaning that
there is no strong indication of melting. These snapshots are from videos
r_C.mp4, v_C.mp4, and g_C.mp4 in the Supplemental Material [32]. Figure 3:
Microscopic structure, with data from a run with an even slower speed of 63.5
mm/s for the exciter and 83.8 mm/s for the shock. As in the other runs,
quadrilaterals were abundant near the shock. The disorder in the postshock was
more than in Fig. 2, and almost as much as in Fig. 1. We note that the speed
of the shock, by itself, does not determine the outcome of whether the
postshock is melted. These snapshots are from videos r_E.mp4, v_E.mp4, and
g_E.mp4 in the Supplemental Material [32].

The shock has a peak density, at a location marked as a solid line in Figs.
1-3. We refer to that point as the shock position. Immediately in front of the
shock position there is a localized high-gradient region, which we analyzed in
Ref. [14] to quantify the finite thickness of the shock.

Our chief result is the finding that microstructures rearrange from hexagonal
to quadrilateral within a shock. Such a quadrilateral structure is seen in
Figs. 1-3, and in the videos in the Supplemental Material [32].

The quadrilateral structure is most often observed at locations behind the
shock position. In some instances, it also extends slightly ahead of the shock
position, as seen in the Glaser-Clark polygon maps. This quadrilateral region
is located between what we call the preshock and postshock. The preshock
region is located at the most positive $x$-coordinate. Here, we distinguish
the quadrilateral region and the postshock not according to a numerical
criterion, but by our qualitative observation of the type of disorder seen in
the local microstructure. In our experiment, the quadrilateral region
typically had a thickness of about 4 mm.

An additional result is a considerable increase in defects in the post-shock
region. The defect fraction is quantified as the number of non-six-fold
Voronoi cells, divided by the total number of cells within a given region.
Values of this defect fraction are printed on each video frame in the Voronoi-
map movies in Supplemental Material [32], and they are summarized in Table I.
In all six runs the defect fraction behind the shock was mainly greater than
the 25% level considered to be a liquid in the melting experiment of Ref.
[40]. For comparison, the defect fraction was only 3% to 14% in the
unperturbed condition, before the shock arrived in the ROI.

Table 1: Runs from the experiment of Ref. [9]. The defect fraction is calculated as the number of non-six-fold Voronoi cells divided by the total number of Voronoi cells, in a spatial region at least 1.25 mm behind the point of highest density. Values for this defect fraction are reported as a typical range over the course of time; values for specific video frames are printed in the videos in the Supplemental Material [32]. The unperturbed conditions were obtained early in each run before the microparticles in the ROI were distrubed noticeably, while the conditions behind the shock are for video frames when the point of maximum density had reached at least the farthest 70% of the ROI. The Exciter speed and Mach number data are the same as in Fig. 5 of Ref. [9]. For comparison, we note that Feng et al. [40] found in their melting experiment that a liquid typically has a defect fraction greater than approximately 0.20. Run |  $\mathrm{Exciter}\,\mathrm{speed}$, mm/s |  $\mathrm{Exciter}\,\mathrm{speed}/a\omega_{p}$ |  Exciter Mach Number |  Defect fraction: unperturbed |  Defect fraction: behind the shock
---|---|---|---|---|---
A |  101.6 |  4.9 |  5.36 |  0.05 - 0.08 |  0.41 – 0.47
B |  101.6 |  4.9 |  5.36 |  0.02 – 0.06 |  0.47 – 0.59
C |  88.9 |  4.3 |  4.69 |  0.07 – 0.08 |  0.17 – 0.35
D |  76.2 |  3.7 |  4.02 |  0.05 – 0.08 |  0.26 – 0.37
E |  63.5 |  3.1 |  3.35 |  0.03 – 0.10 |  0.33 – 0.42
F |  50.8 |  2.5 |  2.68 |  0.10 – 0.14 |  0.41 – 0.48

Glaser-Clark polygon maps are more useful than Voronoi analysis, for the
purpose of identifying a quadrilateral structure. A limitation of Voronoi
polygons is that they are most likely to have 5, 6 or 7 sides, while very
seldomly do they have 4 sides. The Voronoi mapping struggles to map onto a
quadrilateral arrangement particles, doing so only by making some of the
polygon sides extremely short, as can be seen in the magnified view in Fig. 4.
The Voronoi cells in Fig. 4 (a) might at a glance appear to be quadrilateral,
but a close inspection shows that they are mostly 6-sided, along with some
that are 5 or 7-sided, but always with some sides that are extremely short.

Figure 4: Magnification of Fig. 1, showing details in the microscopic
structure at a location near $x=16.4$ mm and $y=7.7$ mm. Under the compressed
condition near the shock, six-fold Voronoi cells (a) have sides that are very
unequal in length, with some sides that are so short that it is not very
meaningful to say that the coordination is six-fold. Under these
nonequilibrium conditions, we find that the polygons in (b) better indicate
the nearly rectangular arrangement of particles.

##  IV Discussion

An analysis method that had not previously been applied to 2D shocks allowed
us to find that within a shock, a hexagonal microstructure rearranged into a
quadrilateral microstructure that had not be reported before. This analysis
method was the Glaser-Clark polygon analysis, which we compared to the more
common Voronoi analysis method. The data that we analyzed were from runs of
the 2D dusty plasma shock experiment of Ref. [9]. An advantage of the Glaser-
Clark polygons that we found was that it easily reveals quadrilateral
structures, while Voronoi polygons are not well suited for this purpose, as
they are mostly 5, 6 or 7-sided. Another advantage of Glaser-Clark polygons,
which was reported previously by Ruhunusiri and Goree [27], is the gradation
of defects, classified by the number of sides of the polygon increases as the
disorder in the microstructure becomes greater.

The quadrilateral structure that we observed is generally not otherwise seen,
in a single layer of identical particles interacting isotropically. A
triangular lattice with six-fold symmetry is the only stable lattice in this
case. As an example, a square lattice configuration is not stable in
equilibrium, as the total interparticle potential energy can be diminished by
alternate rows slipping so that the lattice becomes triangular.

We believe that it may be significant that the quadrilateral structure was
observed at locations where the conditions were strongly nonequilibrium.
Shocks are certainly nonequilibrium, especially near the location of greatest
compression, and it was near this location where we observed quadrilateral
structures.

There are previous dusty plasma experiments, done under nearly 2D conditions,
where quadrilateral microstructures have been observed under steady
conditions. These experiments, which did not involve shocks, fall in two
categories: those with an externally imposed pattern, and those with out-of-
plane particles. An externally imposed pattern was applied to a magnetized
plasma, when a metal mesh with an 0.56 mm square pattern was inserted as a
boundary in the Magnetized Dusty Plasma Experiment (MDPX) device [41, 42, 43],
causing the same square pattern to appear in the arrangement of dust
particles. This effect can be attributed to the external influence of the
mesh, which persisted throughout the volume due to the magnetization of
electrons and ions. For unmagnetized plasmas, we note two experiments that
involved out-of-plane positions of some particles. In both of these
experiments authors attributed their observations of quadrilateral
microstructures to non-isotropic interparticle interactions: ion wakes beneath
the upper particles that disturbed the lower particles. In 2019 Huang et al.
[23] reported an experiment with a binary mixture, and they observed
distinctive regions with four-fold symmetry rather than the six-fold symmetry
that is common in a 2D crystal with monodisperse microparticles. In 2022 Singh
et al. [25] used monodisperse particles and attained a quadrilateral
microstructure by causing the microparticle layer to buckle, meaning that the
horizontal repulsion between microparticles overcame the confining forces in
the vertical direction.

Shocks in a 2D layer have been studied previously in dusty plasma experiments
[5, 44] and in molecular dynamics simulations [45, 7, 10, 15, 16, 12, 13]. The
latter mimicked dusty plasmas by tracking the motion of identical particles
interacting with a repulsive Yukawa potential, but without capturing all the
effects in the experiment such as out-of-plane displacements and finite
dispersion of particle size. In none of the experiment or simulation papers
that we reviewed did the authors comment upon a quadrilateral microstructure.
Nevertheless, we see indications of quadrilateral shapes in both the
experiment in Fig. 1 of Ref. [44] and in Fig. 1(a) of the Supplemental
Material of Ref. [45]. These quadrilateral structures appear in the vicinity
of the peak density of a blast-wave shock in the experiment of Ref. [44] but
on the pre-shock side in the simulation of a steadily-driven shock in Ref.
[45].

Besides our primary result of finding a quadrilateral microstructure, we also
discussed the production of defects and melting. The defect fraction was
considerably higher behind the shock as compared to the unperturbed condition
in six runs that were analyzed. In fact, behind the shock the defect fraction
was generally so high as to be consistent with melting. We did not observe a
distinctive trend for the defect fraction to increase with exciter speed. We
note that all of the exciter speeds in the experiment were at Mach numbers
considerably higher than in recent simulations [13] and [12] where the onset
of shock-induced melting was studied. Some factors that could affect the
defect fraction behind the shock include the presence, in the earlier
unperturbed structure, of domain walls and other defects, which vary from one
experimental run to another. One more factor that can vary from one run to
another is the direction of the lattice alignment in the unperturbed
condition. These conditions for shocks in a solid differ from those of
gasdynamic shocks, where the unshocked region has no microscopic structure or
favored directions.

###### Acknowledgements.

This work was supported by the Army Research Office under MURI Grant No.
W911NF-18-1-0240, the United States Department of Energy under Grant No. DE-
SC0014566 and NASA/JPL RSA Nos. 1689926, and the National Science Foundation
under Grant No. PHY-1740379.

## References

  * Goree [2008] J. Goree, _Fundamentals of Dusty Plasmas_ (Wiley-VCH, 2008), chap. 6, pp. 157–206.
  * Melzer [2019] A. Melzer, _Physics of Dusty Plasmas_ , vol. 962 of _Lecture Notes in Physics_ (Springer International Publishing, Cham, 2019).
  * Samsonov et al. [1999] D. Samsonov, J. Goree, Z. W. Ma, A. Bhattacharjee, H. M. Thomas, and G. E. Morfill, Phys. Rev. Lett. 83, 3649 (1999).
  * Samsonov et al. [2000] D. Samsonov, J. Goree, H. M. Thomas, and G. E. Morfill, Phys. Rev. E 61, 5557 (2000).
  * Samsonov et al. [2004] D. Samsonov, S. K. Zhdanov, R. A. Quinn, S. I. Popel, and G. E. Morfill, Phys. Rev. Lett. 92, 255004 (2004).
  * Couëdel et al. [2012] L. Couëdel, D. Samsonov, C. Durniak, S. Zhdanov, H. M. Thomas, G. E. Morfill, and C. Arnas, Phys. Rev. Lett. 109, 175001 (2012).
  * Lin et al. [2019] W. Lin, M. S. Murillo, and Y. Feng, Phys. Rev. E 100, 043203 (2019).
  * Kananovich and Goree [2020a] A. Kananovich and J. Goree, Phys. Plasmas 27, 113704 (2020a).
  * Kananovich and Goree [2020b] A. Kananovich and J. Goree, Phys. Rev. E 101, 043211 (2020b).
  * Lin et al. [2020] W. Lin, M. S. Murillo, and Y. Feng, Phys. Rev. E 101, 013203 (2020).
  * Tao and Duan [2020] L.-L. Tao and W.-S. Duan, Chin. J. Phys. 68, 950 (2020).
  * Sun and Feng [2021] T. Sun and Y. Feng, Phys. Plasmas 28, 063702 (2021).
  * Qiu et al. [2021] P. Qiu, T. Sun, and Y. Feng, Phys. Plasmas 28, 113702 (2021).
  * Kananovich and Goree [2021] A. Kananovich and J. Goree, Phys. Rev. E 104, 055201 (2021).
  * Sun et al. [2021] T. Sun, M. S. Murillo, and Y. Feng, Phys. Plasmas 28, 103703 (2021).
  * Ding et al. [2021] X. Ding, S. Lu, T. Sun, M. S. Murillo, and Y. Feng, Phys. Rev. E 103, 013202 (2021).
  * Gou et al. [2021] X.-Q. Gou, K.-H. An, and W.-S. Duan, Braz. J. Phys. 51, 1346 (2021).
  * Li and Duan [2021] Z.-Z. Li and W.-S. Duan, Phys. Plasmas 28 (2021).
  * Tao et al. [2021] L.-L. Tao, L. Wei, B. Liu, H. Zhang, and W.-S. Duan, Pramana-J. Phys. 95 (2021).
  * Zhang et al. [2022] X. Zhang, J. Xia, and X. Liu, Phys. Rev. B 105 (2022).
  * Samsonov and Morfill [2008] D. Samsonov and G. Morfill, IEEE Trans. Plasma Sci. 36, 1020 (2008).
  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky, _Principles of Condensed Matter Physics_ (Cambridge University Press, New York, 1995), pp. 50–53.
  * Huang et al. [2019] H. Huang, A. V. Ivlev, V. Nosenko, Y.-F. Lin, and C.-R. Du, Phys. Plasmas 26, 013702 (2019).
  * Zampetaki et al. [2020] A. V. Zampetaki, H. Huang, C.-R. Du, H. Löwen, and A. V. Ivlev, Phys. Rev. E 102, 043204 (2020).
  * Singh et al. [2022] S. Singh, P. Bandyopadhyay, K. Kumar, and A. Sen, Phys. Rev. Lett. 129, 115003 (2022).
  * Voronoi [1908] G. Voronoi, J. Reine Angew. Math. 1908, 198 (1908).
  * Ruhunusiri et al. [2011] W. D. S. Ruhunusiri, J. Goree, Y. Feng, and B. Liu, Phys. Rev. E 83, 66402 (2011).
  * Glaser and Clark [1990] M. A. Glaser and N. A. Clark, Phys. Rev. A 41, 4585 (1990).
  * Glaser and Clark [1992] M. A. Glaser and N. A. Clark, _Melting and Liquid Structure in two Dimensions_ (Wiley, 1992), pp. 543–709.
  * Ruhunusiri [2014] W. D. S. Ruhunusiri, Ph.D. thesis, the University of Iowa (2014).
  * Chaubey et al. [2021] N. Chaubey, J. Goree, S. J. Lanham, and M. J. Kushner, Phys. Plasmas 28, 103702 (2021).
  * Note [1] Note1, the Supplemental Material for this preprint will be published, along with the article itself, in a future journal publication.
  * Feng et al. [2007] Y. Feng, J. Goree, and B. Liu, Rev. Sci. Instrum. 78, 053704 (2007).
  * Feng et al. [2011] Y. Feng, J. Goree, and B. Liu, Rev. Sci. Instrum. 82, 053707 (2011).
  * Feng et al. [2016] Y. Feng, J. Goree, Z. Haralson, C.-S. Wong, A. Kananovich, and W. Li, J. Plasma Phys. 82, 615820303 (2016).
  * Schneider et al. [2012] C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, Nat. Methods 9, 671 (2012).
  * Murray and Van Winkle [1987] C. A. Murray and D. H. Van Winkle, Phys. Rev. Lett. 58, 1200 (1987).
  * Radzvilavčius [2012] A. Radzvilavčius, Phys. Rev. E 86, 051111 (2012).
  * Aldakul et al. [2020] Y. K. Aldakul, Z. A. Moldabekov, and T. S. Ramazanov, Phys. Rev. E 102, 033205 (2020).
  * Feng et al. [2008] Y. Feng, J. Goree, and B. Liu, Phys. Rev. Lett. 100, 205007 (2008).
  * Thomas et al. [2015a] J. Thomas, Edward, U. Konopka, B. Lynch, S. Adams, S. LeBlanc, R. L. Merlino, and M. Rosenberg, Phys. Plasmas 22 (2015a), 113708.
  * Thomas et al. [2015b] J. Thomas, Edward, B. Lynch, U. Konopka, R. L. Merlino, and M. Rosenberg, Phys. Plasmas 22 (2015b), 030701.
  * Hall et al. [2018] T. Hall, J. Thomas, Edward, K. Avinash, R. Merlino, and M. Rosenberg, Phys. Plasmas 25 (2018), 103702.
  * Oxtoby et al. [2013] N. P. Oxtoby, E. J. Griffith, C. Durniak, J. F. Ralph, and D. Samsonov, Phys. Rev. Lett. 111, 015002 (2013).
  * Marciante and Murillo [2017] M. Marciante and M. S. Murillo, Phys. Rev. Lett. 118, 025001 (2017).


SI_final.pdf

# Microscopic theory of colour in lutetium hydride

Sun-Woo Kim,${}^{1,\;{\color[rgb]{1,0,0}*}}$ Lewis J. Conway,1,2 Chris J.
Pickard,1,2

G. Lucian Pascut,3 and Bartomeu
Monserrat${}^{1,4,\;{\color[rgb]{1,0,0}\dagger}}$

_1 Department of Materials Science and Metallurgy, University of Cambridge,
27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom_

_2 Advanced Institute for Materials Research, Tohoku University,
2-1-1 Katahira, Aoba, Sendai 980-8577, Japan_

_3 MANSiD Research Center and Faculty of Forestry,
Stefan Cel Mare University (USV), Suceava 720229, Romania_

_4 Cavendish Laboratory, University of Cambridge,
J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom_

∗ email<EMAIL_ADDRESS>
† email<EMAIL_ADDRESS>
##  I Abstract

Nitrogen-doped lutetium hydride has recently been proposed as a near-ambient-
conditions superconductor. Interestingly, the sample transforms from blue to
pink to red as a function of pressure, but only the pink phase is claimed to
be superconducting. Subsequent experimental studies have failed to reproduce
the superconductivity, but have observed pressure-driven colour changes
including blue, pink, red, violet, and orange. However, discrepancies exist
among these experiments regarding the sequence and pressure at which these
colour changes occur. Given the claimed relationship between colour and
superconductivity, understanding colour changes in nitrogen-doped lutetium
hydride may hold the key to clarifying the possible superconductivity in this
compound. Here, we present a full microscopic theory of colour in lutetium
hydride, revealing that hydrogen-deficient LuH2 is the only phase which
exhibits colour changes under pressure consistent with experimental reports,
with a sequence blue-violet-pink-red-orange. The concentration of hydrogen
vacancies controls the precise sequence and pressure of colour changes,
rationalising seemingly contradictory experiments. Nitrogen doping also
modifies the colour of LuH2 but it plays a secondary role compared to hydrogen
vacancies. Therefore, we propose hydrogen-deficient LuH2 as the key phase for
exploring the superconductivity claim in the lutetium-hydrogen system.
Finally, we find no phonon-mediated superconductivity near room temperature in
the pink phase.

##  II Introduction

The proposal by Ashcroft that hydrogen-rich compounds could host high
temperature phonon-mediated superconductivity under pressure
hydrogen_superconductivity ; ashcroft_superconducting_hydrides  has stimulated
a profusion of theoretical proposals for high pressure superconducting
hydrides H3S-duan ; LaH10-theory-hemley1 ; LaH10-theory-ma ; LaH10-theory-
hemley2 ; pickard-hydride-review ; pickard-sc-high-throughput ; pickett-
hydride-review  and the subsequent experimental discovery of some of these
H3S-eremets ; LaH10-hemley ; LaH10-eremets . This new class of hydride
superconductors has re-ignited the search for superconductivity at ambient
conditions, and Dasenbrock-Gammon and co-workers have recently reported
superconductivity in nitrogen-doped lutetium hydride with a maximum critical
temperature of $294$ K at a moderate pressure of $10$ kbar Lu-H-
N_superconductivity_nature . Interestingly, superconductivity is reported to
coincide with drastic colour changes in the reflectivity of the sample:
increasing pressure transforms a non-superconducting blue metal to a
superconducting pink metal at $3$ kbar, and a further transformation to a non-
superconducting red metal above $30$ kbar.

This remarkable report has sparked a growing number of experimental exp-
colour-change-CPL ; exp-reflectivity-scibul ; exp-colour-res-mag-wen1 ; exp-
colour-res-mag-wen2 ; exp-colour-pressure-medium-arxiv ; exp-lu-h-n-no-sc-
arxiv ; exp-elec-mag-arxiv ; exp-grain-size-arxiv ; exp-luh3-phase-transition-
arxiv ; luh-sc-high-pressure ; salke2023evidence  and theoretical xie-
structure-prediction-arxiv ; structure-prediction-binary-arxiv ; duan-
structure-prediction-arxiv ; zurek-structure-prediction-arxiv ; pickard-
structure-prediction-arxiv ; ndoping-pressure-luh3-arxiv  investigations, most
of which have so far failed at reproducing or explaining near-ambient
superconductivity. On the experimental front, most measurements of resistivity
and magnetic susceptibility find no superconductivity near ambient conditions,
with the exception of a recent work in which resistivity changes consistent
with high temperature superconductivity are reported salke2023evidence .
Puzzlingly, multiple studies report pressure-driven colour changes, but these
include a wide range of seemingly incompatible colour sequences and pressure
conditions: blue-to-pink at $3$ kbar and pink-to-red at $30$ kbar in the
original report Lu-H-N_superconductivity_nature ; blue-to-pink at $22$ kbar
and pink-to-red at $40$ kbar exp-colour-change-CPL ; blue-to-violet upon
contact with a diamond culet, violet-to-red at $30$ kbar and red-to-orange at
$120$ kbar exp-reflectivity-scibul ; blue-to-violet at $94$ kbar exp-
reflectivity-scibul ; blue-to-violet at $120$ kbar, violet-to-pink-to-red
gradually between $160$ kbar and $350$ kbar and red persisting up to $420$
kbar exp-colour-res-mag-wen1 ; exp-colour-res-mag-wen2 ; blue-to-violet-to-
pink-to-red with transition pressures differing by up to $60$ kbar depending
on the pressure medium used exp-colour-pressure-medium-arxiv ; and persistent
blue colour up to $65$ kbar exp-lu-h-n-no-sc-arxiv . Growing evidence suggests
that the colour changes are significantly affected by the initial compression
procedure exp-reflectivity-scibul  and by the pressure medium used in the
diamond anvil cell exp-colour-pressure-medium-arxiv .

On the theoretical front there have been multiple reports of structure
searches in the Lu-H binary and the Lu-H-N ternary systems xie-structure-
prediction-arxiv ; structure-prediction-binary-arxiv ; duan-structure-
prediction-arxiv ; zurek-structure-prediction-arxiv ; pickard-structure-
prediction-arxiv . Most studies only report metastable ternary structures, but
Ferreira and co-workers report a ternary Lu4N2H5 stable structure pickard-
structure-prediction-arxiv . The roles of pressure and nitrogen doping
ndoping-pressure-luh3-arxiv  and of quantum and thermal ionic vibrations
anharmonic-luh3-arxiv  in stabilising the cubic LuH3 structure have also been
studied. None of the predicted stable or metastable structures are found to be
phonon-mediated superconductors near room temperature.

Given the claimed association between superconductivity and colour changes in
the original superconductivity report, understanding colour changes in
nitrogen-doped lutetium hydride holds the key to clarifying the possible
superconductivity in this compound. However, experimental reports provide an
inconsistent picture regarding pressure-driven colour changes, and there are
no theoretical studies yet. In this work, we provide a full microscopic theory
of colour in lutetium hydride.

##  III Results

###  III.1 LuH2 under ambient conditions

Figure 1:  Structure, reflectivity, and colour of LuH2. a. Crystal structure
of $Fm\overline{3}m$ LuH2. b. Reflectivity of LuH2 calculated using semilocal
density functional theory (DFT), DFT corrected with a Hubbard $U$ term
(DFT$+U$), and DFT including electron-phonon coupling (DFT$+$EPC). The
experimental reflectivity is taken from Ref. exp-reflectivity-scibul . c.
Colour and photorealistic rendering of LuH2 calculated using DFT, DFT$+U$, and
DFT$+$EPC. The photorealistic rendering is shown as LuH2 surrounding a grey
ball with an opening in the centre.

Lutetium hydride under ambient conditions crystallises in the LuH2
stoichiometry with cubic space group $Fm\overline{3}m$. As shown in Fig. 1,
LuH2 adopts the fluorite structure with the lutetium atoms occupying the sites
of an fcc lattice, and the hydrogen atoms occupying the tetrahedral
interstitial sites.

LuH2 is a metal whose reflectivity endows it with a blue appearance. We
demonstrate the validity of our computational approach by reporting the
calculated colour of LuH2 at ambient conditions in Fig. 1. We show results
using three distinct computational models to explore the potential role of
electron correlation due to the presence of lutetium $5d$ electrons and the
potential role of strong electron-phonon coupling due to the presence of
hydrogen.

The first model we consider uses semilocal density functional theory (DFT) in
the generalised gradient approximation, labelled DFT in Fig. 1. This model
provides a basic description of the electronic structure without a detailed
treatment of electron correlation and without the inclusion of electron-phonon
effects. The calculated reflectivity is large in the infrared region above
$800$ nm, is strongly suppressed in the red part of the visible spectrum with
a calculated minimum at $710$ nm, and increases gradually towards the blue
part of the visible spectrum. The overall shape of the reflectivity is
consistent with that observed experimentally exp-reflectivity-scibul  and
directly leads to the blue colour of LuH2.

Lutetium has an electronic configuration with a partially filled $5d$ shell,
which suggests that electronic correlation beyond that captured by standard
DFT may contribute to the electronic properties of LuH2. To explore the
possible role of electron correlation, we repeat our calculations using DFT
corrected with a Hubbard $U$ term, labelled as DFT$+U$ in Fig. 1, which
captures static correlation. The reflectivity curve has a similar shape to
that obtained at the DFT level, but the minimum of the reflectivity has a
lower value and occurs at a slightly shorter wavelength of $690$ nm. Combined
with a slighly larger reflectivity in the blue part of the visible spectrum,
we obtain a slightly brighter blue colour for LuH2 using the DFT$+U$ model.
These results indicate that static electron correlation arising from lutetium
only plays a minor role in LuH2. We have also performed dynamical mean field
theory calculations that capture dynamical correlation and also find that they
can be neglected. We rationalise these results by noting that $5d$ orbitals
have a large electron bandwidth spanning multiple eV (see Supplementary Fig.
S1) which implies that the spatial extent of the orbitals is large and the
corresponding local correlations weak.

Hydrogen is the lightest of all elements, and as such it exhibits significant
nuclear motion, even at zero temperature, due to quantum zero-point effects.
This significant nuclear motion can lead to strong electron-phonon coupling,
and this is indeed the prime motivation behind the proposal that hydrogen-rich
compounds could be high temperature phonon-mediated superconductors. To
explore the possible role of electron-phonon interactions in LuH2, we repeat
our calculations including contributions from both zero-point quantum nuclear
motion at $0$ K and thermal nuclear motion at finite temperature, labelled
DFT$+$EPC in Fig. 1. The reflectivity curve has a similar shape to those
obtained with DFT and DFT$+U$, but exhibits a lower value in the infrared
region and a larger value in the red region. We again obtain a blue colour,
indicating that electron-phonon coupling does not significantly modify the
reflectivity of LuH2.

Overall, we find that electronic correlation and electron-phonon interactions
make a small contribution, and that the main features of the reflectivity
curve and the resulting blue colour of LuH2 are correctly captured by
semilocal DFT. Therefore, our subsequent discussion neglects electron
correlation and electron-phonon interactions, but further details about these
contributions are included in the Supplementary Information.

###  III.2 Lutetium hydride colour changes under pressure

Figure 2:  Pressure dependence of the convex hull diagram, colour, and
reflectivity of the lutetium-hydrogen binary system. a. Convex hull diagrams
as a function of pressure for compositions in the range LuH0 to LuH3. Dark
blue circles indicate thermodynamically stable structures and red circles
indicate metastable structures. The dashed vertical lines for selected
stoichiometries are shown for guidance only. b. Colour of LuHn compounds as a
function of pressure. c. Reflectivity of LuHn compounds as a function of
pressure.

To build a microscopic theory of colour in lutetium hydride, we first explore
the pressure-driven colour changes in the lutetium-hydrogen system. We have
performed extensive structure searches for stoichiometries ranging from LuH0
to LuH3 at multiple pressures. The results are summarised in the convex hull
diagrams depicted in Fig. 2a. At $0$ kbar, the only thermodynamically stable
structures are LuH2 ($Fm\overline{3}m$ space group) and LuH3
($P\overline{3}c1$ space group; not cubic). There are multiple metastable
structures in the entire composition space from LuH0 to LuH3 that are within
$60$ meV/atom of the convex hull. Increasing pressure leads to multiple
additional stable structures with stoichiometries intermediate between LuH2
and LuH3. We highlight that substoichiometric LuH2-δ structures are close to
the convex hull at both $0$ and $400$ kbar (within $11$ and $42$ meV/atom,
respectively, up to $\delta=0.25$). We have also checked their dynamical
stability (see Supplementary Fig. S14), so we expect that they can be
accessible experimentally.

We show the calculated colour as a function of pressure for the most stable
structure at each composition between LuH0 and LuH3 in Fig. 2b. Consistently
with experimental observations, we find that pure lutetium has a silvery white
colour exp-elec-mag-arxiv  (see also Supplementary Fig. S17) and, as described
above in Fig. 1, LuH2 has a blue colour. Across the entire composition space,
the only compositions that exhibit a blue colour at ambient pressure occur for
stoichiometries close to LuH2. Similarly, the only stoichiometries that
exhibit a violet colour at high pressure are those close to LuH2.
Specifically, we find this trend is present for substoichiometric LuH2-δ, but
not present for suprastoichimoetric LuH2+δ. We also note that LuH3 has a grey-
green colour at ambient pressure that becomes grey with increasing pressure.

Figure 2c shows reflectivity curves for selected stoichiometries in the range
LuH0 to LuH3 at multiple pressures. We note that only the reflectivity of LuH2
has a minimum in the red part of the visible spectrum leading to an overall
blue colour, as already discussed in Fig. 1 above. The reflectivities of all
other compositions show relatively flat curves across the visible spectrum,
which lead to colours with various tones of grey.

Overall, the results depicted in Fig. 2 show that the only compounds in the
lutetium-hydrogen binary space that are blue at ambient conditions and violet
at high pressure have stoichiometries close to LuH2 with a moderate amount of
hydrogen vacancies. This conclusion still holds when structures in the ternary
lutetium-hydrogen-nitrogen system are considered. To demonstrate this, we
perform extensive crystal structure searches in the full Lu-H-N ternary space
(see Supplementary Fig. S12), and note that our structure searches are the
only ones of all those published that identify a stable ternary compound at
ambient pressure pickard-structure-prediction-arxiv . We have calculated the
colour of this stable compound and also the colour of multiple other ternary
compounds that are not on the convex hull but whose simulated X-ray
diffraction data is consistent with experimental reports. We find that none of
these structures exhibit colours that are consistent with experiment (see
Supplementary Fig. S13).

These observations allow us to conclude that the colour changes in the
lutetium-hydrogen-nitrogen system are dominated by the LuH2 stoichiometry. In
particular, we discard the LuH3 composition proposed by Dasenbrock-Gammon and
co-workers to explain high temperature superconductivity Lu-H-
N_superconductivity_nature  as this structure has a grey-green colour at all
pressures. LuH2 has also been identified as the relevant stoichiometry by
comparing the calculated equation of state zurek-structure-prediction-arxiv
and X-ray diffraction patterns zurek-structure-prediction-arxiv ; structure-
prediction-binary-arxiv ; pickard-structure-prediction-arxiv ; xie-structure-
prediction-arxiv  to experiment, and we also note a recent experimental work
that uses LuH2-based samples and that has successfully reproduced the Raman
spectrum and colour sequence with pressure reported in the original work exp-
colour-res-mag-wen1 .

###  III.3 Hydrogen deficient LuH2-δ

The calculated colour changes from blue to violet in LuH2 are also observed in
multiple experiments exp-reflectivity-scibul ; exp-colour-res-mag-wen1 ; exp-
colour-res-mag-wen2 ; exp-colour-pressure-medium-arxiv , but they occur at
different pressures in different experiments, ranging from $0$ to at least
$190$ kbar, but possibly higher as some experiments only observe a blue phase.
Furthermore, some experimental observations reveal additional colour changes
with increasing pressure, which include pink Lu-H-N_superconductivity_nature ;
exp-colour-change-CPL ; exp-colour-pressure-medium-arxiv , red Lu-H-
N_superconductivity_nature ; exp-colour-change-CPL ; exp-reflectivity-scibul ;
exp-colour-pressure-medium-arxiv , and orange exp-reflectivity-scibul . The
pink colour is particularly important as it is associated with the
superconducting phase in the original report Lu-H-N_superconductivity_nature .
Based on these observations, we next explore the pressure-driven colour
changes of substoichiometric LuH2 in more detail.

Figure 3:  Pressure dependence of the reflectivity and colour of pure and
hydrogen-deficient lutetium dihydrides. a-c. Reflectivity as a function of
pressure (in kbar) for a LuH2, b LuH1.875 and c LuH1.750. d-f. Colour and
photorealistic rendering of d LuH2, e LuH1.875 and f LuH1.750 as a function of
pressure.

We show the pressure evolution of the reflectivity and colour of LuH2 and
hydrogen-deficient LuH1.875 and LuH1.750 in Fig. 3. At ambient conditions LuH2
has the reflectivity described in Fig. 1 and repeated in Fig. 3 with a minimum
in the red part of the spectrum that leads to an overall blue colour. With
increasing pressure, the reflectivity minimum shifts towards shorter
wavelengths, and the reflectivity from the red part of the spectrum increases,
in agreement with experiment exp-reflectivity-scibul . This leads to a gradual
colour change from blue to violet with increasing pressure. LuH2 undergoes a
structural phase transition at $732$ kbar to a phase of space group $P4/nmm$,
which has a grey colour with an orange-red hue (reflectivity and colour shown
in Supplementary Fig. S7).

At ambient conditions, LuH1.875 exhibits a reflectivity with a shape similar
to that of LuH2 but with the minimum occurring at somewhat shorter wavelenghts
of $600$ nm. The resulting colour is still blue. Similar to LuH2, increasing
pressure leads to an overall shift of the reflectivity minimum to shorter
wavelengths and to an increase in the reflectivity in the red part of the
spectrum. As a result we observe a blue-to-violet colour change at a pressure
of about $100$ kbar, significantly lower than the corresponding colour change
in pure LuH2 and in the experimentally observed pressure range. Increasing
pressure further leads to a gradual transition to pink (peaking at about $300$
kbar), followed by red (peaking at about $500$ kbar) and tending towards
orange approaching $1000$ kbar. Therefore, hydrogen-deficient LuH2 exhibits a
sequence of colour changes that includes all colours reported experimentally.

The sequence and pressure of colour changes in hydrogen-deficient LuH2 is
strongly dependent on the concentration of hydrogen vacancies. Figure 3 also
depicts the pressure evolution of the reflectivity and colour changes of
LuH1.750 with a higher concentration of hydrogen vacancies. In this case, the
reflectivity minimum occurs at a wavelength of $550$ nm at ambient conditions,
giving a pink colour. Increasing pressure suppresses the reflectivity in the
blue region, turning the colour from pink towards orange at lower pressures
than those necessary for LuH1.875.

These results suggest that the seemingly contradictory experimental
observations of colour changes in lutetium hydride are likely due to varying
hydrogen vacancy concentrations in LuH2. In particular, Dasenbrock-Gammon and
co-workers observe a pink phase starting with a pressure of $3$ kbar Lu-H-
N_superconductivity_nature , significantly lower than the pressure reported in
multiple subsequent experiments. Our results suggest that this is due to a
higher concentration of hydrogen vacancies in the original work compared to
subsequent studies.

The concentration of nitrogen dopants also induces colour changes in LuH2 (see
Supplementary Fig. S10), but the colour changes driven by nitrogen doping only
play a secondary role compared to hydrogen vacancies. We have also tested the
role that hydrogen vacancies and nitrogen doping have on the reflectivity and
colour of cubic $Fm\overline{3}m$ LuH3 (see Supplementary Fig. S11), as this
phase has been tentatively identified as the parent phase responsible for the
superconductivity claim Lu-H-N_superconductivity_nature . Our results show
that the LuH3 phase retains a grey colour under nitrogen doping, further
discarding it as a relevant phase.

###  III.4 Absence of phonon-mediated superconductivity in the lutetium-
hydrogen-nitrogen system

Most experimental reports since the original announcement of near-ambient
conditions superconductivity in the lutetium-hydrogen-nitrogen system have
been unable to confirm this claim. Similarly, no calculation of stable and
metastable phases in the lutetium-hydrogen-nitrogen system has predicted a
high superconducting critical temperature within a phonon-mediated framework.

Table 1:  Calculated superconducting properties of hydrogen-deficient LuH2 at various pressures. $\lambda$ is the total electron-phonon coupling parameter computed from the Eliashberg function, $\omega_{\text{log}}$ is the logarithmic averaged frequency, $T_{c}^{\mathrm{AD}}$ is the superconducting critical temperature estimated from the semiempirical Allen-Dynes formula TC_AD_formula , and $T_{c}^{\mathrm{E}}$ is the superconducting critical temperature estimated from the isotropic Eliashberg equation. A standard value of $\mu^{*}=0.125$ is used for the Morel-Anderson Coulomb pseudopotential.  Structure | Pressure | $\lambda$ | $\omega_{\text{log}}$ | $T_{c}^{\mathrm{AD}}$ | $T_{c}^{\mathrm{E}}$
---|---|---|---|---|---
LuH1.875 | 20 kbar | 0.333 | 257.516 | 0.107 | 0.190
100 kbar | 0.326 | 257.997 | 0.085 | 0.161
400 kbar | 0.360 | 243.566 | 0.220 | 0.328
LuH1.750 | 20 kbar | 0.277 | 231.482 | 0.008 | 0.035
100 kbar | 0.276 | 224.351 | 0.007 | 0.032
400 kbar | 0.299 | 211.951 | 0.023 | 0.062

Our results suggest that hydrogen-deficient LuH2 is responsible for the colour
changes observed experimentally. Given the claim by Dasenbrock-Gammon and co-
workers that it is the pink phase that is a room temperature superconductor,
we calculate the superconducting critical temperature of LuH1.875 under
pressure, which exhibits the pink phase. We find no room-temperature
superconductivity, with a calculated critical temperature of the order of
$0.1$ K (Table I; see also details in Supplementary Note 7).

##  IV Discussion

Dasenbrock-Gammon and co-workers report superconductivity in nitrogen-doped
lutetium hydride over the pressure range $3$-$30$ kbar Lu-H-
N_superconductivity_nature . Importantly, the claimed pressure-driven
transition to and from the superconducting phase occurs simultaneously with
drastic colour changes in the sample, which is pink in the claimed
superconducting phase, compared to blue below $3$ kbar and red above $30$
kbar. Additionally, they attribute the claimed superconductivity to a cubic
LuH3 phase with some unknown concentration of nitrogen dopants and some
unknown concentration of hydrogen vacancies.

Our calculations show that the cubic LuH3 phase is not consistent with the
colour changes observed experimentally. Additionally, our results suggest that
the only phase that is consistent with the observed colour changes is
hydrogen-deficient cubic LuH2, and that the concentration of hydrogen
vacancies and nitrogen dopants controls the colour at each pressure. Finally,
we also show that hydrogen-deficient LuH2 is unlikely to be a high temperature
phonon-mediated superconductor.

Our work presents a compelling demonstration of how the first principles
prediction of the colour of a material can be exploited to effectively
identify the microscopic characteristics of the corresponding experimental
samples. Given that colour is readily accessible experimentally, while other
structural characterisation techniques can be challenging to implement
(particularly under pressure), our work provides a promising new avenue for
identifying composition and structure of complex samples. It would be
interesting to further explore the applicability of this method to study the
colour of other compounds, including strongly correlated materials Iron_SC_1 ;
Iron_SC_2  and magnetic materials AFM_1 ; AFM_2 .

##  V Methods

Electronic structure calculations. - We perform density functional theory
(DFT) calculations using the Vienna $ab$ $initio$ simulation package (vasp)
VASP1 ; VASP2  implementing the projector-augmented wave (PAW) method PAW . We
treat the $4f$ states of lutetium as valence by employing PAW pseudopotentials
with 25 valence electrons ($4f^{14}5s^{2}5p^{6}5d^{1}6s^{2}$). For the
exchange-correlation energy, we use the generalized-gradient approximation
functional of Perdew-Burke-Ernzerhof modified for solids (PBEsol) PBEsol .
Converged results are obtained with a kinetic energy cutoff for the plane wave
basis of $400$ eV and a $\mathbf{k}$-point grid of size $40\times 40\times 40$
for the LuH2 primitive cell and commensurate grids for other cell sizes and
shapes (see convergence tests in Supplementary Fig. S20). The geometry of the
structures is optimised until all forces are below $0.01$ eV/Å and the
pressure is below $1$ kbar. We also perform select calculations using DFT
corrected with a Hubbard $U$ term, and for these we use a value of $U=3$ eV.
We also perform select calculations using dynamical mean field theory with the
edmft code edmft-1 ; edmft-2 , which implements density functional theory with
embedded dynamical mean field theory (DFT+eDFMT). For the DFT part we have
used the wien2k code wien2k .

Reflectivity and colour. - Our reflectivity calculations follow the
methodology described in Ref. colour-marzari . We calculate the complex
dielectric function within the independent-particle approximation as
implemented in vasp. In the optical limit ($\mathbf{q}\rightarrow 0$), the
dielectric function $\varepsilon(\mathbf{q},\omega)$ is given by the sum of an
intraband Drude-like term $\varepsilon^{\mathrm{intra}}(\mathbf{q},\omega)$
due to the electrons at the Fermi surface and an interband term
$\varepsilon^{\mathrm{inter}}(\mathbf{q},\omega)$ describing vertical
transitions between valence and conduction bands. The explicit form of each
term is given by dielectric_fct1 ; dielectric_fct2 :

$\varepsilon^{\mathrm{intra}}(\mathbf{q},\omega)=-\frac{\omega^{2}_{D}(\mathbf{\hat{q})}}{\omega(\omega+i\gamma)},$
(1)

where the independent particle approximation Drude plasma frequency is

$\omega^{2}_{D}(\mathbf{\hat{q})}=\frac{4\pi}{V}\sum_{\mathbf{k}}\sum_{n}\left|\bra{\psi_{n\textbf{k}}}\mathbf{\hat{q}}\cdot\mathbf{v}\ket{\psi_{n\mathbf{k}}}\right|^{2}\left(-\frac{\partial
f_{n\mathbf{k}}}{\partial E_{n\mathbf{k}}}\right),$ (2)

and

$\varepsilon^{\mathrm{inter}}(\mathbf{q},\omega)=1-\frac{4\pi}{V}\sum_{\mathbf{k}}\sum_{n;n\neq
n^{\prime}}\sum_{n^{\prime}}\frac{\left|\bra{\psi_{n\textbf{k}}}\mathbf{\hat{q}}\cdot\mathbf{v}\ket{\psi_{n^{\prime}\mathbf{k}}}\right|^{2}}{(E_{n\mathbf{k}}-E_{n^{\prime}\mathbf{k}})^{2}}\frac{f_{n\mathbf{k}}-f_{n^{\prime}\mathbf{k}}}{\omega+E_{n\mathbf{k}}-E_{n^{\prime}\mathbf{k}}+i\eta}.$
(3)

In these equations, $V$ is the volume of the system,
$|\psi_{n\mathbf{k}}\rangle$ is an electronic state with associated energy
$E_{n\mathbf{k}}$ and labelled with quantum numbers $(n,\mathbf{k})$,
$\mathbf{v}$ is the velocity operator, and $f_{n\mathbf{k}}$ is the Fermi-
Dirac distribution. We use the empirical broadening parameters
$\gamma=\eta=0.1$ eV. To obtain the reflectivity, we average the dielectric
function and the Drude plasma frequency as

$\varepsilon(\omega)=\frac{\varepsilon(\mathbf{\hat{x}},\omega)+\varepsilon(\mathbf{\hat{y}},\omega)+\varepsilon(\mathbf{\hat{z}},\omega)}{3}\hskip
28.45274pt\text{and}\hskip
28.45274pt\omega^{2}_{D}=\frac{\omega^{2}_{D}(\mathbf{\hat{x}})+\omega^{2}_{D}(\mathbf{\hat{y}})+\omega^{2}_{D}(\mathbf{\hat{z}})}{3}.$
(4)

Using the relation $\varepsilon(\omega)=[n(\omega)+ik(\omega)]^{2}$ with the
refractive index $n(\omega)$ and the optical extinction coefficient
$k(\omega)$, we compute the reflectivity at normal incidence by assuming a
vacuum-material interface as

$R(\omega)=\frac{[n(\omega)-1]^{2}+k(\omega)^{2}}{[n(\omega)+1]^{2}+k(\omega)^{2}}.$
(5)

We follow the method described in Ref. colour-marzari  to obtain the colour
from the reflectivity and we assign names to the calculated colours based on
the online tool ArtyClick artyclick . We note that the intensity of the colour
can change slightly using different exchange-correlation functionals (see
Supplementary Fig. S19).

We use the Mitsuba 3 renderer for the photorealistic rendering mitsuba3 . The
photorealistic images presented in the main text are rendered by assuming an
ideal bulk system with a clean surface. We have also considered the effects of
surface roughness on the colour and find no significant changes (see
Supplementary Fig. S18). We note that colour perception is subjective and that
colour appearance depends on many effects including surface roughness and
thickness of the sample. For this reason, we include all calculated dielectric
functions as Supplementary Data, which allows readers to independently explore
the resulting colours using their own setup.

Electron-phonon coupling calculations. - We use the finite displacement method
in conjunction with nondiagonal supercells nondiagonal_supercells  to
calculate the phonon frequencies $\omega_{\mathbf{q}\nu}$ and eigenvectors
$\mathbf{e}_{\mathbf{q}\nu}$ of a phonon mode labelled by wave vector
$\mathbf{q}$ and branch $\nu$. The electronic structure parameters are the
same as those reported above, and we use a $4\times 4\times 4$ coarse
$\mathbf{q}$-point grid to construct the matrix of force constants.
Representative phonon dispersions are reported in Supplementary Note 6. We
evaluate the imaginary part of the dielectric function at temperature $T$
renormalized by electron-phonon coupling using the Williams-Lax theory
Williams1951 ; Lax1952 :

$\varepsilon_{2}(\omega;T)=\frac{1}{\mathcal{Z}}\sum_{\mathbf{s}}\left\langle\Phi_{\mathbf{s}}(\mathbf{u})\left|\varepsilon_{2}(\omega;\mathbf{u})\right|\Phi_{\mathbf{s}}(\mathbf{u})\right\rangle
e^{-E_{\mathbf{s}}/k_{\mathrm{B}}T},$ (6)

where $\mathcal{Z}$ is the partition function,
$|\Phi_{\mathbf{s}}(\mathbf{u})\rangle$ is a harmonic eigenstate $\mathbf{s}$
of energy $E_{\mathbf{s}}$, $\mathbf{u}=\\{u_{\mathbf{q}\nu}\\}$ is a vector
containing all atomic positions expressed in terms of normal mode amplitudes
$u_{\mathbf{q}\nu}$, and $k_{\mathrm{B}}$ is Boltzmann’s constant. We evaluate
Eq. (6) by Monte Carlo integration accelerated with thermal lines
thermal_lines : we generate atomic configurations in which the atoms are
distributed according to the harmonic nuclear wave function in which every
normal mode has an amplitude of
$\left(\frac{1}{2\omega_{\mathbf{q}\nu}}\left[1+2n_{\mathrm{B}}(\omega_{\mathbf{q}\nu},T)\right]\right)^{1/2}$,
where $n_{\mathrm{B}}(\omega,T)$ is the Bose-Einstein factor. We note that the
two terms in the Bose-Einstein factor imply that the electron-phonon
renormalised dielectric function includes the effects of both quantum zero-
point nuclear vibrations at $T=0$ K and thermal nuclear vibrations at finite
temperature. We build the electron-phonon renormalised reflectivity using the
electron-phonon renormalised dielectric function.

For the superconductivity calculations we evaluate the electronic and
vibrational properties using Quantum Espresso quantum_espresso  together with
GBRV ultrasoft pseudopotentials gbrv_pseudos  and the PBE exchange-correlation
functional pbe . We use a plane-wave cutoff of $50$ Ry, a $4\times 4\times 4$
$\mathbf{k}$-point grid, and a Gaussian smearing of $0.02$ Ry for geometry
optimisations and electronic structure calculations. We use phonon
$\mathbf{q}$-grids of size $4\times 4\times 4$ and dense electron
$\mathbf{k}$-grids of size $8\times 8\times 8$ for the superconductivity
calculations. In the calculation of the Eliashberg spectral function,

$\alpha^{2}F(\omega)=\frac{1}{2}\sum_{\nu}\int_{\text{BZ}}\frac{d\mathbf{q}}{\Omega_{\text{BZ}}}\lambda_{\mathbf{q}\nu}\delta{(\omega-\omega_{\mathbf{q}\nu})},$
(7)

the integral is calculated by a sum over the $\mathbf{q}$-grids and the Dirac-
delta functions are replaced by Gaussians with a width of 0.1 THz, where
$\lambda_{\mathbf{q}\nu}$ is the electron-phonon coupling parameter for
vibrational mode $\nu$ at $\mathbf{q}$. These parameters provide converged
values of $\lambda$ and $\omega_{\text{log}}$ derived from moments of
$\alpha^{2}F(\omega)$. We estimate the superconducting critical temperature
$T^{\text{AD}}_{c}$ using the McMillan-Allen-Dynes formula TC_AD_formula ,

$k_{\mathrm{B}}T_{c}^{\mathrm{AD}}\,=\,\frac{\omega_{\mathrm{log}}}{1.2}\,\exp\Biggl{[}-\,\frac{1.04\,(1+\lambda)}{\lambda(1-0.62\,\mu^{*})-\mu^{*}}\,\Biggr{]},$
(8)

where $\mu^{*}=0.125$. We also calculate $T^{\text{E}}_{c}$ from a numeric
solution to the isotropic Eliashberg equations. The values of $T_{c}$ do not
change significantly for a range of $\mu^{*}$ values from 0.06 to 0.18.

Structure searches - In the study of the lutetium-hydrogen binary system, we
perform structure searches using the same Ephemeral Data-Derived Potential
(EDDP) Pickard2022  as used by some of us in Ferreira and co-workers pickard-
structure-prediction-arxiv . To generate the initial binary structures, we
remove hydrogen atoms in cubic LuH3. Starting with a $2\times 2\times 2$
supercell of the pristine cubic LuH3 structure, we enumerate all possible
symmetrically inequivalent defect structures with hydrogen vacancies using the
disorder code code_disorder_1 ; code_disorder_2 . This results in 55,066
unique structures for compositions in the range LuH0 to LuH3, reducing the
original count from 17,777,214 through symmetry considerations. These
structures are then optimised using the potential. The low energy structures
resulting from the machine learned potential relaxation are carried forward
for subsequent DFT calculations.

We further use the ab initio random structure searching (AIRSS) method
Pickard2006 ; Pickard_AIRSS  to explore lutetium-hydrogen-nitrogen ternary
compounds with cubic symmetry matching that reported in experimental X-ray
diffraction data. We generate 1,000 cubic structures by randomly replacing all
Lu and H sites with H, N, or a vacancy, which supplements the 200,000 randomly
generated Lu-H-N structures sampled in previous work pickard-structure-
prediction-arxiv . We then relax these additional structures using DFT.

##  VI Data availability

The data that support the findings of this study are available within the
paper and Supplementary Information. In particular, the dataset for dielectric
functions and structure files can be found at
https://github.com/monserratlab/Lu_hydrides_colour.

##  VII Code availability

The VASP code used in this study is a commercial electronic structure modeling
software, available from https://www.vasp.at. The Quantum Espresso code used
in this research is open source: https://www.quantum-espresso.org/.

##  VIII References

## References

  * (1) Ashcroft, N. W.  Metallic hydrogen: A high-temperature superconductor?  _Phys. Rev. Lett._ 21, 1748–1749 (1968).  URL http://link.aps.org/doi/10.1103/PhysRevLett.21.1748.
  * (2) Ashcroft, N. W.  Hydrogen dominant metallic alloys: High temperature superconductors?  _Phys. Rev. Lett._ 92, 187002 (2004).  URL https://link.aps.org/doi/10.1103/PhysRevLett.92.187002.
  * (3) Duan, D. _et al._ Pressure-induced metallization of dense (H2S)2H2 with high-Tc superconductivity.  _Scientific Reports_ 4, 6968 (2014).  URL https://www.nature.com/articles/srep06968.
  * (4) Liu, H., Naumov, I. I., Hoffmann, R., Ashcroft, N. W. & Hemley, R. J.  Potential high-Tc superconducting lanthanum and yttrium hydrides at high pressure.  _Proceedings of the National Academy of Sciences_ 114, 6990–6995 (2017).  URL https://www.pnas.org/doi/abs/10.1073/pnas.1704505114.
  * (5) Peng, F. _et al._ Hydrogen clathrate structures in rare earth hydrides at high pressures: Possible route to room-temperature superconductivity.  _Phys. Rev. Lett._ 119, 107001 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevLett.119.107001.
  * (6) Liu, H. _et al._ Dynamics and superconductivity in compressed lanthanum superhydride.  _Phys. Rev. B_ 98, 100102 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.100102.
  * (7) Pickard, C. J., Errea, I. & Eremets, M. I.  Superconducting hydrides under pressure.  _Annual Review of Condensed Matter Physics_ 11, 57 (2020).  URL https://doi.org/10.1146/annurev-conmatphys-031218-013413.
  * (8) Shipley, A. M., Hutcheon, M. J., Needs, R. J. & Pickard, C. J.  High-throughput discovery of high-temperature conventional superconductors.  _Phys. Rev. B_ 104, 054501 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.104.054501.
  * (9) Pickett, W. E.  Colloquium: Room temperature superconductivity: The roles of theory and materials design.  _Rev. Mod. Phys._ 95, 021001 (2023).  URL https://link.aps.org/doi/10.1103/RevModPhys.95.021001.
  * (10) Drozdov, A. P., Eremets, M. I., Troyan, I. A., Ksenofontov, V. & Shylin, S. I.  Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system.  _Nature_ 525, 73 (2015).  URL https://doi.org/10.1038/nature14964.
  * (11) Somayazulu, M. _et al._ Evidence for superconductivity above 260 K in lanthanum superhydride at megabar pressures.  _Phys. Rev. Lett._ 122, 027001 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.027001.
  * (12) Drozdov, A. P. _et al._ Superconductivity at 250 K in lanthanum hydride under high pressures.  _Nature_ 569, 528 (2019).  URL https://doi.org/10.1038/s41586-019-1201-8.
  * (13) Dasenbrock-Gammon, N. _et al._ Evidence of near-ambient superconductivity in a N-doped lutetium hydride.  _Nature_ 615, 244 (2023).  URL https://doi.org/10.1038/s41586-023-05742-0.
  * (14) Shan, P. _et al._ Pressure-induced color change in the lutetium dihydride LuH2.  _Chinese Physics Letters_ 40, 046101 (2023).  URL https://dx.doi.org/10.1088/0256-307X/40/4/046101.
  * (15) Zhao, X. _et al._ Pressure tuning of optical reflectivity in LuH2.  _Science Bulletin_ 68, 883–886 (2023).  URL https://www.sciencedirect.com/science/article/pii/S2095927323002499.
  * (16) Ming, X. _et al._ Absence of near-ambient superconductivity in LuH2±xNy.  _Nature_ 620, 72–77 (2023).  URL https://doi.org/10.1038/s41586-023-06162-w.
  * (17) Zhang, Y.-J. _et al._ Pressure induced color change and evolution of metallic behavior in nitrogen-doped lutetium hydride.  _Science China Physics, Mechanics & Astronomy_ 66, 287411 (2023).  URL https://doi.org/10.1007/s11433-023-2109-4.
  * (18) Xing, X. _et al._ Observation of non-superconducting phase changes in nitrogen doped lutetium hydrides.  _Nature Communications_ 14, 5991 (2023).  URL https://doi.org/10.1038/s41467-023-41777-7.
  * (19) Cai, S. _et al._ No evidence of superconductivity in the compressed sample prepared from the lutetium foil and H2/N2 gas mixture.  _Matter and Radiation at Extremes_ 8, 048001 (2023).  URL https://doi.org/10.1063/5.0153447.
  * (20) Zhang, S. _et al._ Electronic and magnetic properties of Lu and LuH2.  _AIP Advances_ 13, 065117 (2023).  URL https://doi.org/10.1063/5.0153011.
  * (21) Wang, N. _et al._ Percolation-induced resistivity drop in lutetium dihydride with controllable electrical conductivity over six orders of magnitude.  _Science China Physics, Mechanics & Astronomy_ 66, 297412 (2023).  URL https://doi.org/10.1007/s11433-023-2171-8.
  * (22) Moulding, O. _et al._ Trigonal to cubic structural transition in possibly N-doped LuH3 measured by Raman and X-ray diffraction.  _arXiv:2304.04310_ (2023).  URL https://arxiv.org/abs/2304.04310.
  * (23) Li, Z. _et al._ Superconductivity above 70 K observed in lutetium polyhydrides.  _Science China Physics, Mechanics and Astronomy_ 66, 267411 (2023).  URL https://doi.org/10.1007/s11433-023-2101-9.
  * (24) Salke, N. P., Mark, A. C., Ahart, M. & Hemley, R. J.  Evidence for near ambient superconductivity in the Lu-N-H system.  _arXiv:2306.06301_ (2023).  URL https://arxiv.org/abs/2306.06301.
  * (25) Xie, F. _et al._ Lu–H–N phase diagram from first-principles calculations.  _Chinese Physics Letters_ 40, 057401 (2023).  URL https://dx.doi.org/10.1088/0256-307X/40/5/057401.
  * (26) Liu, M. _et al._ Parent structures of near-ambient nitrogen-doped lutetium hydride superconductor.  _Phys. Rev. B_ 108, L020102 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.108.L020102.
  * (27) Huo, Z. _et al._ First-principles study on the conventional superconductivity of N-doped fcc-LuH3.  _Matter and Radiation at Extremes_ 8, 038402 (2023).  URL https://doi.org/10.1063/5.0151844.
  * (28) Hilleke, K. P. _et al._ Structure, stability, and superconductivity of N-doped lutetium hydrides at kbar pressures.  _Phys. Rev. B_ 108, 014511 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.108.014511.
  * (29) Ferreira, P. P. _et al._ Search for ambient superconductivity in the Lu-N-H system.  _Nature Communications_ 14, 5367 (2023).  URL https://doi.org/10.1038/s41467-023-41005-2.
  * (30) Sun, Y., Zhang, F., Wu, S., Antropov, V. & Ho, K.-M.  Effect of nitrogen doping and pressure on the stability of LuH3.  _Phys. Rev. B_ 108, L020101 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.108.L020101.
  * (31) Lucrezi, R., Ferreira, P. P., Aichhorn, M. & Heil, C.  Temperature and quantum anharmonic lattice effects in lutetium trihydride: stability and superconductivity.  _arXiv:2304.06685_ (2023).  URL https://arxiv.org/abs/2304.06685.
  * (32) Tian, Y. C. _et al._ Ultrafast dynamics evidence of high temperature superconductivity in single unit cell FeSe on SrTiO3.  _Phys. Rev. Lett._ 116, 107001 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevLett.116.107001.
  * (33) Wu, Q. _et al._ Ultrafast quasiparticle dynamics and electron-phonon coupling in (Li0.84Fe0.16)OHFe0.98Se.  _Chinese Physics Letters_ 37, 097802 (2020).  URL https://dx.doi.org/10.1088/0256-307X/37/9/097802.
  * (34) Zhao, J., Bragas, A. V., Lockwood, D. J. & Merlin, R.  Magnon squeezing in an antiferromagnet: Reducing the spin noise below the standard quantum limit.  _Phys. Rev. Lett._ 93, 107203 (2004).  URL https://link.aps.org/doi/10.1103/PhysRevLett.93.107203.
  * (35) Zhao, J., Bragas, A. V., Merlin, R. & Lockwood, D. J.  Magnon squeezing in antiferromagnetic MnF2 and FeF2.  _Phys. Rev. B_ 73, 184434 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevB.73.184434.
  * (36) Kresse, G. & Furthmüller, J.  Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set.  _Computational Materials Science_ 6, 15–50 (1996).  URL https://www.sciencedirect.com/science/article/pii/0927025696000080.
  * (37) Kresse, G. & Furthmüller, J.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  _Phys. Rev. B_ 54, 11169–11186 (1996).  URL https://link.aps.org/doi/10.1103/PhysRevB.54.11169.
  * (38) Blöchl, P. E.  Projector augmented-wave method.  _Phys. Rev. B_ 50, 17953–17979 (1994).  URL https://link.aps.org/doi/10.1103/PhysRevB.50.17953.
  * (39) Perdew, J. P. _et al._ Restoring the density-gradient expansion for exchange in solids and surfaces.  _Phys. Rev. Lett._ 100, 136406 (2008).  URL https://link.aps.org/doi/10.1103/PhysRevLett.100.136406.
  * (40) Haule, K., Yee, C.-H. & Kim, K.  Dynamical mean-field theory within the full-potential methods: Electronic structure of CeIrIn5, CeCoIn5, and CeRhIn5.  _Phys. Rev. B_ 81, 195107 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevB.81.195107.
  * (41) Haule, K.  Structural predictions for correlated electron materials using the functional dynamical mean field theory approach.  _Journal of the Physical Society of Japan_ 87, 041005 (2018).  URL https://doi.org/10.7566/JPSJ.87.041005.
  * (42) Blaha, P. _et al._ WIEN2k: An APW+lo program for calculating the properties of solids.  _The Journal of Chemical Physics_ 152, 074101 (2020).  URL https://doi.org/10.1063/1.5143061.
  * (43) Prandini, G., Rignanese, G.-M. & Marzari, N.  Photorealistic modelling of metals from first principles.  _npj Computational Materials_ 5, 129 (2019).  URL https://doi.org/10.1038/s41524-019-0266-0.
  * (44) Marini, A., Onida, G. & Del Sole, R.  Plane-wave DFT-LDA calculation of the electronic structure and absorption spectrum of copper.  _Phys. Rev. B_ 64, 195125 (2001).  URL https://link.aps.org/doi/10.1103/PhysRevB.64.195125.
  * (45) Harl, J.  _The linear response function in density functional theory: optical spectra and improved description of the electron correlation_.  Ph.D. thesis, Wien (2008).  URL https://doi.org/10.25365/thesis.2622.  Universität Wien, Dissertation.
  * (46) https://colors.artyclick.com/color-hue-finder.
  * (47) Jakob, W. _et al._ Mitsuba 3 renderer (2022).  Https://mitsuba-renderer.org.
  * (48) Lloyd-Williams, J. H. & Monserrat, B.  Lattice dynamics and electron-phonon coupling calculations using nondiagonal supercells.  _Phys. Rev. B_ 92, 184301 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevB.92.184301.
  * (49) Williams, F.  Theoretical low temperature spectra of the thallium activated potassium chloride phosphor.  _Phys. Rev._ 82, 281–282 (1951).  URL https://journals.aps.org/pr/pdf/10.1103/PhysRev.82.281.2.
  * (50) Lax, M.  The Franck-Condon principle and its application to crystals.  _J. Chem. Phys._ 20, 1752–1760 (1952).  URL https://aip.scitation.org/doi/pdf/10.1063/1.1700283?class=pdf.
  * (51) Monserrat, B.  Vibrational averages along thermal lines.  _Phys. Rev. B_ 93, 014302 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.93.014302.
  * (52) Giannozzi, P. _et al._ Advanced capabilities for materials modelling with Quantum ESPRESSO.  _Journal of Physics: Condensed Matter_ 29, 465901 (2017).  URL https://dx.doi.org/10.1088/1361-648X/aa8f79.
  * (53) Garrity, K. F., Bennett, J. W., Rabe, K. M. & Vanderbilt, D.  Pseudopotentials for high-throughput DFT calculations.  _Computational Materials Science_ 81, 446–452 (2014).  URL https://www.sciencedirect.com/science/article/pii/S0927025613005077.
  * (54) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Phys. Rev. Lett._ 77, 3865–3868 (1996).  URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3865.
  * (55) Allen, P. B. & Dynes, R. C.  Transition temperature of strong-coupled superconductors reanalyzed.  _Phys. Rev. B_ 12, 905–922 (1975).  URL https://link.aps.org/doi/10.1103/PhysRevB.12.905.
  * (56) Pickard, C. J.  Ephemeral data derived potentials for random structure search.  _Physical Review B_ 106, 014102 (2022).
  * (57) Lian, J.-C., Wu, H.-Y., Huang, W.-Q., Hu, W. & Huang, G.-F.  Algorithm for generating irreducible site-occupancy configurations.  _Phys. Rev. B_ 102, 134209 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.134209.
  * (58) Lian, J.-C. _et al._ Highly efficient tree search algorithm for irreducible site-occupancy configurations.  _Phys. Rev. B_ 105, 014201 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.105.014201.
  * (59) Pickard, C. J. & Needs, R. J.  High-pressure phases of silane.  _Phys. Rev. Lett._ 97, 045504 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevLett.97.045504.
  * (60) Pickard, C. J. & Needs, R. J.  Ab initio random structure searching.  _Journal of Physics: Condensed Matter_ 23, 053201 (2011).  URL https://dx.doi.org/10.1088/0953-8984/23/5/053201.

##  IX Acknowledgments

###### Acknowledgements.

S.-W.K. and B.M. are supported by a UKRI Future Leaders Fellowship
[MR/V023926/1]. B.M. also acknowledges support from the Gianna Angelopoulos
Programme for Science, Technology, and Innovation, and from the Winton
Programme for the Physics of Sustainability. G.L.P. acknowledges funding from
the Ministry of Research, Innovation, and Digitalisation within Program
1-Development of National Research and Development System, Subprogram
1.2-Institutional Performance-RDI Excellence Funding Projects, under contract
no. 10PFE/2021. The computational resources were provided by the Cambridge
Tier-2 system operated by the University of Cambridge Research Computing
Service and funded by EPSRC [EP/P020259/1], by the UK National Supercomputing
Service ARCHER2, for which access was obtained via the UKCP consortium and
funded by EPSRC [EP/X035891/1], and by the SCARF cluster of the STFC
Scientific Computing Department. For the purpose of open access, the authors
have applied a Creative Commons Attribution (CC BY) licence to any Author
Accepted Manuscript version arising from this submission.

##  X Author contributions

B.M. and S.-W.K. conceived the study and planned and supervised the research.
S.-W.K. performed the DFT calculations, S.-W.K. and G.L.P. performed the DMFT
calculations, L.J.C. and C.J.P. performed the structure searches, and all
authors contributed to the analysis. B.M. and S.-W.K. wrote the manuscript
with input from all authors.

##  XI Competing interests

The authors declare no competing interests.

in 1,…,0 See pages \x of SI_final.pdf