Photorealistic modelling of metals from first principles

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.


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 [1,2]), for optoelectronics devices (e.g. in ultra-thin films for transparent conductive electrodes [3,4]), and also for microscopy and optical data storage [5]. 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 [6]. 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. [7] computed the optical properties of several elemental metals whereas, more recently, Werner et al. [8] 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. [9] 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. [10] 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 AuAl 2 , AuGa 2 and AuIn 2 was first computed in Ref. [11] and, afterwards, Keast et al. [12] studied the influence of alloying concentrations on the reflectivity and colour of the intermetallic AuAl 2 by considering the ternary compounds having the Au 1−x Pt x Al 2 composition, with x = 0.0, 0.5, 0.75, 1.0; equivalent computational results were obtained independently by Kecik [13]. Calculated and experimental [14,15] 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 Au 0.5 Cu 0.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 [16], 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.

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.
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) [17] within the generalized gradient approximation (GGA) and relying on the PBE exchange-correlation functional [18]. 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 [19][20][21]), albeit at a largely increased computational cost.
So, while in the present work we just rely on the Kohn-Sham (KS) PBE bands [22], 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 electronhole 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 [23]. In the optical regime the momentum q transferred by the photon is negligible so that we can consider the optical limit, q → 0, of the expression for the IPA dielectric function ε(q, ω). In general, the dielectric function still depends on the direction q = q/|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 Drudelike term ε intra (q, ω) due to the conduction electrons at the Fermi surface and an interband term ε inter (q, ω) due to vertical transitions between occupied and unoccupied bands, so that ε(q, ω) = ε inter (q, ω) + ε intra (q, ω). Using the solutions of the one-particle Schrödinger equation for periodic systems, H KS |ψ nk = E nk |ψ nk (where H KS is the KS Hamiltonian from DFT), the explicit expression of the IPA dielectric function can be written as [23][24][25] where we have defined the IPA Drude plasma frequency as 3) can be found in Ref. [26].
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 so that we can drop the dependence on the directionq. Similarly we also define a corresponding average IPA Drude plasma frequency as ω 2 D = [ω 2 D (x) + ω 2 D (ŷ) + ω 2 D (ẑ)]/3. In order to compute the reflectivity from the knowledge of the dielectric function we first introduce the refractive index n(ω) and the extinction coefficient k(ω) that are defined from the equation [n(ω) + ik(ω)] 2 = ε(ω). The reflectivity at normal incidence and assuming a vacuum-material interface is then simply linked to n(ω) and k(ω) through the Fresnel equations of classical electromagnetism (see for example Ref. [27]): 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) [28] 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(λ) for all the wavelengths λ 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 [29]. In particular, according to the CIE 1931 standard colorimetric observer, the tristimulus values (X, Y , Z) which define the CIE-XY Z colour space completely describe a colour stimulus and are given by the following integrals over the visible range wherex(λ),ȳ(λ) andz(λ) 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(λ) 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(λ) = 1 for all visible wavelengths.
In practice, it is more convenient to work within the CIELAB colour space rather than in the CIE-XY Z 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 [29]. 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 are the CIELAB coordinates of two objects, their colour difference is simply given by Typically, a difference ∆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 [30], 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 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) [31,32], 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 α and β, we use the Bruggeman model [33] to estimate the optical properties of the alloy. Within the Bruggeman model the dielectric function of the mixture, that we indicate as ε Br (ω), is given by the following expression in terms of the dielectric functions ε α (ω) and ε β (ω) of the single phases, and where x α and x β (with x α + x β = 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.
Elemental metals 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 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 [34] and Ag [35] but not for Au, for which G 0 W 0 gives very similar results to PBE [36]. For this latter case, the quasi-particle self-consistent GW (QSGW ) [37,38] approach is required for the occupied 5d bands of gold to be lowered in energy by the right amount [36]. 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 [25] 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 [34,35] (i.e. they do not act as a simple scissor operator). Nonetheless, these approximated simu- (ω)) in PBE is also not at the correct position (similar to the case of semiconductors for which the PBE band gap is systematically underestimated [20]). On the other hand, the shape of ε(ω) for noble metals is reasonably well reproduced.

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).
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 [11,12] 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 AuAl 2 , bluish AuGa 2 and yellow PtAl 2 (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 [39]. 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 ∼ 5 eV below the Fermi level and these do not contribute to the characteristic colours of these compounds [40]. 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 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 [41].
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 ele-

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 [26]. 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).

Workflow
All DFT calculations are performed with the Quantum ESPRESSO distribution [43], which is based on the plane-wave pseudopotential method for the numerical solution of the KS equations. We use Shirley's interpolation method [44,45] 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

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

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. Figure 1: Schematic representation of the workflow, named ColourWorkflow, designed to simulate the reflectivity and colour of a metallic material giving as input its crystal structure.   Experimental data are taken from Ref. [14] for AuAl 2 , AuGa 2 and AuIn 2 , from Ref. [12] for PtAl 2 , from Ref. [61] for NiSi 2 and from Ref. [39] for CoSi 2 and PdIn.

FIGURES
The two vertical dashed lines show the limits of the visible range. Figure 5: Comparison between the simulated rendering of a metallic surface (left panel) and real samples (right panel) of the intermetallic compounds AuAl 2 (top), AuGa 2 (center) and PtAl 2 (bottom). Images of the AuAl 2 and AuGa 2 samples are adapted from Ref. [62] while image of the PtAl 2 sample is adapted from Ref. [15].   [64] (right panel). The results of the simulations are obtained using the Bruggeman model in which the two phases of the system are assumed to be elemental Ag and elemental Cu (dashed lines).    Table 2: Computed values of the IPA Drude plasma frequency ω D (in eV) compared to previous simulations (J. Harl [25]) and experiments (Exp.). The experimental values with no explicit reference are extracted from the data reported in Ref. [25]. 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.  Table 3: Colour differences in CIELAB space between simulated colours (present work) and experimental colours [12,14,39,61] derived from reflectivity data, ∆E exp , and between simulated colours (present work) and previously published simulations [11,12], ∆E sim .