Daniel Mourad
email: daniel.mourad@chemie.uni-hamburg.de

Antoine Guille

Tangi Aubert

Edouard Brainis

Zeger Hens
email: zeger.hens@ugent.be

Random-Alloying Induced Signatures in the Absorption Spectra of Colloidal Quantum Dots

The optical spectra of homogeneously alloyed colloidal CdSe x S 1-x quantum dots (QDs) for sulfur-rich compositions show additional transitions between the first and second absorption features of pure CdSe and CdS. Using tight-binding calculations for a large (N = 50) number of stochastic realizations for each composition of the alloyed QDs, we trace back this effect to transitions from deeper valence-band states to the lowest conduction-band level. These transitions are energetically degenerate with the main absorption lines in the pure systems, but are shifted in-between the lowest two main transitions in alloyed QDs. Moreover, the configurational disorder upon alloying distorts the symmetry properties of electron and hole wave functions,

which results in a broad band of allowed transitions from these deeper valence-band states to the lowest conduction-band level. In agreement with experimental findings, the corresponding absorption feature shows a blueshift relative to the first exciton transition with decreasing sulfur to selenium ratio. Many-body calculations using the configuration interaction scheme reproduce the effect and additionally yield a bowing parameter of the optical gap in very good agreement with experimental findings from UV-Vis spectra.

Tuning materials properties to suit application needs has always been a general aim of materials science. In the case of bulk semiconductors for example, alloying, i. e., the formation of solid solutions of two or more compounds, has been widely used to adjust the band gap and concomitantly tune the color emitted by semiconductor LEDs and lasers. Opposed to this, quantum confinement makes electronic energy levels of semiconductor nanocrystals size and shape dependent. As compared to bulk semiconductors, the quasi-continuous energy bands reduce to a set of discrete, atomic-like energy levels with an energy gap between the highest occupied and the lowest unoccupied state that increases when reducing the nanocrystal dimensions. 1 As a result, opto-electronic properties can be tuned by the dimensions of the crystallite and quantum wells, rods and dots -where size-quantization is active in one, two or three dimensions, respectively -are nowadays widely studied for, e.g., lighting and display applications, lasing, photodetection and solar energy conversion. 2 Both for self-assembled and colloidal quantum dots (QDs) made using molecular beam epitaxy, metal-organic vapor phase epitaxy or wet chemical techniques, respectively, alloying has been used as an additional parameter along with size to tune opto-electronic properties.

In the case of colloidal QDs, this involves ternary and quaternary solid solutions of, e. g. , zinc and cadmium chalcogenides, 3,4 lead chalcogenides, [5][6][7] lead and tin chalcogenides 8 or indium pnictides. 9 In particular, (Cd,Zn)(Se,S) alloyed quantum dots have shown a wide tunability range of their emission over the entire visible spectrum. 10 An important asset of alloyed QDs is that their properties are not simply the arithmetic mean of the pure compounds. Band gaps show bowing as a function of composition, 4,7,11,12 even showing minimum or maximum values, and immiscible bulk compounds may show solid solution nanocrystals such as (Pb,Sn)S. 8 In addition, the random positioning of two types of atoms on the crystal lattice positions breaks the nanocrystal symmetry. Theoretical studies predict that this randomness should lead to a relaxation of the selection rules of interband electronic transitions and to a strong heterogeneity in the fine structure splitting. 13 Although this could provide additional ways to tune the opto-electronic properties of colloidal QDs by alloying, next to the continuous tuning of the band-gap transition, experimental results corroborating these theoretical predictions are still missing.

In the present work, we study the influence of compositional disorder on the interband transitions in homogeneously alloyed CdSe x S 1-x QDs. We demonstrate that addition of selenium to CdS QDs leads to an additional feature in the absorbance spectrum, that shifts to higher energy relative to the fundamental band-gap transition with increasing selenium content. To interpret this feature, we make use of tight-binding (TB) calculations, where for each composition the energy level spectrum and the oscillator strength of interband transitions is determined for a large (N = 50) number of stochastic realizations of the alloyed QDs. We first demonstrate that these calculations yield energy level spectra for CdSe and CdS in line with existing literature and that they reproduce the experimentally observed bowing of the CdSe x S 1-x optical gap (which is an excited-state property) when a many-body approach is used. Next, we show that the additional absorption feature is due to an increased oscillator strength for a band of transitions that brings electrons from deeper valence-band states to the first conduction-band level. In line with the experimental findings, this band shifts deeper in the valence band with increasing selenium content. This result shows that the computational approach used here to analyze randomly alloyed colloidal QDs yields results in line with experimental findings. Hence, its predictions could be used to investigate more deeply the influence of random alloying on the opto-electronic properties of quantum dots and to obtain easily verifiable fingerprints of random alloying -which is often an experimental challenge -in the QD absorption spectrum.

Background: Tight-Binding Modeling of Alloyed Quantum Dots

Due to the large number of different microstates, the exact simulation of substitutional disorder effects in the framework of a multiband TB model is numerically very expensive, especially when many-body effects come into play. But even on a single-particle level, the Hamiltonian matrix dimension scales directly with the number of employed bands n, resulting in a memory requirement of the order n 2 and a computation time scaling of the order n 3 (see, e.g., Reference 14). Therefore, we employ a TB parametrization that allows couplings up to second nearest neighbor unit cells in real space and hence a reliable reproduction of the dispersion throughout the whole Brillouin zone (BZ) (due to higher non-vanishing Fourier coefficients in k-space) within a small basis set of effective Wannier orbitals. 15 In the case of zincblende semiconductors, such as the alloyed CdSe x S 1-x QDs studied here, the optical properties are dominated by contributions from an s like conduction band (CB) edge and three p like valence bands (VBs) per spin direction, namely the heavy hole (HH), light hole (LH) and split-off (SO) band. Consequently, we use eight Wannier orbitals |Rα as localized TB basis states, where α ∈ {s, p x , p y , p z } × {↑, ↓} and R ranges over the unit cells of the underlying fcc lattice (see Figure 1a). [START_REF]This approach is commonly called effective bond-orbital model (EBOM) in the literature and sometimes notationally differentiated from ETBMs that resolve the atomic basis sites. As the alternative usage of atomic orbitals is only a means to an end, namely[END_REF] Further details on the band structure parametrization (see Figure 1b), including the material parameters, and the modeling procedure for the stochastical QD realizations (see Figure 1c) can be found in the Methods section.

To emulate the experimental results, the QD diameter has been set to d ≈ 5 a, where a is the conventional lattice constant. This corresponds to approximately 3 nm, with a slight dependence on composition due to Vegard's law. Diagonalization of the TB Hamiltonian then yields the energies and expansion coefficients c i Rα of the eigenvectors, such that

|ψ i = Rα c i Rα |Rα .
To cover the relevant energy/wavelength range between the two main Figure 1: Calculation scheme for stochastic exact-disorder analysis.

absorption features of the experimentally recorded spectrum, we calculated the energies e i , h i and corresponding eigenstates |ψ i for 4 electron and 40 hole states per spin direction for each configuration and composition. 17

Finally, to simulate the excitonic spectrum and thus the optical properties of alloyed CdSe x S 1-x QDs, configuration interaction calculations can be used (see Figure 1d). Here, the eigenstates |Ψ of the interacting many-body system are expanded in a basis of Slater determinants constructed from the single-particle TB eigenstates of each individual realization of an alloyed QD. It should be noted that the meaning of "configuration interaction" varies between different communities and even authors. In our case, where very large systems are treated, we use an effectively screened Coulomb interaction as common in the solid state physics community (see Supporting Information and, e.g., References 18, 19) and expand the many-body wave function in a small basis set. This differs from the common quantum chemistry use of the term, where the bare Coulomb interaction is used alongside a much larger basis set of single-particle wave functions. From this expansion, the optical transition probability of the stochastic ensemble for converting a photon into an electron-hole pair is obtained from the initial (final) many-particle state of the n-th microstate |Ψ n i (|Ψ n f ) with energy E n i (E n f ) using Fermi's golden rule:

