# Anisotropic electron polaron mobility in V\_2O\_5: from ab initio to kinetic Monte Carlo

Hala Houmsi, Benoit Sklenard, François Triozon, Marc Guillaumont, Jing Li

► To cite this version:

Hala Houmsi, Benoit Sklenard, François Triozon, Marc Guillaumont, Jing Li. Anisotropic electron polaron mobility in V\_2O\_5: from ab initio to kinetic Monte Carlo. SISPAD 2025 - Conference on Simulation of Semiconductor Processes and Devices, Sep 2025, Grenoble, France. 10.1109/SISPAD66650.2025.11186381. cea-05533764

HAL Id: cea-05533764

https://cea.hal.science/cea-05533764v1

Submitted on 2 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Anisotropic electron polaron mobility in V₂O₅: from ab initio to Kinetic Monte Carlo

Hala Houmsi*†‡, Benoit Sklénard†, François Triozon†, Marc Guillaumont*, Jing Li†

*Lynred, 364 Avenue de Valence, 38113, Veurey-Voroize, France

†Univ. Grenoble Alpes, CEA, Leti, F-38000, Grenoble, France

‡hala.houmsi@cea.fr

**Abstract**—In this work, we investigate electron polaron mobility in pristine orthorhombic V₂O₅, providing insights into its anisotropic transport behavior: All hopping pathways in this material are investigated by using DFT+U. The results obtained are employed to calculate the corresponding transfer rate coefficients and the resulting mobility through Landau-Zener theory and Kinetic Monte Carlo (KMC) simulations. The KMC simulations allow us to obtain temperature and direction-dependent mobility, and the results confirm that hopping is favorable along the [010] direction. An evaluation of the assumptions and computational steps used to relate Landau-Zener theory to charge carrier mobility is presented, accompanied by a comparison with experimental data reported in the literature.

**Index Terms**—Polaron hopping, Landau-Zener rate, KMC, anisotropic mobility, ab initio.

## I. INTRODUCTION

V₂O₅ is an orthorhombic layered semiconductor (space group Pmmn) composed of weakly van der Waals bonded layers. Due to its high electron affinity and rich chemical properties, V₂O₅ has gained interest in energy storage, catalysis, and optoelectronics [1]. In an intermediate temperature range, its charge transport is governed by small polaron hopping, where charge carriers localize and move among vanadium sites via thermally activated hopping, leading to a strong temperature dependence of mobility [2]–[4]. In addition, it is characterized by an anisotropic transport behavior [5]–[8].

First-principles calculations enable the determination of transfer rates for each hopping path using Landau-Zener (LZ) equation [9], [10]. This approach is based on the Marcus–Emin–Holstein–Austin–Mott (MEHAM) formalism and has been applied to study small polaron hopping in various materials: TiO₂ [11], Fe₂O₃ [12], FePO₄ [13], BiVO₄ [14] and V₂O₅ [15]. The polaron transport mobility can then be evaluated by Kinetic Monte Carlo (KMC) simulations. In Ref. [15], a simplified expression for the transfer rate, based solely on activation energy was used, and the analysis was limited to four polaron hopping paths.

The aim of this work is to investigate all possible hopping pathways along with providing a critical examination of the assumptions and computational steps involved in connecting LZ theory to charge carrier mobility. Building on our previous work [16], this study refines the methodology to improve the accuracy of parameter extraction for LZ theory. This yields

hopping rates that determine the direction-dependent polaron mobility in V₂O₅.

## II. METHODOLOGY

### A. Computational framework

All hopping pathways in V₂O₅ are computed within Density Functional Theory (DFT) framework. To account for the strong electron correlations in the vanadium 3d orbitals, we apply the Hubbard onsite potential U in the DFT calculation (DFT+U) [17].

We employ the Projector Augmented Wave method (PAW) [18] with plane-waves basis set, as implemented in the Vienna Ab initio Simulation Package (VASP) [19]. DFT calculations were conducted on a 1×4×3 supercell containing 168 atoms, with a cut off energy of 400 eV. Due to the large cell size, the k-point mesh was restricted to the Γ-point. In the electronic self consistent loop, the total energy is converged within 10⁻⁸ eV. For structural relaxation, the ionic positions were fully optimized, achieving forces on atoms below 0.01 eV/Å. We chose the Generalized Gradient Approximation (GGA) functional developed by Perdew, Burke, and Ernzerhof (PBE) [20].

