Surface Functionalization of 2D MXenes: Trends in Distribution, Composition, and Electronic Properties

Using a multiscale computational scheme, we study the trends in distribution and composition of the surface functional groups −O, −OH, and −F on two-dimensional (2D) transition metal carbides and nitrides (MXenes). We consider Ti2N, Ti4N3, Nb2C, Nb4C3, Ti2C, and Ti3C2 to explore MXenes with different chemistry and different number of atomic layers. Using a combination of cluster expansion, Monte Carlo, and density functional theory methods, we study the distribution and composition of functional groups at experimentally relevant conditions. We show that mixtures of functional groups are favorable on all studied MXene surfaces. The distribution of functional groups appears to be largely independent of the type of metal, carbon, or nitrogen species and/or number of atomic layers in the MXene. We further show that some properties (e.g., the work function) strongly depend on the surface composition, while others, for example, the electric conductivity, exhibit only a weak dependence.


Cluster Expansion
We start by generating a cluster expansion, which is constructed based on the assumption that energy of the system can be expanded on a series of clusters (or nearest neighbours depicted on Fig. S1 (a)) and their effective interactions (ECI) Fig. S2. The model Hamiltonian is then obtained by fitting ECI to the set of energies calculated by DFT method. For each system, we produce a DFT set of structures with binaries and ternaries of functional groups within a full range of O, OH, and F concentrations. We used set of 140 structures for Ti 2 N, 80 structures for Ti 4 N 3 , 71 structures for Ti 2 C, 106 structures for Ti 3 C 2 , 132 structures for Nb 2 C, and 119 structures for Nb 4 C 3 .
To evaluate the accuracy of each cluster expansion, we used cross validation scores. Cross validation scores were obtained within a range of 7-15 meV which gives rise to an error at maximum of 10 %. As an example, the scattered plot of DFT energies vs. CE energies as a function of O concentration is shown for Ti 2 N and Nb 2 C in Fig. S1(b). In both regressions, the predicted energies are reproducing the actual calculated energies well, with no outliers. In a low mixing energy region, the structures with no presence of O are mostly common, whereas in a higher mixing energy region, mostly structures with sufficient amount of O are present.

Monte Carlo simulations and special quasi-ordered structures
For MC sampling we used a 40x40x1 conventional supercell with 3200 functional group sites. Sampling was carried out using canonical ensemble at the temperature decreased from 2000 K to 300 K with a 100 K intervals. The calculations are carried out with 200000 steps, where the number of steps was tested with respect to the convergence of energy. Figure S2: Effective cluster interactions for Ti 2 N as a function of clusters radius.
Best representative structures were acquired using a special quasi-ordered structures method (SQoS) [10,5]. Cluster correlations from larger equilibrated structures were used to generate smaller representative structures by minimizing the objective function. Minimization is done by canonical Monte Carlo simulation with 100000 steps for the 4x4x1 supercells. The calculated structures are mimicking distribution of functional groups of the larger structures.

Free energy calculations
In order to compare the stability of MXene sheets with different surface terminations, we determine Gibbs free energy of formation, which is obtained by finding Gibbs free energy for all the constituents. Here we define it for the terminated sheet MXT with respect to the bare, unterminated sheet MX as: where n i are the number of termination atoms of type i and µ i are their chemical potentials. G is the Gibbs free energy of the system and • refers to standard conditions: room temperature and in solution with pressure p = 1 atm. Consequently, the free energy of the sheet should include the vibrational contributions as well as the interaction with the solution. We assume that the two contributions do not depend on each other and thus where E(MXT) is the DFT total energy, ∆ vib F is phonon contribution to free energy at room temperature evaluated in vacuum, and ∆ sol E is the solvation energy evaluated at T = 0 using implicit solvation models. (pV term can be ignored at p = 1 atm.) We use vibrational and solvent contributions for calculating formation energies, since the order of magnitude is comparable, we use the values from calculated earlier Ti 2 C and Ti 3 C 2 [3]. Note, that following framework is adopted from the theoretical concepts proposed by Todorova and Neugebauer [9]. We determine the chemical potentials of H, O, and F in accordance to the experimental conditions. Chemical potential of O is determined via chemical potential of H and water: is solvated water. We use the experimental Gibbs free energy of formation where ∆ f G(H 2 O) = −237.14 kJ/mol = −2.458 eV [2]. The chemical potential of H + ions, would depend on the electron chemical potential. The formation energies of H + and F − ion are written as where ∆ f G • ( +/− ) are the energy of solvated ions. Data are taken from NIST-JANAF thermochemical tables and the hydration energies from Refs. [6,7]. pH is directly related to the H + concentration in the solution: where c 0 = 55.55 mol/l is the concentration of H 2 O molecules in water. From this and by using Eq. S4 we obtain Similarly, for F chemical potential, the Eq. S4 can be rewritten.
and finally T · [ln(55.55) + ln 10 · pH] (S10) Eventually, we have chemical potentials of all species connected to pH, electron chemical potential µ e , and temperature T. Experimental conditions suggest that T would be constant since the synthesis happens at room temperature (T=298 K). However, pH and electron chemical potential may vary. From the dependence of H chemical potential on the electron chemical potential, we can obtain a computational standard hydrogen electrode (SHE) potential equal to 4.7 eV (pH=0). Further, to present our results according to experimental conditions, we vary pH and open circuit potential (OPC), where OPC or (U-U SHE ) would be a difference between negative electron chemical potential and SHE potential (−µ e − U SHE ). Thus, we can vary µ e and pH in equations (S7, S10), to calculate Gibbs free energies of the systems and their minimum energy compositions of functional groups.

DFT calculations
All Density functional theory calculations were performed using the Vienna ab initio simulation package (VASP) [4] together with projector augmented plane wave method (PAW) [1]. Perdew-Burke-Ernzerhof exchange-correlation functional for solids (PBEsol) has been used for all calculations [8], which was selected based on benchmarkings performed in Ref. [3]. The optimal plane-wave cutoff energy was chosen as 550 eV according to the convergence tests. The k-points set of 16x16x1 was chosen as optimal for all unit cell calculations and set of 4x4x1 k-points was used for calculation of best representative structures (4x4x1 size of supercell). we evaluate the charges associated with each atom, using a Bader charge analysis. The averaged number of excess electrons, defined as the difference between the Bader charge and the number of valence electrons of the corresponding species, is shown in the Fig. S4 as a function of O-OH concentration. The averaged number of electrons are associated with O and OH linearly depend on the amount of OH in the system [ Fig. S4(a,c)], whereas the values are overall rather similar and we found no clear correlations with the maximum mixing energy or substrate and the number of atomic layers dependent features. On the other hand, charges of N and C atoms differ depending on the surrounding metallic species. In all Nb-based systems C does not change its charge significantly within a whole range of O-OH concentration, while the N and C atoms in Ti-based systems gain more electrons with an increase of OH content, additionally the carbon atoms in Ti-based systems take more electrons, than N atoms. Difference in the charge distributions might be one of the reason of different functional groups composition in Ti-based systems with carbon and nitrogen.