I(ω) = 2π 1 N N n=1 i,f | Ψ n f |H n D |Ψ n i | 2 δ(E n i -E n f -ω). (1)
Here, is the reduced Planck constant, ω the absorption frequency and H D is the light-matter coupling Hamiltonian in the dipole approximation, of which the one-particle representation is constituted by the MEs of the dipole operator:

d eh ij = e 0 ψ e i |r|ψ h j . (2)
In this expression, e 0 is the electron charge and r is the position operator. The right-hand side of Equation ( 1) is proportional to the experimentally accessible absorption spectrum. To make a proper comparison possible, a finite broadening which emulates the energy resolution of the UV-Vis spectra has additionally been incorporated (see Figure 1e).

Results

Experimental Absorption Spectra of CdSe x S 1-x Quantum Dots A set of CdSe x S 1-x QDs has been synthesized with an average diameter d of 3 nm and a composition ranging from x = 0 to x = 1. As illustrated by the transmission electron microscopy (TEM) image in Figure 2a, these QDs are quasi-spherical and monodisperse.

Moreover, the synthesis method we have used gives QDs that truly consist of solid solutions

as it has been previously demonstrated by Raman spectroscopy. 20 As exemplified by the set of absorption spectra of ≈ 3 nm CdSe x S 1-x QDs shown in Figure 2b, the homogeneous alloying of the anions in these nanocrystals enables their band gap to be tuned through the visible spectrum while keeping their size constant.

Interestingly, the absorption spectra of CdSe x S 1-x also contain an additional feature that is best visible for QDs with a high Se content (x = 0.1-0.5). Figure 2c represents the Holes entire set of absorption spectra of CdSe x S 1-x QDs as a 2D wavelength/composition color map. One sees that the addition of Se to CdS QDs appears to split the absorption peak related to the first exciton transition, a splitting that is more pronounced with increasing Se content until it tends to disappear at x ≈ 0.5. Importantly, as shown in the Supporting Information (Figure S1) this splitting persists when changing the QD diameter -although it becomes hard to discern for larger alloyed QDs -and a comparable effect has been reported by Swafford et al., 4 Ouyang et al. 21 and Zou et al. 22 Since all these authors use different syntheses to form CdSe x S 1-x QDs, our conjecture is that this additional feature is intrinsic to CdSe x S 1-x QDs, and not due to the unintended formation of, e.g., bimodal size distributions.

ψ h 1 ψ h 2 ψ h 3 ψ h 4 ψ h 1 ψ h 2 ψ h 3 ψ h 4 ψ h 1 ψ h 2 ψ h 3 ψ h 4 x = 0 x = 0.5 x = 1
To elucidate the origin of this additional transition and its link to the formation of random alloys, we therefore compared the experimental results with theoretical calculations within the previously described framework of an ETBM combined with configuration interaction.

Calculated Single-Particle States

As outlined in the background section, the spectrum of exciton states of QDs is obtained by calculating first the single-particle states with an ETBM appropriate for modeling alloyed QDs, followed by determining the many-particle exciton states by configuration interaction using Slater determinants of the single particle states as a basis. Although the experimental data should be best compared with the spectrum of exciton states obtained after configuration interaction is taken into account, we will first describe the single particle states obtained by diagonalization of the TB Hamiltonian. This approach anticipates the result that the energetically lowest exciton states can be approximately assigned to transitions between two single-particle levels, shifted by the respective electron-hole Coulomb interaction, when the broadening of Equation 1 is reduced. 23 The composition-resolved single-particle spectra e i (x), h j (x) including their stochastically determined relative frequency can be found in the Supporting Information. Figure 3 represents an example of the corresponding probability densities of the first four electron and hole states for pure CdS (x = 0), pure CdSe (x = 1) and a CdSe 0.5 S 0.5 alloy (x = 0.5). One sees that in the case of CdS and CdSe, the electron states can be readily classified as s and p x , p y and p z like and s, s, p x,y,z and p ± = (1/ √ 2)(p x ± ip y ) like, respectively. In the case of the hole states on the other hand, a strict classification by their nodal structure is not possible even for the pure QDs. Hole band mixing effects render them slightly anisotropic, making them show traits of both s and p like behaviour. If we nevertheless roughly classify them by the number of nodes of the envelopes of the wave function/probability density, we can identify groups of two s like and two p ± like hole states in the pure QDs. Note that without the Coulomb interaction, the hole level ordering as depicted in Figure 3 is not identical for pure CdSe and CdS QDs, as the energetically highest holes have p character for x = 1 but s character for x = 0. Due to the small hole level spacing however, the single-particle level order is often sensitive to the input parameter set. 24 Moreover, when the Coulomb interaction is subsequently included, the lowest transition is in both cases bright and involves an s like electron and an s like hole. The larger wavefunction overlap renders the respective Coulomb attraction larger for s-s than for s-p combinations, a difference exceeding the corresponding hole level splitting. This emphasizes the general importance of the simulation of excited-state properties rather than staying in the singleparticle picture when discussing electronic optical properties of QDs.

The coarseness and faceting of the isosurfaces that is visible in the pure QDs has two origins. First, the non-sphericity of the potential landscape as induced by the QDs surface facets and the crystal lattice results in a real, physical effect on the wave functions that is reproduced by real-space lattice models. Second, the interpolation of the wave functions from the fcc lattice to a cartesian grid leads to a coarsening artefact due to the visualization procedure. For the alloyed QDs, strong disorder-induced distortions of the isosurfaces can be identified and the symmetry character of the wave functions is visibly affected.

When we count the twofold Kramers degeneracy for each single-particle level in by introduction of an additional index {+, -}, the electron ground state ψ e 1+ , ψ e 1-is twofold degenerate, while the hole ground state levels ψ h 1+ , ψ h 1-, ψ h 2+ , ψ h 2-as well as the first excited hole states ψ h 3+ , ψ h 3-, ψ h 4+ , ψ h 4-show a fourfold degeneracy in both the CdSe and the CdS QDs. This agrees with results from other theoretical approaches for pure zincblende CdSe QDs, 25,26 while similar QDs with a wurtzite structure will have a slightly different level structure due to the lower symmetry of the system. 18,27 In total, the lowest excitonic state in the pure unalloyed QDs thus shows an eightfold degeneracy, which will subsequently be lifted on a finer energy scale by electron-hole exchange effects and distortions from a quasi-spherical shape. 25 Due to the configurational disorder in the alloyed QDs, however, all point group symmetries of the potential landscape are lost and the only remaining symmetry is the timereversal invariance, which leads to the essential Kramers doublets of the single-particle levels also for the alloyed QDs. As the symmetry character of the wave functions can in principle be deduced from symmetries of the corresponding Hamiltonian by group theoretical considerations, the wave functions lose their spatial symmetry properties in the disordered QDs Information of Reference 28. Some disparities to the experimental results can be spotted on the low energy/large wavelength side of the Se-rich part of the spectrum, where the splitting between the first two main absorption bands is slightly too large in theory. Also, an additional small peak just below 500 nm develops when Coulomb effects are switched on, which is not identifiable in the UV-Vis data. However, a more sophisticated treatment of the temperature dependence, the work function of the QDs (which is in reality finite and should in principle be concentration dependent as well) and the inclusion of size dispersion effects will affect the calculated optical spectra. This may further improve the quantitative agreement between calculated and experimental spectra, especially for small QDs with large confinement energies as used here.

Calculated Optical Spectra

Bowing of the Optical Gap

The energetically lowest bright transition is an excited-state property that defines the optical gap E opt g (x). The bowing coefficient b of an A x B 1-x alloy is obtained by fitting the composition dependent gap E g (x) to

E g (x) = x E A g + (1 -x) E B g -x(1 -x) b, (3)
where E A g and E B g are the energy gaps of the pure A and B constituents, respectively.

Overall, we observe a very good agreement between theory and experiment. When Equa-

This presumed weak size-dependence is also in agreement with theoretical findings for (Cd,Zn)Se QDs, 24 but has to be separately analyzed for each size in each material system under consideration due to the compensatory effects of Coulomb interaction and quantum confinement on the optical bowing.

