Photorealistic modelling of metals from first principles

TL;DR

This paper presents a first-principles computational framework to accurately predict the optical properties and colours of various metals and alloys, validated against experimental data and capable of photorealistic rendering.

Contribution

It introduces a comprehensive, first-principles method for predicting metallic reflectivity and colour, including complex alloys, with validation against experimental results.

Findings

Accurate prediction of metal reflectivity and colour from first principles.

Successful validation with experimental data and photorealistic rendering.

Applicable to elemental metals and various metallic compounds.

Abstract

The colours of metals have attracted the attention of humanity since ancient times, and coloured metals, in particular gold compounds, have been employed for tools and objects symbolizing the aesthetics of power. In this work we develop a comprehensive framework to obtain the reflectivity and colour of metals, and show that the trends in optical properties and the colours can be predicted by straightforward first-principles techniques based on standard approximations. We apply this to predict reflectivity and colour of several elemental metals and of different types of metallic compounds (intermetallics, solid solutions and heterogeneous alloys), considering mainly binary alloys based on noble metals. We validate the numerical approach through an extensive comparison with experimental data and the photorealistic rendering of known coloured metals.

Figures

Tables

Table 1. Table 1: Different types of compounds and corresponding simulation method used in this work for the first-principles simulation of these systems.

Binary alloys	Computational method
Intermetallic compounds	Primitive cell
Solid solutions	Supercell (SQS)
Heterogeneous alloys	Bruggeman model

Table 2. Table 2: Computed values of the IPA Drude plasma frequency ω D subscript 𝜔 D \omega_{\text{D}} (in eV) compared to previous simulations (J. Harl  Harl2008 ) and experiments (Exp.). The experimental values with no explicit reference are extracted from the data reported in Ref.  Harl2008 . For transition metals there are no experimental data available because, due to the presence of interband transitions even at vanishingly small frequencies, the Drude plasma frequency cannot be extracted by fitting experimental optical data to the Drude model, even at very low energies.

Element	This work	J. Harl Harl2008	Exp.
Cu	8.8	9.1	8.8, 8.9
Ag	8.9	9.2	8.9 $\pm$ 0.2 Yang2015 , 8.9
Au	8.6	9.0	8.45 Olmon2012 , 8.7
Li	6.4	6.5	6.4
Na	6.0	5.9	5.7
Ca	4.1	4.3	5.7
Al	12.5	12.6	12.3 , 12.5
Rh	9.6	10.1
Pd	7.0	7.4
Pt	8.4	8.8

Table 3. Table 3: Colour differences in CIELAB space between simulated colours (present work) and experimental colours  Vishnubhatla1967 ; Keast2013 ; Steinemann1997 ; Amiotti1990 derived from reflectivity data, Δ ​ E exp Δ subscript 𝐸 exp \Delta E_{\text{exp}} , and between simulated colours (present work) and previously published simulations  Keast2013 ; Keast2011 , Δ ​ E sim Δ subscript 𝐸 sim \Delta E_{\text{sim}} .

Compound	$\Delta E_{\text{exp}}$	$\Delta E_{\text{sim}}$
AuAl2	11 Vishnubhatla1967	8 Keast2013 , 4 Keast2011
AuGa2	2 Vishnubhatla1967	1 Keast2011
AuIn2	4 Vishnubhatla1967	1 Keast2011
PtAl2	12 Keast2013	2 Keast2013
CoSi2	3 Steinemann1997
NiSi2	5 Amiotti1990
PdIn	10 Steinemann1997

Equations

$${\epsilon }^{\text{inter}}(\hat{\mathbf{q}},\omega )$$

$${\epsilon }^{\text{intra}}(\hat{\mathbf{q}},\omega )$$

$${\omega }_{\text{D}}^2(\hat{\mathbf{q}})=\frac{4\pi }{V}\sum _{\mathbf{k}}\sum _n|⟨{\psi }_{n\mathbf{k}}|\hat{\mathbf{q}}\cdot \mathbf{v}|{\psi }_{n\mathbf{k}}⟩{|}^2\left(-\frac{\partial f_{n\mathbf{k}}}{\partial E_{n\mathbf{k}}}\right).$$

$$\epsilon (\omega )=\frac{\epsilon (\hat{\mathbf{x}},\omega )+\epsilon (\hat{\mathbf{y}},\omega )+\epsilon (\hat{\mathbf{z}},\omega )}{3},$$

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

$$X$$

$$Y$$

$$Z$$

$$\Delta E=\sqrt{(L_1^{\ast }-L_2^{\ast }{)}^2+(a_1^{\ast }-a_2^{\ast })+(b_1^{\ast }-b_2^{\ast }{)}^2}.$$

