Computer Simulation of the Incorporation of V 2 + , V 3 + , V 4 + , V 5 + and Mo 3 + , Mo 4 + , Mo 5 + , Mo 6 + Dopants in LiNbO 3

: The doping of LiNbO 3 with V 2 + , V 3 + , V 4 + and V 5 + as well as Mo 3 + , Mo 4 + , Mo 5 + and Mo 6 + ions is of interest in enhancing its photorefractive properties. In this paper, possible incorporation mechanisms for these ions in LiNbO 3 are modelled, using a new set of interaction potentials ﬁtted to the oxides VO, V 2 O 3 , VO 2 , V 2 O 5 and to LiMoO 2

Previous work on vanadium and molybdenum doped lithium niobate has included experimental studies of how its photorefractive properties are enhanced by doping with molybdenum ions [16,17] where it is suggested that the Mo 6+ ion dopes at the Nb 5+ site.Another study looks at LiNbO 3 co-doped with Mg and V, concluding that some of the vanadium dopes at the Nb site in the 5+ charge state, but that V 4+  Li , V 3+  Li and V 2+ Li defects are also observed [18].Finally, another recent publication [19] has looked at the photorefractive response of Zn and Mo co-doped LiNbO 3 in the visible region, and concluded that the presence of Mo 6+ ions helps promote fast response and multi-wavelength holographic storage, which is attributed to their occupation of regular niobium sites in the lattice.
In a Density Functional Theory (DFT) study [20], vanadium doping was modelled, and it was concluded that vanadium substitutes at the Li + site as V 4+ , but that it dopes at the Nb site as a neutral defect as the Fermi level is increased.In another DFT study [21], molybdenum doping was modelled and it was concluded that the most stable configuration involves doping at the Nb 5+ site, in agreement with the previously mentioned experimental studies [16,17].It is noted that in the DFT

Interatomic Potentials
The interatomic potentials used in this work consist of Buckingham potentials, supplemented by an electrostatic term, as given below: This expression shows that for each pair of ions it is necessary to determine three parameters: A ij , ρ ij and C ij , which are constants for each interaction, q i , q j represent the charges of the ions i and j, and r ij is the interatomic distance.The parameters are determined by empirical fitting, and formal charges are used for q i and q j .The procedure by which potentials were obtained for LiNbO 3 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 LiNbO 3 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 V 2+ , V 3+ , V 4+ and V 5+ as well as Mo 3+ , Mo 4+ , Mo 5+ and Mo 6+ doped lithium niobate, with the aim of establishing the optimal doping site and charge compensation scheme for both sets of ions.

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.

Derivation of Interatomic Potential Parameters
It was necessary to derive potential parameters for the dopant oxide structures: VO, V 2 O 3, VO 2 and V 2 O 5 as well as LiMoO 2 , Li 2 MoO 3 , Li 3 MoO 4 and Li 2 MoO 4 .For interactions, a new set of potentials was derived empirically by fitting to the observed structures as shown in Table 1.The O 2− -O 2− potential was obtained by Sanders et al. [29] and uses the shell model for O [30], which is a representation of ionic polarisability, in which each ion is represented by a core and a shell, coupled by a harmonic spring, and the Li-O potential was taken from [22].In all cases, the dopant-oxide potentials were obtained by fitting to parent oxide structures.[32] , VO 2 [33] and V 2 O 5 [34] oxides as well as LiMoO 2 [35], Li 2 MoO 3 [36], Li 3 MoO 4 [37] and Li 2 MoO 4 [38] lithium molybdate structures, using the potentials in Table 1.It is seen that the experimental and calculated lattice parameters differ by less than 1%, confirming that the potentials can be used in further simulations of defect properties.The calculations were carried at 0 K (the default for the modelling code and used in most other theoretical studies) and at 293 K for comparison with room temperature results.In this way, we can see how the structure and energies vary with temperature.