Finally, it should be noted that the here obtained excellent agreement between theoretical and experimental bowing supposedly slightly exceeds the predictive power of our TB calculations. First, the spatial resolution is inherently limited to the underlying crystal lattice. Second, several model assumptions (see Methods section for a detailed discussion) were made to incorporate the effect of substitutional disorder into the TB calculation scheme.

A closer look at Figure 5 reveals that the experimental and theoretical data points do in fact show some discrepancy. However, the influence of this dispersion on the bowing parameter is mostly compensated by the small deviation of the curve from a parabola and the slightly larger theoretical result for the x = 0 (pure CdS) end point of the curve. Although the numerical solution of the eigenvalue problem is essentially exact, the complexity of the underlying theory renders any attempt to evaluate error bars unreasonable. Nevertheless, we cautiously claim that the good agreement between theory and experiment a posteriori confirms the validity of our model assumptions. Figure 6: Ensemble-averaged transition probabilities as given by the modulus squared dipole matrix elements |d eh ij | 2 between the lowest electron state and the highest 40 hole states for each value of the Se fraction x. The energetically lowest main absorption feature can be attributed to hole levels h j<13 . In the alloyed cases 0 < x < 1 a broad band of non-vanishing transitions develops for deeper hole states, identifiable by the "ridge" spanning across the right half of the bar plot.

Understanding the Additional Transition in CdSe

x S 1-x Quantum

Dots

The results obtained from theory allow us to further analyze the origin of the additional absorption channel in a single-particle representation. Figure 6 shows the calculated ensembleaveraged transition probabilities as given by the modulus squared dipole MEs |d eh ij | 2 between the lowest electron state and the highest 40 hole states for each value of the Se fraction x (see Method section for further details on the calculation). The corresponding transitions cover the energy range between the first two main absorption features on the red (low energy) side of the spectrum. In the pure limits x = 0 and x = 1, those can be attributed to allowed transitions between levels e 1 -h j and e 2,3,4 -h j , where e 2 , e 3 , e 4 are quasi-degenerate. The light polarization has been arbitrarily chosen to the [100] direction in the disordered system under consideration. Please note again that the single-particle indices do not assign a strict 15 (for x = 0). Depicted are isosurfaces of 50% of the maximum value and |ψ(x, y, z = const.)| 2 for a (001) cut plane through the center of the QD, respectively. A replica of Figure 6 is included as a guide to the eye. energetic hierarchy, as the Coulomb interaction affects the levels in a different way, e.g., the x = 1 ground state transition is not dark in the interacting case, as already discussed above.

In contrast to previous work from the literature, which have exhaustively analyzed the nature of the lowest transitions for unalloyed II-VI QDs, 18,25,27,30 we focus here on the band of deeper valence band states, i. e. , the levels h j with j > 10. In pure CdS and CdSe (x = 0 and x = 1), most of the corresponding transitions from the electronic ground state are either dipole forbidden or show very low oscillator strengths, except for three transitions, at h 13-15 for x = 0 and at h 34 for x = 1. 31 These hole states show a relatively large probability density near the center of the QD, which leads to a considerable overlap with the electron ground state. For the sake of completeness, they are exemplarily depicted in Figure 7.

Although a classification of these highly excited states according to their nodal structure is ambiguous due to the band mixing effects, their non-vanishing overlap with the even-parity s-like electron ground state indicates an even-parity character at least for the envelope ξ h of these hole states, as the dipole MEs in envelope function approximation (EFA) are directly proportional to the overlap integral ξ e i |ξ h j . 32 It should be emphasized, however, that our analysis goes beyond the (multiband) EFA, as the here employed TB approach also probes the orbital character of the underlying effective Wannier orbitals, as the corresponding site amplitudes enter the dipole MEs (see Equation 5).

Opposite from the pure QDs, Figure 6 shows that all transitions acquire a non-zero oscillator strength in the case of alloyed CdSe x S 1-x QDs. As an example, Figure 8 shows the corresponding modulus squared electron (top row) and hole wave functions (bottom row)

for the electron ground state |ψ e 1 | 2 and the excited hole state |ψ h 8 | 2 for Se fractions x = 0 (pure CdS), x = 0.5 (arbitrary configuration) and x = 1 (pure CdSe). To interpret it again in the simple EFA picture, in the pure cases, the hole wave function envelope has uneven parity, resulting in a zero overlap with the even parity electron ground state and thus a vanishing oscillator strength for the corresponding electronic transition. For the disordered alloy however, both the electron and the hole state lose their well-defined symmetry, which leads to a non-zero overlap and a concomitant brightening of the electronic transition.

Energetically, the transitions which form the clearly identifiable non-vanishing block with concentrations 0 < x < 1 and hole levels h 1 to h 12 form the first main transition feature.

Looking at Figure 6, one notices that an additional broad band of non-vanishing transitions develops for h j with j > 12 that form a clearly identifiable ridge that moves from the h 13-15

region for x = 0 (CdS) to h 34 for x = 1 (CdSe). Figure 4 shows that only the part 0.5 ≥ x > 0 of this ridge is truly visible in the optical spectra as the additional transition that appears upon addition of Se to CdS QDs. At the Se-rich side, the energy at which the transition occurs merges with the continuous second absorption band that contains transitions into the excited electronic states ψ e 2 , ψ e 3 , ψ e 4 . Interestingly, this ridge already appears between the two partially allowed deep-hole transitions of pure CdS and CdSe when only single stochastic realizations of the alloyed QDs are studied for each composition and a quasi-Gaussian shape including the spectral broadening quickly develops for low realization numbers N (see Supporting Information).

We already showed that the additional feature visibly develops from transitions from hole states that are initially also allowed in the pure CdS and CdSe material, but energetically merge with the first (for x = 0) and second (for x = 1) main transitions. Therefore, it seems natural to interpret the additional line as a transition from deep valence band states into the lowest electron state that is initially hidden in the main absorption bands due to the high degeneracy of the unperturbed hole levels. When alloying is introduced, all degeneracies besides the time-reversal Kramers symmetry are lifted and the feature is energetically shifted in-between the two lowest main absorption lines. Also, it is subsequently spectrally broadened due to the inherent disorder, which distorts the wave functions and thus relaxes corresponding selection rules. No virtual-crystal like mean-field calculation can reproduce this effect, as it does not only lack broadening effects, but also cannot reproduce the nonlin-ear nature of alloying-induced energy shifts without additional free parameters. 24 Dynamical mean-field approaches like the coherent potential approximation (CPA) can reproduce nonlinear shifts and broadening effects, but cannot simulate the lifting of level degeneracies.

Furthermore, the CPA is exact in the limit of infinite coordination numbers, making it in reverse not suitable for the application to nanosystems with only few lattice sites. It should also be kept in mind that the high geometrical frustration of the underlying fcc crystal lattice will in general always connect alloying and disorder effects, as corresponding binary superstructures are non-trivial. 33

Discussion

The absorption spectrum provides a quick and convenient way of analyzing the properties of a dispersion of colloidal QDs. In particular the position and the width of the first exciton peak can be directly translated into an estimate of the mean diameter and the size dispersion of the QDs. The finding that in the case of CdSe x S 1-x alloyed QDs, the absorption spectrum contains additional features related to homogeneous alloying further strengthens the position of absorption spectroscopy as the first method of choice to analyze sample quality. Since this additional feature is closely linked to the breaking of selection rules of electronic transitions induced by random alloying, we expect that other alloyed quantum dots may show a similar fingerprint of random alloying in their absorption spectrum.

Next to the band gap tuning resulting from homogeneous alloying, the formation of core/shell QDs with graded, alloyed interfaces has proven a crucial step towards non-blinking QDs. In these systems, the strong reduction of Auger recombination rates of trion or multiexciton states is linked to the smoothening of the confinement potential due to a gradual change in composition between core and shell. 34 The close correspondence found here between the experimental and simulated absorption spectra of CdSe x S 1-x alloyed QDs indicates that the ETBM used here yields an accurate description of the opto-electronic properties of QDs with compositional disorder. This indicates that also the electronic states of graded core/shell QDs can be determined and Auger rates be calculated within the framework of an atomistic disorder model, in an effort to go beyond the continuum-model based calculations in Reference 34. Recently, changes in the surface composition and termination 35 were also found as a means to tune the luminescence behavior of colloidal QDs. However, the optical properties of surface-composition tuned QDs can most likely not reliably be treated within our calculation scheme, as the ETBM is not the method of choice for the predictive simulation of surface and interface effects on the nanometer scale.