$$(1-x_{\beta })\frac{{\epsilon }_{\alpha }(\omega )-{\epsilon }_{\text{Br}}(\omega )}{{\epsilon }_{\alpha }(\omega )+2{\epsilon }_{\text{Br}}(\omega )}+x_{\beta }\frac{{\epsilon }_{\beta }(\omega )-{\epsilon }_{\text{Br}}(\omega )}{{\epsilon }_{\beta }(\omega )+2{\epsilon }_{\text{Br}}(\omega )}=0,$$

Full text

I Introduction

Optical properties of metals are important for novel technological applications where the optical response needs to be engineered for specific purposes, such as for plasmonic devices (e.g. in spectrally-selective coatings Bilokur2017 ; Guo2014 ), for optoelectronics devices (e.g. in ultra-thin films for transparent conductive electrodes Yun2017 ; Ren2015 ), and also for microscopy and optical data storage Hatwar1997 . Also, the colours of metals (which are related to the optical properties within the visible range of the electromagnetic spectrum) play a significant role in the jewellery industry, decoration and dentistry. For these applications, the most used materials are metallic alloys based on gold or other coinage or precious metals, such as silver, copper, palladium and platinum. In particular, gold and copper are the few elemental metals that show a characteristic colour, due to the presence of a drop in the reflectivity curve inside the visible range; reflectivities of nearly all other metals are instead generally high and flat for all visible frequencies, making them appear shiny and silvery white. Moreover, gold alloys and intermetallics are known to show a broad spectrum of colours (yellow, red, purple, to name a few), which can be tuned by varying the alloying elements and concentrations in the material Cretu1999 . Since in the jewellery industry there is the demand, due to market and fashion trends, for new precious-metal alloys with specific colours, it is also of great interest there the search and identification of novel alloys with novel optical properties.

Generally speaking, the common route followed by researchers and manufacturers in order to identify any type of novel materials is through trial-and-error experiments which, however, have the drawback to be time-consuming and, if dealing with precious-metal-based systems, expensive. In order to streamline this process, an alternative route that can help in guiding the search for new promising candidate systems is computational modelling, so that the physical properties under investigation are assessed through computer simulations rather than by real experiments. Here, we show and discuss how it is possible to perform photorealistic simulations of metals by means of first-principles methods and, as a consequence, predict the colour of novel metallic alloys. Previously published studies about first-principles simulation of optical properties of both elemental metals and alloys already point towards the feasibility of this approach. Indeed, in 1988 Maksimov et al. Maksimov1988 computed the optical properties of several elemental metals whereas, more recently, Werner et al. Werner2009 performed a similar study on 17 elemental metals; both studies found qualitative agreement with experimental results. For compounds, on the other hand, Blaber et al. Blaber2009 calculated the optical properties of several intermetallic compounds, with a particular focus on alkali-noble intermetallics, for new possible candidates as plasmonic materials while, in another work, Keast et al. Keast2014 computed the density of states and dielectric function of gold intermetallics compounds and gold binary alloys. Regarding the simulation of specific coloured intermetallic compounds, the reflectivity and colour of the three well-known coloured gold intermetallics AuAl2, AuGa2 and AuIn2 was first computed in Ref. Keast2011 and, afterwards, Keast et al. Keast2013 studied the influence of alloying concentrations on the reflectivity and colour of the intermetallic AuAl2 by considering the ternary compounds having the Au1-xPtxAl2 composition, with $x=0.0,0.5,0.75,1.0$; equivalent computational results were obtained independently by Kecik Kecik2013 . Calculated and experimental Vishnubhatla1967 ; Furrer2014 reflectivity curves and colours for these compounds showed good agreement and the trends in colour as a function of the composition were well reproduced. In addition, the effect of disorder on the optical properties of Au0.5Cu0.5 was studied by comparing the dielectric function of the random solid solution, simulated using the supercell approach, with that of the ordered intermetallic compound DeSilva2015 , and the main spectral differences between the two different types of compounds were captured by the simulations.

In this work, first we establish a general computational approach that can be used for the photorealistic simulation of metals, showing how the reflectivity and colour of metallic crystals can be estimated by means of first-principles techniques. We then demonstrate through a systematic study on elemental metals and extensive comparisons with experimental data that the theoretical and numerical approximations adopted are able to reproduce the correct behaviour of the reflectivity curves and to capture the main differences in optical properties across the periodic table. Finally, we perform a similar study on metal alloys by considering different types of compounds, i.e. ordered intermetallics, disordered solid solutions and heterogeneous alloys. In particular, we show through a comparison with experimental results that, if the appropriate methods are used for the simulation of the different types of compounds, (i) the simulated colours of known coloured intermetallics are in qualitative and most often in quantitative agreement with experiments and that (ii) one can reproduce the main colour trends in noble-metal-based binary alloys.

II Results

II.1 Computational approach

The computational workflow that allows one to obtain the reflectivity and colour of a given metal from an initial crystal structure, schematically depicted in Fig. 1, can be divided into four main computational steps: (i) evaluation of the electronic structure, (ii) calculation of the dielectric function, (iii) calculation of the reflectivity and colour and (iv) photorealistic rendering of the material.

