Unravelling local environments in mixed TiO2–SiO2 thin films by XPS and ab initio calculations

Abstract Mixed TixSi1−xO2 oxide can exhibit a partial phase separation of the TiO2 and SiO2 phases at the atomic level. The quantification of TiO2–SiO2 mixing in the amorphous material is complicated and was so far done mostly by infrared spectroscopy. We developed a new approach to the fitting of X-ray photoelectron spectroscopy data for the quantification of partial phase separation in amorphous TixSi1−xO2 thin films deposited by plasma enhanced chemical vapour deposition. Several fitting constraints reducing the total number of degrees of freedom in the fits and thus the fit uncertainty were obtained by using core electron binding energies predicted by density functional theory calculations on TixSi1−xO2 amorphous supercells. Consequently, a decomposition of the O 1s peak into TiO2, SiO2 and mixed components was possible. The component areas ratios were compared with the ratios predicted by older theoretical models based on the atomic environment statistics and we also developed several new models corresponding to more realistic atomic structure and partial mixing. Based on the comparison we conclude that the studied films are mostly disordered, with only a moderate phase separation.


Introduction
The mixed TiO 2 -SiO 2 oxides have attracted considerable attention in the area of photocatalysis as it has been shown that they are more active than pure TiO 2 [1][2][3]. The mixed oxides serve as base components for complex nanomaterials with applications in photoluminescence [4], electrochemical sensors [5] or catalytic hydrogenation [6] dehydrogenation [7] and hydrodeoxygenation [8]. There are also multiple possible optical application for those coatings [9], with demonstrated use in waveguides [10,11], laser mirrors [12] and rugate filters [13].
The atomic structure of mixed TiO 2 -SiO 2 films has been studied quite extensively as it influences both the optical and catalytic activity. Often debated topic is the dominant coordination number for Ti and Si atoms. In the crystalline TiO 2 and SiO 2 , the Ti cation is present in the 6coordinated octahedral position and the Si cation in the 4-coordinated tetrahedral position, with the only exception being some high pressure polymorphs such as stishovite. The majority of reports on TiO 2 -SiO 2 shows that Si atoms are also almost always 4-coordinated [14], whereas for Ti the coordination depends on composition. Ti evolves from mostly 4-coordinated at low Ti concentration towards dominantly 6-coordinated at higher concentration -nevertheless, 5-coordinated Ti atoms were also reported [15][16][17][18][19]. There are however some contradictory reports mentioning dominantly 6-coordinated Ti atoms even at low TiO 2 concentrations [20,21]. It has been also established, both experimentally and theoretically, that there are no direct Ti-Ti, Si-Si or Ti-Si bonds, and the Ti and Si atoms are always interconnected through the bridging oxygen atoms [16,22,23].
An extensively studied topic for mixed TiO 2 -SiO 2 films is whether the TiO 2 and SiO 2 phases are well mixed at the atomic scale or whether they consist of more than one separated phases. Multiple reports show that TiO 2 -SiO 2 phase separation can be easily induced by annealing [24,25], especially for samples prepared by the sol-gel method and high TiO 2 concentrations [26][27][28]. However, it was also shown that it is possible to prepare homogeneous mixtures of the TiO 2 and SiO 2 phases in the whole composition range [22,17,29]. The thorough ab initio report on atomic structure of a-Ti x Si 1−x O 2 by Landmann et al. [16] reported no phase separation at x as high as 0.8, suggesting that the atomically mixed phases can be formed even at high Ti concentrations, the possibility of phase separation was not ruled out though and more recent ab initio reports agree that the TiO 2 and SiO 2 mixing (in bulk) is not energetically preferable [30,31].
Larouche et al. [22] proposed a simple statistical model to quantify the probability of occurrences of Si-O-Si, Ti-O-Si and Ti-O-Ti groups under the assumption of 4-coordinated Si and Ti atoms and 2-coordinated O atoms for the two cases of perfectly random mixing and separate phases. Busani et al. [29] calculated the same probabilities numerically using Monte Carlo simulations.
The most often used method for the determination of the mixing in TiO 2 -SiO 2 films is the infrared (IR) absorption spectroscopy [17,22]. Absorption around 950 cm −1 is attributed to the Ti-O-Si stretching vibrations [32] and by comparing the relative intensities of the Si-O-Si and Ti-O-Si vibration modes in the absorbance measurements with the statistical Larouche's model some effort was dedicated to quantify the mixing [17,22]. Unfortunately, there are some difficulties with this approach. It is unclear whether the Larouche's statistical model works well for high TiO 2 concentrations when the Ti atoms are mostly 6-coordinated and the O atoms mostly 3-coordinated. Moreover, the existing attempts [17,22] show only the evolution of the Si-O-Si and Ti-O-Si components and the Ti-O-Ti component was not quantified at all, possibly due to the overlap with rocking Si-O-Si mode, rather broad nature of the Ti-O-Ti peak, weak activity, and its position in the far IR range, partially extending beyond the measurement range. Further problem arises from the unknown infrared absorption cross section (oscillator strength) of the Ti-O-Si vibrations. While it can be estimated from the pure phases for the Ti-O-Ti and Si-O-Si modes, this is not possible for the Ti-O-Si mode and hence its normalization factor is unknown. Another issue is related to the large refractive index difference between TiO 2 and SiO 2 . Especially at high TiO 2 concentrations the effects such as multiple reflections in the films might be significant due to the high refractive index and hence a simple extraction of absorbance might not be enough to properly quantify the absorption, a proper fitting of the spectra might be needed. All those factors increase the uncertainty in the mixing determination.
X-ray photoelectron-spectroscopy (XPS) is a suitable method for obtaining the atomic composition of TiO 2 -SiO 2 films [22,33] and it can be also used to extract information about atomic environments; in this case to detect the presence of the mixed phase [25,[34][35][36]. The analysis of TiO 2 -SiO 2 mixing by a careful evaluation of the XPS spectra (specifically the O 1s and Ti/Si 2p) does not suffer from some of the shortcomings of the IR analysis. Specifically, we can easily observe all the three components (TiO 2 , mixed Ti-O-Si and SiO 2 ) and their intensity depends on the concentration only. However, while in IR spectra the peaks for components are well separated, this is not the case in XPS. A careful fitting of the XPS peaks for different compositions of TiO 2 -SiO 2 films is necessary because even the sub-peaks belonging to the same component (e. g. TiO 2 ) shift with the composition [35]. Another downside of XPS is that, contrary to the IR spectroscopy, it is just a surface technique. This introduces additional challenges; for example some reports show different atomic concentration at the surface and in the bulk [26], or surface contamination [25]. The work of Orignac et al. [36] is a rare attempt to quantitatively analyse the mixed TiO 2 -SiO 2 films (with x between 0.05 and 0.3), by fitting the XPS O 1s peak using the three components' contributions. The results for these sol-gel films were somewhat surprising because the TiO 2 component had a larger area than the mixed component in the whole deposited x range. This differs significantly from the ratios expected from the Larouche's model for random mixing, and suggests a significant phase separation. Several recent reports highlighted the possibility of using ab initio calculations in order to achieve a better understanding and help with the analysis of XPS measurements [37][38][39].
In this paper, we develop a more quantitative approach to the determination of the phase separation in mixed TiO 2 -SiO 2 films. At first, the XPS spectra of Ti x Si 1−x O 2 are modelled from the first principles using the density functional theory (DFT) in order to check the fitting assumptions (e. g. whether the different types of environments have indeed significant differences in mean energy) and suggest some approaches for the fitting procedure. In the following section, an experimental approach to the fitting of XPS data is introduced and discussed with application to a series of TiO 2 -SiO 2 films deposited by plasma enhanced chemical vapour deposition (PECVD). The last section discusses several possible generalizations of the Larouche's statistical model for the case of a more realistic structure and also for the case of a partial phase separation.

DFT calculations
Structural models of amorphous Ti x Si 1−x O 2 were produced for seven different mixture ratios: x = 0, 0.1875, 0.34375, 0.5, 0.65625, 0.8125, and 1. The size of the amorphous cell was 96 atoms, with 64 oxygen atoms and the remaining 32 atoms were titanium or silicon, depending on the particular composition (0, 6, 11, 16, 21, and 26 Ti atoms, respectively). The structures were initiated with mass densities linearly decreasing from 3.8 g/cm 3 for TiO 2 to 2.2 g/cm 3 for SiO 2 and were generated using the ab initio simulated annealing approach described elsewhere [30] followed by a structural relaxation at 0 K. These steps were done using Vienna Ab initio Simulation Package (VASP) [40,41] with GGA-PBE [42] PAW-enabled pseudopotentials [43], a plane wave cut-off energy of 400 eV and a × × 3 3 3 k-mesh. The resulting models of amorphous structures were used to evaluate the core electron binding energies (BE) for the O 1s, Ti 2p and Si 2p core levels. The linearized augmented plane wave method (LAPW) as implemented in the Wien2k density functional theory package [44] was employed. We used a standard PBE exchange-correlation functional [42], with a × × 3 3 3 k-grid and the R K · mt max matrix size parameter was 7.5. The Slater transition state method was used to obtain the core electron binding energies. In this method, half of an electron is removed from the corresponding core state (compensated by a background charge) and the binding energy is calculated as a difference between the Fermi energy and the eigenvalue of this state [45]. We note that our calculations do not produce energies directly comparable to the experimental ones. The Fermi level is always at the valence band maximum (VBM) in the calculations due to the perfect stoichiometry of the simulated cells. This is not the case in experiment, where, for example, the Fermi level for TiO 2 was reported to lie just below the conduction band minimum due to n-doping [46,47]. Hence all the calculated BEs are aligned to the VBM and we denote them the modified binding energies (MBEs) in the following text to make this distinction clear. More details about the ab initio methodology can be found in the Supplemental Material.
The chemical environment corresponding to each atom was estimated by the coordination number analysis. To obtain the coordination number for each atom, a fixed radius r min was selected using the values of the first minimum in the pair distribution function [16,30]. The coordination number (and the specific environment) were obtained as the number and types of neighbour atoms inside this radius.

Preparation of Ti x Si 1−x O 2 Films
Thin Ti x Si 1−x O 2 films were prepared by PECVD in a low pressure plasma reactor consisting of an inductively coupled plasma (ICP) source (13.56 MHz) and a diffusion chamber with the substrate holder, mounted below the source. The reactor details are described by Granier et al. [48] and the PECVD of Ti x Si 1−x O 2 films is discussed in the previous paper [33]. The deposition proceeded from hexamethyldisiloxane (C 6 H 18 OSi 2 , HMDSO, 98%, Sigma Aldrich Ltd.) and titanium isopropoxide (Ti(-OC 3 H 7 ) 4 , TTIP, 99.999%, Sigma Aldrich Ltd.) vapours mixed with oxygen. Oxygen was introduced at the top of the ICP source, while both precursors were injected into the diffusion chamber via a dispersal ring, 15 cm in diameter, located 8 cm above the substrate. The TTIP vapours were kept in a heated container (60°C) and introduced through a heated line (63°C) with an MKS mass flowmeter (80°C). Another MKS mass flowmeter was used to control the amount of the HMDSO vapours. HMDSO was kept in a container heated to 35°C and the vapours were injected into the diffusion chamber through a line heated to 40°C. Both precursors were introduced into the chamber at the same time. The film composition was varied by changing the flow rates of HMDSO and TTIP (Table 1) whereas the flow rate of O 2 was constant, 16.4 sccm. The deposition times were selected to keep the film thickness around 250 nm. Crystalline (1 0 0)-oriented silicon wafers with thickness of (381 ± 25) µm were used as substrates for the deposited films. The substrate holder was at the floating potential.

Characterization of films
XPS spectra were recorded using a Kratos Axis Ultra operating at 1486.6 eV (Al source). The pass energy to measure the region of interest was fixed to 20 eV allowing an resolution in energy equal to 0.1 eV. Due to the insulating character of the samples a charge neutralization was used when acquiring the spectra. The C 1s signal from adventitious carbon was used afterwards for energy calibration, setting its position at 284.7 eV.
The atomic composition determined from XPS measurements was used to calculated the titanium-to-silicon ratio, i. e. the titanium stoichiometry x in the Ti x Si 1−x O 2 films, for varied alkoxide flow rates. The results are summarized in Table 1 and revealed that it was possible to prepare films with a controlled x from pure SiO 2 till TiO 2 .
After removal of a Tougaard-type background, an approximate Voigt lineshape (70% Gaussian, 30% Lorentzian) in the product form [49] was used for peak fitting. The fitting of highly-resolved atomic signals of Ti 2p, Si 2p and O 1s involved a multisample fitting procedure combined with some constrains of the fitting parameters. This advanced procedure was developed to unravel the local chemical environment of Ti x Si 1−x O 2 material prepared by PECVD and it is described in a later section together with its results and implications.
The films were previously characterized by optical spectroscopic methods [30] and their basic characterization concerning morphology, crystallinity and chemical structure was carried out by scanning electron microscopy, X-ray diffraction and infrared spectroscopy, respectively [33]. The TiO 2 film has a polycrystalline anatase structure, all samples containing Si are amorphous. The film density, given in Table 1, was determined from the film thickness measured by ellipsometry and the film mass.

Ab initio modelling
Ab initio modelling techniques were employed to test if a decomposition of the O 1s peak into just the TiO 2 , mixed and SiO 2 parts in a real three dimensional structure can be justified and to reveal any trends which could be used to reduce the uncertainty during the experimental fits.
As can be seen in Fig. 1, the DFT results are in a good qualitative agreement with the normalized experimental spectra. Although the qualitative agreement is excellent, a more quantitative test is needed to verify the aforementioned assumptions. With the information from analysis of coordination numbers and first neighbours we can easily sort the binding energies into components corresponding to the various atomic environments. The amorphous structure, represented by our model at least to the degree possible by the limited supercell size, contains a variety of local environments. An oxygen atom can have two, three or four neighbours with different number of Ti or Si atoms -see Fig. 2a for the environments found in structures with = x 0.5. This finegrained decomposition of the spectra is interesting from the fundamental point of view, but it has to be simplified for the experimental data evaluation. Fitting this many overlapping, poorly resolved peaks would lead to highly correlated components with huge uncertaintiesessentially, the results would be random and hence useless.
Therefore, we simply group the different atomic environments into the three simplified components. The Ti 2 , Ti 3 and Ti 4 local O environments into the TiO 2 component, the TiSi, TiSi 2 , SiTi 2 environments into the mixed Ti-O-Si component and the Si 2 environment into the SiO 2 component (subscripts denote the total numbers of neighbours of each Table 1 The list of Ti x Si 1−x O 2 films prepared at varied flow rates of HMDSO and TTIP, Q HMDSO and Q TTIP . The titanium stoichiometry x in the films was determined by XPS and the density was determined from the film thickness measured by ellipsometry and film mass. kind; not all are always present due to limited statistics of our rather small supercells). Fig. 2b shows the O 1s peak for one selected concentration = x 0.5 and highlights that the initial assumption, i.e. that overall peak shape can be fitted by three separate components, seems reasonable. Now we take a look at the structural models and try to explain some of the differences between the calculated and experimental data. While a perfect agreement cannot be expected due to the different mixing in the calculations and the experiment, there are a few noteworthy features.
The already mentioned absolute difference in the energies is unfortunately a problem of the DFT method (or the Slater transition state approach), and is also highly dependent on a particular choice of the xc functional. When replacing the GGA-PBE xc potential with local density approximation (LDA), the absolute values shift by~4 eV; however, the relative shifts stay the same.
It can also be noticed that the calculated broadened MBEs do not shift as much with compositional changes as the experimental ones. This is due to the fundamental differences between the calculated and experimental values, specifically as the experimental values are given with respect to the Fermi level, a shift of Fermi level with respect to VBM between TiO 2 and SiO 2 can cause this discrepancy.
Due to the perfect composition and stoichiometry of the calculated structures, the spectral features associated with the OH-groups such as the shoulder around 531.5 eV are missing. The presence of the OH peak, which is hidden in the pure SiO 2 such as seen in [25], can also explain why the experimental peak of pure SiO 2 is broader than of TiO 2 , although they seem quite similar in the calculated spectra.
To do a more thorough analysis of the fitting constraints, mean peak positions for the components of the O 1s peak and also for the Ti 2p and Si 2p peaks were calculated as a mean value of the individual binding energies. The results are shown in Fig. 3a, c and     Full widths at half maximum (FWHMs) of the components were obtained as the standard deviation of the MBEs and are shown in Fig. 3b. The FWHM of the mixed peak is bigger than that of the pure components due to the more diverse specific environment. Obviously, there are no clear trends in the peak widths. Due to the nature of DFT calculations, this is only the part of the peak broadening associated with the structural disorder as other contributions, such as lifetime broadening, instrumental broadening or temperature (phonon) broadening, are not easily accessible by ab initio methods (at least not by the Slater transition state approach). The question is whether the structural disorder is comparable, larger or smaller as the rest of the contributions. In the first approximation, the FWHM B can be estimated as Surprisingly, the trends are reversed here, but it must be noted that the experimental TiO 2 is actually polycrystalline anatase and hence its FWHM is expected to be somewhat narrower than in the amorphous material. On the contrary, the Si-O-Si part in SiO 2 can not be easily distinguished from the impurity component, so the resulting peak is broader. Additionally, there can also be small differences due to different conductivity of TiO 2 and SiO 2 . Using those values, nevertheless, the total broadening from all the other contributions has FWHM in the range between 0.9 eV and 1.3 eV. When compared with the calculated values in Fig. 3b, it can be seen seen that the disorder broadening is the most significant for the mixed component where it is the dominant broadening factor, while it is less significant for the TiO 2 and SiO 2 components.

XPS data analysis
As indicated by the DFT results, it is reasonable to fit the O 1s peak with the three (Gaussian-Lorentzian) components in order to compare their ratios and hence to deduce the overall mixing. Unfortunately, while the mean position of the components are in fact separated, the energy differences are not very big, e. g. as small as 0.6 eV for the TiO 2 and mixed component and the components are not resolved well enough. Moreover, as can be seen from peak shapes in Fig. 1, the O 1s peak can be fitted with only two components for many compositions (if the areas, positions and FWHMs are left free). The situation is also further complicated by the oxygen engagement in other bonds (e. g. OH).
To reduce the uncertainty, the degree of freedom during fitting must be reduced. The areas under the peaks need to stay as free parameters in order for the method to work. It is also not possible to fix the energy of the components completely even when the MBE seems constant (as for example the mixed components of the O 1s peak), since the VBM to E F distance changes significantly between TiO 2 and SiO 2 . We must therefore utilize relative shifts. The first choice is to fix the relative shift between the SiO 2 component of the O 1s peak and Si 2p peak, which is quite constant as seen in Fig. 4b. It would be convenient to fix also the same for the TiO 2 component of the O 1s peak and Ti 2p 3/2 peak, but although the difference (as shown in Fig. 4a) is constant for x 0.5, it decreases significantly for < x 0.5. This originates form the downshift of MBE of the TiO 2 O 1s component as seen in the Fig. 3a. However, besides the large uncertainties, this is likely an artefact of the DFT calculations as already discussed in the previous section, since this corresponds to the right shoulder in the spectra and this feature is seen neither in our experimental data nor in any other reported spectra [22,25,50]. Hence it seems reasonable to neglect those two points and assume the distance between the TiO 2 component of the O 1s peak and Ti 2p 3/2 peak to be constant as well. The same can be then reasoned for the energy difference between TiO 2 and mixed components of the O 1s peak. This three fixed energy shifts are the basis of our fitting approach. The shift between O 1s TiO 2 and mixed components was 0.65 eV, obtained as average value from the DFT data. The shift between the TiO 2 component of the O 1s peak and Ti 2p 3/2 peak was 71.26 eV, and the shift between the SiO 2 component of the O 1s peak and Si 2p peak was 429.28 eV. These values were extracted from the pure TiO 2 and SiO 2 experimental samples.
Our preliminary fits also indicated, that it is necessary to impose some other constraints regarding the FWHMs. For some fitting parameters, when left completely free, non-physically small values were obtained together with unrealistic fluctuations between the compositions. They cannot be completely fixed either, since such fixing leads to bad fits. The calculated disorder-dependent contribution to the FWHMs are smaller then the rest of contributions to the peak FWHM (as estimated in previous section) for the TiO 2 and SiO 2 parts, and comparable for the mixed part. We can therefore choose one composition-independent parameter expressing the disorder part B i of FWHM of each component i and one sample-dependent parameter B j , where the total FWHM for the given component and composition (sample) is calculated as This allows for some freedom compensating the fluctuations and at the same time it reduces the number of free FWHM-related parameters in the multisample fit of n samples from n 3 to + n 3. In order to facilitate compositionally independent (global) parameters in the fits, all the spectra were fitted simultaneously.
The last difficulty is a proper treatment of oxygen atoms bonded to atoms other than Ti and Si, such as O in OH and CO groups, which we will together denote impurities. They are even worse resolved than the three main components, and fitting them independently was out of the question. Instead, we linked their positions and areas to the main components. The peak shapes of the pure TiO 2 and SiO 2 samples were decomposed into the MO 2 (M = Ti or Si) and impurity components ( Fig. S1 in Supplemental Material). Area ratios and relative energy shifts between MO 2 and impurity components from these fits were fixed in the subsequent multisample fit. Both impurity components had one common global disorder-dependent FWHM parameter, similar as for the other three components, and followed the same scheme for calculation of the total FWHM (Eq. 1).
This approach was at first tested for a reference bi-layer sample constituted of SiO 2 (2 nm) on TiO 2 (20 nm). A very good agreement between the experimental and calculated spectra was obtained using only the TiO 2 and SiO 2 components as expected (Fig. S2 in Supplemental material). It indicated that the chosen fitting procedure works very well for the sample containing simultaneously both SiO 2 and TiO 2 materials. The same approach fails for the PECVD Ti x Si 1−x O 2 = x 0.57 sample as shown in the same figure, hence indicating that a reasonable fit without the mixed component is impossible in such case.
Example fit of the O 1s peak from the aforementioned multisample fitting approach including all the constraints is shown in Fig. 5 and displays a nice agreement with the measured peak shape. The full multisample fit including the Ti 2p and Si 2p peaks is shown in Fig. S3 in Supplemental Material. Fig. 6 shows the evolution of the peak positions for the O 1s components as well as for the Ti 2p 3/2 and Si 2p peaks. The positions decrease significantly in the SiO 2 -rich region and are mostly constant for the x 0.3. Such significant changes are not seen in the calculated MBEs and hence they most likely originate in the shift of the Fermi level with respect to the VBM. This interpretation is further supported by the fact that the shifts are the same for O 1s, Ti 2p and Si 2p peaks, and, as will be shown in the next section, especially the Si environment is almost constant over the whole composition range and therefore no shifts due to local environment changes could be expected. In fact the nonlinear trend is very similar to the trends in the evolution of the band gap [30] and hence we can conclude that the Fermi level is just below the minimum of the conduction band for majority of the compositions; only for the SiO 2 -rich compositions it shifts towards the middle of the band gap as expected for SiO 2 .
It is plausible that with a high resolution measurements of the valence bands it would be possible to align the spectra according to VBM and then fix the component position to absolute energy values to simplify the fitting procedure.
The proportions of the O 1s components as a function of the Ti metal fraction, x, in the films are plotted in Fig. 7. The percentage in oxygen atoms engaged in a TiO 2 -like framework increases with x whereas those belonging to a SiO 2 -like environment decreases as expected. Surprisingly, the fraction of the O 1s mixed component is significantly smaller than predicted by the Larouche's mixing model [22] and shows some asymmetry. The implications are thoroughly discussed in the following section.

Phase separation modelling
Larouche et al. [22] proposed a simple statistical model for the O 1s components. Under the assumption of 2-coordinated O atoms, 4-coordinated Ti and Si atoms, and a perfect random mixing, the expected fractions are x x x , 2 (1 ) 2 and x (1 ) 2 for the TiO 2 , mixed and SiO 2 components, respectively. For completely separated phases, the fractions are the straight lines x and x 1 , with the mixed component absent.
The experimentally determined fractions of the three components are plotted as functions of composition x in Fig. 8 alongside the model curves for phase separation and random mixing. It is evident that real PECVD film structures are far from either of these theoretical border cases. Interpretation of the results thus requires a model incorporating a partial phase separation.
A complete model would include not only the evolution of the atomic network topology with changing x, as studied by Landmann et al. [16], but also its realistic dependence on phase separation, and possibly additional parameters and peculiarities. This is beyond the scope of this work. Nevertheless, even relatively simple approaches allow progressing beyond the Larouche's model. Two important aspects to consider are an improvement of statistics of O environments based on ab initio calculations and a physical order-disorder model describing partially separated phases.

Oxygen environment statistics
The assumptions of Larouche's model are actually not valid in the Ti x Si 1−x O 2 case [16], as also illustrated in Fig. 9. While Si is approximately 4-coordinated over the majority of compositional range, there is a significant amount of 5 and 6-coordinated Ti atoms and 3-coordinated O atoms. It was already argued [17] that this can be the cause of differences between the model and measured components. Is this explanation supported by our calculations?
We start by analysing the effect of changes in the average coordination numbers in the case of random mixing. Two factors come into play here. First, in the Larouche's model the probabilities that a   , respectively. However, as Ti atoms have larger average coordination number, O atoms tend to have more Ti neighbours than Si neighbours. Hence, looking at a single random M-O bond, the probability of finding Ti-O is larger than for Si-O in more than half of the composition range. Here C t denotes the average coordination number of the atom = t Ti, Si and O. Note that the denominators can be simplified using the relation which follows from the stoichiometry. Second, a 3-coordinated oxygen has lower probabilities of pure environments than 2-coordinated simply due to having more neighbours.
The fraction of a pure component MO 2 is , where the f O2 and f O3 are the fractions of 2 and 3-coordinated O atoms, respectively. If we disregard 4-coordinated O, which occurs only rarely, they can be written simply as Si and mixed for random mixing. The average coordination numbers can be taken from our DFT structural models or the extensive Landmann's review (Fig. 9). Looking at the Landmann's data, C Ti is almost perfectly linear Si is close to 4 in the majority or the compositional range and can be approximated by + C x 4 0.5 Si 2 , and C O is then obtained from the relation (4). Our DFT data are less smooth than Landmann's due to the smaller size of the cells. Nevertheless, the trends are similar and the main difference is that all coordination numbers are somewhat smaller due to a lower mass density [30]. The average coordination number for Ti ends at 5.6 and the highest average coordination number for Si is 4.16 at = x 0.8125. Hence, Si is essentially 4-coordinated in the entire range.
The result of feeding our coordination number trends into the model can be seen in Fig. 10. The final curves are surprisingly similar to the Larouche's model. A decrease in the SiO 2 component and corresponding increase of the mixed component are visible, whereas the TiO 2 component is largely unchanged. We note that the differences increase when Landmann's data are used instead, but only moderately. The overall shapes match nicely the data points from the DFT model, confirming that the DFT cells are well mixed. Indeed, it would be surprising to see phase separation in the DFT models because the cooling was quite fast. There are small discrepancies in the ratios of the 2 and 3coordinated parts of the TiO 2 and mixed component. However, they can be explained by either the small cell size or not taking into account in the model that some groups might be more energetically preferable. We can conclude that the simple Larouche's statistical model is a good approximation even though its initial assumptions are not justified.
Since the DFT cells correspond well to random mixing, they can serve directly as models of perfectly mixed structures. Instead of utilizing just the average coordination numbers, we can incorporate the complete structural data by taking the calculated fractions f x ( ) t of all possible O environments types = t Si 2 , TiSi, Ti 2 , Ti 3 , Ti 2 Si, … found in the simulated structures. They can be classified further to individual coordination polyhedra types [16]; nevertheless, the chemical identity of the neighbours is sufficient here.

Ti
Si and mixed are then obtained by sorting the environments t into three groups, Ti -only Ti neighbours, Si -only Si, and mixed -all combinations of different metal atoms. The individual   [30], empty points are from DFT calculations by Landmann et al. [16].
contributions in each group are then summed, for instance Ti  Ti  Ti  Ti   2  3 4 for the titania component. Of course, in this case the resulting dependencies are very similar to the approximate model because the contribution of disregarded environments is negligible.
As a final remark we note that the presented splitting of O environments into simple pure and mixed components is unambiguous for 2-coordinated O, but for higher coordination numbers another scheme is possible. It considers all M-O-M configurations that can be found in the environment of a specific O atom and weights them accordingly. For instance, an O atom in Ti 2 Si environment would contribute 2/3 to the mixed component and 1/3 to the TiO 2 component. One serious disadvantage of this scheme is that the components are then combinations of more individual environment types. This results in rather odd component line shapes if we plot them as combinations of individual peaks as obtained from DFT, bringing additional uncertainty to the peak fitting constraints. Nevertheless the results of one approach can be in principle recalculated to the other using the nearest neighbour statistics.

Simple partial phase separation models
The simplest partial phase separation model can be constructed by considering the material as a combination of mixed and separated materials. For each environment, the compositional dependence for separated phases (chemically ordered structures) is the straight line constructed from the fractions in pure materials. For a material composed of mixed and separated phases with fractions q and q 1 , the dependence can be expressed Formula (5) is not a realistic description of an order-disorder evolution. However, it allows modelling and quantification of partial phase separation using simulation results without having to produce structures covering the entire two-dimensional space of composition x and disorder q.
It is not necessary to use the detailed O environment data f t as the input. The average coordination numbers can also serve as the starting point, as well as the Larouche's model. In the latter case there are only three environments, Ti 2 , Si 2 and TiSi, with parabolic dependencies providing a simple analytical model for partial phase separation. Expressions (6) can be also derived in an alternative manner which leads to an interesting interpretation of the parameter q. We imagine traversing -M-O-M-O-M-chains in the atomic network. In a random mixture the metal atoms M we encounter are randomly Ti or Si, while for separate phases we see long chains of one type of metal atom. The sequences of encountered atoms can be constructed as Markov chains [51]. If we propose that the probability of switching between Ti and Si is reduced by a factor q (0, 1] with respect to random mixture (independently on composition) then the transition matrix is The limit values = q 1 and q 0 correspond to perfect random mixture and complete phase separation, respectively. In other words, chain segments or 'clusters' composed of pure phases are, on average, q 1/ larger than in the random mixture. The transition matrix (7) gives the stationary distributions of O environments (6).
By fitting this simplest model to the experimental dependencies we already improve the agreement, as illustrated in Fig. 8 (curve Larouche/ Markov). The expressions (6) are symmetrical with respect to swapping the metals and simultaneously x with x 1 . Hence, the model is not capable of capturing the evident asymmetry in experimental XPS data. The asymmetry is probably in part related to the larger average coordination numbers of O in the Ti-rich region compared to the Si-rich region. Nonetheless, even curves calculated from the ab initio structural data exhibit only a slight asymmetry, so this is unlikely the full explanation; see the fit with a model obtained by analysing oxygen neighbours in DFT cells in Fig. 11 (curve DFT neighbours). The theoretical compositional dependence is asymmetrical, but only very slightly.

Statistical model
For a more physical approach to phase separation, we return to regular atomic networks with 2-coordinated O and describe a model based on the concept of quenched disorder in a regular solid solution. This type of model is known under many names, such as Ising or lattice gas model [52][53][54][55][56][57]. While separated phases are energetically preferred, the non-equilibrium deposition process produces an atomic network with a relatively large disorder (large mixing), corresponding to a high temperature. The network is frozen in this state as the energy barriers impeding further relaxation are too high compared to k T B at room temperature T (k B denotes the Boltzmann constant). The disorder then depends on = E k T / B , where T is the effective temperature and E is the interchange energy which is related to the energy of formation [30] E f of a perfect random mixture as follows (see Supplemental Information): An explicit solution is known for one-dimensional (1D) chains that are equivalent to the 1D Ising model [55], leading to expressions with the same overall form as (6). However, q is no longer a constant; it is now a function of both and x (see Supplemental Information): As a result, the mixed curve is flatter around = x 0.5 compared to the parabolic dependencies obtained with constant q. The fit of the experimental dependencies with this model is also illustrated in Fig. 11.
The 1D model can be considered realistic when the material is actually formed mostly by linear chains or the phase separation is relatively low (small ) and the topology of the atomic network plays little role because there is no long-range order yet. For larger , the 1D model exhibits a gradual progression to separated phases. This is typical for such 1D systems with short-range interactions and generally differs from the behaviour of higher-dimensional systems where a phase transition between disordered and phase-separated states is observed [55].
To obtain x ( , ) dependencies that capture this phenomenon at least qualitatively, we performed a simple numerical simulation of twodimensional (2D) atomic networks using the Metropolis-Hastings algorithm [58][59][60]. For simplicity, a square geometry in which both metals are 4-coordinated was used. By annealing atomic networks of varying composition x to different finite temperatures, corresponding to different values, we sampled the function in a mesh of x ( , ) points. Values for arbitrary x and could then obtained by interpolation, permitting fitting the experimental data with the 2D model.
The resulting theoretical mixed curve plotted in Fig. 11 is even flatter in the centre than for the 1D model. This effect is even more marked for smaller degree of disorder, when the function becomes almost constant over a wide range of compositions (not shown). We can expect that for 3D structures the dependence would have the same overall character because the Ising model behaviour changes qualitatively only between 1D and higher dimensions, despite the various effects influencing it in detail: different topology of the atomic network, including mean coordination numbers and ring size (that are on average larger and smaller, respectively, in realistic atomic networks [16]), 3D structures being generally closer to a mean field model than for 2D [55], etc.
For both the 1D and 2D models, the most notable discrepancy is nevertheless still the asymmetry, which they cannot describe because they are inherently symmetrical with respect to the two metals. We also note that even though the theoretical mixed curves have dissimilar shapes for the individual models, the sums of squared residuals differ by about 10% between all the various fits. Hence, the differences are not significant.
Considering how the model was constructed, one is of course curious what is the thermodynamic temperature T corresponding to the observed disorder. An estimate of E is necessary, which can be obtained by fitting the formation energy dependencies on composition [30] by formula (9). For the amorphous phase, which is probably the most relevant, we obtain E 0.17 eV, giving T 2400 K. Although this temperature may seem high at the first sight, it is consistent with similar approaches also considering only configurational entropy. For example, the maximum temperature for spinodal decomposition in Ti Al x x 1 N system was theoretically predicted to be 5000-8000 K [61,62] while experimentally the decomposition starts at around 1300 K [63]. Later theoretical reports suggested that this discrepancy is mostly due to neglecting other entropic terms (e. g. related to lattice vibration or anharmonic effects [64]). Consequently, also the here predicted temperature for phase separation is actually expected to be significantly lower in reality.
It should be also noted that the phase separation transition occurs around of unity, depending on the dimension and composition. The fitted value is somewhat smaller, but relatively close. Disorder close to critical may arise naturally when the deposition process only heats the material locally and the energy dissipates as is the case of PECVD, because decreasing temperature also means falling probabilities of overcoming potential barriers and the progress of phase separation beyond the critical point is very slow in such case.

Conclusions
We quantified from XPS measurements the phase separation in 10 amorphous thin Ti x Si 1−x O 2 films prepared by PECVD for a series of compositions x covering the entire interval x [0, 1]. The analysis was based on structural models constructed using ab initio DFT calculations, which provided O 1s, Ti 2p and Si 2p binding energies of each individual atom in simulated amorphous cells. The DFT results showed that the O 1s peak can be decomposed to three simple components, corresponding to the TiO 2 , SiO 2 and mixed environment of the O atoms. The necessity for a mixed component was further demonstrated by the impossibility to fit O 1s peaks by just the two peak shapes obtained for pure materials. The DFT results also justified several constraints for BEs and FWHMs of the components that were utilized in the XPS data fitting, in particular constant relative shifts between pairs of components, O 1s TiO 2 and mixed (0.65 eV), O 1s SiO 2 and Si 2p (429.28 eV), and O 1s TiO 2 and Ti 2p 3/2 (71.26 eV). Using a multisample analysis of all samples, we then obtained dependencies of the TiO 2 , SiO 2 and mixed components on x. The dependencies did not lie close to either of the previously published simple models of random mixing and phase separation. They also exhibited a notable asymmetry, with the mixed component being more prevalent in Ti-rich materials. Since it had been previously suggested that the asymmetry can arise from changes in average coordination numbers with x, we utilized the O environment statistics from ab initio results to construct more realistic models. However, the resulting theoretical curves were still only slightly asymmetrical, and we have to conclude the changes in coordination number alone do not explain the observed asymmetry. Finally, we presented several models suitable for fitting the functional dependencies of partial phase separation obtained from XPS data. The models range from a simple Markov chain-like model, which does not require any detailed information about the material, to models utilizing the coordination number and environment statistics and also a statistical phase separation model. They all agree that the studied PECVD are mostly disordered, with only a moderate phase separation tendency. Since modelling and methods form a large part of this work, we should also comment on their applicability to other materials and analyses. Although concrete DFT calculations need to be done anew in each case, the overall procedure is quite general. We utilize theoretical BEs from DFT to derive constraints which then provide guidance for XPS spectra fitting. Importantly, this approach remains practical even when DFT modelling of spectra is not feasible, such as for complex amorphous materials. The presented phase separation models can be adopted more directly. Individual models vary in complexity, assumptions and required information about the material -and this decides their suitability in specific cases. Nonetheless, all enable the characterization of disorder by combining XPS analyses for a range of material compositions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.