The markedly distorted probability distributions found for randomly alloyed QDs as compared to the single compounds indicates that also more detailed opto-electronic properties, such as biexciton addition energies, the electron-phonon coupling strength or the fine structure splitting (FSS) of the lowest exciton transition, may be significantly affected by random alloying. Regarding future tasks, not only the simulation of surface effects, but also the analysis of excited-state properties on a much smaller energy scale than in the present work (e.g., the composition dependence of the FSS) 36 might show the limitations of the underlying ETBM approach for very small QD systems. Consequently, we have recently begun to complement our calculation scheme using conceptually more sophisticated, but computationally far more expensive, methods like DFT-LDA for structural/vibrational properties 37,38 and pseudopotential calculations for electronic properties. 39,40 As part of ongoing research, these calculations will supplement the TB results not only for CdSe x S 1-x but also for further QD material systems. Due to the high computational effort, however, their combination with the stochastic framework as presented in this paper will likely be limited to smaller basis sets for the CI expansion, pre-selected compositions and smaller ensemble sizes than in the ETBM at least in the near future.

Conclusions

In summary, we have shown that an additional transition appears in absorption spectra of alloyed CdSe x S 1-x QDs. In order to investigate this effect, we carried out theoretical calculations with an empirical tight-binding model, augmented by the configuration interaction scheme for the many-body properties. Together with a stochastic approach used to describe the distribution of elements in the crystal lattice, this model has proven itself very efficient to reproduce experimental data, such as the occurrence of additional features in the absorption spectra and the bowing of the optical gap.

In addition, we found that the additional feature present in the S-rich absorption spectrum is unambiguously found in homogeneously alloyed CdSe x S 1-x QDs in experiment as well as in theory. It could therefore be used as a signature confirming that a given synthesis method gives homogeneous CdSe x S 1-x alloyed QDs and not QDs with a gradient composition. In that view, our theoretical calculation method can also be implemented for other alloyed QD systems. For example, in the case of (Cd,Zn)Se QDs, 11,24 the difference of reactivity between Cd and Zn precursors makes it difficult to create homogeneous alloys.

A comparison between theoretical and experimental absorption spectra may therefore give useful indications on the degree of homogeneity and disorder in their structure. However, such an analysis will require a thoroughly performed simulation and, if possible, measurement of the effect of different composition profiles on the optical properties for each material system, which is beyond the scope of the present study. As indicated in the present paper, it might in many cases also be sufficient to perform calculations on the single-particle level in order to qualitatively identify corresponding features, which does significantly limit the computational and conceptional effort.

This work also gives evidence that the selection rules of the optical transitions are relaxed because of the local substitutional disorder that is introduced in the structure of alloyed QDs.

As a result, by comparison to pure QDs, the symmetry character of the charge carriers' wave functions is affected in alloyed QDs and, in the case of originally forbidden transitions, the overlap of the electrons and holes states can become non-negligible. Together with the obligatory lifting of symmetry-related level degeneracies, alloying thus gives large tunability of the electronic and optical properties, but also affects those characteristics in a qualitative way.

Methods

Synthesis of Alloyed CdSe x S 1-x QDs. Zincblende CdSe x S 1-x alloyed QDs were synthesized according to a previously reported procedure. 20 This method makes use of a heterogeneous dispersion of black Se powder in 1-octadecene (Se-ODE), and a homogenous S-ODE solution, where the S is dissolved in the ODE by heating. For the synthesis of the Cd(Se,S) alloyed QDs, a mixture of these two precursors was injected under air in a reaction mixture made of CdO dissolved in ODE using myristic acid and kept at 270 °C. After 5 minutes of reaction at 260 °C, the reaction was quenched, the QDs were purified by precipitation and centrifugation, and the ligands were exchanged for oleic acid. ∆ so , the conventional lattice constant a, and additionally to the X-point energies of the conduction band X 1c , the HH and LH bands X 5v (degenerate at X), and the SO band X 3v (see Figure 1b). The input material parameters for zincblende CdSe and CdS including their sources are explicitly given in Table 1. The temperature dependence has been incorporated merely by the use of room temperature band gaps in the calculations, as the corresponding input data for the other parameters like band curvatures or BZ edge energies is only sparsely available for the material system under consideration.

Characterization of Alloyed

Tight-binding simulation of random alloys. The TBME for the pure constituents can be used for the calculation of disordered CdSe x S 1-x QDs in a stochastic approach. We simulate the alloying as a random process, i.e., neither are clustering effects included, nor do we distinguish whether some realizations would in theory be energetically more favorable than others. Assuming uncorrelated substitutional disorder, each unit cell of the QD then contains the corresponding CdSe or CdS matrix elements (ME), with probabilities x or 1 -x, and intersite (hopping) ME between CdSe and CdS unit cells are approximated by the arithmetic average. The resolution on the scale of unit cells allows for an unambiguous assignment of the site-diagonal TBME to the CdSe or CdS input band structure, whereas atomic TB parameters would need an averaging procedure, 46 as they depend on the binary host material. 47 To validate this approach and the input parameters first, TB calculations on finite supercells with periodic boundary conditions have been performed to emulate a bulk alloy. 48 Here, supercells with 4000 unit cells and 30 different configurations per composition were sufficient to obtain convergence (see Supporting Information for a convergence test).

When Equation (3) is applied, this yields a bulk bowing parameter of b ≈ (0.46 ± 0.03) eV. This value corresponds to the d → ∞ bulk limit of the QD and is in agreement with the findings of Swafford et al. in Reference 4 considering their experimental accuracy.

To use this empirical tight-binding method (ETBM) for simulating the properties of an ensemble of alloyed QDs, where each single QD shows a different microscopic distribution of Se and S atoms, N = 50 randomized configurations are generated for various compositions in steps of ∆x = 0.1 (see Figure 1c). With this number of configurations, all energies, e. g., the spectral positions of the bright transitions, converged on the energy scale relevant for the discussion in this paper (accuracy around 10 meV) regardless of the composition. In order to concentrate on the effect of alloying and configurational disorder on the electronic and optical properties of the QDs, the following model assumptions were made:

1. Similar to References 26,49-51, the hopping probability at the surface is set to zero and the surrounding medium is shifted far above the electron energies calculated here, which corresponds to a perfect surface passivation. The influence of the height of the potential barrier, dangling bonds and surface relaxations have extensively been discussed in the literature on different levels of sophistication, [52][53][54] while the simulation of the influence of ligand-related states goes beyond the scope of the present model.

2. Due to the usage of fcc-centered Wannier orbitals, all effects on a smaller length scale like local bond length modifications due to lattice distortions are absorbed into the hopping ME between unit cells of different binary materials. As we found the results in the case of CdSe and CdS were not sensitive to the exact treatment of such hoppings, bond length changes and angle distortions are neglected. TB approaches for QDs with slightly larger lattice mismatch obtained good results without explicitly taking angular distortions into account. 55 Pure core QDs show practically no strain. 56 In core/shell structures, however, strain and relaxation effects and resulting internal fields will likely have a larger impact on the electronic properties and should be incorporated on an appropriate level. [57][58][59] 3. The energetic offset between CdSe and CdS is incorporated via the bulk interface valence band offset ∆E vb .

Calculation of