To stabilize a polaron on a specific vanadium atom, we applied a 10% distortion (expansion) to the neutral bulk structure in order to generate an initial structure. This is done by selectively altering the bond distances between a chosen vanadium atom and its nearest oxygen neighbors within a cut-off radius of 3 Å. Then, we added an electron and relaxed the atomic positions to obtain a polaron state. The same procedure was applied to stabilize the polarons located on neighboring atoms. To determine polaron hopping barriers, intermediate atomic structures (images) along the polaron hopping path were generated using Linear Interpolation (LI) between the initial and final atomic structures. The coordinates of all atoms along the polaron hoping path from atom A to atom B are defined as:

$$R(x) = R_A + x(R_B - R_A) \quad \text{for } 0 \le x \le 1, \quad (1)$$

where x is the reduced Reaction Coordinate RC, R_A the initial and R_B the final coordinates. This approach was selected over the Nudged Elastic Band (NEB) method, which allows structural relaxation of images along the migration path. In
Ref. [15], authors show that NEB gave similar barriers to LI in both [100] and [010] directions in V₂O₅. Additionally, LI offers a more practical balance between computational efficiency and accuracy, particularly given the need to evaluate multiple polaron hopping barriers. DFT+U calculations were performed for each image along the hopping path to access the variation of total energy. The activation energy, denoted by Eₐ, was extracted from the total energy barriers in Fig. 1.

The Hubbard U parameter can either be fixed from first principles [21]–[23] or semi-empirically by fitting experimentally verifiable properties such as the electronic gap [24], or the oxidation enthalpies [25]. It is important to mention that the chosen U influences the extracted LZ parameters as it impacts the degree of polaron localisation in the structure [16]. While various methods exist to compute U, it ultimately remains a tunable parameter within a certain range, introducing an inherent degree of uncertainty in the results. The parameter U was set to 3.5 eV, following the methodology proposed by Falleta et al. [23] for systems with polaronic defects. The procedure used to extract this value is described in detail in our previous work [26].

<figure><img src="image_2.png" /><figcaption>Fig. 1. The polaron migration barriers in different directions: Total energy barriers (dots) computed using DFT+U. Solid lines are included as guides for the eye.</figcaption></figure>

Based on a two site single particle model Hamiltonian, the electronic coupling *t*, was determined from the splitting of the two single particle polaron states at the transition for a neutral system, which has been examined as the most accurate scheme [26]. Extracting a coupling strength free from charging energy, avoids its overestimation and the consequent prediction of overestimated transfer probabilities in the following steps. In fact, the activation energy and coupling strength are the key parameters to compute the transfer rate (Table I).

## B. Transfer rate and mobility

Within the general framework of the MEHAM model, the hopping process is typically described as a one-dimensional motion coupled to a single effective phonon mode $ν_{eff}$. Incorporating the full vibrational spectrum is more accurate

but computationally demanding. Based on Marcus theory, the polaron transfer along the 1D RC is expressed using a total energy two-site Hamiltonian model under harmonic approximation [26], [27]. The transfer rate relies on the LZ term that expresses the probability of diabatic transfer between two harmonic states [9], [10], [28]:

$$k_{\text{et}} = \frac{2P}{1+P} \nu_{\text{eff}} \exp \left( -\frac{E_a}{k_B T} \right) \quad (2)$$

where $k_B$ is the Boltzmann constant and T the temperature. $P = 1 - \exp[-2\pi\gamma]$ represents the probability of an adiabatic transition (remaining on the low energy branch on the potential energy surface), in which the exponential represents the LZ term. The diabacity parameter reads:

$$2\pi\gamma = \frac{\pi^{3/2}t^2}{h\nu_{\text{eff}}\sqrt{\lambda k_B T}} \quad (3)$$

where *h* is the Planck constant and *λ* the reorganization energy. According to Marcus theory, the harmonic approximation and the presence of symmetric barriers in V₂O₅ allow us to simplify the reorganization energy as follows: $λ = 4(Eₐ + t)$. The effective phonon mode $ν_{eff}$ is computed using the approach detailed in Ref. [29]:

$$\nu_{\text{eff}} = \sqrt{\frac{8(E_a + t)}{\sum_{\alpha} m_{\alpha} |\mathbf{r}_{\alpha}^{\text{f}} - \mathbf{r}_{\alpha}^{\text{i}}|^2}} \quad (4)$$

