I. INTRODUCTION

In advanced semiconductor technologies, copper (Cu) and tungsten (W) have long been employed as the primary interconnect materials in logic and memory devices ^[1]. However, as interconnect dimensions shrink below the electron mean free path (MFP)---typically 39.9 nm for Cu and 19.1 nm for W---enhanced electron scattering at surfaces and interfaces leads to a significant increase in resistivity ^[2-4]. This resistivity rise, combined with increased RC delay and power consumption, has made interconnect resistance a major bottleneck in chip performance ^[2] Accordingly, alternative materials to replace Cu in logic devices and W in 3D NAND word lines are under active investigation. Among them, molybdenum (Mo) has emerged as a promising candidate for advanced interconnects, particularly at the 2 nm technology node. When deposited via atomic layer deposition (ALD), Mo demonstrates exceptionally lower resistivity (186 n$\Omega$·m at a 24 nm thickness) and possesses shorter MFP (${\sim}$11.2 nm) as compared to Cu and W, making it less sensitive to scaling effects ^[2,4,5]. Moreover, due to its high melting point and excellent electromigration resistance even above 10${}^{8}$ A cm${}^{-2}$, Mo can be integrated directly on dielectrics, thereby eliminating the need for diffusion barrier layers such as Ta/TaN and TiN/WN, which are essential in Cu and W-based interconnects respectively to prevent metal diffusion into adjacent dielectric materials ^[2]. As illustrated in Fig. 1, Mo demonstrates lower resistivity than W across a wide range of thicknesses and maintains comparable performance to Ru at nanoscale dimensions ^[6-8].

This highlights Mo's potential to overcome the resistivity scaling challenges faced by conventional interconnect metals. These attributes underscore Mo's compatibility and potential for next generation scaled semiconductor applications.

To enable practical integration of Mo via ALD, the selection of an appropriate precursor is critical. Since the ALD process relies on sequential, self-limiting surface reaction, which requires precursors to have sufficient volatility and thermal stability to achieve uniform film growth. Commonly studied Mo precursors include molybdenum dichloride dioxide (MoO${}_{2}$Cl${}_{2}$), molybdenum pentachloride (MoCl${}_{5}$), and molybdenum hexacarbonyl (Mo(CO)${}_{6}$) ^[9,10]. Among these, Mo(CO)${}_{6}$ enables Mo and MoS${}_{2}$ film deposition at low temperatures (155--175 ${{}^\circ}$C), reducing the thermal budget by over 60\% compared to MoCl${}_{5}$ (440--470 ${{}^\circ}$C), which is highly advantageous for thermally sensitive substrates ^[11]. Although solid at room temperature, Mo(CO)${}_{6}$ exhibits high volatility and a substantial vapor pressure (0.1-0.27 torr at 25-30 ${{}^\circ}$C), enabling easy sublimation and gas-phase delivery without extensive heating ^[12,13]. These characteristics make Mo(CO)${}_{6}$ an ideal precursor for conformal Mo and MoS${}_{2}$ film growth under thermally constrained conditions and within complex 3D structures. Due to its promising properties, numerous experimental studies have been undertaken to characterize the vapor pressure of Mo(CO)${}_{6}$, including static manometric methods using an MKS Baratron gauge ^[14], paraffin oil bath--based techniques for enhanced thermal precision ^[15], and effusion measurements fitted using the Clarke--Glew equation with various orifice sizes ^[16]. However, computational studies to predict the vapor pressure of Mo(CO)${}_{6}$ still remain limited. This underscores the need of conducting computational investigations to complement experimental results and provide a more comprehensive understanding of its thermodynamic properties.

Concurrently, extensive research studies have also been conducted using Mo(CO)${}_{6}$ as a precursor, aiming to evaluate its film characteristics and assess its suitability across a range of thin-film applications such as Mo, MoO${}_{3}$, and MoS${}_{2}$ ^[17-20]. Joo et al. performed a plasma enhanced atomic layer deposition (PEALD) to deposit nanocrystalline $\gamma$-Mo${}_{2}$N using Mo(CO)${}_{6}$/NH${}_{3}$-plasma as reactants ^[18]. Taijun et al. demonstrated a PEALD process utilizing Mo(CO)${}_{6}$/O${}_{2}$-plasma to deposit high-dielectric ($\kappa \approx 17$) MoO${}_{3}$ oxide films suitable for MoS${}_{2}$ field-effect transistors ^[19]. Sachin et al. further studied the surface reaction mechanism and growth chemistry of MoO${}_{3}$ thin films using thermal ALD with Mo(CO)${}_{6\ }$and ozone as reactants; however, their study revealed a key limitation---residual carbon contamination---resulting from the incomplete decomposition of the CO ligands at lower substrate temperature (398 K) ^[20]. This limitation reveals incomplete understanding of how CO ligands are preferentially desorbed, suggesting the need for further simulation-based investigations to address this gap.