d eh ij = e 0 RR αα c i,e Rα * c j,h R α e p • R δ RR δ αα . (5)
Even in the non-alloyed CdSe and CdS QDs, the dipole selection rules cannot be obtained from total angular momentum selection rules, as valence band mixing effects will hinder a corresponding classification especially for the hole states. Furthermore, the small size of the QDs leads to significant deviations from a spherically symmetric potential due to surface faceting and the underlying crystal structure. Such effects can also be incorporated into multiband effective-mass/k • p models. 25 However, even an explicit numerical calculation of the dipole MEs then cannot simulate the influence of the alloying-induced disorder, as the envelope function approximation lacks a lattice-resolved discretization of the potential landscape and thus only allows the incorporation of alloy effects on a mean-field level. In contrast, our approach incorporates the influence of the visible wave function distortions on the selection rules. Nevertheless, it should be mentioned that d eh ij as calculated in Equation 5neglects the spatial dependence of the wave function inside the unit cells, the so-called "orbital part". If desired, it can be calculated using (orthogonalized) Slater orbitals. 55,62 This part is of course automatically included in more sophisticated, but computationally also much more costful pseudopotential treatments like Reference 40. As it is not necessary to explain the optical feature examined in the present paper and its relative contribution was found to be small for other QD systems, 55 the orbital part has been neglected.

Simulation of Coulomb Interaction and Screening. The interaction between the confined charge carriers is incorporated via the excitonic Coulomb interaction Hamiltonian

H Coul = - ijkl V eh ijkl a † i,e a † j,h a k,h a l,e . (6)
Here, V eh ijkl are the Coulomb MEs and a † i,e (a † i,h )/a i,e (a i,h ) denote the creation/annihilation operators for electrons (holes) in the states |ψ e i (|ψ h i ). The calculation has been performed on the same level of approximation as in References 24,63, leading to (29) The error range here only takes the deviation from a parabolic behavior into account.

Figure 2 :

 2 Figure 2: Experimental results for the 3 nm sized alloyed CdSe x S 1-x QDs. (a) TEM image for Se fraction CdSe 0.5 S 0.5 . (b) UV-Vis absorption spectra for several Se fractions. (c) Colormap representation of absorption spectra for the entire set of the CdSe x S 1-x QDs. Solid lines mark the first two main transitions. The additional, alloying-induced transition is marked with a dashed line.

Figure 3 :

 3 Figure 3: Modulus squared electron and hole envelope wave functions |ψ i | 2 = |c i Rα | 2 of the first four electron and hole states for Se fractions x = 0, x = 0.5 (one arbitrary configuration) and x = 1. Depicted are isosurfaces of 50% of the maximum value.

Figure 4 :

 4 Figure 4: Optical spectra as obtained from the stochastic ETBM. (a) Including Coulomb interaction. (b) Coulomb effects artificially switched off.

Figure 4

 4 Figure4represents the optical properties of CdSe x S 1-x QDs as a function of wavelength and composition as calculated with or without Coulomb interactions according to the procedure described in the background and the methods section. In order to facilitate the comparison between experiment and theory, a constant blueshift has been applied to the theoretical spectra, which matches the calculated wavelength of the first exciton transition for CdSe (x = 1) QDs to the experimental value. This shift absorbs influences from model assumptions (infinite work function) as well as input parameter ambiguities (e. g., size dispersion, temperature dependence of band curvature parameters). On the other hand, since we kept the shift energetically constant over the whole composition range, it does neither alter the nature of the transition channels nor the bowing discussed further below. Aside from this cosmetic correction, no free parameters are involved in the present calculation of optical properties, as all input parameters are directly connected to experimentally or theoretically accessible properties. The colormap in Figure4aincludes the Coulomb interaction between

  Figure 5 brings together the experimental and theoretical values for E opt g

Figure 5 :

 5 Figure 5: Experimental and theoretical results for the composition dependence of the optical gap E opt g (x) of d ≈ 3 nm CdSe x S 1-x QDs including least-squares parabolic fits with bowing parameter b.

  tion 3 is applied to E opt g (x), a least-squares fit yields a theoretical bowing parameter of b = (0.48 ± 0.04) eV, while the experimental data results in a value amounting to b = (0.47 ± 0.03) eV. 29 These values agree with the findings of Swafford et al. in Reference 4 within their range of accuracy. Interestingly, our TB calculations yield a d → ∞ bulk limit of b = (0.46 ± 0.03) eV (see Methods section). Together with the also size-independent results from Swafford et al., this suggests that the size dependence of the bowing for CdSe x S 1-xQDs might be negligible in practice. For effective mass or k • p approaches (also with diameters larger than 3 nm), we therefore recommend a constant room temperature bowing coefficient of b Cd(Se,S) ≈ 0.47 eV.

Figure 7 :

 7 Figure 7: Modulus squared electron (top row) and hole (bottom row) envelope wave functions |ψ i | 2 = |c iRα | 2 for the electron ground state ψ e 1 and the excited hole state ψ h 34 (for x = 1) and ψ h 15 (for x = 0). Depicted are isosurfaces of 50% of the maximum value and |ψ(x, y, z = const.)| 2 for a (001) cut plane through the center of the QD, respectively. A replica of Figure6is included as a guide to the eye.

Figure 8 :

 8 Figure 8: Modulus squared electron (top row) and hole wave functions (bottom row) for the electron ground state |ψ e 1 | 2 and the excited hole state |ψ h 8 | 2 for Se fractions x = 0 (pure CdS), x = 0.5 (arbitrary configuration) and x = 1 (pure CdSe). Depicted are isosurfaces of 50% of the maximum value and |ψ(x, y, z = const.)| 2 for a (001) cut plane through the center of the QD, respectively.

  CdSe x S 1-x QDs. The size and morphology of the QDs was determined by TEM using a Cs-corrected JEOL 2200 FS microscope. The zincblende crystal structure has been confirmed by X-Ray diffraction (XRD) analyses (see also FigureS3in the Supporting Information of Reference 20). The actual composition of the QDs was determined by X-ray fluorescence (XRF) analysis using a Rigaku NEX CG EDXRF spectrometer and/or by energy dispersive X-ray spectroscopy (EDX) jointly with the TEM analyses. Absorption spectra were taken using a Perking Elmer Lambda 950 spectrometer. Tight-binding parametrization scheme. Using the parametrization scheme from Reference 41, the resulting TB matrix elements (TBME) Rα| H |R α contain couplings up to second nearest neighbors and fit the band structure to the bulk band gap E g , the conduction band effective mass m e , the three Luttinger parameters γ 1 , γ 2 , γ 3 , the spin-orbit splitting

  Dipole Matrix Elements. Following References 51,60,61, the position operator in a TB model can be approximated by r ≈ Rα |Rα R Rα|. With the light polarization vector e p , the dipole MEs d eh ij = e p • d eh ij are then given as

2 0

 2 Rα ) * (c j,h R α ) * c k,h R α c l,e Rα e 4πε 0 ε r |R -R | ,(7)with ε 0 as the vacuum permettivity. The R = R contribution is treated semi-analytically following Reference 64. The calculation of the concentration-, size-and distance-dependentcharge screening ε r = ε QD ∞ (x, d, |R -R |) (ingeneralization of an approach first used by Franceschetti et al. for unalloyed QDs 18 ) is given in the Supporting Information. It should be noted that alternatives to this rather simple approach (which is based on the generalized Penn model 65 ) exist in the literature,

( 21 )

 21 Ouyang, J.; Vincent, M.; Kingston, D.; Descours, P.; Boivineau, T.; Zaman, M. B.; Wu, X.; Yu, K. J. Phys. Chem. C 2009, 113, 5193. (22) Zou, Y.; Li, D.; Yang, D. Nanoscale Res Lett 2010, 5, 966. (23) The finite broadening as used for Figure 4 will merge bright transitions to higher hole states into the visible peaks. (24) Mourad, D.; Czycholl, G. Eur. Phys. J. B 2010, 78, 497. (25) Efros, A. L.; Rosen, M.; Kuno, M.; Nirmal, M.; Norris, D. J.; Bawendi, M. Phys. Rev. B 1996, 54, 4843. (26) Schulz, S.; Czycholl, G. Phys. Rev. B 2005, 72, 165317. (27) Korkusinski, M.; Voznyy, O.; Hawrylak, P. Phys. Rev. B 2010, 82, 245304. (28) Jasieniak, J.; Califano, M.; Watkins, S. E. ACS Nano 2011, 5, 5888.