The quantity we consider for the first-principles simulation of optical properties is the complex, wavevector- and frequency-dependent, dielectric function $\varepsilon(\mathbf{q},\omega)=\varepsilon_{1}(\mathbf{q},\omega)+i\varepsilon_{2}(\mathbf{q},\omega)$. In fact, the knowledge of the dielectric function gives then access to all the optical constants measurable by optical experiments, such as absorption coefficient and reflectivity.

Throughout this work the electronic structure is computed using density-functional theory (DFT) Kohn1964dft within the generalized gradient approximation (GGA) and relying on the PBE exchange-correlation functional Perdew1996 . We emphasize here that the electronic structure could alternatively be obtained with more accurate techniques; for example, the accuracy of the band structures could be improved by computing quasi-particle corrections on top of PBE results (typically at the $GW$ level Hybertsen1986 ; Onida2002 ; Reining2017 ), albeit at a largely increased computational cost. So, while in the present work we just rely on the Kohn-Sham (KS) PBE bands Kohn1965 , the use of conceptually and quantitatively correct $GW$ bands would take place seamlessly inside this workflow. Subsequently, we calculate the dielectric function within the independent particle approximation (IPA), which amounts to neglecting (i) effects related to electron-hole interactions (excitonic effects), since these are effectively screened by the conduction electrons and (ii) effects related to the rapidly varying microscopic electric fields inside the material (local-field effects) since these are typically small in homogeneous systems such as bulk metals Marini2001 . In the optical regime the momentum $\mathbf{q}$ transferred by the photon is negligible so that we can consider the optical limit, $\mathbf{q}\to\mathbf{0}$, of the expression for the IPA dielectric function $\varepsilon(\mathbf{q},\omega)$. In general, the dielectric function still depends on the direction $\hat{\mathbf{q}}=\mathbf{q}/|\mathbf{q}|$ of the perturbing electric field and only for crystals with cubic symmetry it is the same in every direction. In the optical limit it is convenient to divide the evaluation of the IPA dielectric function of metals into two separate contributions, an intraband Drude-like term $\varepsilon^{\text{intra}}(\hat{\mathbf{q}},\omega)$ due to the conduction electrons at the Fermi surface and an interband term $\varepsilon^{\text{inter}}(\hat{\mathbf{q}},\omega)$ due to vertical transitions between occupied and unoccupied bands, so that $\varepsilon(\hat{\mathbf{q}},\omega)=\varepsilon^{\text{inter}}(\hat{\mathbf{q}},\omega)+\varepsilon^{\text{intra}}(\hat{\mathbf{q}},\omega)$. Using the solutions of the one-particle Schrödinger equation for periodic systems, $H^{\text{KS}}\ket{\psi_{n\mathbf{k}}}=E_{n\mathbf{k}}\ket{\psi_{n\mathbf{k}}}$ (where $H^{\text{KS}}$ is the KS Hamiltonian from DFT), the explicit expression of the IPA dielectric function can be written as Wooten1972 ; Marini2001 ; Harl2008

[TABLE]

where we have defined the IPA Drude plasma frequency as

[TABLE]

In the expressions above $\mathbf{v}=-i\,[\mathbf{r},H^{\text{KS}}]$ is the velocity operator, $f_{n\mathbf{k}}$ is the Fermi-Dirac occupation function of the KS Bloch state $\ket{\psi_{n\mathbf{k}}}$ identified by band index $n$ and wavevector $\mathbf{k}$ within the Brillouin zone (BZ) and $V$ is the volume of the crystal. Instead $\eta$ is an infinitesimal broadening introduced to perform the adiabatic switch of the perturbation within linear-response theory and, in practical calculations, it is used as an empirical broadening which accounts for scattering processes, always present in real materials, and/or for finite experimental resolution. Similarly, $\gamma$ is an empirical broadening representing dissipation effects of the conduction electrons (see the Methods section for more details on the parameters effectively used in the simulations). A more extensive discussion on the first-principles theory of optical properties and the derivation of the expression of the IPA dielectric function in the optical limit (Eq. 1, Eq. 2 and Eq. 3) can be found in Ref. Prandini2019 .

As the most typical experimental situation is to have polycrystalline materials in which grains have random orientations, in the following we always deal with the dielectric function averaged over the three Cartesian directions

[TABLE]

so that we can drop the dependence on the direction $\hat{\mathbf{q}}$. Similarly we also define a corresponding average IPA Drude plasma frequency as $\omega^{2}_{\text{D}}=[\omega^{2}_{\text{D}}(\hat{\mathbf{x}})+\omega^{2}_{\text{D}}(\hat{\mathbf{y}})+\omega^{2}_{\text{D}}(\hat{\mathbf{z}})]/3$.

