## 2 Computational Methods  
All calculations were performed within the framework of density functional theory, as implemented in the VASP code [24]. The generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) was used as the exchange correlation functional. The plane wave cut-off energy was set to 520 eV. The electronic structures were calculated using 41 \(\times\) 41 \(\times\) 1 Monkhorst-Pack \(k\)-point sampling for Mo\({}_{2}\)CT\({}_{2}\), where T represents surface functional groups (O, OH, and F). The 1 \(\times\) 1 \(\times\) 1 unit cell of a bulk crystal was selected on the basis of previous calculations [23]. The monolayers and multilayers were constructed by adding a layer into the 1 \(\times\) 1 \(\times\) 5 supercell in VESTA and were optimized when a layer was added. The optimized lattice constants and bond lengths are shown in Table 1. The electronic structures were calculated along the highly symmetric k-points line \(G\)-\(A\)-\(H\)-\(K\)-\(G\)-\(M\)-\(L\)-\(H\) in the first Brillouin zone.  
The optical properties of MXenes have a high dependence on the energy structure, specifically on the bandgap energy and the type of bandgap (direct or indirect). In addition, the optical properties are related to the dielectric function, which is described by the optical frequency (\(\omega\)), and the real and imaginary part as follows:  
\[\varepsilon(\omega)=\varepsilon_{1}(\omega)+i\varepsilon_{2}(\omega)\]  
The imaginary part is investigated based on the Fermi distribution and the crystal wave function [25; 26]:  
\[\varepsilon_{2}(\omega)=\frac{4\pi^{2}e^{2}}{\Omega m^{2}\omega^{2}}\sum_{kmn ^{\prime}\sigma}\langle k,n,\sigma,|p_{i}|,k,n^{\prime}\rangle\sigma\times \langle k,n^{\prime},\sigma,|p_{i}|,k,n,\sigma\rangle k_{kn}(1-f_{km^{\prime}} )\delta(e_{kn^{\prime}}-e_{kn}-\hbar\omega)\]Here \(e\), \(m\), and \(\Omega\) are the electron charge, mass, and crystal volume, respectively. The crystal wave function \(\mid kn\sigma\)\(>\) is proportional to the eigenvalue \(n\)th, the crystal momentum \(k\), and the spin \(\sigma\). By using the Kramers-Kronig transformation [27; 28; 29], the real part \(\varepsilon_{1}(\omega)\) can be extracted from the imaginary part \(\varepsilon_{2}(\omega)\). Once the dielectric function \(\varepsilon(\omega)\) has been determined, based on the standard optical relationships, all other optical parameters, such as the refractivity, reflectivity, and absorption or optical conductivity, can be estimated.  
The Boltzmann's transport equation (BTE) within the constant relaxation time (\(\tau\)) approximation is employed for calculating the thermoelectric (TE) properties of Mo\({}_{2}\)C-MXenes based on the electronic structure results. The Boltztrap code is a well-established computational simulation tool for TE properties [30]. The dimensionless figure of merit (ZT) describes the TE performance of a material, and is given as:  
\[ZT=\frac{S^{2}\sigma}{\kappa}T\]  
where \(S\) is the Seebeck coefficient; \(\sigma\) is the electrical conductivity; \(\kappa\) is the thermal conductivity; and \(T\) is the absolute temperature [31].  
It should be noted that the ideal TE material is defined as having a high \(ZT\) value, in which the specifications include a high power factor with a high Seebeck coefficient (\(S\)), high electrical conductivity (\(\sigma\)), and low thermal conductivity (\(\kappa\)). On the other hand, these quantities depend on \(\tau\) and are difficult to track down solely through band structure analysis. However, all investigated Mo\({}_{2}\)C- have the same hexagonal structure, and they initially all belong to the same material family as MAX phases. Therefore, the electron-scattering mechanisms in Mo\({}_{2}\)C- can be predicted to be similar, and the comparison of thermoelectric results is meaningful. The calculation results of the thermoelectric properties are presented and compared at different temperatures in the range of 300-600 K for all monolayer and multilayer structures.