Thermodynamic Stability and Structural Insights for CH3NH3Pb1−xSixI3, CH3NH3Pb1−xGexI3, and CH3NH3Pb1−xSnxI3 Hybrid Perovskite Alloys: A Statistical Approach from First Principles Calculations

The recent reaching of 20% of conversion efficiency by solar cells based on metal hybrid perovskites (MHP), e.g., the methylammonium (MA) lead iodide, CH3NH3PbI3 (MAPbI3), has excited the scientific community devoted to the photovoltaic materials. However, the toxicity of Pb is a hindrance for large scale commercial of MHP and motivates the search of another congener eco-friendly metal. Here, we employed first-principles calculations via density functional theory combined with the generalized quasichemical approximation to investigate the structural, thermodynamic, and ordering properties of MAPb1−xSixI3, MAPb1−xGexI3, and MAPb1−xSnxI3 alloys as pseudo-cubic structures. The inclusion of a smaller second metal, as Si and Ge, strongly affects the structural properties, reducing the cavity volume occupied by the organic cation and limitating the free orientation under high temperature effects. Unstable and metaestable phases are observed at room temperature for MAPb1−xSixI3, whereas MAPb1−xGexI3 is energetically favored for Pb-rich in ordered phases even at very low temperatures. Conversely, the high miscibility of Pb and Sn into MAPb1−xSnxI3 yields an alloy energetically favored as a pseudo-cubic random alloy with tunable properties at room temperature.

Justified by the imminent scarcity of energy sources based on conventional fossil fuels, the recent rise of metal halide perovskites (MHP defined by ABX 3 ) as alternative of low cost photovoltaic material has excited the community centered around silicon, which has been considered the principal element in solar cells [1][2][3][4][5] . MHP based on the use of lead iodide ( − PbI 3 ) and methylammonium ( = + + CH NH MA 3 3 ), i.e., MAPbI 3 6-8 , reached remarkable 20% 9 of efficiency in lighting conversion devices, which has put it as background for improvements of its photovoltaic performance [10][11][12][13][14] . However, a deeper comprehension of the chemical and structural properties correlated with the optical efficiency is needed. Additionally, for a large scale commercialization of solar cells based on MHP, combining thermodynamic stability and high photovoltaic performance is the key point for the viability of those devices [15][16][17] .
Experiments have revealed the MAPbI 3 stability in different structural motifs into a relative short range of temperatures. For example, below 163 K the orthorhombic (Amm2 space group, = .
a 8 87, = . b 12 67), and above 328 K 19 MAPbI 3 has been suggested as pseudo-cubic (P4mm space group, = . a 6 31) 20 . Additionally, the thermodynamic stability of MHP has been investigated aiming their obstacles against the inclement weather, such as UV, moisture, heat, and oxygen, which is crucial for MHP durability of photovoltaic cells [21][22][23][24] . Experiments of differential thermal analysis has indicated the decomposition of MAPbI 3 tetragonal phase in CH 3 NH 3 PbI 3 (s) → PbI 2 (s) + CH 3 NH 2 (g) + HI(g), in order that for temperatures from 403 K the perovskite gradually starts to be decomposed 25 . It is reported through X-ray diffraction that even after the MAPbI 3 systems be submitted under temperature of 443 K the sample keeps as 69% of MAPbI 3 and 31% of PbI 2 25 . For MASnI 3 , for instance, X-ray diffraction experiments revealed the presence of tetragonal structure at 423 K, and at room temperature by considering an MASn x I 3 for . ≤ ≤ .
x 0 9 1 4 as relative quantities between MA:Sn 2+ (in 1:x) used throughout the synthesis process, the perovskite adopts a pseudo-cubic structure for some x values 26 . However, the thermal decomposition starts only at 473 K, which is a higher than for MAPbI 3 .
From the last years the mixtures (alloys) MAPbI 3 -based perovskites has provided a new perspective to stabilize and tune MHP properties from their composition through several different ways, such as: (i) changing the MA + organic cation by another keeping the charge balance [27][28][29] ; (ii) replacing Pb 2+ atoms by another cation, e.g., Sn 2+ or Ge 2+30-33 ; or (iii) varying the halogen [34][35][36] . This approach brought a tremendous progress in the development of MAPbI 3 -based for photovoltaic devices, especially for the MAPb 1−x B x I 3 alloys, from which the toxicity of Pb can be suppressed through the use of another congener eco-friendly metal (e.g. B = Sn or Ge) 37,38 . Based on that, those MHP alloys open an enhancement field for the photovoltaic performance by chemical control of the thermodynamic stability and optical properties 20,[39][40][41] . Even though MASnI 3 has been investigated as an alternative for lead-free perovskite, its low power conversion efficiency 31 and low oxidation resistence 30 are some motivatory hindrances to workaround through the use of alloys. For instance, the MAPb x Sn 1−x I 3 stable alloy was recently investigated by Hao et al. 31 , who showed experimentally the control of the band gap of the MAPbI 3 (1.55 eV) for compositions towards MASnI 3 pure (1.30 eV), as the lower band gaps in 1.17 eV and 1.24 eV for MAPb 0.5 Sn 0.5 I 3 and MAPb 0.75 Sn 0.25 I 3 . Furthermore, the study revealed that the MAPb 0.5 Sn 0.5 I 3 alloy adopts a pseudo-cubic structure, while in so far as the content of Pb increases the structure adopts a tetragonal configuration, i.e., gradually reaching the stable phase of MAPbI 3 at room temperature. In others studies focused on the optical and eletrochemical properties 30,42 , it was found an increase for the incident photon wavelength for MAPb 0.5 Sn 0.5 I 3 , which was red-shifted to 1060 nm, corresponding to the 260 nm displacement with respect to the MAPbI 3 pure. Beyond that, since a large band gap of 1.90 eV for MAGeI 3 has been found 43 , the MAPb x Ge 1−x I 3 alloy as tetragonal structure has also been investigated through theoretical calculations 44 . The alloys presented narrower band gaps than their pure perovskites counterparts, so that MAGe 0.75 Pb 0.25 I 3 composition has presented the highest theoretical efficiency of about 24%. However, this study is restricted to few configurations and a deeper understanding of the structural stability is still needed.
As first attempt to determine the stability of a hibrid perovskite from a specific composition, the Goldschmidt's tolerance factor (t) is a geometric parameter initially used to predict the ability to form a 3D perovskite 45 , which empirically lie into . < < . t 0 80 11 range 7,46 . The t is part of an empirical relation given by where R A is the effective radii of organic cation, R B the radii for bivalent metal cation, R X for halide anion. However, the Goldschmidt's tolerance factor is limited to predict the perovskite alloys formation, since that parameters as the miscibility between the different metals involved within the crystal, i.e., concerning the octahedral inner sites occupied by Pb or a second metal B, as well as temperature relative to the thermodynamic favoring associated to the alloy stability, are crucial features for the comprehension of their electronic and atomic properties in dependence with the composition 47 . Furthermore, the Pb/B ratio for the metal size creates crystalline distortions (combined with the different magnitude for the spin-orbit contributions) which gives important insights for electronic characterization of those systems [48][49][50] . As such, a theoretical study for perovskite alloys needs a proper statistical approach relative to the configurational sampling constituting the statistical ensemble, which is required to calculate the average of thermodynamic and structural properties.
Here, we have performed first-principle calculations based on the Density Functional Theory (DFT) to investigate possible perovskite MAPbI 3 -based alloys. The generalized quasi-chemical approximation (GQCA) was used as statistical method, from which thermodynamic properties and averages of the structural parameters can be calculated for a wide chemical range at arbitrary temperatures. Thus, an improved picture on the perovskite alloys, since Si, Ge, and Sn metals present different relative atomic size with respect Pb, were studied in a pseudo-cubic MAPbI 3 structure, considering their local impact on the structure for different direction within the crystal.