Defect Calculations
In this section, calculated energies for dopant ions in LiNbO 3 are reported.The divalent, trivalent, tetravalent, pentavalent and hexavalent dopants can substitute at Li and Nb sites in the LiNbO 3 matrix with charge compensation taking place in a number of ways.The proposed schemes described in the following subsections are written as solid state reactions using the Kroger-Vink notation [39].This notation appears in the tables in Sections 3.2.1-3.2.5 where the dot/bullet (•) means a net positive charge and the dash/prime ( ) means a net negative charge.

Divalent Dopants
The substitution of the divalent dopant V 2+ in the Li + and Nb 5+ host sites requires a charge-compensating defect, which can involve Li and Nb vacancies, Nb Li anti-sites, interstitial oxygen, self-compensation and oxygen vacancies.The modes of substitution considered for divalent cations are shown in Table 3.The solution energies for the divalent (V 2+ ) dopant with different charge-compensating mechanisms were evaluated and plotted as a function of the reaction schemes.Based on the lowest energy value, it seems that the incorporation of a divalent (V 2+ ) ion is energetically favourable at the lithium and niobium sites, taking into account the first in relation to the c axis.In schemes (i) and (iv), the energy difference in eV is small at both temperatures in the first neighbours, indicating that it can be incorporated at the lithium site compensated by a lithium vacancy as well as by self-compensation as shown in Figure 1.This can be attributed to the similarity between the ionic radius of V 2+ , which is 0.79 Å, and those of the Li + site, which varies between 0.59 and 0.74 Å, and the Nb 5+ site, which varies between 0.32 and 0.71 Å [40].

Divalent Dopants
The substitution of the divalent dopant V 2+ in the Li + and Nb 5+ host sites requires a charge-compensating defect, which can involve Li and Nb vacancies, NbLi anti-sites, interstitial oxygen, self-compensation and oxygen vacancies.The modes of substitution considered for divalent cations are shown in Table 3.
Table 3. Types of defects considered due to M = V 2+ incorporation in LiNbO3.

Charge Compensation Reaction
Li + Lithium Vacancies (i) Lithium Vacancies and Anti-site The solution energies for the divalent (V 2+ ) dopant with different charge-compensating mechanisms were evaluated and plotted as a function of the reaction schemes.Based on the lowest energy value, it seems that the incorporation of a divalent (V 2+ ) ion is energetically favourable at the lithium and niobium sites, taking into account the first in relation to the c axis.In schemes (i) and (iv), the energy difference in eV is small at both temperatures in the first neighbours, indicating that it can be incorporated at the lithium site compensated by a lithium vacancy as well as by self-compensation as shown in Figure 1.This can be attributed to the similarity between the ionic radius of V 2+ , which is 0.79 Å, and those of the Li + site, which varies between 0.59 and 0.74 Å , and the Nb 5+ site, which varies between 0.32 and 0.71 Å [40].

Trivalent Dopants
As with the divalent ion V 2+ , the trivalent V 3+ and Mo 3+ dopants can be incorporated at the lithium and niobium sites in the LiNbO 3 matrix through various schemes as shown in Tables 4 and 5.When these ions are substituted at Li and Nb sites, the extra positive charge can, as noted earlier, be compensated by the creation of vacancies, interstitials, anti-site defects or self-compensation.According to Figures 2 and 3 for the first and second neighbours with respect to the c axis, the trivalent V 3+ and Mo 3+ ions prefer to occupy both the Li and Nb sites according to scheme (iv) which is also observed in other trivalent ions [23][24][25].This can be attributed to the similarity between the ionic radius of V 3+ which is 0.64 Å and Mo 3+ which is 0.67 Å [40] and that of Li + and Nb 5 + .The ionic radius of Li + varies between 0.59 Å and 0.74 Å and Nb 5+ varies from 0.32 Å to 0. 66 Å [40].All these ionic radii are in relation to the coordination sphere with oxygen atoms.

Tetravalent Dopants
Like other divalent and trivalent cations, tetravalent V 4+ and M 4+ dopant ions can also substitute at either the Li + or Nb 5+ sites.When these ions substitute at the Li + and Nb 5+ site charge compensation is required, and various schemes involving vacancies, interstitials, anti-sites and self-compensation are adopted, as shown in Tables 6 and 7.
Table 6.Types of defects considered due to M = V 4+ incorporation in LiNbO3.