In this study, we applied Density functional theory (DFT) and ab-initio thermodynamics, incorporating vibrational contributions to enhance the accuracy of Gibbs free energy--based vapor pressure predictions for Mo(CO)${}_{6}$, and its initial surface interaction with $\beta$-cristobalite SiO${}_{2}(111)$ surface. To gain further insight into the mechanism and energy requirement for complete removal of the first CO bond from Mo(CO)${}_{6}$, bond dissociation energies (BDE) of Mo-CO at 298 K was calculated. Consequently, the precursor-adsorption chemistry and two CO desorption pathways were studied under the experimentally relevant condition of 435 K and $1.31 \times 10^{-8}$ atm reported by Taijun et al. ^[19]; the desorption pathway with the lowest computed Gibbs free energy was identified as the most thermodynamically preferred route. Given that experimental determination of precursor vapor pressure generally requires high-purity sample synthesis, precise thermal-pressure control, and sophisticated equipment, we instead rely on DFT-based ab-initio thermodynamic simulation approach---a faster and cost-effective computational alternative---providing reliable, predictive insights into precursor behavior.

Fig. 1. Resistivity of metal thin films as a function of film thickness. Mo data are from ^[7] for as-deposited PVD films on SiO$_2$; Ru and Cu data are from ^[8] for PVD and annealed films; W data are from ^[4].

II. CALCULATIONAL DETAILS

Fig. 2 shows Mo(CO)${}_{6}$ structures in solid and gaseous phases. Mo, C, and O atoms are shown in purple, brown, and red, respectively. Solid Mo(CO)${}_{6}$ (space group Pnma) has a unit cell of $6.44\times 11.37 \times 11.97$ Å$^{3}$ with 52 atoms (Mo: 4, C: 24, O: 24), and gaseous Mo(CO)${}_{6}$ (point group $O_{h}$) is modeled in a $15 \times 15 \times 15$ Å$^{3}$ cubic cell.

Fig. 2. Structure of Mo(CO)$_6$: (a) solid and (b) gaseous structures. Solid Mo(CO)$_6$ (space group $Pnma$) has a unit cell containing 52 atoms, while gaseous Mo(CO)$_6$ (point group $O_h$) consists of one Mo atom and six CO ligands. The purple, brown, red colors correspond to molybdenum, carbon, and oxygen atoms, respectively.

All DFT calculations were carried out using the Vienna Ab-initio Simulation Package (VASP) code ^[21-23]. The electron wave functions ^[24] were treated with the projector augmented wave (PAW) method developed by Kresse and Joubert, and the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) was employed for the exchange-correlation energy ^[25]. K-point sampling was performed using $3 \times 2 \times 2$ Monkhorst-Pack grid for solid Mo(CO)${}_{6}$ and $1 \times 1 \times 1$ grid for gaseous phase ^[26]. The calculational parameters for the solid phase were obtained from Materials Project database ^[27]. Convergence criteria of $10^{-6}$ eV and $10^{-2}$ eV/Å were applied for electronic and ionic optimizations, respectively ^[28-30].

The Gibbs free energy ($G$) was computed using DFT methods in conjunction with the equations detailed in reference ^[9]. The expression for Gibbs free energy is given by

where $H$ is enthalpy, $T$ is temperature, $S$ is entropy, $U$ is internal energy, $P$ is pressure, and $V$ is volume. For a solid phase, $PV = 0$, while for a gas phase, $PV = k_{B} T$. Here, ${{E}}_{\rm{elec}}$ represents the electronic energy from DFT at 0 K, $ZPE$ is the zero-point energy, $\int{{{C}}_{\rm{p}}{dT}}$ accounts for the enthalpy change due to temperature, and $S_{\rm t}$, $S_{\rm r}$ and $S_{\rm v}$ represent translational, rotational, and vibrational entropies, respectively. $ZPE$, $\int{{{C}}_{\rm{p}}{dT}}$, and $S_{\rm v}$ were calculated using VASPKIT ^[31-33]. The terms $S_{\rm t}$, $S_{\rm r}$, and $S_{\rm v}$ were derived using the Sackur-Tetrode equation, the rigid-rotor model, and harmonic oscillator approximation, respectively. The vibrational degrees of freedom for both solid and molecular phases are defined as $3N-6$, where $N$ is the number of atoms. To avoid any anomalous corrections, vibrational frequencies below $50$ cm$^{-1}$ were adjusted to $50$ cm$^{-1}$. Given that translational and rotational entropies are negligible in the solid phase, only vibrational contributions were considered in the analysis. In this study, the vapor pressure of Mo(CO)${}_{6}$ was calculated by equating the Gibbs free energies of the solid and gaseous phases.