Table 1 :

 1 Material parameters including references.

			CdSe	Ref. CdS	Ref.
	E g	(eV)	1.68	42	2.46	42
	m c	(m 0 ) *	0.12	43	0.14	42
	γ 1		3.33	43	4.11	42
	γ 2		1.11	43	0.77	42
	γ 3		1.45	43	1.53	42
	∆ so	(eV)	0.41	43	0.078	42
	X c 1	(eV)	2.94	44	5.24	42
	X v 5	(eV)	-1.98	44 -1.9	42
	X v 3	(eV)	-4.28	44 -4.8	42
	∆E vb (eV)	0.4	45	0	45
	a	( Å)	6.078	43	5.825	42

* Free electron rest mass.

  66,67 if a more sophisticated treatment of screening effects is needed. Consistent with theoretical results for similar systems, 27 additional electron-hole exchange matrix elements were found to be considerably smaller and thus do not contribute on the energy scale discussed in the present paper.to build a proper Wannier basis by suitable linear combinations, we refrain from such a distinction. Modern ab initio band structure suites (VASP, Quantum ESPRESSO) also offer postprocessing interfaces for a direct calculation of Wannier states with the Wannier90 algorithm.68 (17) To prevent confusion, we do not count the obligatory Kramers degeneracy of each single-particle state in, if not explicitly mentioned otherwise. (18) Franceschetti, A.; Fu, H.; Wang, L. W.; Zunger, A. Phys. Rev. B 1999, 60, 1819.(19) Bester, G. J. Phys.: Condens. Matter 2009, 21, 023202. (20) Aubert, T.; Cirillo, M.; Flamee, S.; Van Deun, R.; Lange, H.; Thomsen, C.; Hens, Z.

Chem. Mater. 2013, 25, 2388.

Acknowledgement

D. Mourad would like to thank G. Czycholl, G. Bester and J. Jackson for helpful discussions, thankfully acknowledges a grant of the North-German Supercomputing Alliance HLRN for computation time and also acknowledges the developers of the free Python numpy and matplotlib modules. Z. Hens acknowledges BelSPo (IAP 7.35, photonics@be) and the FWO-Vlaanderen (research project G.0794.12) for funding.

Supporting Information Available

The Supporting Information contains the calculated CdSe x S 1-x bulk bowing including convergency behaviour with respect to the supercell size, the calculated spectral distribution of single-particle electron and hole levels and additional transition probabilities for a single stochastical realization. Furthermore, it contains the determination of the size-dependent screening constant and the screening radius. This material is available free of charge via the Internet at http://pubs.acs.org/.

Notes and References

(1) Rogach, A. L.; Kornowski, A.; Gao, M.; Eychmüller, A.; Weller, H. J. Phys. Chem. B 1999, 103, 3065.

(2) Talapin, D. V.; Lee, J.-S.; Kovalenko, M. V.; Shevchenko, E. V. Chem. Rev. 2010, 110,

389.

(3) Bailey, R. E.; Nie, S. J. Am. Chem. Soc. 2003, 125, 7100.


Oscar Laurent

Jean-François Moyen

Jörn-Frederik Wotzlaw

Jana Björnsen

Olivier Bachmann

Early Earth zircons formed in residual granitic melts produced by tonalite differentiation

The oldest geological materials on Earth are Hadean (>4 Ga) detrital zircon grains. Their chemistry and apparently low Ti-in-zircon temperatures (≤700°C) are considered to be inconsistent with crystallization in a magma of the tonalite-trondhjemite-granodiorite (TTG) suite, although these are the dominant Archean (4.0-2.5 Ga old) silicic rocks. Using a new dataset of trace element contents in zircons from Paleoarchean Barberton TTGs (South Africa) and thermodynamic modelling, we show that these zircons have crystallized at near-solidus conditions from a compositionally uniform granitic melt. This melt is residual from the upper crustal crystallization of a less evolved (tonalitic) parent and thereby shows major and trace element compositions different from bulk TTG rocks. A global compilation reveals that most Hadean detrital and Archean TTG-hosted grains share a peculiar zircon trace element signature that is distinct from the chemical trends defined by Phanerozoic zircons. Our model shows that the low Ti contents of early Earth zircons reflect crystallization at higher temperatures (720-800°C) than initially inferred, due to lower modelled TiO2 activity in the melt relative to previous estimates. We therefore propose that near-solidus zircon crystallization from a chemically evolved melt in a TTG-like magmatic environment was the dominant zircon-forming process on the early Earth.

INTRODUCTION