where $m_α$, $\mathbf{r}_α^f$ and $\mathbf{r}_α^i$ are the mass, the position in the final image and initial image of atom α, respectively.

With these transition rates, we employed KMC to evaluate the mobility of polarons in the V₂O₅ hopping network, as described in Ref. [14]. A lattice model that includes all polaron hopping events described in Table I is constructed. At each iteration, a transition from the previous polaron site A to a randomly selected neighboring site B is performed, with a probability weighted by the transfer rate of the transition. The associated time increment is computed, and the simulation proceeds from site B. This approach enables the propagation of a single polaron through the lattice, assuming the absence of significant polaron-polaron interaction. 10⁵ simulations were conducted, with each simulation consisting of 10³ hopping events. This allows the calculation of the polaron's mean squared displacement L²(t) as a function of time (see Fig. 2). The diffusion coefficient along each directions follows:

$$D_i = \frac{L_i^2}{2t}; \quad i = x, y, z \quad (5)$$

where *t* is the time of the simulation. The mobility is computed using Einstein-Smoluchowski equation:

$$\mu_i = \frac{qD_i}{k_B T}; \quad i = x, y, z \quad (6)$$

where *q* is the absolute value of the carrier charge.
<table><caption>TABLE I
LANDAU-ZENER PARAMETERS FOR DIFFERENT JUMPS IN V₂O₅ AT
T = 300 K: Eₐ IS THE ACTIVATION ENERGY, t THE ELECTRONIC
COUPLING, νₐ₇₅ THE ATTEMPT FREQUENCY OF A BARRIER CROSSING, P
THE PROBABILITY OF AN ADIABATIC JUMP, AND kₐ₇₅ THE TRANSFER RATE.</caption><thead><tr><th>Jump</th><th>E<sub>a</sub> [meV]</th><th>t [meV]</th><th>hν<sub>eff</sub> [meV]</th><th>P</th><th>k<sub>et</sub> [Hz]</th></tr></thead><tbody><tr><td>[100]</td><td>96</td><td>150</td><td>245</td><td>0.96</td><td>1.42 × 10<sup>12</sup></td></tr><tr><td>[010]</td><td>224</td><td>105</td><td>208</td><td>0.80</td><td>7.69 × 10<sup>9</sup></td></tr><tr><td>[001]</td><td>269</td><td>20</td><td>139</td><td>0.09</td><td>1.65 × 10<sup>8</sup></td></tr><tr><td>[110]</td><td>214</td><td>7.6</td><td>181</td><td>0.01</td><td>2.57 × 10<sup>8</sup></td></tr><tr><td>[101]</td><td>271</td><td>22</td><td>145</td><td>0.10</td><td>1.81 × 10<sup>8</sup></td></tr><tr><td>[110]</td><td>241</td><td>99</td><td>207</td><td>0.75</td><td>3.85 × 10<sup>9</sup></td></tr><tr><td>[111]</td><td>258</td><td>8</td><td>148</td><td>0.01</td><td>4.71 × 10<sup>7</sup></td></tr><tr><td>[011]</td><td>275</td><td>19</td><td>154</td><td>0.07</td><td>1.20 × 10<sup>8</sup></td></tr></tbody></table>

### III. RESULTS

The activation energy, coupling strength, attempt frequency, probability of an adiabatic jump and transfer rate for each jump at T = 300 K are summarized in Table I. We notice the presence of both adiabatic and diabatic jumps in V₂O₅ (high and low P values). The energy barriers of the four first jumps in the table are in excellent agreement with those computed in Ref. [15].

Fig. 3 displays the three paths with the highest transfer probabilities in the x-y plane indicated by orange, blue and gold arrows. The corresponding potential energy surfaces are shown in the same colour in Fig. 1. The hopping along x direction, involves two consecutive jumps, the [100] jump (1.42 × 10¹² Hz, orange arrow) and the [110] (2.57 × 10⁸ Hz, gray arrow) jump. The later one limits the polaron hopping in x direction. This behaviour is reflected in Fig. 2 where the fluctuations of L²ₓ observed at short times is related to fast "back and forth" oscillations along the orange arrow. The hopping along z direction, i.e. inter-layer hopping is weak, as indicated by the [001] jump (1.65 × 10⁸ Hz). The remaining two favorable jumps are [010] (7.69 × 10⁹ Hz, gold arrow) and [110] (3.85 × 10⁹ Hz, equivalent to blue arrow), resulting in an intra-layer transport along y direction. The qualitative y-direction favorable transport is recovered and agrees with previous findings from both experimental [2], [30] and theoretical studies [15].