In order to compute the reflectivity from the knowledge of the dielectric function we first introduce the refractive index $n(\omega)$ and the extinction coefficient $k(\omega)$ that are defined from the equation $[n(\omega)+ik(\omega)]^{2}=\varepsilon(\omega)$. The reflectivity at normal incidence and assuming a vacuum-material interface is then simply linked to $n(\omega)$ and $k(\omega)$ through the Fresnel equations of classical electromagnetism (see for example Ref. Griffiths2007 ):

[TABLE]

Eventually, we relate the reflectivity of a material to its perceived colour using the standard colour spaces introduced by the Commission Internationale de l’Eclairage (CIE) cie_webpage for quantitative measures of colour. For this purpose trichromatic theory gives the rigorous mathematical framework that permits to estimate the colour of an opaque material (e.g. a metal) by knowing its reflectivity $R(\lambda)$ for all the wavelengths $\lambda$ in the visible range (i.e. in the range [380, 780] nm), and to condense this information into three numbers, i.e. the colour coordinates Schanda2007 . In particular, according to the CIE 1931 standard colorimetric observer, the tristimulus values ($X$, $Y$, $Z$) which define the CIE-$XYZ$ colour space completely describe a colour stimulus and are given by the following integrals over the visible range

[TABLE]

where $\mathrm{bar}{\text{x}}(\lambda)$, $\mathrm{bar}{\text{y}}(\lambda)$ and $\mathrm{bar}{\text{z}}(\lambda)$ are the three so-called colour-matching functions and describe the chromatic response of the observer, being related to the sensitivity of the three different colour-sensitive photoreceptors present in the human eye. $S(\lambda)$ is instead the spectral power distribution of one of the standard CIE illuminant (throughout this work the D65 illuminant is used which corresponds to average daylight) while the constant $k$ is chosen so that $Y=100$ for objects for which $R(\lambda)=1$ for all visible wavelengths.

In practice, it is more convenient to work within the CIELAB colour space rather than in the CIE-$XYZ$ colour space, which is defined by three coordinates ($L^{*}$, $a^{*}$, $b^{*}$) that are easily computed from the knowledge of the tristimulus values ($X$, $Y$, $Z$) through a coordinate transformation Schanda2007 . Indeed, since the CIELAB colour space is nearly uniform, euclidean distances can be used to approximately represent the perceived magnitude of colour differences between two objects in the same external conditions. Therefore, if ($L_{1}^{*}$, $a_{1}^{*}$, $b_{1}^{*}$) and ($L_{2}^{*}$, $a_{2}^{*}$, $b_{2}^{*}$) are the CIELAB coordinates of two objects, their colour difference is simply given by

[TABLE]

Typically, a difference $\Delta E>1-2$ can be perceived by the human eye. In addition, we use photorealistic rendering, which is based on the solution of the light-transport equation Kajiya1986 , to simulate the actual appearance of an object made of a material with specified optical constants in the visible range within a realistic 3D scene.

Our goal is to apply the computational approach described above to study metals in their crystalline form. From the point of view of first-principles calculations, elemental crystals are the easiest and most computationally efficient systems to simulate since they are periodic and their primitive cell, which typically consists of only a few atoms, can simply be taken as the simulation cell. For multi-component systems instead (in this work we focus on binary alloys), we distinguish different types of compounds according to their atomic configuration and microstructure. In particular, we consider the following three limiting cases: (i) perfectly ordered phases, i.e. pure intermetallic compounds, (ii) perfectly disordered phases, i.e. pure solid solutions, and (iii) heterogeneous alloys, i.e. alloys consisting of a mixture of two different phases. Since we are exclusively interested in the study of the intrinsic bulk colours of metals, we neglect the influence on the optical properties of defects (e.g. vacancies, dislocations, etc.) and any type of surface effects.

We use different simulation methods in order to properly model the reflectivity and colour of these three different types of compounds, as summarized in Table 1. As for the case of elemental crystals, pure intermetallic compounds are periodic systems and are simply simulated in their primitive cell. On the other hand, we use the supercell approach, based on the use of special quasi-random structures (SQS) Zunger1990 ; Wei1990 , to take into account effects related to disorder in the simulation of the optical properties of solid solutions (see Supplementary Discussion 1 for a comparison between the SQS supercell approach and the virtual-crystal approximation). Instead, for heterogeneous alloys made of two phases $\alpha$ and $\beta$, we use the Bruggeman model Niklasson1981 to estimate the optical properties of the alloy. Within the Bruggeman model the dielectric function of the mixture, that we indicate as $\varepsilon_{\text{Br}}(\omega)$, is given by the following expression

[TABLE]

in terms of the dielectric functions $\varepsilon_{\alpha}(\omega)$ and $\varepsilon_{\beta}(\omega)$ of the single phases, and where $x_{\alpha}$ and $x_{\beta}$ (with $x_{\alpha}+x_{\beta}=1$) are the fractions of the two phases present in the material. The dielectric function of the single phases can be obtained with the methods of Table 1 for intermetallic compounds and solid solutions.