Charge Compensation Reaction
Li + Lithium Vacancies (i)

Tetravalent Dopants
Like other divalent and trivalent cations, tetravalent V 4+ and M 4+ dopant ions can also substitute at either the Li + or Nb 5+ sites.When these ions substitute at the Li + and Nb 5+ site charge compensation is required, and various schemes involving vacancies, interstitials, anti-sites and self-compensation are adopted, as shown in Tables 6 and 7.The results obtained from these calculations are given in Figures 4 and 5.By inspecting these figures, it can be seen that the tetravalent cation V 4+ prefers to be incorporated at the Li + and Nb 5+ sites through scheme (iv), while the Mo 4+ ion prefers to be incorporated at the niobium site compensated by an oxygen vacancy according to scheme (ix).Similar to the divalent and trivalent dopants, this preference is related to the proximity with the ionic radii of Li + and Nb 5+ .
The results obtained from these calculations are given in Figures 4 and 5.By inspecting these figures, it can be seen that the tetravalent cation V 4+ prefers to be incorporated at the Li + and Nb 5+ sites through scheme (iv), while the Mo 4+ ion prefers to be incorporated at the niobium site compensated by an oxygen vacancy according to scheme (ix).Similar to the divalent and trivalent dopants, this preference is related to the proximity with the ionic radii of Li + and Nb 5+ .

Pentavalent Dopants
For the pentavalent dopants V 5+ and Mo 5+ , no charge compensation is required for the substitution at the Nb 5+ host site, but it is required when the substitution is at the Li + host site, as shown in Tables 8 and 9.

Charge Compensation Reaction
Li + Lithium Vacancies (i)

Pentavalent Dopants
For the pentavalent dopants V 5+ and Mo 5+ , no charge compensation is required for the substitution at the Nb 5+ host site, but it is required when the substitution is at the Li + host site, as shown in Tables 8 and 9.
Table 8.Types of defects considered due to M = V 5+ incorporation in LiNbO 3 .The solution energies for the pentavalent (V 5+ ) and (Mo 5+ ) dopants with different charge compensation mechanisms were evaluated and plotted as a function of the reaction scheme.Based on the lowest energy value, it seems that the incorporation of pentavalent (V 5+ ) and (Mo 5+ ) ions at an Nb site is energetically more favourable than at an Li site, according to scheme (iv) as shown in Figures 6  and 7 at temperatures 0 K and 293 K.This can be attributed to the similarity between the charge of the V 5+ and Mo 5+ ions and the Nb 5+ host, which can contribute to a small deformation in the lattice and consequently a lower solution energy.Experimental results by Kong et al. [17] and Tian et al. [16] show that substitution occurs at the Nb 5+ site.

Hexavalent Dopants
For the hexavalent dopant Mo 6+ , as with the pentavalent ions, there is no self-compensation mechanism and charge compensation schemes are possible when replacing Li and Nb in the LiNbO3 matrix as shown in Table 10.
Table 10.Types of defects considered due to Mo 6+ incorporation in LiNbO3.

Charge Compensation Reaction
Li + Lithium Vacancies (i) Solution energy per dopant ion(eV)

Hexavalent Dopants
For the hexavalent dopant Mo 6+ , as with the pentavalent ions, there is no self-compensation mechanism and charge compensation schemes are possible when replacing Li and Nb in the LiNbO3 matrix as shown in Table 10.
Table 10.Types of defects considered due to Mo 6+ incorporation in LiNbO3.

Charge Compensation Reaction
Li + Lithium Vacancies (i) The solution energies for the hexavalent (Mo 6+ ) dopants with different charge-compensation mechanisms were evaluated and plotted as a function of the reaction scheme.Based on the lowest