This relationship allows the determination of the pressure $P$ at a given $T$, which corresponds to the vapor pressure. The Gibbs free energy for the gas phase includes the translational entropy contribution $S_{\rm t}$, which accounts for $P$. For an ideal gas, the translational entropy $S_{\rm trans}$ is given by

where $k_{B}$ is the Boltzmann constant, $m$ is the molecular mass, and $h$ is Planck's constant.

For the calculation of BDE of Mo(CO)${}_{6}$ at 298 K (standard room temperature condition), $D_{298}$(Mo(CO)${}_{6}$), we have incorporated the thermodynamic correction by adding the integrated heat capacity term to the BDE at 0 K, $D_{0}$(Mo(CO)${}_{6}$), as expressed in the following equations:

where $C_{P}$ denotes the heat capacity term for each molecule, ${{E}}_0{=}{{E}}_{{elec}}+ {ZPE}$ ^[34].

III. RESULTS AND DISCUSSION

Fig. 3 presents the relative Gibbs free energies of Mo(CO)${}_{6}$ in the solid and gaseous phases as a function of temperature, normalized per Mo(CO)${}_{6}$ molecule, considering four molecules per unit cell. The relative Gibbs free energies were calculated using the solid phase Gibbs free energy as a reference state. The gaseous phase Gibbs free energy is plotted under two distinct pressure conditions (1 atm and 10${}^{-3}$ atm). As temperature increases, the Gibbs free energy of the gaseous phase decreases more rapidly than that of the solid phase due to the gaseous phase's higher entropy. The intersection between solid and gaseous Gibbs free energy curves corresponds to the sublimation temperature of Mo(CO)${}_{6}$, which was calculated to be 385 K at 1 atm and 296 K at 10${}^{-3}$ atm.

Fig. 3. Relative Gibbs free energy of Mo(CO)$_6$ in the solid and gaseous phases as a function of temperature. Red dotted and yellow dashed curves represent the gaseous phase at 1 atm and $10^{-3}$ atm, respectively, while the black solid curve indicates the solid phase reference. Sublimation temperatures are identified by the intersection points: 385 K at 1 atm and 296 K at $10^{-3}$ atm.

Fig. 4. Vapor pressure of Mo(CO)$_6$ as a function of temperature on logarithmic scale. The solid line indicates the sublimation curve predicted by the ab-initio thermodynamics. Vertical dotted lines represent the theoretical sublimation temperature (385 K) and experimental value (424 K). Experimental data are taken from ^[14-16] for comparison.

Fig. 4 shows the vapor pressure of Mo(CO)${}_{6}$ as a function of temperature on logarithmic scale. Vapor pressures were obtained from ab-initio thermodynamic calculations by equating Gibbs free energy of solid and gas phases, as described in the calculation details. The theoretical sublimation temperature at atmospheric pressure (1 atm) is 385 K, whereas the experimentally reported value is approximately 424 K ^[15,16]. The vapor pressure calculated using ab-initio thermodynamics is consistently higher than the experimentally measured values across the examined temperature range. This discrepancy between experimental and theoretical vapor pressures of Mo(CO)${}_{6}$ can be reasonably explained by considering the evaporation coefficient, as reported for structurally analogous volatile oxides such as RuO${}_{4}$ ^[35]. In the case of RuO${}_{4}$, experimental vapor pressure was lower than theoretical predictions and quantitatively attributed to a low evaporation coefficient of approximately 0.027, requiring correction by a factor of ${\sim}$44 to match calorimetric enthalpy data. This behavior originates from limited vaporization kinetics, likely influenced by surface interactions and molecular complexity. Despite this deviation, the theoretical vapor pressure data remain valuable for preliminary screening of potential precursors, offering a cost- and time-efficient alternative to real experimental measurements.