In the following, we apply and validate the computational approach described here, and discuss its limitations, on several elemental metals and binary compounds.

II.2 Elemental metals

Fig. 2 shows the comparison between IPA results and experimental data for the reflectivity curves of 18 elemental metals, focusing on frequencies centered around the visible range (i.e. in the range [1.59, 3.26] eV). Experimentally, we observe high and flat reflectivities along the visible spectrum for the “precious” transition metals (i.e. Rh, Ir and Pd) while we observe flat but slightly lower reflectivities for the other transition metals considered (i.e. V, Nb, Ta, Cr, Mo and W). As a consequence, in terms of CIELAB colour coordinates, metals in the first group have a large CIELAB brightness $L^{*}$ and thus whitish colour, while the others have smaller brightness and thus a more greyish colour (e.g. rhodium has $L^{*}$=90 while vanadium has $L^{*}$=78). An interesting exception among the transition metals is osmium, that shows a reflectivity curve that is low in the low-energy part of the visible spectrum but then suddenly rises in the blue-violet part, thus giving a bluish tint to pure osmium. A similar behaviour is found also in tantalum, but the rise of the reflectivity curve in the blue-violet region is significantly smaller and, consequently, also the bluish tint of the material is less pronounced. Instead, the simple $sp$ metals lithium, potassium and aluminium all have very high and nearly flat reflectivity curves in the visible range (and therefore whitish colour) while in beryllium the intensity of the reflectivity is lower, and comparable to that of the transition metals (and thus having a greyish colour). Interestingly, the reflectivity curve of caesium decreases significantly within the visible range, so that red and yellow radiation is strongly reflected, while all other visible frequencies are absorbed, giving a yellow tint to the material. As clearly shown in Fig. 2, the IPA simulations reproduce these different features of the elemental metals. In contrast, for noble metals, while the characteristic drop in the reflectivity curve in the visible range (for Cu and Au) or in the ultraviolet (for Ag) is also reproduced by the simulations, it happens at smaller energies compared to experiments due to the well-known underestimation of the interband gap between valence $d$ bands and conduction $sp$ bands of PBE band structures. This discrepancy can be corrected using approaches beyond DFT, such as the $GW$ approximation of many-body perturbation theory. By correcting the DFT band energies at the $G_{0}W_{0}$ level, a quantitative agreement with respect to experiments is obtained for the optical spectra of Cu Marini2001a and Ag Marini2002 but not for Au, for which $G_{0}W_{0}$ gives very similar results to PBE Rangel2012 . For this latter case, the quasi-particle self-consistent $GW$ ($QSGW$) vanSchilfgaarde2006 ; Kotani2007 approach is required for the occupied $5d$ bands of gold to be lowered in energy by the right amount Rangel2012 .

A quantitative measure of the accuracy of the simulations can be obtained through the colour difference $\Delta E$ (given in Eq. 9) with respect to experiments. Its average value is found to be $<\Delta E>=6.4$. For a more qualitative visual comparison, Fig. 3 shows the simulated rendering of a metallic surface of elemental gold, osmium and caesium together with the appearance of experimental samples of the same materials. In gold, the shift of the reflectivity edge in the simulations with respect to experiments makes the rendered colour more reddish than the true red-yellow colour of pure gold. On the other hand the bluish colour of osmium and the yellow colour of caesium are well reproduced by the IPA simulations.

Moreover, as shown in Table 2, the IPA results for the Drude plasma frequency are in good agreement both with experiments and with previous simulations Harl2008 performed at the same level of theory for some elemental metals.

From all these results, we conclude that the IPA approach applied on top of PBE band structures predicts the reflectivity and colour of elemental metals surprisingly well. Although the colour is not always in quantitative agreement with experiments, the shape and the main features of the experimental reflectivity curves are reproduced in elemental metals. These results are somewhat surprising because it is known that quasi-particle corrections modify significantly the PBE band structure in metals and the corrections are k-dependent Marini2001a ; Marini2002 (i.e. they do not act as a simple scissor operator). Nonetheless, these approximated simulations manage to capture the correct features of the optical constants. This can intuitively be understood by the fact that the dielectric function is given by the sum of all possible vertical transitions over all the BZ and small differences in the positions and features of the bands (like gradient and curvature) are averaged out in the spectra. In the special case of noble metals, the position of the occupied $d$ bands in PBE is not correct and, since there are no other allowed interband transitions in that energy range, the onset of interband optical absorption (i.e. $\varepsilon_{2}^{\text{inter}}(\omega)$) in PBE is also not at the correct position (similar to the case of semiconductors for which the PBE band gap is systematically underestimated Onida2002 ). On the other hand, the shape of $\varepsilon(\omega)$ for noble metals is reasonably well reproduced.

