*2.1. Interatomic Potentials*

The interatomic potentials used in this work consist of Buckingham potentials, supplemented by an electrostatic term, as given below:

$$V(r\_{ij}) = \frac{q\_i q\_j}{r\_{ij}} + A\_{ij} \exp\left(\frac{-r\_{ij}}{\rho\_{ij}}\right) - \mathbf{C}\_{ij} r\_{ij}^{-6} \tag{1}$$

This expression shows that for each pair of ions it is necessary to determine three parameters: *Aij*, ρ*ij* and *Cij*, which are constants for each interaction, *qi*, *qj* represent the charges of the ions *i* and *j*, and *rij* is the interatomic distance. The parameters are determined by empirical fitting, and formal charges are used for *qi* and *qj*. The procedure by which potentials were obtained for LiNbO3 is explained in the work of Jackson and Valério [22], and the derivation of the potentials for the vanadium and molybdenum dopants is described in Section 3.1 below. The potentials for LiNbO3 have been the subject of recent studies on the doping of the structure with rare earth ions [23,24], doping with Sc, Cr, Fe and In [25], metal co-doping [26] and doping with Hf [27]. These papers show that modelling can predict the energetically optimal locations of the dopant ions and calculate the energy involved in the doping process. This paper extends this procedure to the study of V2+, V3+, V4<sup>+</sup> and V5<sup>+</sup> as well as Mo3<sup>+</sup>, Mo4<sup>+</sup>, Mo5<sup>+</sup> and Mo6<sup>+</sup> doped lithium niobate, with the aim of establishing the optimal doping site and charge compensation scheme for both sets of ions.

#### *2.2. Defect Formation Energies*

The calculation of defect formation energies is carried out using the Mott–Littleton approximation [28], in which the crystal is divided into two regions: region I, which contains the defect, and region II, which extends from the edge of region I to infinity. In region I, the positions of the ions are adjusted until the resulting force is zero. The radius of region I is selected such that the forces in region II are relatively weak and the relaxation can be treated according to the harmonic response to the defect (a dielectric continuum). An interfacial region IIa is introduced to treat interactions between region I and region II.