Hexavalent Dopants
For the hexavalent dopant Mo 6+ , as with the pentavalent ions, there is no self-compensation mechanism and charge compensation schemes are possible when replacing Li and Nb in the LiNbO 3 matrix as shown in Table 10.The solution energies for the hexavalent (Mo 6+ ) dopants with different charge-compensation mechanisms were evaluated and plotted as a function of the reaction scheme.Based on the lowest Crystals 2020, 10, 457 9 of 12 energy value, it seems that the incorporation of hexavalent (Mo 6+ ) ions at an Nb site is energetically more favourable than at an Li site, according to scheme (iv) as shown in Figure 8 at temperatures 0 K and 293 K.This can be attributed to the similarity between the ionic radii of Mo 6+ ions and the Nb 5+ host site (0.32-0.71Å) [40].The ionic radii of Mo 6+ , taking into account the coordination number, vary between 0.42 and 0.67 Å [40], and the small difference between the Mo 6+ dopant ions and Nb 5+ ions can contribute to a small deformation in the lattice and consequently a lower solution energy.This result reveals that global trends of dopant solution energies are controlled by the combination of dopant ion size [40] and its electrostatic interactions, demonstrating that there is a relation between the energetically preferred site and the types of defect mechanisms involved in the doping process.Experimental results from Kong et al. [17] and Zhu et al. [41] show that substitution occurs at the Nb 5+ site.
Crystals 2020, 10, x FOR PEER REVIEW 10 of 13 energy value, it seems that the incorporation of hexavalent (Mo 6+ ) ions at an Nb site is energetically more favourable than at an Li site, according to scheme (iv) as shown in Figure 8 at temperatures 0 K and 293 K.This can be attributed to the similarity between the ionic radii of Mo 6+ ions and the Nb 5+ host site (0.32-0.71Å) [40].The ionic radii of Mo 6+ , taking into account the coordination number, vary between 0.42 and 0.67 Å [40], and the small difference between the Mo 6+ dopant ions and Nb 5+ ions can contribute to a small deformation in the lattice and consequently a lower solution energy.This result reveals that global trends of dopant solution energies are controlled by the combination of dopant ion size [40] and its electrostatic interactions, demonstrating that there is a relation between the energetically preferred site and the types of defect mechanisms involved in the doping process.Experimental results from Kong et al. [17] and Zhu et al. [41] show that substitution occurs at the Nb 5+ site.In all cases, the energy involved in doping was obtained by calculating the solution energy, which includes all terms of the thermodynamic cycle involved in the solution process.For example, the solution energy, Esol, corresponding to the incorporation of V 2+ at the Li + site (second equation in Table 3) is given by: where the Elatt and EDef terms are lattice energies and defect energy.All energies were normalised by the number of dopants, i.e., the solution energy is divided by the number of dopants involved.For example, for scheme (ii) of Table 3, the energy must be divided by five, since five lithium sites are occupied.This is done because the number of dopants varies for each mechanism.Lattice energies, Elatt, required to calculate the solution energies are given in Table 11.In all cases, the energy involved in doping was obtained by calculating the solution energy, which includes all terms of the thermodynamic cycle involved in the solution process.For example, the solution energy, E sol , corresponding to the incorporation of V 2+ at the Li + site (second equation in Table 3) is given by: where the E latt and E Def terms are lattice energies and defect energy.All energies were normalised by the number of dopants, i.e., the solution energy is divided by the number of dopants involved.For example, for scheme (ii) of Table 3, the energy must be divided by five, since five lithium sites are occupied.This is done because the number of dopants varies for each mechanism.Lattice energies, E latt , required to calculate the solution energies are given in Table 11.In this sub-section, the results presented in the last five subsections are summarised.Divalent dopants: the calculations predict that, for V 2+ , self-compensation (simultaneous doping at lithium and niobium sites) and doping at the lithium site with lithium vacancy compensation are most likely.It is noted that V 2+  Li defects have been observed experimentally [18].Trivalent dopants: both V 3+ and Mo 3+ ions are predicted to self-compensate.Experimental data from [18] support V 3+ doping at the lithium site, as with V 2+ .
Tetravalent dopants: here, different behaviour is predicted for vanadium and molybdenum.V 4+ is predicted to self-compensate, while Mo 4+ is predicted to occupy a niobium site with oxygen vacancy charge compensation.Again, [18] suggests that V 4+ can dope at a lithium site.
Pentavalent dopants: both V 5+ and Mo 5+ are predicted to dope at the niobium site (no charge compensation is needed), agreeing with experimental results [16,17].
Hexavalent dopants: Mo 6+ is predicted to dope at the niobium site, with charge compensation by lithium vacancy formation.The occupation of the niobium site is supported by experimental data [16,17,19].