II.3 Alloys

In order to validate the theoretical approach used on binary compounds, we first compare the reflectivity and colour between simulations and experiments for known coloured intermetallic compounds, as previously done for elemental metals. Second, we check the predictive accuracy of the simulations by studying in noble-metal-based alloys both the trends in reflectivity with respect to composition (in Ag-Au and Ag-Cu) and the differences in optical properties among different types of compounds for a given alloy composition (in Au-Cu and Ag-Cu).

Intermetallics.

We first simulate the reflectivity and colour of intermetallic compounds that are experimentally known to be coloured. The compounds studied are the purple AuAl2, blue AuIn2, bluish AuGa2, yellow PtAl2, red PdIn, blue-grey NiSi2 and dark blue CoSi2 Steinemann1997 ; Cretu1999 . All these intermetallics have cubic symmetry: AuAl2, AuGa2, AuIn2, PtAl2, CoSi2 and NiSi2 crystallize in the FCC CaF2 prototype structure (space group Fm$\mathrm{bar}{3}$m) while PdIn crystallizes in the BCC CsCl prototype structure (space group Pm$\mathrm{bar}{3}$m).

As shown in Fig. 4, the experimental shape of the reflectivity curve for the coloured intermetallics is well reproduced by the simulations. The colour differences between simulations and experiments are summarized in Table 3, where the comparison with other first-principles simulations Keast2011 ; Keast2013 is also reported. The agreement with previous simulations is satisfactory and, moreover, we reproduce the true colour of the intermetallic compounds studied (although the CIELAB brightness is typically overestimated by the simulations). For example, the comparison between photorealistic rendering and real material samples clearly shows that the simulations predict the correct colours of purple AuAl2, bluish AuGa2 and yellow PtAl2 (see Fig. 5). The characteristic colours of these highly symmetric intermetallic compounds are due to selective optical absorption in confined regions of the visible spectrum Steinemann1997 . For the gold compounds, the optical absorption inside the visible range is given by transitions from $sp$ conduction states below the Fermi level to unoccupied states above the Fermi level. The bands originating from the $5d$ states of gold, that are problematic in the study of elemental gold, are located at $\sim 5$ eV below the Fermi level and these do not contribute to the characteristic colours of these compounds Keast2015 . This explains the better agreement with experiments found for the gold intermetallics compounds compared to the case of elemental gold.

Au-Ag-Cu.

The Au-Ag-Cu system is an ideal test case for the application of the computational approach described above to alloys since (i) several experimental optical data on this system are available, (ii) its constituent binaries show very different behaviours in terms of phase stability and so different types of compounds are observed and (iii) it is the basis of the most common jewellery and dental alloys in use today. Concerning the phase stability of the constituent binaries, Ag is completely soluble in Au thus Au and Ag form solid solutions for each composition and no long-range order is observed at low temperatures. Also Au and Cu form solid solutions over all concentrations at high temperatures but, for certain composition ranges, ordered intermetallic phases can be obtained at lower temperatures. In particular the known intermetallic compounds are the cubic AuCu3 and Au3Cu (space group Pm$\mathrm{bar}{3}$m), the low-temperature phase AuCu(I) (space group P4/mmm) and the high-temperature phase AuCu(II) (space group Imma). The phase diagram of Ag-Cu instead exhibits eutectic behaviour with a wide miscibility gap and the system tends to segregate in phases of nearly pure Ag and pure Cu at room temperature Okamoto1990 .

We study the effect of composition on the reflectivity of the Ag-Au system and compare experimental data of solid solutions with SQS simulations. Fig. 6 shows that the gradual shift to lower wavelengths of the reflectivity edge of gold by increasing the Ag content is reproduced by the simulations. However, as already discussed above for the case of elemental noble metals, the position of the reflectivity edge in IPA simulations based on PBE band structures does not correspond to the experimental one, but it is instead systematically shifted to longer wavelengths for each atomic concentration $x$ considered. Although the simulations are not in quantitative agreement with experiments, the qualitative trends in reflectivity, and thus in colour, with respect to the alloy composition of Ag-Au are reproduced. Similarly, we simulate the optical properties of Ag1-xCux two-phase alloys by employing the Bruggeman model described above and study also in this system the effect of composition on the reflectivity. The $\alpha$ and $\beta$ phases entering in the expression for the alloy dielectric function $\varepsilon_{\text{Br}}(\omega)$ of Eq. 10 are assumed to be elemental Ag and elemental Cu, respectively. And the dielectric functions of the two constituent elements are taken from the simulations of elemental metals discussed above. As shown in Fig. 7, Ag additions in Cu increase the reflectivity at wavelengths shorter than the reflectivity edge of elemental Cu but do not shift, as it happens in Ag-Au solid solutions, the position of the edge. The Bruggeman model provides the correct trend with composition but the effect on the drop in the reflectivity is less evident because the reflectivity edge of elemental Cu in IPA simulations is less steep than the experimental one. Note that the application of the Bruggeman model to experimental data of the dielectric function of elemental Ag and elemental Cu gives very good agreement with experimental data for the two-phase alloy and validates the use of the model.