The mobility computed in the y direction at 300 K is equal to 5.94 × 10⁻⁴ cm²/(V·s) (see Fig. 4) and is around two times bigger than that reported in Ref. [15] (2.7 × 10⁻⁴ cm²/(V·s)). This discrepancy results from the lack of the adiabatic hopping event in the [110] direction in the lattice model in Ref. [15], which impedes mobility in the x-y plane. Additionally, the omission of coupling strength and the effective frequency for a barrier crossing influence the computed adiabatic transfer probability and thus transfer rates. The experimental mobilities reported in the literature are highly scattered and span over several orders of magnitude as can be seen in Fig. 4. In Ref [3], a mobility of 1.8 × 10⁻² cm²/V·s is obtained in y direction from DC electrical conductivity measurements on a V₂O₅ single crystal. Carrier density is taken equal to spin con-

<figure><img src="image_3.png" /><figcaption>Fig. 2. Variation of L²ₓ, i = x, y, z vs. time showing that the highest diffusion coefficient is obtained in the y direction. The fluctuation along x at short time is related to fast "back and forth" oscillations along the orange arrow.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. V₂O₅ intra-layer hopping showing that the low probability of transfer of some jumps along the [110] direction reduces the mobility in the x direction (red cross on the grey arrow): oxygen atoms are shown in red and vanadium atoms in blue. The charge density isosurface of the electron polaron localized on a vanadium atom is shown in yellow.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Variation of the computed mobilities for different temperatures in the y direction. A macroscopic activation energy of 0.202 eV is extracted from the slope of the curve. Values of mobility measured from experimental works are also reported. In Ref [31], different mobilities are reported. One is measured from Hall effect at 293 K (green square) and several are measured from conductivity and ESR at various temperatures (green crosses).</figcaption></figure>
centration measured by double integration of the Electron Spin Resonance (ESR) spectrum at 300 K (2 × 10¹⁸ spin/cm³). The studied material is characterized by an *n*-type semiconductor behavior due to non-stoichiometry, where oxygen vacancies are compensated by vanadium ions in a lower oxidation state. The presence of such defects was not considered in our simulations, which can justify the discrepancy.

Hall effect measurements were performed in Ref. [31] at 293 K under a magnetic field of 2.2 T on different samples of 99.9% pure V₂O₅ single crystals. This allowed the extraction of a Hall mobility in the *y* direction of 3 × 10⁻² cm²/(V·s). This extracted value can depend on crystal purity, and on the strength of the applied magnetic field. In the same work, those samples were subject to conductivity measurements as function of temperature. Carrier density was again estimated from the concentration of ESR centers (3 × 10¹⁸ to 10¹⁹ spin/cm³). The extracted mobility is around 1.25 × 10⁻³ cm²/(V·s) at 300 K. It is also reported with green crosses for other temperatures in Fig. 4.

The measured mobility in V₂O₅ is quite low and lies between three orders of magnitude [10⁻³, 1] cm²/V·s at 300 K [3], [30], [31]. Most of the discrepancies reported in the literature can be attributed to variations in the sample types, such as thin films fabricated using different methods, single crystals, or differences in purity, as well as the characterization techniques employed and the approximations made during the extraction of mobility values. In the present study, we consider pristine V₂O₅ with a single self-trapped polaron, starting from a defect-free and stoichiometric structure. Our mobility in the *y* direction might be underestimated due to the absence of defects in our model. In real materials, oxygen vacancies, dislocations, and grain boundaries can contribute to transport via defect-mediated hopping, or increase the carrier density [15], leading to a higher experimental mobility compared to the one simulated for the perfect system.

Fig. 4 shows the variation of the computed mobility for different temperatures in the *y* direction. This allows the extraction of a macroscopic activation energy of 0.202 eV. The experimental values reported in the literature, measured from the slope of conductivity versus temperature on various samples, lie in the range [0.17, 0.27] eV [3], [30], [31]. Those values remain reasonably close to the one computed in this work considering the different approximations employed in experiments and simulations.

In perspective, it would be of interest to investigate the sensitivity of the LZ parameters with respect to the U correction. This could enable us to extract the range of variation in charge carrier mobility as a function of the range of *U* values and to subsequently assess its physical plausibility.

## IV. CONCLUSION