cluster expansion and thermodynamic treatment
The structural and thermodynamic behaviour of the perovskite alloys were investigated through a rigorous and systematic statistical description based on the GQCA 51 . In the GQCA, the alloy (mixture) is described as an ensemble of clusters (herein our supercell), statistically and energetically independent of the surrounding atomic configuration. It has been demonstrated that this model successfully describes the physical properties of several 2D and 3D alloys, as well as to 2D sheets [52][53][54][55] . Furthermore, the GQCA method also has been employed in the thermodynamic analysis of perovskite alloys of MAPb(I 1−x Br x ) 3 56 , however, the method was still not employed for perovsksite alloys from the metal perspective.
Within the GQCA size and shape of the clusters play an important role, wherein the supercell model has two advantages: (i) it has a reasonable size for taking into account the local correlation; and (ii) it has an exact counting scheme for the configurational entropy, since no two clusters share the same alloying atom. Based on that, the Fig. 1(a) shows a representation of a MHP as a cubic structure (symmetry group O h ) with the CH 3 + NH 3 cations balancing the − MI 3 anions charges of the octahedrals. We used a supercell with 2 × 2 × 2 expansion of a cubic perovskite by starting from the MAPbI 3 system, from which the alloys are made by replacing the 8 octahedral central sites by Sn, Ge, and Si, named by the letter B in the general case, to build the CH 3 NH 3 Pb 1−x Sn x I 3 , CH 3 NH 3 Pb 1−x Ge x I 3 , and CH 3 NH 3 Pb 1−x Si x I 3 systems, respectively.
Regarding the 8 sites involving the replacement of 2 metal species, as shown in Fig. 1(b), the total number of possible atomic configurations is given by 2 n , where n is the number of sites labeled by 12345678, i.e., resulting on 2 8 = 256 possible configurations for each alloy. However, the 256 atomic configurations can be organized in = J 22 symmetry equivalent classes by considering all the O h space group operations. The Table 1 describes the 22 classes with respect the replacement of the octahedral sites, wherein Pb atoms are labeled by A and the Sn, Ge, and Si atoms by B.
The Fig. 2 shows a representation of the relative positions of the octhedral occupied by Pb (blue) and B (red) of the 22 classes, as well as their respective compositions x and degeneracies g j . Thereby, to describe our statistical ensemble for the perovskite alloys, we considered the set of 9 compositions, as x = 0, 0.125, 0.250, 0.375, 0.500, 0.625, 0.750, 0.875, 1, which were defined by the quantities of both metals involved in the alloy formation. Thus, for a given N as the total number of metals involved (or as the total number of sites occupied aforementioned), = x n n j with n j as the number of Sn, Ge, and Si atoms, and − n n j the number of Pb atoms in the cluster j. Thus, the excess energy of each of those j configurations among the 22 possibilities with internal mixing energy Δε j can be defined by where, E j , E MAPbI 3 , and E MABI 3 are the total energies of the cluster configuration j, the cluster of MAPbI 3 , and the cluster of MABI 3 with B = Pb, Sn, Ge, and Si pure perovskites. As such, the internal energy is calculated by , where x j is the probability distribution for the occurence of a cluster with configuration j. As described elsewhere 51,52,54,55 , the occurence probability x j of equivalence class j can be estimated by the constrained minimization of the Helmholtz free energy, i.e., and average of composition x as calculated by ∑ = = n x x T nx ( , ) j J j j 1 51,52,57 . Thereby, the x x T ( , ) j distribution is given by    3 octahedrals to study perovskite alloys with their n j B atoms (Sn, Ge, and Si). The sequence 12345678 labeling the sites in the cluster can be found in Fig. 1(b), where A is Pb and B are the Sn, Ge, and Si atoms to each alloy, where g j is the degeneracy factor. www.nature.com/scientificreports www.nature.com/scientificreports/ , and η is an adimensional parameter obtained by the average composition constrain, and g j is the degeneracy defined to each j as described in Table 1. The set of probabilities x j permits to calculate any arbitrary property p x T ( , ) for the alloy by where p j is the local property of each cluster class j.
The mixing entropy in ΔS x T ( , ) equation is calculated as