Conclusions
This paper has presented a computational study of VO, V 2 O 3 , VO 2 and V 2 O 5 as well as LiMoO It was found that divalent (V 2+ ), trivalent (V 3+ , Mo 3+ ) and tetravalent (V 4+ ) ions are more favourably incorporated at the Li and Nb sites through the self-compensation mechanism.The tetravalent (Mo 4+ ) ion is more favourably incorporated at the niobium site, compensated by an oxygen vacancy.The pentavalent ions (V 5+ , Mo 5+ ) and hexavalent (Mo 6+ ) ion are more favourably incorporated at the Nb site, and the lowest energy schemes involve, respectively, no charge compensation, and for the Mo 6+ ion, charge compensation with lithium vacancy.This is shown to be consistent with some experimental data, although future calculations involving finite V 5+ and Mo 6+ concentrations will be carried out to investigate this further.
Finally, to summarise, in this paper we have looked in detail at vanadium and molybdenum dopants in various charge states in LiNbO 3 , and through the use of solution energies, identified the energetically favoured sites and charge compensation mechanisms, while comparing the results with available experimental and theoretical work in this field.

Figure 1 .
Figure 1.Bar chart of solution energies vs. solution schemes for divalent dopant (V 2+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 1 .
Figure 1.Bar chart of solution energies vs. solution schemes for divalent dopant (V 2+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

13 Figure 2 .
Figure 2. Bar chart of solution energies vs. solution schemes for trivalent dopant (V 3+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 2 .
Figure 2. Bar chart of solution energies vs. solution schemes for trivalent dopant (V 3+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 2 .
Figure 2. Bar chart of solution energies vs. solution schemes for trivalent dopant (V 3+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 3 .
Figure 3. Bar chart of solution energies vs. solution schemes for trivalent dopant (Mo 3+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 3 .
Figure 3. Bar chart of solution energies vs. solution schemes for trivalent dopant (Mo 3+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 4 .
Figure 4. Bar chart of solution energies vs. solution schemes for tetravalent dopant (V 4+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 4 . 13 Figure 5 .
Figure 4. Bar chart of solution energies vs. solution schemes for tetravalent dopant (V 4+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.Crystals 2020, 10, x FOR PEER REVIEW 8 of 13

Figure 5 .
Figure 5. Bar chart of solution energies vs. solution schemes for tetravalent dopant (Mo 4+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Crystals 2020 , 13 Figure 6 .
Figure 6.Bar chart of solution energies vs. solution schemes for pentavalent dopant (V 5+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 7 .
Figure 7. Bar chart of solution energies vs. solution schemes for pentavalent dopant (Mo 5+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 6 . 13 Figure 6 .
Figure 6.Bar chart of solution energies vs. solution schemes for pentavalent dopant (V 5+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 7 .
Figure 7. Bar chart of solution energies vs. solution schemes for pentavalent dopant (Mo 5+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 7 .
Figure 7. Bar chart of solution energies vs. solution schemes for pentavalent dopant (Mo 5+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 8 .
Figure 8. Bar chart of solution energies vs. solution schemes for hexavalent dopant (Mo 6+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Figure 8 .
Figure 8. Bar chart of solution energies vs. solution schemes for hexavalent dopant (Mo 6+ ) at the Li and Nb sites, considering the first neighbours in relation to the c axis.

Table 2
compares experimental and calculated structures of VO [31], V 2 O 3

Table 8 .
Types of defects considered due to M = V 5+ incorporation in LiNbO3.

Table 9 .
Types of defects considered due to Mo 5+ incorporation in LiNbO3.

Table 11 .
Lattice energies used in the solution energy calculations (eV).

Table 11 .
Lattice energies used in the solution energy calculations (eV).Summary of Results for Vanadium and Molybdenum Dopants in LiNbO 3