Table 1 shows bond length and BDE of Mo(CO)${}_{6}$ that are in good agreement with literature ^[39-43]. The computed BDE of the Mo--CO bond (1.74 eV) is significantly lower than that of the C--O bond (11.87 eV), indicating that the initial step in the adsorption mechanism is likely characterized by dissociative chemisorption, in which a CO ligand detaches from Mo(CO)${}_{6}$ rather than undergoing decomposition into C and O atoms. This facilitates the movement of the Mo center toward the oxygen-rich sites on the SiO${}_{2}(111)$ surface, thereby promoting the formation of stable precursor--substrate interactions.

To simulate realistic surface conditions, this study employed a bare SiO${}_{2}(111)$ surface instead of a hydrogen-terminated one, under the assumption that surface hydrogen species get removed through thermal annealing. While direct experimental evidence specifically addressing thermal dehydroxylation of crystalline SiO${}_{2}(111)$ surfaces remains limited, this assumption is indirectly supported by a previous study demonstrating successful removal of surface hydrogen species from amorphous silica surfaces via high-temperature annealing ^[36]. Moreover, crystalline $\beta$-cristobalite SiO${}_{2}(111)$ closely resembles amorphous silica in structural characteristics such as density, bond lengths, and bond angles ^[37], making it a widely accepted structural substitute in computational studies ^[9,38]. We acknowledge, however, that this dehydroxylated model is an approximation and may not fully represent actual ALD conditions, where hydroxyl coverage is influenced by temperature, exposure time, and other gas-phase interactions. Additional experimental work under realistic ALD conditions would therefore be valuable for substantiating and refining this modeling assumption.

Fig. 5 illustrates two desorption pathways which were considered for the interaction between Mo(CO)${}_{6}$ and the SiO${}_{2}$ surface and the Gibbs free energy for both reaction condition was calculated:

where the dotted line ($\cdots $) indicates that the surface O's are coordinated to Mo in Mo(CO)${}_{5}$. In the first pathway (Fig. 5(a)), bond dissociation of one of the six Mo-CO bonds leads to the loss of metal-to-ligand $\pi$-back-donation, thereby strengthening the C-O bond. Thus, we observed a reduction in bond length from 1.16 Å to 1.14 Å. While in the second pathway (Fig. 5(b)), the desorbed CO ligand subsequently reacts with one of the surface oxygen atoms from the substrate to form CO${}_{2}$, which then desorbs. In order to evaluate the thermodynamic preference between these two reaction pathways, the relative Gibbs free energies of the optimized configurations were compared, resulting in the CO$_{2}$ desorption configuration (Fig. 5(b)) to be 2.52 eV lower than the CO desorption configuration (Fig. 5(a)) under the experimental conditions of 435 K and $1.31 \times 10^{-8}$ atm ^[19]. This substantial energy difference strongly supports that CO${}_{2}$ desorption is more thermodynamically favorable than CO desorption. These findings clearly demonstrate that surface oxygen atoms actively participate in the oxidation of CO ligands, thereby facilitating the formation and desorption of CO${}_{2}$ during the initial decomposition of Mo(CO)${}_{6}$ on the SiO${}_{2}$ surface. This result is consistent with previous experimental studies reporting CO${}_{2}$ as the primary gaseous byproduct during the initial decomposition of Mo(CO)${}_{6}$ on SiO${}_{2}$ ^[19,20].

Fig. 5. Mo(CO)$_6$ reaction on the SiO$_2(111)$ surface: (a) CO desorption and (b) CO$_2$ desorption. Shaded regions indicate the desorbed species. Configuration (b) is 2.52 eV more stable than (a), indicating a thermodynamic preference for CO$_2$ formation. The CO$_2$ desorption shown in (b) occurs under the experimental condition of $T = 435$ K and $P = 1.31 \times 10^{-8}

IV. CONCLUSION

In this study, the thermodynamic properties and reaction pathways of the Mo(CO)${}_{6}$ precursor interacting with the SiO${}_{2}(111)$ surface were analyzed using DFT-based ab-initio thermodynamics. Vapor pressure of Mo(CO)${}_{6}$ was determined to be 385 K at 1 atm. Thermodynamic comparisons of ligand desorption pathways revealed that CO${}_{2}$ formation, induced by surface oxygen, is more favorable than simple CO desorption, highlighting the active role of surface oxygen in promoting precursor decomposition. These results demonstrate that first-principles simulations can serve as a powerful and efficient tool for predicting key physicochemical properties of precursors, offering a cost-effective alternative to complex and resource-intensive experimental procedures.