Summarizing, for Ag-Au solid solutions, where there is a gradual shift of the reflectivity edge by varying alloying additions from elemental Au to elemental Ag, the colour of the alloy changes from red-yellow to yellow, pale greenish-yellow and eventually white of pure Ag. Au-Cu solid solutions show a similar behaviour Rivory1977 and the colour of the alloy changes from red-yellow to reddish and eventually red of pure Cu. Instead, in Ag-Cu two-phase alloys there is no shift of the reflectivity edge but, for all wavelengths in the visible range below the reflectivity edge of elemental Cu, the reflectivity curve rises roughly uniformly so that the colour of Ag-Cu changes from the red of pure Cu to reddish and then directly to whitish and white of pure Ag Cretu1999 .

After considering the effect of composition on the reflectivity of binary alloys, we now study the effect of different types of compounds for a given fixed atomic concentration $x$ directly on the dielectric function of the Au-Cu and Ag-Cu systems. Indeed, for Au-Cu at the composition $x=0.81$, experimental data are available in the literature for the optical absorption of both the solid solution and the intermetallic compound AuCu3 Rivory1977 . Analogously, for Ag-Cu at the composition $x=0.30$, experimental data are available for both a segregated two-phase sample made of a pure Cu phase and a pure Ag phase, and for a metastable solid solution obtained by vapor quenching Rivory1977 . We compare the optical absorption of the Au1-xCux solid solution, at $x=0.81$ in experiments and at $x=0.75$ in simulations, with the optical absorption of the intermetallic compound appearing around the composition $x=0.75$, i.e. the cubic AuCu3 phase. The purpose of this comparison is to study the differences in optical properties between ordered and disordered phases. As shown in Fig. 8, the optical absorption of the intermetallic compound is very similar to the one of the random alloy with the notable exception of the presence of an additional peak at around 3.6 eV, which is missing in $\varepsilon_{2}(\omega)$ for the solid solution. The comparison of the SQS results for the disordered alloy with the simulated results of the intermetallic compound shows that the simulations clearly capture this small difference. Nonetheless, we underline that there is no significant change in the resulting colour between ordered and disordered alloy for this system because the position of the onset of optical absorption is not modified by the presence of long-range order, and thus neither is the colour.

Similarly, Fig. 9 reports the comparison between the optical absorption of the Ag1-xCux two-phase alloy, at $x=0.70$ in experiments and at $x=0.75$ in simulations, with respect to that of the metastable solid solution having the same composition. In the two-phase alloy, where the alloy optical properties are well approximated by a combination of those of pure Cu and pure Ag (Bruggeman model), we observe two onsets of absorption: the first one at $\sim$ 2.1 eV corresponding to the absorption edge of pure Cu and the second one at $\sim$ 4.0 eV corresponding to the absorption edge of pure Ag. The optical absorption of the solid solution instead is very similar to the one of pure Ag but, in addition, we observe the presence of a supplementary broad peak at energies below the onset of absorption of pure Ag due to Cu impurity states. The SQS results for the solid solution and the results of the Bruggeman model applied on the IPA dielectric function of elemental Ag and Cu reproduce the two different trends, although the SQS shows a small blueshift of the peak that follows the absorption edge of pure Ag which is not observed experimentally.

III Discussion

We have shown that the theoretical methods and approximations considered in this paper, i.e. IPA optical spectra computed on top of the DFT-PBE electronic structure, can be employed in systematic studies on the optical properties of metals in order to predict trends in real metallic systems and to help the search for novel materials with specific optical properties, and therefore also colours, by exploring the composition space through the computational screening of materials Prandini2019 . Moreover, this work could help stimulate future studies aiming to achieve the photorealistic simulation of different types of materials by means of first-principles techniques. For example, the systematic validation of the approach performed on elemental metals and binary alloys can be seen as a necessary preliminary step for the photorealistic simulation of more complex metallic alloys having a larger number of constituent elements, such as ternaries, quaternaries, etc., which are more relevant for technological applications (e.g. superalloys and high-entropy alloys).

IV Methods

IV.1 Workflow