Results and Discussion
We discuss the structural parameters, such as lattice parameters (on the orthogonal directions a, b, and c), local M-I distances (d M-I ), angles I-M-I, between the lattice constants (α, β, and γ), and the volume (Å 3 ) of the unit cell for the MAPb 1−x B x I 3 perovskite alloys as a function of the composition and temperature. By taking the Pb atom as reference, the atomic sizes decrease rising in the IV group of the periodic table, as B = Sn, Ge, and Si which are, respectively, 4.08, 17.01, and 32.65% smaller with respect to the Pb atom 66 . These differences in the atomic sizes of the metals correlated with the organic cation occupying the different cavity sizes made by the octahedrals, taking the relative orientations on the a, b, and c directions (as represented in Fig. 1(c,d)), permit a detailed atomistic comprehension for the pure and alloys perovskites in different compositions. Furthermore, a thermodynamic characterization is provided through the mixing internal energy (ΔU), mixing entropy (ΔS), excess of free energy (ΔF), as well as the construction of the − T x phase diagram of the perovskite alloys.
We found that the smaller atomic size for Si and Ge when compared with Pb contributes to decrease the lattice constants in up to 4.91% (relative to the b direction) for both MASiI 3 and MAGeI 3 in comparison with MAPbI 3 . As consequence, their octahedrals are locally more distorted, as can be seen in Table 2 through the differences between the shortest and largest d M-I values on all a, b, and c directions. We found that, in general, throughout the < < < Si Ge Sn Pb sequence for the atomic size the shortest d M-I distances increase while the largest d M-I distances decrease, which is an effect of the competition of the metals into neighbor octahedrals by the I in the vertice between them. The angles values between the lattice constants (α, β, and γ) and the octahedral connection angles, i.e., φ − − M I M , reveals that for MASnI 3 and MAPbI 3 the local distortions are similar, since their atomic sizes for Sn and Pb are similar. However, for MASiI 3 and MAGeI 3 the small metal occuping the octahedral sites promote higher deviations for the α, β, and γ angles with respect to the 90°, by leading also to the decreasing of the φ − − M I M on all directions also as a local distortion effect on the octahedrals. Our unit cell volume results increasing as MAPbI3 in correlation with the metal size, i.e., < < < Si Ge Sn Pb, suggest the same tendency relative to the cavity size where the organic cation is sited. For instance, the relative similarity between the MAPbI 3 and MASnI 3 pseudo-cubic structures also can be seen as a similar effect of the organic cation orientation on the a, b, and c lattice directions, yielding a low structural distortion on the pseudo-cubic motif and a low dependency of the structural parameters on a, b, and c directions with respect to the organic cation orientation. Consequently, the largest and shortest d M-I values are similar on b for MAPbI 3 (3.17 Å) and MASnI 3 (3.13 Å) due to the CH 3 and NH 3 hydrogen, while on a and c the    www.nature.com/scientificreports www.nature.com/scientificreports/ C-N bond axis its slope effects in the cavity are more pronounced on large and short d M-I values. Conversely, as an effect of the small metal size and a smaller cavity volume, the stronger distortion observed for MASiI 3 and MAGeI 3 by comparing with MAPbI 3 indicates a higher dependency relative to the organic cation orientation.
Therefore, we considered the momentary orientation of the organic cation to understand its effects on the − MI 3 inorganic octahedra. As such, Fig. 1(c,d) shows the MA + C-N bond axis as momentarily oriented on a, giving the C-N bond axis sloped in the cavity on b, providing lowest energy configuration for the CH 3 NH 3 group as reported by several atomistic simulation studies 56,69,70 . Thus, it is reasonable to expect that even though the high temperature effects promote the MA + free reorientation in the cavity for MAPbI 3 , while the reorientation may be slightly limited in the MASiI 3 and MAGeI 3 pseudo-cubic structures.
Lattice parameters of the alloyed perovskites. The optimization of synthesis process of pure 4,18,34 and alloy 30,31,41 MHP at room temperature have widely been investigated, especially through self-assembling principles from the chemical precursors for the metal halides. As such, our statistic averages were calculated through GQCA at 300 K from the weighted contribution of each j configuration, providing the average of the structural parameters for the MAPb 1−x B x I 3 alloys as a function of the composition at room temperature.
We calculated the average lattice constants into the supercell on the a, b, and c directions, as well as the angles between them and the volume for the unit cell for each j cluster alloy (Fig. 3). Thus, the results connect the values for the MAPbI 3 ( = x 0) and MABI 3 ( = x 1), B = Si, Ge, and Sn. We found that the lattice parameters for the Fig. 3) alloy follow the Vegard's law 71 on the a and b directions, i.e., linearly decrease as a chemical specie with smaller atomic size is included into the bulk, while for the c direction it is observed a bowing. This result is due to the effects of the organic cation orientation taken as reference, wherein the C-N bond into the small cavity size yields different constraints on the lattice on the different directions. For example, on the plane made by b and c directions, on which the C-N bond of the CH 3 + NH 3 is perpendicular, there is a deviation of the linearity with respect the composition as an effect of greater permissiveness of lattice distance adjustments with respect to the composition. Furthermore, as a consequence of the higher contraction of the lattice parameters as the Si atoms amount increases, we found a crossing over of the lattice parameters on the a and c, wherein the organic cation orientation yeilds lattice distances as < a c and > a c for the compositions < .
For the MAPb 1−x Ge x I 3 alloy, the lattice parameter results were similar with the MAPb 1−x Si x I 3 (panel (b) in Fig. 3). We found that the Vegard's law is followed for all the composition range for the a and b directions. The crossing over between a and c appears from > .
x 0 750, from which lattice parameters are > a c. Similarly, this result is also explained for the gradual contraction of the lattice parameters due to the small size of the Ge, as a consequence of the replacement of the Pb by Ge atoms, by yeilding a decreasing of the cavity size. As such, even though the C-N atoms of the MA + are oriented perpendicular to the plane made by b and c orientations, the   www.nature.com/scientificreports www.nature.com/scientificreports/ crossing over between a and c parameters for MAPb 1−x Ge x I 3 appears for lower quantities of Ge when compared with MAPb 1−x Si x I 3 , which is a consequence of larger Ge size by comparing with Si.
For the MAPb 1−x Sn x I 3 lattice parameters shown into the panel (c) in Fig. 3, since the atomic sizes of the Pb and Sn atoms are similar there is no crossing over between a and c parameters, and the Vegard's law is followed in all composition range connecting linearly the lattice parameter of the MAPbI 3 and MASnI 3 pure perovskites. As such, the linearity connecting the lattice parameters for the Pb-I-Pb, Pb-I-Sn, or Sn-I-Sn combinations are independent of the direction, suggesting that the pseudo-cubic structure for MAPb 1−x Sn x I 3 alloy is quite resistent with respect to the composition.
Lattice angles and volume of the alloyed perovskites. The panels (d), (e), and (f) in Fig. 3 show the lattice angles (α, β, and γ) for the MAPb 1−x Si x I 3 , MAPb 1−x Ge x I 3 , and MAPb 1−x Sn x I 3 alloys as a function of the composition. We found that α and β angles slightly increase between = x 0 and = x 1 for MAPb 1−x Si x I 3 , lying into the interval 90°-92°. Conversely, γ decrease sharply with angle from 90° up to 84°, which is explained by the strong distortion on the pseudo-cubic structure due to the gradual replacement of Pb by Si atoms. Additionally, the volume of the unit cell for the MAPbI 3 and MASiI 3 pure perovskites are linearlly connected as function of the composition, with values lying between 265.79Å 3 and 235.29Å 3 , which describes the constraction of the alloy by correlating with the Vegard's law.
For MAPb 1−x Ge x I 3 alloys, we found that α and β are close to 90° between x = 0 and 0.875, while γ decrease sharply similarly with respect to the MAPb 1−x Si x I 3 , that is between 90° up to 85°. This result shows the effects of the metals size differences, as well as the linear contraction for the volume of the unit cell between MAPbI 3 and MAGeI 3 . This behaviour is also indicated for the α and β kept in 90° for the MAPb 1−x Sn x I 3 due to the similar size by comparing Pb and Sn, while γ lie into a short interval between 89°-90°. The shortest d M-I values (Fig. 4 leftmost) in the alloys are determined by the Si-I, Ge-I, and Sn-I distances, which is an effect of the metal size differences with respect to the size of the Pb. One observes that for the a and b directions that for Pb-rich compositions the largest d Pb-I values are higher than d M-I values, wherein for few quantities of B the shortening of the B-I distance in an particular octahedral results in an elongation for the Pb-I distance relative to the neighbor octahedral. Thanks to these differences for the metal sizes into the clusters j, one observes an increasing of the amplitude for the shortest and largest d M- 3 . This behaviour is explained by the local distortions on the octahedrals as the metal size differences are pronounced, also as  With the results above discussed, we note the important role of the atomic size difference between the metals involved in the perovskite alloy formation. For MAPb 1−x Sn x I 3 , as a case of similar size for the metals, the small local distortions into the octahedral and the linearity correlation between the MAPbI 3 and MASnI 3 pure perovskites show a preference in preserving the pseudo-cubic structure similar to the pure perovskites in the whole range of compositions. Conversely, the MAPb 1−x Si x I 3 and MAPb 1−x Ge x I 3 alloys are examples of large difference between the atomic size of the metals, we found that the composition is an additional variable with respect to the temperature to promotes strong distortions into the phase, reinforcing the necessity of a proper statistical analysis to correlates the thermodynamic stability with the structural motifs for the alloy. thermodynamic parameters and ordering preference. To predict the most favorable local arrangement of metal in the octahedral inner sites, i.e., the − PbI 3 and − BI 3 relative configuration, the alloy excess energies (Δε j ) were calculated in order to determine the composition-dependent cluster probabilities (x j ). Consequently, by knowing x j as dependent of Δε j and the degeneracies g j for each j-configuration, we calculate the mixing free energy ΔF x T ( , ) from the contributions of the interplay between the configurational entropy ΔS x T ( , ) and the internal energy ΔU x T ( , ) through the GQCA. As such, below we provide a thermodynamic discussion to enlighten the preferential ordering correlated to the stability of the MAPb 1−x B x I 3 perovskite alloys. The panel (a) shows all the positive Δε j values for the MAPb 1−x Si x I 3 alloy, which means that at = T 0 K there is a high stability of the MAPbI 3 and MASiI 3 pure perovskites in detriment of the alloy. We found that all the pseudo-cubic configurations strongly distorted between < < x 0 1 lie into Δε j values between 10 and 63 meV/ metal, which is an evidence of the high strain yielded by the difference of the atomic size between Pb and Si. By comparing with the Ge alloy, in panel (b), an energetically favored cluster with Δε = − . 1 28 meV/metal j is observed at = .
x 0 125, which correlates with a tendency to form a long-range ordered alloy depending on the temperature. However, all the distorted pseudo-cubic configurations for x 0 125 1 .
< < present Δε j values between 7 and 30 meV/metal for MAPb 1−x Ge x I 3 . This result suggests that for an MAPb 1−x B x I 3 (with B = Si, Ge, or Sn) perovskite alloy energetivally favorable two stability parameters are correlated: (i) the proportion (composition for the alloy) between the metals occupying the octahedral sites; and the (ii) magnitude of the atomic size difference between the metals involved.
As a consequence of small difference between the atomic size for Pb and Sn in the MAPb 1−x Sn x I 3 , the Δε j values lie in an interval of energies between −9 and 4 meV/metal. Thus, several configurations can be easily favorable when the entropy effects be considered. Therefore, as previously discussed for the structural parameters, such as the lattice parameters, d M-I , and α M-I-M as a function of the a, b, c directions, this results suggest that the replacement of Pb by Sn yields only slight changing in the MAPb 1−x Sn x I 3 structure. Among all the configurations between = x 0 and x 1 = for the short range of Δε j values for the MAPb 1−x Sn x I 3 alloy, additionally to the = .
x 0 250 are represented by j = 3, 4, and 5, which present Δε j in −3.43, −0.88, and −0.82 meV/metal respectively. We observe that the ordering j = 3 as represented in Fig. 5 is the most favored, in which the stability is reached by the stacking of the intercalated − PbI 3 and − SnI 3 octahedral rows.
Perovskite alloys free energies and ordering. Here, we discuss the statistical contributions of the Δε j values for the thermodynamic properties for the alloys under the temperature effects through the GQCA method. The variation in the energy of mixing (ΔU) and entropy of mixing (ΔS) used to calculate the Helmholtz free energies (ΔF in m/metal) for the MAPb 1−x Si x I 3 , MAPb 1−x Ge x I 3 , and MAPb 1−x Sn x I 3 alloys within the GQCA are shown in Fig. 6. In order to verify the entropy effects for the stabilities of the alloys, we analysed these parameters as a function of low and high temperatures, e.g., 100, 300, 500, 700, and 900 K. www.nature.com/scientificreports www.nature.com/scientificreports/ One observes by the ΔS symmetrical curves with temperature around = .
x 0 500, panels (d), (e), and (f), indicating that all the alloys follow an ideal entropy expression at high temperatures, i.e., − . The ΔU curves for MAPb 1−x Si x I 3 -panel (a) -present a positive parabolic behaviour due to the higher stability of the MAPbI 3 and MASiI 3 pure perovskites in comparison with the alloy. Thus, the profile of the ΔU and ΔS curves indicates that the alloy can be stabilized by entropic contributions, consequently by increasing the magnitude of disorder through the insertion of Si atoms, which promotes the contribution of several j configurations. The panel (g) shows a behaviour slightly asymmetric for the ΔF curve around x 0 500 = . , so that for < T 300 K we found Δ > F 0 as an evidence of the instability of the alloy at low temperatures. However, for > T 300 K one observes that Δ < F 0 and the alloy starts to be stable, and for temperatures between < < T 300 K 500 K there are points throughout ΔF with same tangent, indicating the existence of a miscibility gap for an extensive range of temperatures.
For the MAPb 1−x Ge x I 3 alloy, we found that the ΔF -panel (h) -presents points with same tangent for < < T 100 K 500 K, which is a range of lower temperatures for the miscibility gap than for MAPb 1−x Si x I 3 . The ΔF reaches symmetrical curves for temperatures higher than 500 K, which the entropy effects start to be dominant over the small negative ΔU values, panel (b), for few Ge quantities. Conversely, with the increasing of the temperature, the disordering is reached with the weighted contributions of all − PbI 3 and GeI 3 − octahedrals configurations, from which the random configurations for = .
x 0 500 compositions are the most favorable. The ΔU curves profile for the MAPb 1−x Sn x I 3 alloy -panel (c) -show the effect of the favorable ordering for compositions with excess of both Pb and Sn metals, as x = 0.125, and 0.875. Firstly, this yields two regions for Δ < U 0 relative to the orderings as represented in Fig. 5, so that the alloy stabilizes when the − SnI 3 individual octahedrals are completely involved by − PbI 3 octahedrals, as well as for the opposite configuration, i.e., − PbI 3 individual octahedrals completely involved by − SnI 3 . Secondly, the short range of excess energies for MAPb 1−x Sn x I 3 yields a short interval of ΔU variation as a function of the composition and temperature. Thus, for temperatures higher than 100 K the entropy effects are dominant, so that the shape of the ΔF curve becomes more symmetric in order that the contribution of all configurations increases with the temperature, consequently, increasing the disordering of − PbI 3 and − SnI 3 positions in the alloy. As such, it is expected to observe a miscibility gap in MAPb 1−x Sn x I 3 alloy with pseudo-cubic structure only for very low temperatures, since there is no effective x 0 125 ( = j 2) and = .
To investigate the similarity between the GQCA probability x x T ( , ) j and a random alloy, relative to the random contribution of a particular j configuration in a range of temperatures, we present the KL divergence, namely, Fig. 7. The maximum divergence at low temperatures means that in a particular composition the configuration j relative to the x j dominates over the others, and in so far as the temperature increases the divergence goes to zero, i.e., the system starts to behave as a random alloy. In Fig. 7(a-c) are plotted the D KL (  x x j j 0 ) for MAPb 1−x Si x I 3 , MAPb 1−x Ge x I 3 , and MAPb 1−x Sn x I 3 for the compositions at x = 0.125 and 0.875 as a function of the temperature. For MAPb 1−x Si x I 3 at very low temperatures one observes a tendency for phase segregation with the formation of MAPbI 3 and MASiI 3 pure perovskites, as observed through x j plots in Fig. 7(d,g), wherein there is a predominancy of x 1 and x 22 configurations at compositions x = 0.125 and 0.875, respectively. For MAPb 1−x Ge x I 3 , panel (b), clearly it is seen that at = .
x 0 125 the divergence is smaller than for = .
x 0 875 at very low temperatures, so that the ordering given by the = x 1 2 configuration dominates at = .
x 0 125 composition for the alloying at Pb-rich compositions (panel (e)), as well as yields a small contribution at x 0 875 = .
together with the dominant x 22 configuration (panel (h)). Conversely, for MAPb 1−x Sn x I 3 at = .
x 0 875 and at low temperature, panel (c), one observes the high miscibility between Pb and Sn as an effect of the similar metal size. Thereby, the x j plots, panels (f) and (i), show the predominancy of the = x 1 x 0 875) at 0 K, demonstrating the tendency of the system in organizing itself in energetically favored alloyed configurations even at very small temperatures. The observed ordering of atomic distribution, however, does not persist for temperatures above 150 K.
Phase diagram of the alloys. To identify regions of stability and metaestability as a function of the temperature and composition, we built the phase diagram for the alloys at the pseudo-cubic structure, as shown in Fig. 8. We observe for MAPb 1−x Si x I 3 (leftmost), MAPb 1−x Ge x I 3 (middle), and MAPb 1−x Sn x I 3 (rightmost) critical temperatures (T c ) of 527, 440, and 204 K, respectively. Above T c the solid solution are stable for any composition, whereas below them are evidenced the presence of miscibility gaps to each alloy defined by spinodal lines, given by ′ x 1 and ′ x 2 (blue regions), as well as binodal lines from the x 1 and x 2 points defining the − T x metaestability regions. For MAPb 1−x Si x I 3 , due to the ΔF profile observed especially for its formation at room temperature (300 K), as shown in Fig. 6, its phase diagram presents unstable regions from Pb-to Si-rich compositions at low temperatures, yielding two miscibility gaps in dependence of the composition region. For instance, for Pb-rich compositions the first miscibility gap lies between ′ = .
x 0 19 2 , whereas the second one, relative to the Si-rich compositions, lies in the interval of ′ = .
x 0 89 2 . One observes that the first miscibility gap reduces as the temperature increases up to 445 K, from which a solid solution is formed for Pb-rich compositions. However, only from = T 527 K c at = .
x 0 68 the solution solid is stable into all composition interval. Furthermore, from the end of the first miscibility gap up to the beginning of the second one, i.e., between the set of For the MAPb 1−x Ge x I 3 alloy, a stable solid solution is observed in all range of temperatures only for Pb-rich compositions between < < x x 0 1 for = .
x 0 20 1 , while at 300 K the miscibility gap appears between ′ = . . At 400 the intervals for miscibility gap and metaestable phases are shorter than at room temperature. By comparing the MAPb 1−x Si x I 3 and MAPb 1−x Ge x I 3 alloys at compositions Si-, Ge-, and Pb-rich, one observes a behavior very different due to the effect of the Pb/Si and Pb/Ge metal size differences. Even though there is a stability of the MAPb 1−x Ge x I 3 alloy for Pb-rich into all temperatures, the symmetrical-like spinodal line at 300 K yields a stability for a range of Ge-rich compositions. Additionally, metaestable phases are observed into x = 0.20 − 0.31 and x = 0.70 − 0.81 intervals of compositions.
We found that the critical temperatures T c for MAPb 1−x Si x I 3 , MAPb 1−x Ge x I 3 , and MAPb 1−x Sn x I 3 , from which the solid solution at all compositions is stable, correlates with the atomic size difference for the metals involved. For example, the Fig. 8 shows also the phase diagram for the MAPb 1−x Sn x I 3 , in which one observes the effects of small difference between the Pb and Sn atomic size from the high solubility of the metals into MASnI 3 and MAPbI 3 , respectively. Since the critical temperature is = T 204 K c , at 300 K a stable solid solution is observed within all range of compositions, which is in agreement with Hao et al. 31 experiments for the synthesis of Other configs.  www.nature.com/scientificreports www.nature.com/scientificreports/ MAPb 1−x Sn x I 3 who observed a high stability of the pseudo-cubic structure of the MAPb 0.5 Sn 0.5 I 3 alloy, as well as in others compositions. Furthermore, we found a miscibility gap slightly symmetrical for MAPb 1−x Sn x I 3 appearing only at very low temperatures, since the local distortions into the structure are suppressed and the entropic effects are restricted to the configurations of the − PbI 3 and − SnI 3 octahedrals. Therefore, by taking as reference = T 443 K and = T 473 K as experimental temperatures in which the MAPbI 3 and MASnI 3 pure perovskites start to be decomposed 25,26 , our results show that there is a range of temperatures from = T 204 K c in which the MAPb 1−x Sn x I 3 is stable as a random alloy before a possible thermal decomposition. Furthermore, for the others alloys, those results may be as a guide for future synthetic process for the MAPb 1−x Si x I 3 and MAPb 1−x Ge x I 3 alloys, from which it is expected the phase segregations for some range of compositions.