We investigated electron polaron transfer in orthorhombic V₂O₅ at room temperature by outlining a methodology that bridges *ab initio* calculations to mobility via the generalized

Landau-Zener transfer rates and Kinetic Monte Carlo simulations. The anisotropy of the extracted directional mobility is demonstrated and aligns with previous literature. The computed mobility is challenging to compare with experimental data due to the assumption of a pristine system. Experimentally, assigning a definitive value for the electronic mobility of V₂O₅ is further complicated by the dependence of the extracted values on the specific fabrication process, characterization techniques, and underlying assumptions employed. The computed activation energy remains in good agreement with the ones reported in experimental works. A robustness study of the simulation could also be carried out in order to evaluate the influence of *U* correction on the resulting charge carrier mobility.

## ACKNOWLEDGMENT

This project is financed by the French government France 2030 investment program and the french national agency of research (ANR-CIFRE No. 2023/0835). Part of the calculations were performed on computational resources provided by GENCI-IDRIS (Grant No. 2025-A0170912036).

## REFERENCES

[1] P. Hu et al., Chem. Rev., vol. 123, no. 8, mar 2023.

[2] J. Haemers et al., Phys. Status Solidi, vol. 20, no. 1, 1973.

[3] C. Sanchez et al., J. Phys. C: Solid State Phys., vol. 15, no. 35, dec 1982.

[4] P. Watthaisong et al., RSC Adv., vol. 9, no. 34, jun 2019.

[5] A. Yosikawa et al., J. Mater. Sci. Lett., vol. 16, no. 8, jan 1997.

[6] C. Lamsal et al., JOM, vol. 70, no. 4, apr 2018.

[7] P. Falun et al., J. Phys. Chem. C, vol. 128, no. 26, jul 2024.

[8] S. Sucharitakul et al., ACS Appl. Mater Interfaces, vol. 9, no. 28, jul 2017.

[9] B. S. Brunschwig et al., J. Am. Chem. Soc., vol. 102, no. 18, aug 1980.

[10] M. D. Newton et al., Annu. Rev. Phys. Chem., vol. 35, no. Volume 35., oct 1984.

[11] G. Palermo et al., Phys. Rev. B, vol. 110, no. 23, dec 2024.

[12] N. Iordanova et al., The Journal of Chemical Physics, vol. 122, no. 14, Apr. 2005.

[13] T. Chen et al., Physical Chemistry Chemical Physics, vol. 25, no. 12, Mar. 2023.

[14] F. Wu et al., J. Mater. Chem. A, vol. 6, no. 41, oct 2018.

[15] L. Ngamwongwan et al., Phys. Chem. Chem. Phys., vol. 23, no. 19, may 2021.

[16] R. Defrance et al., Solid-State Electron., vol. 198, dec 2022.

[17] V. Anisimov et al., J. Physics: Condens. Matter, vol. 9, no. 4, jan 1997.

[18] P. E. Blöchl et al., Physical Review B, vol. 50, no. 24, Dec. 1994.

[19] J. Hafner et al., Journal of Computational Chemistry, vol. 29, no. 13, 2008.

[20] J. P. Perdew et al., Physical Review Letters, vol. 77, no. 18, Oct. 1996.

[21] M. Cococcioni et al., Physical Review B, vol. 71, no. 3, Jan. 2005.

[22] F. Aryasetiawan et al., Physical Review B, vol. 74, no. 12, Sep. 2006.

[23] S. Falletta et al., Npj Comput. Mater., vol. 8, no. 1, dec 2022.

[24] K. Dhaked et al., International Journal of Modern Physics B, vol. 39, no. 10, Apr. 2025.

[25] L. Wang et al., Physical Review B, vol. 73, no. 19, May 2006.

[26] H. Houmsi et al., Physical Review B, vol. 111, no. 16, Apr. 2025.

[27] R. A. Marcus and N. Sutin, Biochimica et Biophysica Acta (BBA) - Reviews on Bioenergetics, vol. 811, no. 3, Aug. 1985.

[28] Z. C., Proc. R. Soc. Lond. A, vol. 137, no. 833, 1932.

[29] A. Alkauskas et al., Phys. Rev. Lett., vol. 109, no. 26, dec 2012.

[30] J. C. Badot et al., J. Mater. Chem., vol. 14, no. 23, nov 2004.

[31] V. A. Ioffe et al., physica status solidi (b), vol. 40, no. 1, 1970.