All DFT calculations are performed with the Quantum ESPRESSO distribution Giannozzi2009 , which is based on the plane-wave pseudopotential method for the numerical solution of the KS equations. We use Shirley’s interpolation method Shirley1996 ; Prendergast2009 as implemented in the SIMPLE code simple_cpc to evaluate the IPA dielectric function of metals including both interband and intraband contributions. Photorealistic rendering is performed with the Mitsuba renderer mitsuba_webpage . Pseudopotentials and plane-wave cutoffs are chosen according to the results of the standard solid-state pseudopotential (SSSP) protocol sssp_npj in order to have reliable and converged band structures as the starting ingredients for the evaluation of the IPA dielectric function. Since the SIMPLE code supports only norm-conserving pseudopotentials, we use optimized norm-conserving Vanderbilt (ONCV) Hamann2013 pseudopotentials from the SG15 Schlipf2015 and PseudoDojo Dojo2017 PBE pseudopotential libraries for all elements considered (see Supplementary Discussion 2 for more details on the choice of the pseudopotentials from the SSSP database of tests). For the purpose of automation, the sequence of calculations required by the computational approach described in this work is implemented as a workflow within the framework of the AiiDA Pizzi2016 infrastructure for computational science. Thanks to this ColourWorkflow (see Fig. 1), it is possible, giving as input a generic crystal structure, to obtain directly as output the reflectivity and colour of a given material.

In all simulations, relativistic effects are accounted for at the scalar-relativistic level (see Supplementary Discussion 3 for an analysis on the effect of spin-orbit coupling on the optical properties of heavy elements) while the IPA dielectric function is always evaluated by including the non-local contribution of the pseudopotentials in the computation of the velocity matrix elements, as implemented in SIMPLE.

IV.2 Elemental metals

All calculations on elemental metals are performed on the ground-state crystal structures at zero temperature, as provided in Ref. Lejaeghere2014 . The equilibrium volume of each structure corresponds to the reference PBE value obtained by extensively tested all-electron calculations for the equation of state Lejaeghere2016 . If needed, the crystal structures are reduced to the primitive cell using the spglib library spglib_webpage . Spin-polarization is not included in our calculations. In the self-consistent DFT calculations for the evaluation of the ground-state density we use a Monkhorst-Pack grid MonkhorstPack of $24\times 24\times 24$ and a cold smearing Marzari1999 of 0.02 Ry. In the non self-consistent band structure calculations needed for the construction of the Shirley’s basis we use a uniform k-grid of $2\times 2\times 2$ including the seven periodic images of the $\Gamma$-point of the BZ and at least 30 empty conduction bands. From a convergence study on the dielectric function we decide to employ an interpolation k-grid of $64\times 64\times 64$ and $\eta=\gamma=0.1$ eV in SIMPLE for each elemental metal considered, with the exception of elemental aluminium for which, because of a very slow convergence of $\varepsilon^{\text{inter}}(\omega)$ with respect to k-points sampling, the interpolation k-grid used is $80\times 80\times 80$ and $\eta$ is set to 0.2 eV. The Shirley’s basis is constructed setting the threshold for the Gram-Schmidt orthonormalization algorithm equal to 0.0075 a.u. (input variable named $s_{b}$ in SIMPLE).

IV.3 Alloys

For the simulation of all binary compounds considered, we always use as plane-wave cutoff the largest value between the plane-wave cutoffs of the two constituent elements, as taken from Supplementary Table 1. The Shirley’s basis is constructed setting $s_{b}=0.01$ a.u. in SIMPLE and considering a number of empty bands at least equal to the number of occupied bands. We choose the interpolation k-grid to be used in the evaluation of the dielectric function in terms of a k-point density, which is defined as the maximum distance between adjacent k-points along the reciprocal axes (in Å*-1*). For all the seven cubic intermetallic compounds considered we select as k-point density 0.04 Å*-1*. With this choice the number of k-points included in the uniform k-grids is of the order $O(10^{5})$, which correspond to uniform k-grids in the range from $46\times 46\times 46$ up to $56\times 56\times 56$.

All the SQSs used in this work to simulate solid solutions of the systems Ag1-xAux, Au1-xCux and Ag1-xCux are generated with the ATAT package vandewalle2002 ; vandewalle2013 . Since we consider only the simple stoichiometric ratios $x=0.25,0.5,0.75$, we use small FCC SQSs with 16 atoms per cell. The interpolation k-grid is set according to a k-point density of 0.04 Å*-1* (corresponding roughly to 11’000 points in the BZ).

Acknowledgements

The authors warmly thank Fanny Lalire and Frédéric Diologent for several useful discussions and for sharing with us confidential experimental results.

Competing Interests

The authors declare no competing interests.

Author contributions

N. M. and G.-M. R. designed the study; G. P. developed the computational workflow, performed the calculations and wrote the manuscript. All authors discussed and analysed the results and commented on the manuscript.

Funding

This research was supported by Varinor SA (CH 2800 Delémont, Switzerland).

Data availability

The data that support the findings of this study are available from the corresponding authors upon reasonable request. The source code of the ColourWorkflow and the input scripts necessary in order to reproduce the simulations performed for this work are available at https://github.com/giprandini/colour-workflow.

Additional information

Supplementary information is available at npj Computational Materials website.

V References

VI Figures

VII Tables