The only remnants of the earliest felsic crust are Hadean (4.4-4.0 Ga-old) detrital zircon grains found in younger sedimentary rocks at Jack Hills, Western Australia (e.g. [START_REF] Harrison | Hadean Jack Hills zircon geochemistry[END_REF].

Constraining the nature of the magma in which these zircons formed is key to understanding the geodynamic context of the early Earth. Although rocks of the tonalite-trondhjemite-granodiorite (TTG) suite dominate the earliest preserved Archean (4.0-2.5 Ga old) continental nuclei [START_REF] Moyen | The composite Archaean grey gneisses: Petrological significance, and evidence for a non-unique tectonic setting for Archaean crustal growth[END_REF], they have usually been discounted as potential parents for the Jack Hills zircons (e.g. [START_REF] Harrison | Hadean Jack Hills zircon geochemistry[END_REF]. Specifically, the apparent crystallization temperatures of Jack Hills zircons (≤700 °C based on Ti-in-zircon thermometry; [START_REF] Watson | Zircon thermometer reveals minimum melting conditions on earliest Earth[END_REF] are lower than those of TTG magma formation (>850 °C), and their trace element compositions do not seem to be in equilibrium with those of TTG bulk rocks [START_REF] Rollinson | Ophiolitic trondhjemites: a possible analogue for Hadean felsic 'crust[END_REF][START_REF] Reimink | A comparison between zircons from the Acasta Gneiss Complex and the Jack Hills region[END_REF][START_REF] Turner | An andesitic source for Jack Hills zircon supports onset of plate tectonics in the Hadean[END_REF].

However, the study of zircon and zircon-hosted inclusions from a 3.81 Ga-old metatonalite from the Itsaq Gneiss Complex in Greenland showed that TTG magmas are initially zircon-undersaturated, so zircon should only crystallize at near-solidus conditions in residual granitic melts that formed after extensive cooling and crystallization [START_REF] Nutman | A granitic inclusion suite within igneous zircons from a 3.81 Ga tonalite (W. Greenland): Restrictions for Hadean crustal evolution studies using detrital zircons[END_REF]. Moreover, changes in melt chemistry during crystallization may result in zircons forming from a melt showing significantly different trace element compositions than TTG bulk rocks. These possibilities must be systematically addressed to evaluate whether TTG magmas could be plausible parents for Hadean detrital zircons.

For this purpose, we investigated the nature and composition of zircon-forming melts in TTG suites, based on new zircon trace element data from the Paleoarchean (3.46 Ga-old) TTGs of the Barberton area, Kaapvaal Craton, South Africa (Fig. 1). The compositional evolution of magmas in the course of crystallization, from the tonalitic starting compositions to the nearsolidus granitic melts, is well constrained following a recent study [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF]. The significance of the results is evaluated through a global comparison of Hadean and Archean zircons, based on a compilation of trace element compositions in Jack Hills detrital zircons (n = 233) and igneous TTG-hosted zircon (n = 947) in 44 samples from Archean cratons worldwide (Fig. 1A). Zircon compositions altered by the effects of radiation damage and/or post-magmatic (hydro-)thermal events were filtered out based on chemical and isotopic criteria from both newly acquired and compiled datasets. Details about data acquisition, compilation and filtering are provided in the Supplementary Material.

ZIRCON-FORMING MELTS IN ARCHEAN TTG SUITES

The 3.46 Ga Barberton TTGs consist of trondhjemite plutons with minor diorite, tonalite and granite; the plutons are associated with coeval and co-genetic felsic volcanic rocks of the Buck Ridge Volcanic Complex [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF] (Fig. 1B). The diverse igneous lithologies have been proposed to result from the differentiation of a parental liquid of tonalitic composition, whereby the compositional range of the plutonic rocks arises from various degrees of crystal accumulation and residual melt loss along the liquid line of descent (Fig. 1C).

Zircons from nine samples representative of the Barberton TTG suite show overlapping trace element concentrations (Fig. 2A and Supplementary Material, Fig. S1). In particular, Ti and Hf are powerful monitors of the crystallization environment of zircon, as Ti relates to the temperature and TiO2 activity of the host melt [START_REF] Ferry | New thermodynamic models and revised calibrations for the Ti-inzircon and Zr-in-rutile thermometers[END_REF] and Hf is a proxy for melt fractionation [START_REF] Claiborne | Tracking magmatic processes through Zr/Hf ratios in rocks and Hf and Ti zoning in zircons: an example from the Spirit Mountain batholith[END_REF]. Although bulk-rock SiO2 contents vary from ~60 to 76 wt%, zircon Ti and Hf concentrations in the Barberton samples show restricted ranges of 5.1 ± 3.2 ppm and 0.95 ± 0.17 wt.% (2 standard deviations [SD], n = 58) respectively (Fig. 2A).

Similar observations were made from Phanerozoic igneous suites and interpreted as reflecting late crystallization of zircon from a compositionally homogeneous residual melt in all rock types, for instance in the Boggy Plain zoned pluton, southeastern Australia [START_REF] Ickert | Ti in zircon from the Boggy Plain zoned pluton: implications for zircon petrology and Hadean tectonics[END_REF]. The observed distribution of zircon Ti and Hf concentrations in the Paleoarchean Barberton TTG plutons suggests a comparable situation.

The inversion of zircon compositions in the Barberton samples, using the parameterization of zircon/melt partition coefficients as a function of Ti concentration of [START_REF] Claiborne | Zircon as Magma Monitor: Robust, Temperature-Dependent Partition Coefficients from Glass and Zircon Surface and Rim Measurements from Natural Systems[END_REF], further supports this interpretation. The composition of the calculated melts indeed overlaps with that of the high-silica granitic melt inclusions from the BRVC (Fig. 2B), which represent the most differentiated, residual melts from crystallization of the parent tonalitic magma [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF]. Importantly, both the calculated melts and BRVC melt inclusions share the same distinctive differences with average TTGs, namely strikingly higher Th and U contents and a lower middle-to heavy rare earth element ratio (Fig. 2B). These features, which likely reflect the highly evolved nature of the melts and significant amphibole fractionation, explain the apparent disequilibrium between trace element concentrations of zircon and the composition of their host TTG rocks [START_REF] Reimink | A comparison between zircons from the Acasta Gneiss Complex and the Jack Hills region[END_REF].

We simulated the crystallization of the parental liquid of the Barberton TTG suite (as determined by [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF] at 3 kbar, using the mineral and melt models of [START_REF] Holland | Melting of peridotites through to granites: a simple thermodynamic model in the system KNCFMASHTOCr[END_REF] and the PERPLE_X software program (https://www.perplex.ethz.ch/; [START_REF] Connolly | The geodynamic equation of state: what and how[END_REF].

The phase equilibria and thermodynamic data obtained served to model the evolution of melt Zr concentration as well as zircon saturation temperatures and Ti contents (see Supplementary Material for details). The results show a first step (from the liquidus at 1050 °C to ~40 wt.% liquid at 725 °C) of dominant plagioclase (~40 wt.%), subordinate amphibole (~15 wt.%), and minor pyroxene and Fe-Ti oxide crystallization (Fig. 3A). The melt fraction then decreases from 725 °C to the solidus (~700 °C), corresponding to the near-eutectic crystallization of quartz, plagioclase, K-feldspar and biotite. These results match the modal compositions recorded by natural samples [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF] and the experimental results of upper crustal tonalite crystallization [START_REF] Marxer | Crystallization and zircon saturation of calc-alkaline tonalite from the Adamello Batholith at upper crustal conditions: an experimental study[END_REF].

The low initial Zr content of the system (50-200 ppm), combined with buffering of the melt Zr concentration by amphibole crystallization, result in a range of zircon saturation temperatures (considering model uncertainties) between 800 °C and 720 °C (Fig. 3A). Zircon crystallizes in a high-silica residual melt similar to the BRVC melt inclusions in both major (Fig. 1C) and trace elements (Fig. 2B), which agrees with the observation of zircon in the interstitial domains between plagioclase and hornblende phenocrysts in natural samples [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF]. These results are consistent with the inferred zircon crystallization temperatures (700-800°C) and granitic mineralogy of zircon-hosted inclusions in a 3.81 Ga-old tonalite from Greenland [START_REF] Nutman | A granitic inclusion suite within igneous zircons from a 3.81 Ga tonalite (W. Greenland): Restrictions for Hadean crustal evolution studies using detrital zircons[END_REF].

ZIRCON CRYSTALLIZATION CONDITIONS ON THE EARLY EARTH

Most TTG-hosted igneous zircons and many Jack Hills detrital zircons from our compilation compositionally overlap with the Paleoarchean Barberton zircons investigated here (Fig. 2A, 4 and Supplementary Material, Fig. S1). In Figure 4, we compared their Ti and Hf contents to those of zircon from various Phanerozoic igneous systems. Phanerozoic zircons show a negative correlation between Ti and Hf contents (Fig. 4) that reflects crystallization during magma cooling and differentiation [START_REF] Claiborne | Tracking magmatic processes through Zr/Hf ratios in rocks and Hf and Ti zoning in zircons: an example from the Spirit Mountain batholith[END_REF]. Such correlation may be recorded by zircons from single hand samples, as illustrated by one monzonite and one granite from the Miocene Spirit Mountain Batholith, southern Nevada, USA [START_REF] Claiborne | Tracking magmatic processes through Zr/Hf ratios in rocks and Hf and Ti zoning in zircons: an example from the Spirit Mountain batholith[END_REF] (Fig. 4). In turn, distinct parental melt compositions and magma evolution pathways result in distinct zircon Ti-Hf co-variations for different geodynamic settings (Fig. 4).

In contrast, most early Earth zircons define a tight compositional cluster in the Ti-Hf space (Fig. 4) with very low Ti contents (2-9 ppm) [START_REF] Carley | Iceland is not a magmatic analog for the Hadean: evidence from the zircon record[END_REF] at moderate Hf concentrations (0.8-1.2 wt%). Strikingly, the variations in zircon Ti and Hf contents in the >40 early Earth samples taken together are more restricted than in single samples of the Spirit Mountain Batholith (Fig. 4). Although the early Earth zircons partly overlap with the compositionally most evolved, modern convergent margin and post-orogenic zircons, they clearly lack the high Ti (>15 ppm) and/or low Hf (<0.7 wt%) parts of the Phanerozoic trends (Fig. 4). Altogether, this supports the existence of specific conditions of zircon growth for Hadean and Archean zircons. As defined based on the Barberton case, these conditions are the crystallization from a near-solidus, residual granitic melt resulting from differentiation of a tonalitic parent magma, with little (if any) zircon fractionation having occurred beforehand. This condition may reflect the generally lower Zr concentrations of early Earth magmas (<150 ppm) relative to Phanerozoic magmas [START_REF] Keller | Temporal variation in relative zircon abundance throughout Earth evolution: Geochemical Perspective Letters[END_REF], which would delay zircon saturation to near-solidus temperatures. This is confirmed by our crystallization model, in which zircon saturates only at ≤800 °C and accordingly shows the low Ti concentrations (<10 ppm) (Fig. 3A-B) typical of early Earth grains (Fig. 2A, Fig. 4). The modelled cumulative curve of Ti contents in zircon predicts a greater proportion of low-Ti zircons than observed in nature (Fig. 3B). This may reflect that during analysis of natural grains, thin zircon rims crystallized at low temperatures were undersampled relative to large, higher-temperature cores [START_REF] Ickert | Ti in zircon from the Boggy Plain zoned pluton: implications for zircon petrology and Hadean tectonics[END_REF]. However, the observed concentration range is exactly the same (2-9 ppm Ti; Fig. 3B). This is a robust indicator of the model's validity, as this range is controlled by zircon saturation temperature and TiO2 activity, which depend only on the initial Zr content of the liquid and phase relations during crystallization. Importantly, the model indicates that the low Ti concentrations of the Jack Hills zircons can be explained by crystallization at significantly higher temperatures (720-800 °C; Fig. 3B) than have been estimated so far (650-720 °C), because the TiO2 activity of the system was much lower (ca. 0.4-0.5 as modelled here) than was previously assumed for Ti-in-zircon thermometry (up to 1; [START_REF] Watson | Zircon thermometer reveals minimum melting conditions on earliest Earth[END_REF].

CONCLUSIONS AND IMPLICATIONS

We propose that (1) zircon in TTG suites forms at near-solidus temperatures (720-800 °C) from a residual granitic melt resulting from extensive crystallization of a tonalitic parent and, hence, compositionally different from the bulk rocks; and (2) this zircon-forming mechanism is common to all early Earth zircons, including the Jack Hills detrital grains, as shown by their shared peculiar trace element signature that is distinct from that of Phanerozoic zircon. This model reconciles the paradoxical observations that the low-Ti nature of the Jack Hills zircons [START_REF] Watson | Zircon thermometer reveals minimum melting conditions on earliest Earth[END_REF] reflects crystallization in a melt close to the granite minimum, although such granites are scarce in the early Earth [START_REF] Moyen | The composite Archaean grey gneisses: Petrological significance, and evidence for a non-unique tectonic setting for Archaean crustal growth[END_REF]. This liquid may in fact be similar to the most differentiated melt of the TTG suites, which remained trapped in the crystalline framework and/or is not preserved as co-magmatic silicic volcanic rocks [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF].

Whether the Jack Hills zircons formed in voluminous TTG plutons, as observed in Archean cratons, or in lower-volume felsic veins and pods in a dominantly mafic crust (e.g. [START_REF] Rollinson | Ophiolitic trondhjemites: a possible analogue for Hadean felsic 'crust[END_REF] is still uncertain. In fact, the parameters controlling zircon chemistry, namely the petrogenesis of felsic magmas (intracrustal differentiation of hydrous mafic material) and hence their phase relationships upon crystallization, would be comparable in both cases.

Nevertheless, closed-system crystallization as modelled here, which is more likely to occur in small, isolated melt pockets, may also explain the distributions of zircon Ti contents skewed to lower concentrations as compared to crystallization in large, open-system plutons (Fig. 3C).

Interestingly, the compositional range of the Jack Hills zircons extends to lower Ti (1-3 ppm) and/or higher Hf (1.2-1.4 wt%) than that of the TTG-hosted zircons (Fig. 4), which may be explained by the smaller amounts of silicic melt and closed-system conditions, and/or even lower initial Zr contents, in the Hadean. Either way, our results indicate that the trace element composition of Hadean detrital zircons is consistent with a TTG-like magmatic environment. Complex [START_REF] Whitehouse | A rare earth element study of complex zircons from early Arcaean Amîtsoq gneisses, Godthåbsfjord, south-west Greenland[END_REF][START_REF] Nutman | ~3,850 Ma tonalites in the Nuuk region, Greenland: geochemistryand their reworking within an Eoarchaean gneiss complex[END_REF][START_REF] Kielman | A tonalitic analogue to ancient detrital zircon[END_REF]; (4) Uivak Gneiss [START_REF] Vezinet | Hydrothermally-altered mafic crust as source for early Earth TTG: Pb/Hf/O isotope and trace element evidence in zircon from TTG of the Eoarchean Saglek Block[END_REF]; (5) Jiaobei Province [START_REF] Liu | The growth, reworking and metamorphism of early Precambrian crust in the Jiaobei terrane, the North China Craton: Constraints from U-Th-Pb and Lu-Hf isotopic systematics, and REE concentrations of zircon from Archean granitoid gneisses[END_REF][START_REF] Shan | Zircon U-Pb ages, geochemistry, and Nd-Hf isotopes of the TTG gneisses from the Jiaobei terrane: Implications for Neoarchean crustal evolution in the North China Craton[END_REF]; (6) North Pole Dome [START_REF] Kitajima | Relative retention of trace element and oxygen isotope ratios in zircon from Archean rhyolite, Panorama Formation, North Pole Dome[END_REF]; (7) Irkut [START_REF] Turkina | U-Pb (SHRIMP II), Lu-Hf isotope and trace element geochemistry of zircons from high-grade metamorphic rocks of the Irkut terrane, Sharyzhalgay Uplift: Implications for the Neoarchaean evolution of the Siberian Craton[END_REF]; (8) Superior Province [START_REF] Bouvier | Li isotopes and trace elements as a petrogenetic tracer in zircon: insights from Archean TTGs and sanukitoids[END_REF]; (9) Tarim [START_REF] Zong | The generation and evolution of Archean continental crust in the Dunhuang block, northeastern Tarim craton, northwestern China[END_REF]; (10) Jack Hills [START_REF] Cavosie | Correlated microanalysis of zircon: trace element, δ 18 O, and U-Th-Pb isotopic constraints on the igneous origin of complex >3900 Ma detrital grains[END_REF][START_REF] Bell | Eoarchean crustal evolution of the Jack Hills zircon source and loss of Hadean crust[END_REF][START_REF] Turner | An andesitic source for Jack Hills zircon supports onset of plate tectonics in the Hadean[END_REF]. B: Geological map of the Barberton Paleoarchean plutonicvolcanic suite, South Africa [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF]; BRVC = Buck Ridge Volcanic Complex. C: K2O vs. SiO2 composition of bulk rock samples from the Barberton suite, BRVC melt inclusions and tonalite liquid lines of descent (experimental from Marxer and Ulmer, 2019 and modelled here -see text); the gray field represents mixed compositions between melts and plagioclase (Pl) + amphibole (Am) cumulates (adapted from [START_REF] Laurent | Earth's earliest granitoids are crystal-rich magma reservoirs tapped by silicic eruptions[END_REF]. Trace element concentrations of zircon-forming melts (normalized to average TTG from [START_REF] Moyen | The composite Archaean grey gneisses: Petrological significance, and evidence for a non-unique tectonic setting for Archaean crustal growth[END_REF] in the Barberton samples obtained from inversion of zircon compositions, compared to those of the BRVC melt inclusions; the low-, medium-and high-pressure (LP-, MP-, HP-)TTG groups [START_REF] Moyen | The composite Archaean grey gneisses: Petrological significance, and evidence for a non-unique tectonic setting for Archaean crustal growth[END_REF]; and the melts saturating zircon in the crystallization model (see text). All Th(t) and U(t) concentrations (including normalization values) are corrected for radioactive decay. zircons from different geodynamic settings (data from [START_REF] Burnham | Formation of Hadean granites by melting of igneous crust[END_REF][START_REF] Carley | Iceland is not a magmatic analog for the Hadean: evidence from the zircon record[END_REF][START_REF] Grimes | On the occurrence, trace element geochemistry, and crystallization history of zircon from in situ ocean lithosphere[END_REF][START_REF] Karakas | The pace of crustal-scale magma accretion and differentiation beneath silicic caldera volcanoes[END_REF] as kernel density estimates contoured for 25, 50, 75 and 90% data; the arrows approximate the highest point density trends and correspond to the zircon compositional evolution throughout magma differentiation from mafic-intermediate to

Figure 1 .

 1 Figure 1. A: Location of zircon trace element data examined here: Barberton (this study); (1)

Figure 2 .

 2 Figure 2. A. Plot of the distribution of zircon Ti and Hf contents in the Paleoarchean Barberton

Figure 3 .

 3 Figure 3. Results of the tonalite crystallization and zircon saturation model. A. Temperature

Figure 4 .

 4 Figure 4. Titanium vs. Hf plots showing the composition of early Earth zircons and Phanerozoic

GSA Data Repository item

202Xxxx, containing details about analytical and modelling methods, data compilation, filtering, and zircon/melt partitioning, is available online at www.geosociety.org/pubs/ft20XX.htm.

ACKNOWLEDGMENTS

We thank M.P. Silva for help during data acquisition, and T. Johnson, A.P. Nutman and H.

Rollinson for constructive reviews. This study was funded by the Swiss National Science Foundation (grant 200021_178928 to OB) and the Australian Research Council (grant FL160100168 to JFM). This is International Research Project BuCoMo contribution A9.