conclusions
In summary, we have performed first-principles calculations combined with a statistical approach based on cluster expansion to study the stability, effects of disorder, distortions, thermodynamic properties, and phase segregation of the pseudo-cubic phase of MAPb 1−x B x I 3 alloys for B = Si, Ge, and Sn.
Our results indicated that the metal atomic size plays an important role on the pseudo-cubic properties of the pure perovskites, e.g., as the similar local distortions for the MAPbI 3 and MASnI 3 octahedrals since their metals have almost the same atomic size. As such, the MAPb 1−x Sn x I 3 alloy presents lattice parameters in good agreement with the Vegard's law for the whole range of compositions, wherein the alloy adopts a random − PbI 3 and − SnI 3 octahedral configurations. Conversely, MASiI 3 and MAGeI 3 in pseudo-cubic structure are strongly distorted as an effect of their second smaller metal in comparison with Pb, suggesting a higher limitation of the organic cation orientation on the lattice directions for the MAPb 1−x Si x I 3 and MAPb 1−x Ge x I 3 alloys, since the cavity volume is reduced. For those cases, the alloys follow the Vegard's law for some particular lattice directions, whereas the others there is a pronounced bowing throughout the range of compositions.
The thermodynamic results revealed that the MAPb 1−x Ge x I 3 alloy is stable for Pb-rich compositions, i.e., between < < .
x 0 020 at 300 K, by presenting a preference for an ordered configuration in which one − GeI 3 octahedral is surrounded by PbI 3 − octahedrals. Conversely, MAPb 1−x Si x I 3 is not favored into very large range of compositions, and even though has presented an interval of compositions into which the alloy is metaestable (into x = 0.19 − 0.45), it indicated a high tendency for segregation phase in MAPbI 3 and MASiI 3 pure perovskites. Thus, the addition of small metal atoms yields strong local distortions on the octahedrals, resulting in very high critical temperatures for these alloys. As an exemple of miscibility, the MAPb 1−x Sn x I 3 alloy presented a critical temperature lower than the room temperature, at 204 K, which is very lower than the temperature of decomposition for the MAPbI 3 and MASnI 3 pure perovskites. Thus, the alloy is favored as a random alloy in all compositions, revealing that there is a safe range of temperatures in which the MAPb 1−x Sn x I 3 alloy properties can be tuned before the material be thermally decomposed.
Therefore, beyond the temperature as variable, the correlation between composition and atomic size, relative to the second metal in MAPbI 3 -based alloys, is a crucial element to promotes the phase stability. The thermodynamic characterization of these alloys for intermediate Pb-based compositions showed the importance of the planning relative to the experimental synthesis conditions, such as temperature and composition, aiming the structural motifs correlated with their performance into solar cells devices. theoretical Approach and computational Details. In this study, to calculate the total energy of the configurations of the alloy in all the range of compositions, we employed spin-polarized calculations based on DFT 73,74 within the semilocal Perdew-Burke-Ernzerhof 75 (PBE) formulation for the exchange-correlation energy functional. The projected augmented wave 76,77 (PAW) method as implemented in the Vienna ab initio simulation package (VASP), version 5.4.1. 78,79 was used to solve the Kohn-Sham (KS) equations, in which the scalar-relativistic approximation is considered to the core states, as well the spin-orbit coupling (SOC) interactions. However, SOC is an important relativistic phenomenon in Pb-based perovskites [63][64][65] , especially occurring within non-spherical atomic orbitals and affecting the directionality of the metal bonds 7 , so that we included SOC interactions also for the valence states in all our calculations.
For total energy calculations, we employed a plane-waves expansion with kinetic energy cutoff of 500, by integrating over the Brillouin zone calculated considering a Monkhorst-Pack k-mesh of 4 × 4 × 4. The total energy convergence to 1.0 × 10 −5 eV with Gaussian broadening parameter of 50 for all calculations. Finally, the equilibrium of the Hellmann-Feynman forces on every atom were reached to forces smaller than 0.010 eV/Å.