Understanding the Structural and Electronic Properties of Photoactive Tungsten Oxide Nanoparticles from Density Functional Theory and GW Approaches

ABSTRACT


Introduction
Heterogenous photocatalytic water splitting represents one of the most attractive technologies, since only photons with the appropriate energy are required to drive the solar-to-chemical energy conversion. 1 The first step of this process is the light absorption by the photoactive semiconductor to form electron-hole pairs, which have to fast and efficiently migrate toward the surface of the catalyst, where water molecules are reduced/oxidized by the electrons/holes to form H 2 /O 2 .A visible-light-driven photocatalyst must exhibit an energy gap below 3.0 eV.The design of such photocatalyst is fundamental to guarantee an optimal visible-light absorption.In principle, this requires a deep knowledge of structure-property relationships, which allow one monitoring the optical and electronic properties by engineering of size and morphology of the nanoparticles (NPs).As a matter of fact, the design of NPs with diameters below 5 nm via facet engineering 2 is one of the most promising avenues to improve the photocatalytic activity.At this size the effects of the quantum confinement facilitate the transfer of charge carriers. 3,4This implies that the size and morphology of NPs are two important factors in determining their electronic structure and, hence, their potential use as photocatalysts. 57][8] In this context, setting up a reliable but still affordable methodology capable to deliver accurate information on the structure, the stability and the band-edge offsets of realistic nanocrystals is crucial for the design of efficient NPs for photocatalytic applications. 9A proper modelling of NPs usually implies to tackle two main issues: (i) the construction of suitable structural models at the atomistic scale that are able to capture both the effects of morphology and size, which is often done by following top-down or bottom-up approaches; and (ii) the setup of a high precise and unbiased approach which can be employed to determine the electronic properties of NPs at the nanoscale level within an affordable computational cost.TiO 2 NPs constitute by far the most studied photocatalytic material, 10,11 although the need for UV radiation makes it inappropriate for large scale water splitting applications.In this context, WO 3 derived nanostructures have emerged as promising photocatalytic materials due to their Earth abundancy, low cost, harmlessness, suitable band edge energies for water splitting reactions, and high chemical stability, even in the presence of strong oxidizing conditions, making them relevant for many technological applications. 12More importantly, the experimental band gap energy (E g ) of bulk WO 3 is in the 2.5-3.0 eV range, 13,14 sensibly smaller than the corresponding value for bulk TiO 2 which is in the 3.0-3.2eV range. 15WO 3 can, thus, act as visible light-driven photocatalyst.However, a systematic and fundamental knowledge of its photoactivity as a function of morphology and size is still lacking.1][22] There are also studies regarding very small (WO 3 ) 3 nanoclusters (NCs) 23,24 but till the date, the modelling of realistic WO 3 NPs is still missing.To fill this gap, we have followed the successfully computational strategy set up by some of us for TiO 2 [25][26][27] and ZnO 28 nanostructures, which permitted us to design different morphologies for a given composition or to vary the composition for a certain given morphology.
Concerning the choice of the proper level of theory, GW formalism (G represents the timeordered one-b y G ' f wh W depicts the screened Coulomb potential) has attracted the attention of both Physics and Chemistry communities during the last years due to its excellent performance in terms of accuracy and computational cost, which is inherent of its N 3 scalability (being N the number of atoms of the system under study). 29Notably, this perturbative method confirmed as appropriate method to study the electronic properties of medium size TiO 2 [30][31][32] and small transition metal nanostructures. 33Nonetheless, due to the large number of atoms in our realistic WO 3 nanostructures (>1500), for the largest systems only standard DFT methods (i.e.Local Density Approximation (LDA) or Generalized Gradient Approximation (GGA)) are applicable, which clearly translates into a larger dependence on the choice of the exchange-correlation potential with a possible significant loss of accuracy.Here we overcome such limitation by establishing a combined DFT-GW protocol, relying on the correlation between the DFT and the GW results within a representative subset of nanoparticles, spanning all the examined morphologies.This enables us to predict the properties of a complete set of WO 3 nanostructures from few cluster units to realistic NPs (up to 1680 atoms) at a higher level of accuracy.
In this work, we discuss the merits and limitations of our combined GW-DFT method and its implications in the design of WO 3 NPs for photocatalysis.We first address the evolution of the relative stability of the WO 3 NPs with respect to their size by considering different structural families and the most likely morphology at the lowest dimension scale.Next, we benchmark the electronic structure properties of these particles on the basis of the total variational energy ff (ΔS to evaluate the suitability of our GW based methodology.Finally, we use our DFT-GW approach to assess the evolution of the NPs energetic levels with respect to their size with the ultimate goal of pointing out the optimal combinations of size and morphology to enhance the photoactivity of the material.

Definition of the structural models
The present WO 3 NPs are all stochiometric and have been designed by using the top-down approach. 34Our study involves NPs sizing from ⁓0.4 ⁓5 w g h z g of the smallest synthetized NPs found in the literature (0.7-5 nm), [35][36][37][38][39] exhibiting cubic, rectangular, octahedral, spherical, hexagonal nanosheet and cuboid nanowires morphologies; as depicted in Figure 1.(002) surfaces. 40,41More in detail, the rectangular NPs are built following the structural features previously reported, with the (002) surface being the preferential one with respect to the (200) and (020) surface. 41Octahedral WO 3 NPs are designed by cutting the m-WO 3 along the (111) directions. 42For the particular case of the spherical NPs, these nanostructures are obtained by imposing symmetrical cuts (non-symmetric nanostructures were found to be unstable upon relaxation) to the monoclinic bulk phase along several directions, thus resulting in a polyhedral ' h -k ' h .R g g h h x g h h (001 f conserved according to previously reported works. 43,44On the other hand, the hexagonal nanowires expose the (110) surfaces as found in experiments, thus resulting in a cuboid prismatic-shaped structure. 43,44 the present study we also include the exposed (100) and (010) surfaces with the aim of getting stochiometric and symmetric nanowires with a section involving the lowest possible area, thus leading to a prismatic nanostructure with a hexagonal moiety.Note that cuboid prismatic 1D nanostructures with the exposed (110) surface have also been tested but all were found to be unstable upon geometry optimization.Once the nanostructures exposing the desired surfaces have been built, we systemically modify the number oxygen surface atoms, keeping the same coordination for all superficial chemical units, to obtain a quasi-stochiometric symmetric nanostructures.Finally, once this symmetric structure is relaxed, we symmetrically add or remove certain oxygen surface atoms to preserve the stoichiometry of the system, prior to the final relaxation of the nanostructure.It is noteworthy to mention that we have selected those relaxed nanostructures which exhibited the larger stability and non-metallic behaviour among different candidates, thus reducing the randomness in the construction of our models.In addition, we are aware that the surface geometries of the pristine nanostructures may be slightly modified by the presence of the solvent molecules around the NP surface with a fraction being probably dissociated.Nonetheless, the investigation of the interaction between the solvent and the NPs is beyond the scope of the present paper, and it will be the motivation for prospective works.
Table 1.Number of WO 3 units (n) and geometrical parameters for all NPs studied in this work, as they are depicted in the schemes on the left side of the Table .The R factor corresponds to the ratio between the largest and the lowest dimensions of the NP.For each one of the structures shown in Figure 1, geometry optimization was carried out in the framework of Density Functional Theory (DFT) using the Perdew-Burke-Ernzerhof (PBE) functional, 45 including the Periodic Boundary Conditions (PBC) and employing a combination of Gaussian and plane wave basis set as implemented in CP2K package. 46To avoid interactions between periodically repeated images of the NPs, the size of the unit cell was set to assure a minimum vacuum distance of 15 Å between neighbouring nanostructures.The calculations were carried out with a Double-zeta Valence Polarized (DZVP) and an energy cut-off for the auxiliary plane wave basis equal to 500 Ry, whereas the core electrons were described by means of Goedecker-Teter-Hutter (GTH) pseudopotentials. 47

Electronic structure calculations
The electronic properties of optimized structures were studied using an all-electron approach and employing a numerical atom-centered (NAO) orbital basis set.These calculations were carried out using the Fritz Haber Institute ab initio molecular simulations (FHI-AIMS) program package. 48A light grid and tier-1 basis set is used, which lead to a similar quality to a valence . 27 Due to the presence of heavy atoms such as W (Z=74), relativistic effects were accounted for through the zeroth order regular approximation (ZORA) 49,50 to the Dirac equation, which has been shown to be accurate yet computationally affordable.We noted that relativistic effects were especially relevant for medium size NPs (see Figs. S1 and S2).For the sake of coherence with the level of theory employed for the geometry relaxation, electronic structure calculations were also carried out by using the GGA PBE and the hybrid PBE0 51 at the PBE structure.
The lack of available experimental data to benchmark the electronic properties of small WO 3 NPs makes compulsory the adoption of highly accurate computational methods, which still allow to handle systems with hundreds of atoms, as the ones considered here.For that reason, semiempirical approaches such as PBE+U have not been considered here in view of the absence of an experimental reference to properly set the values U parameter, which is responsible of preventing the excessive delocalization of the 5d W orbitals in the NP.In this context, the electron affinity (EA) and ionization potentials (IP) energies 52 of the optimized NPs have been estimated by following three different approximation of increasing level of theory: i) the Kohn-Sham (KS) DFT energies; 53,54 ii) the quasi-particle (QP) energies as obtained from many body perturbation theory (MBPT) based GW techniques, 55,56 and, finally, iii h ΔS h wh h h best possible estimated for a given density functional as it involves the difference of two variational energies.The GW technique is especially well suited for extended systems where ΔSCF become cumbersome due to the presence of charged states.
In the case of the KS energies, the IP/EA energies were taken directly from the values of the KS Highest/Lowest Occupied/Unoccupied Molecular Orbitals (HOMO/LUMO) energy levels with respect the vacuum level (VL).It is important to stress, however, that the equivalence between the fundamental gap (IP-EA) and the HOMO-LUMO KS gap strictly holds for the exact exchange and correlation functional 57 and that, in general, the frontier orbital E g is significantly different, often smaller, than the experimental fundamental gap. 58To predict more reliable E g values for semiconductors, the most employed approach relies on the use of hybrid density functionals.These include a (moderate) fraction of Hartree-Fock exchange in the GGA or LDA standard exchange correlation (xc) density functional.By using hybrids functionals, one can correct the sizeable (30-40%) underestimation of E g with respect to experiment 59,60 by standard LDA or GGA functionals. 61,62However, E g appears to be extremely sensitive to the amount of Hartree-Fock exchange employed, and, in practice, one should always benchmark the calculated quantities against experimental data.For instance, for the monoclinic bulk WO 3 , the experimental gap is reproduced by including a 11% of non-local Hartree-Fock exchange (see Figure S3), which goes in line with the percentage reported in previous studies (15%), 17 meaning that PBE0 (25%) is significantly overestimating this quantity by a ⁓1.1 V.
The second analysis consists of the GW method, which allows us accessing to the quasiparticle (QP) energies by perturbating the KS energy levels.In particular, two different approaches were used here: the one-shot GW approximation (G 0 W 0 ) above mentioned and the more accurate eigenvalue partially self-consistent GW (evGW), where the self-gy (∑ w h partial self-consistency in the eigenvalues, while keeping the KS MO coefficients unchanged.
The so-called G 0 W 0 approach has been applied to a large dataset of molecular systems 63 with excellent agreement to experiments and also successfully applied to TiO 2 NPs. 31,32The GW calculations provide QP energies related to physical quantities, such as the electron injection/extraction energies, that can be directly compared with photoemission spectroscopy data. 64Furthermore, GW techniques beyond G 0 W 0 also present the advantage of providing QP energies which are less dependent on the amount of non-local exchange fraction in the starting DFT functional with respect to the KS energies. 65For both GW approaches, we relied on a two pole fitting type for the analytical continuation of the self-energy, as it is implemented in the routines of FHI-AIMS. 66h ΔS h h P EA g w by employing the following formulas: where E(N), E(N-1) and E(N+1) correspond to the total energies for the neutral, +1 and -1 charged systems calculated at their equilibrium ground states (neutral) geometries.As it was the case of the one-shot G 0 W 0 approximation, ΔS h may exhibit an important dependence of the initial KS energy orbitals, which depend on the density functional employed.

Assessing the structural properties
The first property that we address is the energetic stability of WO 3 NPs as a function of their size, as summarized in Figure 2. Clearly, the octahedral NPs are less stable when compared to the other morphologies, thus confirming the experimentally assessed metastability of such NP. 42cordingly, we have decided to discard this family of NPs in the following analysis.Regarding the rest of NPs belonging to the monoclinic bulk phase, a crossover around 50 units of WO 3 is found implying that, below this number, rectangular WO 3 NPs emerge as the most stable ones, whereas cubic WO 3 NPs (a⁓1 b 1 g y f f n > 50.
Interestingly, cubic NPs reached an asymptotic behaviour in terms of energy for a relatively small number of WO 3 units (n⁓250 g h f h y f h gy z NP (a ⁓2.1 b y b w h y g h gh gy ll with respect to the bulk phase.Above the n = 50 crossover, the range of energies for the spherical NPs are found in the middle between the cubic and rectangular NP energies, while below this, the stability of these three families (cubic, rectangular, spherical) was found very similar.One can, thus, infers that the effects of the morphology in the NP stability are especially relevant only for medium and large size NPs.Moving to the derived hexagonal bulk phase NPs, while nanosheet NPs followed similar trends with respect to the monoclinic NPs, resulting in similar values for the asymptotic energies, nanowire NPs energies converged faster in terms of number of units (n ⁓125 .h f y however, about 0.15 eV per unit less stable compared to the rest of NPs, as one may expect due to their 1D morphology, implying larger surface effects.In particular, we predicted that the stability of the nanowires as a function of their lengths already converges for ratios R around 11.
Notably, we found a second crossover at n = 165 (which corresponds to R ratios of 9.8 and 14.2 for the nanosheet and nanowire WO 3 NPs, respectively).Nanowire and nanosheet structures are stable below and above such crossover.

Benchmarking the electronic structure methods: the effect of Hartree-Fock exchange
For the sake of illustration, we have taken the cubic (WO 3 ) 8 NP as reference system and we have estimated the IP and EA by using different level of approximations as the standard PBE and the hybrid PBE0 and PBEh50 density functionals, including 0%, 25% and 50% of Hartree-Fock exchange, respectively; the data are displayed in Figure 3.These data still show a certain degree of dependence, larger at ΔSCF level, on the fraction of non-local exchange included in the DFT functional.Finally, the values obtained from the partial self-consistent GW method are 9.2 (evGW@PBE), 9.1 (evGW@PBE0), and 8.9 eV (evGW@PBEh50), interestingly indicating that evGW exhibits a smaller dependence on the employed density functional than the G 0 W 0 results.This does not necessarily imply, however, a more accurate description; especially in view of the discrepancies with ΔSCF method and the lack of experimental data to contrast these results.For instance, b et al. 33 found that the G 0 W 0 method provided a description of the experimental IP/EAs in small transition metal molecules, better than that arising from other self-consistent GW techniques.In any case, the matching between evGW@P Eh50 ΔS @P Eh50 k b b b b y misleading since such a large amount of non-local exchange will largely overestimate the band gap of bulk WO 3 (⁓3.2V e Figure S3).It is also worth noting that similar E g trends are reported for the hexagonal (WO 3 ) 12 and spherical (WO 3 ) 19 NPs as shown in Figures S4 and S5 reported in the Supporting Information (SI).

Establishing the DFT-GW protocol on a representative NPs subset
Unfortunately, the computational burden for calculations with hybrid functionals and GW methods makes the study of realistic WO 3 NPs unaffordable.As a matter of fact, the electronic structure of large NPs composed by thousands of atoms can be figured out only at the standard DFT level.The idea that we develop in this section is, thus, trying to estimate the electronic properties (E g , EA, IP) at a higher level of theory (GW and/or ΔSCF) by performing affordable GGA calculations.smallest NPs depicted in Figure 1, which have been calculated at the G 0 W 0 @PBE vs KS@PBE.
To this end, we investigated the correlation between the three electronic quantities of interest obtained from the PBE KS orbital energies with the corresponding G 0 W 0 @PBE values (Figure 4).This seems a logical choice since GW allows one to go from the orbital energies to the QPs which have a well-defined physical meaning.The comparison, presented in Figure 4, is carried out for the subset of eleven smaller NPs which covers all possible morphologies: two cubic, one rectangular, two spherical, two nanosheet and four nanowire NPs.This strategy allows to discard possible deviations arising from the use of a given morphology.Moreover, it also permits to compare the properties of the largest NPs with the experimental band edge energies of the monoclinic bulk phase. 68We rely on the one shot G 0 W 0 approximation that has been shown to provide excellent agreement with the experimental data for a large number of tested molecules 63 and is still accessible for the size of the NPs composing the selected subset.These correlations showed a nearly perfect linear dependence between the EA and E g values estimated by both levels of theory, whereas IP quantities were slightly diverging from the linear behaviour.Similar trends were found in the correlations between ΔSCF@PBE and KS@PBE (see Figure S7).

Predicting the electronic properties over the full set on NPs
Figure 5 presents the estimated DFT-GW values of the IP, EA and E g for all NPs studied in this work, as obtained by applying the grey line correlations of Figure 4 to the quantities calculated at PBE level.The corresponding quantities estimated at ΔSCF@PBE0 and ΔSCF@PBE are reported in Figure S9.Note that the systematic deviation of IPs, EAs and on their corresponding gaps predicted by GW ΔS h y ( gure S8).Thus, the agreement in the trends followed by the IP/EA energies with the NP size estimated from either G 0 W 0 @P E ΔS @P E ΔS @P E0 k b ( gure S9).Therefore, even if any of the three methods can be used to investigate trends with size and since PBE0 clearly overestimates the E g in WO 3 , we find it convenient to henceforth focus on the G 0 W 0 @PBE estimates only.Note that in the following discussion we will make use of the terms IP/EA and conduction/valence band (CB/VB) edges to depict the frontier energy levels of the NPs and bulk materials, respectively.Interestingly, the plots in Figure 5 indicate that the IP of the WO 3 NPs do not so strongly vary with the size or the morphology of the NP as compared with their EAs, lying in a range 4-5.5 V vs SHE.Indeed, both estimates are in a rather fair agreement and appear to significantly overestimate the experimental VB edge (i.e., 3.1 V vs SHE). 68On the other hand, the EAs exhibit a significant drop in energy when going below the 100 units and, what is more interesting, they show a large dependence on the NP morphology.While cubic, rectangular and nanowire NPs display an identical behaviour with respect to the NPs size (leading to converged EA values of 0.9 V vs SHE), nanosheet NPs present EA values lower than those of the rest of NP families (with converged energies about -0.7 V vs SHE).Comparing these quantities to the monoclinic WO 3 CB energy of 0.5 V vs SHE 68 and the H 2 /H 2 O reduction potential (0V vs SHE), it appears that the electron injection driving forces with respect to this potential are significantly reduced when decreasing the NP size, with the exception of the nanosheet NPs.Remarkably, the stronger impact of the quantum confinement effects in the EA with respect the IP energies has been observed also experimentally observed, 38 where a change of 1.02 eV in the CB edge was found when moving from bulk phase to quantum dots, whereas this shift only amounted by 0.57 eV for the VB.To get a deeper understanding of the different behaviour of IP and EA energies with the NP size and morphology, we have plotted, in Figure 6, the frontier MOs (HOMO/LUMO) of a representative medium size NP of each family owing a comparable number of chemical units composed by 500-600 atoms.For all the families the HOMO KS orbital is highly localized over a few full coordinated oxygen atoms, essentially dominated by the O(2p) orbitals.As a result, the dependence of the HOMO localization with the NP sizes is almost negligible, as shown by the similar IP energies calculated for all series of NPs.On the contrary, LUMO KS orbitals are fully delocalized on W atoms and mostly dominated by the 5d orbitals of W atoms located at the NP surfaces, thus implying a stronger impact of the size and morphology of the NP.Regarding the evolution of E g , which seems to stabilise at around 150 WO 3 units, all families of NPs present a significant reduction of the gap as the NP size increases.This gap decreasing ranked from 7.5 eV for the smallest NPs to converged gap values for the monoclinic NPs about 2.7-3.5 eV, which are in excellent agreement with the measured (2.62 eV) 69 and calculated gaps at G 0 W 0 @LDA (2.92 eV) 70 and at the PBE0 (3.72 eV, see Figure S3) levels of theory for the monoclinic WO 3 bulk.
To sum up, we can set to 150 the critical number of chemical units where the quantum confinement effects start to be important for WO 3 NPs.Nevertheless, a significant increase in E g values can be detrimental for the photoactivity of the material, since small NPs featuring an E g larger than 3 eV will absorb light in the UV region, thus covering a narrow portion of the solar spectrum.In this viewpoint, by taking as a reference the experimental value of the E g (2.62 eV), we can state that those NPs which exhibit a E g lower than 0.4 eV with respect to their asymptotic value will still work in the VIS region.By looking at the data plotted in Figure 5, we can conclude that all NPs which possess a number of chemical units larger than ca. 100 to satisfy this condition ((E g ) n -(E g )  < 0.4 eV).On overall, these results indicate that the NPs presenting a number of chemical units lying in the interval of 100 < n < 150 are likely to display an excellent photocatalytic behavior.This corresponds to sizes of a= 1.

Conclusions
In this work we have established a solid theoretical protocol, based on combining DFT and stateof-the-art GW calculations to investigate the intertwined role played by the size and morphology of realistic WO 3 NPs on the relevant opto-electronic properties which drive its photocatalytic behaviour.The first addressed issue was the relative stability for the different NP morphologies as the size of the NP increases.We found that cubic/nanosheet NPs are energetically favoured for medium and large size monoclinic/hexagonal phase NPs, while, for small NPs (n < 50/125), rectangular/nanowire NPs displayed a larger stability.In all the cases the NPs energies reached an asymptotic dependency with their size for a few hundred of units.Then, we investigated their optoelectronic properties, showing that the proper description of the electron injection/extraction energies can be achieved only if the QP correction to the KS orbital energies is applied through a GW approach or by invoking the ΔSCF method, although the choice of the exchange-correlation remains an open issue.Nevertheless, for a representative subset of NPs a linear correlation emerges between KS orbital energies and either G 0 W 0 @PBE or ΔSCF@PBE0 calculated values of IP and EA indicating that the error in the KS orbital energies is systematic.Next, by making use of these correlations we were able to obtain EA/IP energies of G 0 W 0 @PBE quality for all the set of NPs studied in this work.The high localization of the highest occupied levels on few oxygen atoms of the NP was found to be at the origin of the almost independent behaviour of the IP energies with respect the size and morphology of the NP.On the other hand, a large delocalization over the whole NP of the LUMO on the 5d W orbitals induced a large dependence of EAs upon size and morphology, with a sizeable drop in magnitudes when moving to small NP sizes (n < 100).The present results indicate that the photocatalytic efficiency of WO 3 nanostructures for water oxidation can be substantially improved by quantum confinement effects when moving towards nanocrystals below about 150 units of WO 3 and choosing the proper morphology.This is the first time that the structure and electronic properties of a complete set of realisticsize WO 3 NPs, spanning all phases and morphologies experimentally reported, are systematically investigated by high level DFT and large-scale GW calculations.We argue that these results may serve as a basis to encourage the experimental community working in photocatalysis to dedicate more research efforts in the synthesis and characterization of small-sized WO 3 NPs which could benefit from the quantum confinement effects.From a theoretical and computational viewpoint, the WO 3 nanocluster models developed in this work can be used to study charge generation and dynamics within the NPs, as well as investigate the mechanism of the water oxidation reaction, addressing finite size and surface effects; whereas our DFT-GW based method can be extended to the study of other relevant materials for photocatalysis, such as ZnO or SnO 2 NPs.
ASSOCIATED CONTENT

Figure 1 .
Figure 1.The full set of WO 3 NPs with different sizes and morphologies investigated here.

Figure 2 .
Figure 2. Ground state energies per chemical unit as function of the number of WO 3 units for the

Figure 3 .
Figure 3. Difference between IP and EA for the monoclinic cubic (WO 3 ) 8 NPs as estimated by

Figure 4 .
Figure 4. Correlations between the IPs, EAs and E g (from the left to the right) for the eleven

Figure 5 .
Figure 5. Ground EA, IP, E g energies as function of their number of chemical units for the

Figure 6 .
Figure 6.Isodensity plots of the HOMO (yellow) and LUMO (green) for the (WO 3 ) 125 cubic, GW protocol, clearly show the urgent need to have fine control over the size of the synthetized NPs to enhance the photoactivity of the material and open the door to exploit the smallest NPs for wide gap absorption applications.