Multiplicity of Zn coordination sites at cubic spinel ferrites: magnetism and influence of the Zn d band

Zinc substituted nickel ferrite (ZnxNi1 − xFe2O4) is investigated under density functional theory (DFT) within the DFT + U approximation for x ≤ 0.50, with particular interest in understanding the effect of Zn on the net magnetization. Using as a reference ZnFe2O4, the localization of the Zn d band is proved to have a large impact on the preference for Zn to occupy either tetrahedral (ZnA) or octahedral (ZnB) coordination sites, which in ZnFe2O4 is equivalent to the relative stability of the direct and inverse spinel forms. This affects the lattice volume, with ZnA favoring larger lattice expansions. Additional important consequences emerge on the magnetism of the system, as ZnA and ZnB alter the balance of atoms at the magnetic sublattices in a different way: while ZnA enhances the global magnetization by reducing the minority spin contribution, the opposite occurs for ZnB. On the other hand, the dominant magnetic exchange interactions are not significantly altered by Zn independently of its distribution, while the magnetic anisotropy of soft NiFe2O4 is further weakened. Our simulations support the presence of a significant ratio of Zn atoms at octahedral positions at ZnxNi1 − xFe2O4, mainly as the Zn concentration increases, putting limits to the ability to increase the magnetization of NiFe2O4 by Zn substitution.


Introduction
Nanoparticles based on cubic spinel ferrites are key elements for numerous magnetic applications, including rare-earth-free permanent magnets, switching and recording devices, color imaging, magnetic refrigeration or diverse biomedical solutions, from drug delivery to detoxification of fluids [1][2][3][4][5]. Their stability under extreme conditions, together with their robust magnetism, low eddy current losses, abundance and low cost confer spinel ferrites an added value. Among these applications, composite ferrite nanostructures combining soft and hard magnets have gained recent interest as one of the best prospective candidates to replace rare-earth permanent magnets, which are critical raw materials. However, to approach the record efficiency of rareearths, it is still necessary to improve the magnetic energy product of the nanostructures, increasing both the magnetization at remanence of the soft magnet and the coercivity provided by the hard component.
NiFe 2 O 4 (NFO) and ZnFe 2 O 4 (ZFO) are widely used cubic spinel ferrites that behave as soft magnets. They show opposite distribution of divalent and trivalent cations: while NFO is an inverse spinel, with Ni 2? (Fe 3? ) at octahedral (tetrahedral) positions, ZFO is usually classified as a direct spinel, where Zn 2? cations occupy tetrahedral sites. The magnetic properties of cubic spinel ferrites are governed by the dominant antiferromagnetic exchange interaction between tetrahedral (A) and octahedral (B) cations, as described by the Néel model [6]. The resulting collinear ferrimagnetic order, that we label c0, is characterized by the parallel alignment of the magnetic moments of all B atoms and the opposite orientation of the A sublattice magnetization. Thus, in general, the net magnetization corresponds to the difference between the B sublattice contribution minus that from the A sublattice. This balance alters when nonmagnetic elements occupy the A sites, as occurs in ZFO: then, the competition of magnetic exchange interactions between B cations determines the global magnetic order, leading to a rich magnetic phase diagram involving frustration [7][8][9].
In the case of NFO, the largest local magnetic moments correspond to Fe 3? cations at A sites; thus, a route to enhance the net magnetization is to remove these cations introducing a partial substitution of Ni by Zn atoms that are expected to occupy A positions. This method has opened large expectations to improve the performance of nanostructured permanent magnets where NFO is the soft component. For moderate Zn concentrations, these Ni-Zn ferrites usually show higher values of the saturation magnetization as compared to pure NFO [10,11]. But the method fails as the Zn concentration increases, and the reasons for such failure are not clear. One possibility could be the occupation of substitutional B sites by Zn, as the specific cation distribution at nanostructures depends on the processing parameters and preparation conditions, particularly for mixtures of both oxides [1,[12][13][14][15][16][17][18][19][20]. Supporting this idea, contrary to the robust tendency of Ni to avoid tetrahedral positions [21], the occupation of octahedral sites by Zn has been reported in pure ZFO [9,19], particularly in thin films [22][23][24][25][26][27][28][29] and nanoparticles [30][31][32][33], and also in other Zn substituted cubic spinel ferrites [34][35][36]. Although many studies discard the presence of Zn atoms at octahedral positions in NFO, the precise identification is difficult in the often disordered nanostructures, and in fact, a mixed spinel with Zn atoms both at A and B sites has sometimes been reported [11,30,37,38] and even claimed to favor the mixing of NFO and ZFO [39]. Another possibility behind the failure to promote a high magnetization lies in the emergence of a low temperature phase rooted in the Yafet-Kittel canting of the local magnetic moments at B sites [18,40]. Both effects are detrimental to enhance the magnetization, but their consequences are different: Magnetic canting manifests only at low temperatures and evolves to the usual ferrimagnetic collinear order as the temperature increases [40], while partial Zn inversion is an intrinsic feature of the nanostructure.
From a theoretical point of view, the electronic structure of Zn doped NFO has been considered only scarcely. Early density functional theory (DFT) calculations described the full range of concentrations at Ni x Zn 1−x Fe 2 O 4 based on LDA or GGA exchange correlation functionals that fail to describe the cubic spinels as insulators [41][42][43]. Later simulations focused on the solubility of NFO/ZFO mixtures included a local U term to correct this description [39]. The value of U was taken from previous calculations of the corresponding monoxides, and in particular, a large value over 7 eV was used for Zn d electrons. This is contrast with previous DFT works devoted to ZFO, where a moderate or even null U (Zn d ) is claimed to properly capture the experimental features [8,44]. In this work, we revisit the substitution of Ni by Zn at NFO, addressing in detail the tendency of Zn to occupy either A or B sites, and the subsequent magnetic and electronic properties as described by DFT?U. As we will show, special attention should be deserved to the previously unexplored influence of U(Zn d ) on the final results, not only at ZNFO but also in the description of ZFO. We conclude that, while the relative position of the Zn d band plays a minor role in the electronic features of cubic ferrites containing Zn, it is crucial to determine the relative stability between the inverse, direct and mixed spinel forms, which has consequences on the magnetism of the final system.
The paper is organized as follows: after introducing the theoretical procedure, we present our results for the pure NFO and ZFO systems. Then, we consider two different representative concentrations of the Zn doped NFO that can be modeled with moderately large unit cells, namely Zn x Ni 1−x Fe 2 O 4 with x =0.25, 0.50. We explore different Zn distributions, always allowing Zn atoms to occupy A and B sites. For all structures, we determine in detail the structural, electronic and magnetic properties, also including the spin-orbit coupling to evaluate the magnetic anisotropy. We end with a discussion of our results leading to the final conclusion that an important number of octahedral Zn positions needs to be considered in Zn substituted NFO, limiting its ability to enhance the magnetization of the soft magnet.

Theoretical method
We have performed ab initio simulations within the DFT as implemented in the Vienna Ab initio Simulation Package (VASP) [45,46]. The Perdew-Burke-Ernzerhof (PBE) parametrization of the Generalized Gradient approximation (GGA) and the PBE revised for solids (PBEsol) has been considered as exchange correlation functionals, adding an effective Hubbard U parameter (GGA?U) at the d valence band (VB) orbitals following the Dudarev approach. The adequate value of U for Fe and Ni at the systems under consideration has already been thoroughly studied [47][48][49][50], and we adhere to the usual description based on intermediate U values around 3-5 eV. As we will show later, setting U(Ni d , Fe d )=4 eV is the best choice to properly describe the structural, electronic and magnetic properties of NFO under our theoretical conditions and also to reproduce the ZFO experimental band gap around 2 eV [51,52]. Regarding U (Znd), previous DFT simulations of ZFO [8,9,44,[53][54][55] have considered different values, from 0 to 5 eV, but without performing a detailed study of its influence on the obtained properties. We will address such study in the next sections.
The core electrons are described with the Projector-Augmented Wave method (PAW) using an energy cutoff of 500 eV. The local charges and magnetic moments are determined based on a post-DFT Bader analysis of the real space charge distribution [56]. The Brillouin zone is sampled with a (1091097) k-mesh for a unit cell of 28 atoms as shown in Fig. 1. All structures considered are allowed to fully relax, including the lattice vectors and the atomic positions, until the forces on the ions are below 0.01 eV/A. This way we determine both the equilibrium volume and the optimized positions. Zn x Ni 1−x Fe 2 O 4 (ZNFO) is modeled replacing one (x=0.25) or two (x=0.50) Ni atoms of this unit cell by Zn. Substitutional Zn is allowed to occupy both tetrahedral and octahedral cation sites, and different Zn distributions are considered, as detailed in the next sections.
The magnetic configuration of each system is explored comparing different reference magnetic orders that serve to extract exchange interactions mapping to a simple Heisenberg Hamiltonian [57]. The spin-orbit coupling is included only to evaluate the magnetic anisotropy energy (MAE) that is determined directly from ab initio total energy differences of situations where the net magnetization is oriented along different high symmetry directions.

NiFe 2 O 4
As mentioned in the introduction, NFO has been largely studied based on DFT, and the optimized conditions to reproduce the experimental features are well known. In Fig. 2a, we show the evolution of the spin-resolved total density of states (DOS) of our 28 atoms unit cell as a function of the effective U value applied to the Ni and Fe d orbitals that we always assign to be same. It is evident that the lack of inclusion of a Hubbard term fails to describe the oxide as an insulator, the experimental gap of 1.64 eV [58] being obtained for U=4-5 eV. Figure 2b) provides the detailed atomic projections of the DOS, evidencing the large overlap of the Ni and O states, and the localization of the Fe d orbitals at the bottom VB and top conduction band. Figure 3 shows the evolution of several structural and electronic features of NFO as a function of U, using the experimental cubic lattice parameter a= 8.34 Å [59] and under PBEsol. We explore different magnetic configurations detailed in section A of the Appendix, following the procedure explained in [57]. The total energies in the top panel evidence that, no matter the value of U, the ground state magnetic order always corresponds to c0, with moderate variations in the energy barriers between configurations. Neither the net magnetization (intermediate panel) has a significant dependence on U for each selected configuration, even though the local magnetic moments (not shown) are enhanced as U increases, reflecting the impact of U on the localization of the d states. Regarding the atomic structure, the bottom panel provides the internal pressure that goes to zero at the equilibrium volume. The graph evidences two features: first, the (nonlinear) dependence of the lattice parameter on U is different at each configuration; this alters the energy barrier between different magnetic orders, that is, the magnetic exchange constants, pointing to the crucial role of U (Ni d , Fe d ) on the description of the ferrite. Second, restricting to the equilibrium magnetic order c0, our model provides a slightly compressed structure (negative pressure) that better approaches the experimental solution for U values≥4 eV. Further relaxation of the structure at U=4 eV under c0 serves us to determine our reference theoretical cubic lattice parameter to be a=8.30 Å for NFO.
It is worth to remark that the same study based on the PBE exchange correlation functional instead of PBEsol provides a significantly expanded lattice (a= 8.43 Å ). This is expected, as PBEsol was conceived to correct the bias of GGA parametrizations like PBE toward the description of free-atom energies, thus reducing their tendency to overestimate lattice volumes. Nevertheless, the electronic structure and the robust stability of the c0 magnetic order coincide under both functional forms.

ZnFe 2 O 4
ZFO crystallizes in the direct spinel structure, with a cubic lattice parameter of 8.52 Å according to the most accurate experimental determination [7], though values around 8.45 Å are commonly reported [60][61][62]. The lack of magnetic atoms at A sites triggers the competition between different B-B magnetic exchange interactions, leading to a complex magnetic phase diagram where spin-glass behavior, magnetic frustration and the coexistence of long-range and short-range orders have been proposed [7,54]. Though recent works support the dominant antiferromagnetic nature of the exchange interactions [8,9], anomalous ferromagnetism has often been measured, its origin being mainly attributed to the partial inversion of ZFO [53,63].
Previous ab initio calculations of direct ZFO have determined the ground state magnetic configuration, identified as the antiferromagnetic AF2 order represented in Fig. 4 [54]. Alternative solutions, such as the AF1 also shown in the figure, are very close in energy, justifying the complex magnetic phase diagram and low magnetic ordering temperature. The degree of inversion is defined by the parameter λ that describes the number of B sites occupied by Zn as follows: (Zn 1−λ Fe λ ) [Zn λ Fe 2− λ ]O 4 , where round and square brackets denote A and B sites, respectively. All previous theoretical results claim the stability of the direct form over the inverse one [8,9,53]. Partial inversion favors the local emergence of the ferrimagnetic order typical of NFO, c0, that becomes the ground state under full inversion. The case λ=1 is also characterized by the sparse distribution of Zn B atoms, together with the slight contraction below 0.02 Å of the lattice parameter as compared to the ZFO direct spinel form [53]. This contraction is expected: at the direct form, the large Zn atoms dictate the volume of the entire A sublattice, while at the inverse form, the presence of Fe atoms at both A and B sites conditions the volume of both sublattices. Experiments confirm a volume reduction under inversion even higher than theoretical predictions, with a cubic lattice parameter a=8.38 Å for λ=0.5 [33]. Contrary to the large effect on the magnetic order, the band gap and local magnetic moments of ZFO are hardly affected by inversion.
Here, we consider the two extreme situations, λ=0 and λ=1, and explore the influence of U(Zn d ) on the structural, electronic and magnetic properties of ZFO, Figure 1 General structure of a cubic spinel ferrite, A refers to tetrahedrally coordinated sites and B refers to octahedrally coordinated sites.
initially based on PBEsol. We restrict to the cubic lattice, as we have checked that full lattice relaxation leads to a negligible distortion from cubic symmetry, but the atomic positions are fully relaxed. To take into account the magnetic configurations of the direct form, a large cubic unit cell of 56 atoms is used. Figure 5 shows the spin-resolved DOS of the direct and inverse spinel structures, each one at the corresponding equilibrium magnetic configuration, under the two limiting U(Zn d ) of 0 eV and 10 eV. We will focus on these extreme U values, after considering intermediate ones that confirm a gradual evolution of the properties of ZFO as a function of U(Zn d ). The figure evidences that the Zn d states locate far from the top VB even under U=0 eV, gradually shifting to the core region as U increases. Consequently, the Zn d band plays a minor role in the determination of the gap and the distribution of states at the top VB. The moderate variation of the DOS upon inversion, that we obtain for all U(Zn d ), is also in good agreement with previous calculations [9,53].
Next, we explore the energy barrier to ferromagnetism of the different magnetic orders, c0 for inverse ZFO and AF1, AF2 for the direct spinel, based on the experimental lattice parameter of 8.52 Å . The results are shown in Fig. 6. It is evident that, no matter the U value, the magnetic energy barriers keep essentially unaffected, the ground state magnetic configuration corresponding to the inverse c0 and direct AF2 orders. We have checked that the equilibrium volume is only moderately dependent on the magnetic order, so that the magnetic energy differences are not expected to change upon relaxation of the structures. This low influence of U(Zn d ) on the magnetic properties is also reflected on the local magnetic moments, that keep constant across all U values for each magnetic configuration. The Fe, Zn and O magnetic moments in μ B are, respectively, 4.21, 0.00 and 0.11 for the direct AF2 structure, and 4.15, 0.05 and 0.02 for the inverse c0 one. These values compare well with both previous calculations and the experimental evidence [53,61].
On the contrary, the figure reflects that the relative stability of the direct structure over the inverse one gradually reduces as U increases. This result might be further affected by the different equilibrium volumes of both forms, as we are using the experimental lattice of the direct structure. To explore this aspect, we determine the equilibrium lattice parameter of the inverse and direct structures; the results are shown in Fig. 7. From the figure, it is clear that the stability of the inverse structure is reinforced when both forms are compared at their respective equilibrium volumes, an effect enhanced the larger the value of U (Zn d ). Independently of the chosen U(Zn d ), the low energy difference between structures with Zn at A or B sites explains the experimental evidences of partial inversion at ZFO.
From these results, we can conclude that U(Zn d ) has a minor impact on the description of the electronic and magnetic properties of ZFO, and particularly on any measured observable that could be used as a reference to choose its most adequate value. However, it is crucial in determining the preference for Zn to occupy tetrahedral or octahedral coordination sites, an essential point at our subsequent study of ZNFO. Additional guidance can be obtained at former calculations of the electronic structure of wurtzite ZnO, that have evidenced the need to correct the standard DFT description of the Zn d band by shifting it to the core region [65]. In our model, this would be equivalent to the use of the largest U(Zn d ). Though Zn only adopts tetrahedral coordination at ZnO, the similar charge of Zn at the inverse and direct spinel forms, with differences under 0.05 e, discards any dependence of U(Zn d ) on the coordination environment. All this invites to set U=10 eV in our description of Zn ferrites. Nevertheless, as the presence of Zn at A or B positions is key to determine the system magnetization, and taking into account the lack of any experimental data confirming the best choice of U(Zn d ), we will keep it as a free parameter in our study of ZNFO.
As a final comment, Fig. 7 reveals that our approach largely underestimates the experimental ZFO lattice parameter, particularly under large U (Zn d ). The origin of this discrepancy can be traced back to our choice of the exchange correlation functional, as we will demonstrate now. We have performed additional calculations based on the PBE functional that, as evidenced for NFO, is expected to favor larger volumes. Table 1 shows the equilibrium lattice parameters of inverse and direct ZFO under both PBE and PBEsol, proving the better agreement of PBE with the experimental information. These results can be understood as PBEsol is particularly adequate to describe densely packed solids, and the large size of the Zn atoms reduces the compactness of the structure. But the important aspect is that the  electronic properties and the magnetic ground state are univocally described by both exchange correlation functionals, similarly to the case of NFO. Noticeably, also, both functionals provide a similar tendency of volume reduction as we go from λ=0 to λ=1, more pronounced under U(Zn d )=10 eV, though always lower than experimentally observed. Another relevant point for our study of ZNFO is that the trend to favor inversion as a function of U(Zn d ) is not altered by the choice of either PBE or PBEsol. This is further proved in Table 2, that compiles the excess energy of the inverse structure over the direct one for the different exchange correlation functionals and under different choices of the lattice parameter: either the experimental one or the theoretical value minimizing the energy, in this case also allowing full relaxation of the atomic positions. The positive energy in the table indicates that in all situations, the direct spinel is the ground state. However, the inverse form becomes closer in energy as U(Zn d ) increases, particularly when the structure is allowed to relax. All this confirms that the choice of either PBE or PBEsol does not alter our description of the relative energy trends when comparing the inverse and direct structures. In the following, and as we are interested in the description of Zn impurities in NFO, a material better described by PBEsol, we will keep this choice of the exchange correlation functional.

Zn substituted nickel ferrite (ZNFO)
We address now the properties of Zn x Ni 1−x Fe 2 O 4 with x=0.25, 0.50. In principle, Zn atoms enter as substitutional cations of Ni that occupies octahedral (B) sites. However, the previous study of ZFO has evidenced that Zn prefers to fill tetrahedral (A) positions at the spinel ferrite structure. Both situations will then be considered. In the following, we present the results for each Zn concentration separately, starting by x=0.25.

x=0.25
Within our unit cell of 28 atoms, this concentration is equivalent to the replacement of one Ni atom by Zn. Figure 8 shows the resulting spin-resolved DOS of this unit cell at the magnetic ground state of NFO, c0, with Zn either at A (ZnA) or B (ZnB) sites and for the two extreme values of U(Zn d ) considered here. Similarly to the case of ZFO, the choice of the largest U (Zn d ) shifts the Zn d band close to the core region. But even for U=0 eV, the position of the Zn d states lies at the bottom VB, and the features at the upper VB are not affected by the Zn d states. It is also remarkable the low influence of the A or B position of Zn on the global features of the DOS, as occurred at ZFO.
When non-magnetic Zn cations enter NFO replacing Ni, the magnetic order may be altered. We have considered the possible couplings between the different magnetic sublattices shown in Table A1, analyzing their relative stability as a function of U(Zn d ). Combining all these configurations with different U (Zn d ) values imply a considerable number of  simulations. To reduce the computational cost, first we have checked the smooth dependence of the system features on U(Zn d ), restricting to the unrelaxed structures. The results are shown at section B of the Appendix. From this evidence, we will limit the detailed analysis of the relaxed structures to only a few U(Zn d ) values. Relaxation does not modify the cubic symmetry for Zn A , but it introduces a slight deformation when Zn occupies B sites. In both cases, there is also an enlargement of the unit cell volume with respect to NFO, of 0.3% in the case of Zn A and 0.2% for Zn B , which is consistent with the experimental evidence [10]. This enlargement of the unit cell can be ascribed to the larger size of Zn ions as compared to Ni. In addition, and similarly to the situation at ZFO, also at Zn 0.25 Ni 0.75 Fe 2 O 4 setting, U (Zn d )=0 eV stabilizes a volume larger by 0.15 Å 3 / atom than using U(Zn d )=10 eV. Figure 9 shows the relative stability of the different magnetic orders after relaxing the structure at each magnetic configuration for different U(Zn d ). We can see that U(Zn d ) has a negligible effect on the magnetic energy barriers at a fixed Zn site, with c0 always as the ground state. Also, the presence of Zn does not seem to alter significantly the hierarchy between magnetic orders observed at pure NFO. However, an important conclusion from Fig. 9 is the low energy difference between Zn A and Zn B , particularly at c0. Moreover, and as already observed at ZFO, the choice of U(Zn d ) influences which position becomes more stable: Zn A (Zn B ) is favored at the lowest (highest) U(Zn d ).
This study of the energetics enables a detailed analysis of the magnetic exchange interactions. Following the procedure explained in Ref. [57], we map our ab initio total energy differences to a Heisenberg model of the form H ex =Σ ij J ij S i ·S j , where J ij are the exchange coupling constants between magnetic sublattices i and j, and Si represents the magnetic moment of sublattice i. The procedure requires the simulation of different magnetic couplings, sometimes far away from equilibrium. As we will see, if such configurations do not lie in a local energy minimum, it may become impossible to converge the selfconsistent procedure and thus, to extract the corresponding J ij . We follow the convention that positive (negative) J ij denotes antiferro-(ferro-)magnetic interactions. As our system is formed by three different magnetic sublattices, Ni, Fe A and Fe B , we end up with six different exchange coupling constants, Figure 7 Evolution of the total energy as a function of the cubic lattice parameter for the inverse (blue) and direct (red) forms of ZFO at their respective magnetic ground states for U(Zn d )=0 eV (left panel) and 10 eV (right panel). At fixed U(Zn d ), the energy zero is set at the structure with minimum energy. In all cases, the atomic positions have been fully relaxed. The vertical dotted red (blue) lines indicate the equilibrium lattice parameter of the direct (inverse) structure.  Table 3. We will refer to this set of six J ij as Js from now on. The influence on the Js of the Zn distribution for different U(Zn d ) is represented in Fig. 10. We could not converge values for J NA and J NN at Zn A under U (Zn d )=10 eV. As could be inferred from the robust stability of c0 shown in Fig. 9, the dominance of the antiferromagnetic exchange interactions between cations at octahedral and tetrahedral sites (J AB , J NA ) is maintained, no matter the Zn position or the value of U(Zn d ). The presence of Zn induces some changes of sign of the weakest interactions that further depend on U(Zn d ), but their low values make negligible their influence on the magnetic order. Under U(Zn d )= 10 eV, apart from the slight reduction in the dominant exchange constants, J AB and J NA , the most relevant effect is the noticeable increase in the presence of Zn B of J BB , that starts to compete with J AB , J NA . The mutual antiferromagnetic alignment of Fe B cations due to the reduction in the exchange coupling between A and B sublattices has already been suggested at ZFO [53].
Restricting now to the ground state magnetic order c0, we have found no significant influence of U(Zn d ) neither on the Bader charges nor in the local magnetic moments. This becomes evident at Table 4, that compiles the average magnetic moment of each atomic species for different U(Zn d ), both for Zn A and Zn B and for different Zn concentrations. As a reference, ZFO and NFO are also shown. Some values are marked with a plus/minus sign, indicating that each half of the corresponding atoms of the sublattice takes a different orientation (leading to a null mean value). Focusing on case x=0.25, only the Zn position, either Zn A or Zn B , introduces some minor differences in the local magnetic moments. Variations are negligible for the magnetic cations and significant at the hybridized moments acquired by the non-magnetic atoms, but they always remain in the 10 -1 µ B . On the contrary, the net magnetization of the entire system is significantly modified depending on the Zn site, because the balance between magnetic A and B sublattices is altered when a non-magnetic Zn atom replaces the contribution of a magnetic sublattice site. As a result, Zn A doubles the magnetization of NFO, while Zn B reduces it by 25%.
Finally, we have determined the magnetic anisotropy energy computing the total energy of configurations with the magnetization oriented along the and Zn B , respectively, which is close to our accuracy limit in the determination of the MAE. In all cases, the easy axis at [110] is maintained. This is in good agreement with experimental results, where the coercivity of the material reduces with the increase in Zn content [10,30].

x=0.50
Substituting two Zn cations at our unit cell of 28 atoms is equivalent to a concentration x=0.50. Now, there are more possible atomic arrangements, with both Zn cations at A sites (Zn A ), both at B sites (Zn B ) or each one at a different coordination site (Zn AB ). In addition, each situation admits different Zn-Zn interatomic distances and different Ni distributions. Addressing a detailed study of diverse magnetic configurations under all these possibilities is a computational challenge out of the scope of this study. Instead, in the Supplementary section C, we analyze the most stable conditions for configurations Zn A , Zn B and Zn AB . The main conclusion is that several Zn distributions may coexist, though a set of three representative atomic arrangements can be identified, each corresponding to one of the three configurations (Zn A , Zn B , Zn AB ). These structures provide the minimum energy of their corresponding configuration independently of U(Zn d ). They also provide the representative magnetization, as the local moments are not affected by the interatomic Zn-Zn distance neither by the particular Ni distribution. More delicate is the limitation to a reduced set of structures in order to determine the local exchange interactions, we will come to this point later. In the following, we restrict our analysis to the three representative distributions selected in Supplementary section C.
According to this choice, Fig. 11 shows the relative stability of Zn A , Zn B and Zn AB , after full structural relaxation and at the magnetic ground state, c0. The figure is just a simplification of Figure C4, and thus, it again evidences that different Zn distributions may coexist, and that U(Zn d ) has a large influence on the identification of the most stable one. Similarly to x= 0.25, Zn B is favored over Zn A under large U(Zn d ). Further, independently of U(Zn d ), the simultaneous presence of Zn atoms at both A and B sites cannot be discarded, even providing the largest stability under U(Zn d )=0 eV.
The structural features after relaxation keep trends similar to those observed at x=0.25 and ZFO. Always, U(Zn d )=0 eV leads to higher volumes than U(Zn d )= 10 eV (about 0.3 Å 3 /atom above), while Zn A tends to slightly expand the lattice over Zn B (by 0.05 Å 3 / atom). The volume of the unit cell is enlarged by 0.4% with respect to pure NFO, again consistent with the experimental evidence [10], and as expected from the increase in the number of Zn atoms in the structure.
Regarding the energy barrier between different magnetic configurations, here, the trends are similar to those observed at lower Zn concentrations: c0 is a robust magnetic ground state, and the presence of Zn cations at B sites reinforces its stability. Figure 12 shows the magnetic exchange constants, obtained following the procedure explained in the previous subsection. It is evident the dominance of the antiferromagnetic exchange between cations at A and B sites, supporting the stability of c0. Further, there is a strong dependence of the Js on the Zn distribution, though in general, the hierarchy between the relative J ij strengths is preserved. The only exception is the enhancement and sign reversal of J BB when Zn occupies only B sites. This result is affected by the specific structure used to extract the Js (2ZnB-1, see Appendix), where Fe B atoms are not nearest neighbors. This result puts a word of caution before extracting conclusions about the detailed balance of exchange interactions without a more thorough consideration of different atomic distributions. However, the general trends evidenced in Fig. 12 suggest that, still at x=0.50, the antiferromagnetic coupling between A and B sublattices is the dominant magnetic interaction that guarantees the stability of c0. Last, it is also noticeable the tendency to increase the Js under U(Zn d )=10 eV, not observed at x=0.25, and that here manifests as more relevant when there are Zn atoms at B positions.
The magnetization and local magnetic moments at c0 are in the lower part of Table 4. In general, the net magnetization is again enlarged under Zn A , as more Fe atoms are forced to move to octahedral positions, with a record enhancement of 4.0 μB/f.u. over pure NFO. However, the individual magnetic moments do not vary noticeably depending on the Zn position. There is only a slight reduction in the moments of the B cations under Zn B , and more dispersed values for the non-magnetic atoms depending on the Zn distribution.
Finally, our results for the MAE indicate a further reduction as the Zn content increases, with values between 0.010-0.014 meV/Fe depending on the Zn positions. This is again consistent with the experimental evidence of a reduction in the coercive field as the Zn content increases.  Table 3 Magnetic exchange interactions between magnetic sublattices at ZNFO. As J XY = J YX , for clarity, we only show the upper part of the Table  Discussion The main purpose of the present study is to understand from ab initio simulations the failure of Zn substitution to enhance the magnetization of NFO. In the previous sections, we have proved the reliability of our model, that reproduces the measured experimental features induced by Zn at NFO: a moderate volume expansion that increases with the Zn concentration, the stability of the ferrimagnetic order and the reduction in coercivity. Further, our description supports the feasibility to tune the magnetization of NFO by substitutional Zn without altering its global electronic and magnetic properties, an important aspect for applications. However, the results also evidence the main drawback of this approach, that arises from the strong tendency of Zn to remain at the B site of the Ni atom it substitutes. We show that this occurs even at low Zn concentrations, even though the effect enhances with the Zn content.
The evolution with the Zn content can be understood in the framework of the itinerant electron model for magnetic oxides (IEO model) that describes the magnetic exchange on the basis of the hopping process of itinerant O2p electrons taking into account the incomplete ionic character of the cation-anion bonds in actual oxides [66]. The properties of cubic spinel ferrites can be accurately described within this model [66]. In the case of NFO, it has been proved [67] that the B sites are mainly occupied by Fe 2? and Ni 2? cations, while Fe 3? goes to A sites, justifying the robustness of the c0 magnetic order as obtained here. When Zn atoms are incorporated replacing Ni, there is a depletion of Ni 2? at B sites that reduces the contribution of the majority spin population. Only if the Zn 2? cations enter A sites displacing Fe 3? to the B lattice, this reduction will be compensated by an enhancement of the majority spin owing to the larger magnetic moment of Fe as compared to Ni. We recall that Fe 2? , Fe 3? and Ni 2? have all 3d electron numbers n d ≥5 and thus, within the IEO model they are expected to align their moments mutually parallel when placed at the same magnetic sublattice. But beyond these features, the model further provides a clue to understand why Zn atoms tend to occupy B sites as their concentration increases. Similarly to the situation observed for Zn x Mn 1-x Fe 2 O 4 at high x [68], the replacement of Fe 3? cations by Zn 2? increases the Pauli repulsion energy at the A sublattice, and there will be a point where the A sublattice cannot accept more Zn 2? cations. This establishes a clear limit to the maximum Zn concentration that enables enhancement of the magnetization at NFO by Zn substitution.
All these results justify the diversity of experimental evidences claiming the presence of Zn atoms both at octahedral and tetrahedral sites [10,[30][31][32][33]64]. Identification of the Zn coordination environment at nanostructures is a complex task that requires advanced techniques often based on synchrotron radiation light sources. Our simulations invite to carefully revise the amount of Zn B cations in previously reported ZNFO samples with low magnetization.
A remaining important aspect that naturally emerges is how to revert the tendency of Zn to occupy octahedral positions, thus enabling large magnetizations at ZNFO. Though the final response requires additional research, some clues can be extracted from the results presented here. On one hand, Zn A and Zn B induce different volume expansions, larger by 0.05Å 3 /atom in the case of Zn A , which implies increasing differences of lattice volume as the Zn content augments. Consequently, it is possible that the growth of slightly expanded structures could favor the dominance of Zn A . Figure 7 supports this idea, as the stability of the inverse ZFO spinel form deteriorates under expansion. Other possibilites emerge from the alteration of the local exchange interactions brought by Zn, particularly the enhancement of J BB in the presence of Zn B . The subsequent new balance of interactions does not revert the stability of c0 over the rest of magnetic orders considered (c1, c2, c3), but the barrier to ferromagnetism is higher than for Zn A , as evidenced in Fig. 9. This suggests that the growth in the presence of magnetic fields might be used to promote occupation of Zn A sites, as the excess energy of Zn B over Zn A would be enhanced. Further ideas could be devised by introducing at ZNFO additional substitutional cations with more tendency than Zn toward octahedral coordination. However, care needs to be taken not to alter the desired NFO electronic and magnetic features as we increase the number of chemical components in the system.

Conclusions
We have addressed the DFT?U simulation of Zn substituted NFO, Zn x Ni 1−x Fe 2 O 4 , for x=0.25, 0.50, providing a detailed description of the structural, electronic and magnetic properties. Special attention has been paid to the position of the Zn d states at the VB, tuned by U(Zn d ). As also observed at ZFO, this U value has a negligible effect on the local magnetic moments and the electronic structure, but a large one on the relative stability of the A or B site for Zn substitution. This is an important result, as the magnetic interactions and net magnetization are largely dependent on the Zn distribution, the most the higher the Zn concentration. Further, no matter U (Zn d ), for all Zn contents there is evidence of the low energy barrier between both Zn coordination sites, supporting the coexistence of Zn A and Zn B in the cubic ferrite structure. Consequently, our results indicate that the partial substitution of Ni by Zn atoms as a route to enhance the net magnetization is severely limited by the simultaneous presence of Zn A and Zn B . As more detailed conclusions, we also observe that while the ferrimagnetic ground state c0 is preserved up to x=0.50 for all atomic arrangements, differences emerge regarding the competition between local exchange interactions. The dominance of the antiferromagnetic exchange between A and B magnetic sublattices seems robust, but the particular atomic distribution may alter significantly the strength and sign of other exchange interactions. Finally, the magnetocrystalline anisotropy is weakened with the increase in the Zn content, ensuring ZNFO as a robust soft magnet, in good agreement with the experimental evidence.
Our study also provides an important remark regarding the value of U(Zn d ) in cubic ferrites. Though U(Zn d ) has no relevant influence on the determination of the band gap of ZFO, neither it modifies the magnetism of the system (including both local magnetic moments and long-range magnetic order), it has a large effect on the relative stability of the inverse form, that is significantly enhanced. In addition, the larger U (Zn d ), the larger the contraction of the inverse lattice over the direct one. Though a final experimental proof still lacks, high U(Zn d ) values bringing the Zn d band close to the core levels seem to describe more accurately the measured lattice volume and to provide a better framework for the experimentally observed coexistence of Zn octahedral and tetrahedral coordination sites. 130957B-C53. Part of the simulations has been performed in the supercomputer Finis Terrae of the CESGA.

Author contributions
CT-C carried most calculations, contributed to the analysis of data and conclusions, prepared the figures and wrote the manuscript. RR did most calculations concerning ZFO. SG designed and supervised the work, did some calculations and contributed to the analysis of data, extraction of conclusions and manuscript preparation.

Funding
Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licen ses/by/4.0/.

Data availability
All data can be made available to interest readers upon request. The used code is the licensed package VASP.

Declarations
Competing interests The authors declare no conflict of interest for the research presented here.  Table A1, considering Zn both at tetrahedral and octahedral sites. The results prove the smooth shift of the Zn d band with U(Zn d ), and the independence of the global electronic features from the position of the Zn d states, as described in Fig. 9 for the relaxed structures. Figure 14 shows the evolution with U(Zn d ) of the total energy and internal pressure of these systems, the upper panel being equivalent to Fig. 9, but referred to the unrelaxed structures. Again, the figure reflects a smooth dependence on U(Zn d ). It is also evident the robust stability of the ground state magnetic configuration of NFO, c0, no matter the Zn position or the U(Zn d ) value. The most notable difference with respect to the relaxed structures is the relative stability of Zn A and Zn B .

Appendix A. Magnetic configurations at NFO and ZNFO
The internal pressure at the lower panel of Fig. 14 a measure of how comfortable is the system at the lattice parameter used in the calculation (that of NFO in our case). It goes to zero at the equilibrium volume. Figure 14 reflects the same results shown in the main text for the relaxed structures: The lattice contracts as U(Zn d ) increases, something already observed at ZFO. In addition, there is a different dependence with U(Zn d ) of the equilibrium volume for Zn A and Zn B , which correlates with the relative stabilities observed in the upper panel. This justifies why the relative stability of Zn A and Zn B varies with U(Zn d ) when the NFO lattice is kept fixed, while they become almost degenerate when the structure is allowed to relax. ? -- Figure 13 Evolution with U(Zn d ) of the spin-resolved DOS of unrelaxed ZNFO under x=0. 25. Majority (minority) spin contributions are shown as positive (negative) values. a) and b) correspond to Zn at octahedral and tetrahedral sites, respectively, while each column refers to a different magnetic order. The projection of the Zn states is represented by the red line..
Appendix C. Zn distribution at ZNFO under x=0.50 Figure 15 shows the different atomic arrangements considered for ZNFO at x=0.50, named according to the distribution of Zn atoms at A and B sites. The corresponding relative distance between Zn-Zn closest neighbors is also provided. Modifying the distribution of Zn at B sites necessary alters that of Ni, and different combinations have been taken into account. Our calculations do not pretend to cover exhaustively all possibilities, but to identify the most stable structures, and the existence or not of a competition between configurations Zn A , Zn B and Zn AB that may affect the magnetic properties of the system. The relative stability of each of these relaxed structures is shown in Fig. 16 as a function of U(Zn d ).
The most relevant feature is the strong dependence of the atomic distribution on U(Zn d ), that affects mainly the hierarchy between configurations Zn A and Zn B . But in general, the energy differences are low, evidencing that all these structures may be found at real samples. The results also indicate a tendency to favor sparse Zn distributions, even though the interatomic distance between Zn neighbors is not determinant. Other factors, as the distribution of Ni atoms, also influence the final energy of the system.
It is known that the arrangement of Ni atoms does not have relevant consequences on the magnetism of NFO. Here, we have confirmed that, for each selected configuration (Zn A , Zn B or Zn AB ), nor the local moments neither the net magnetization of the different structures in Figure C3 show any variation. In addition, both values of U(Zn d ) identify the same structure as the most stable for each configuration, allowing us to restrict our detailed study of the magnetic properties to three representative cases, namely 2ZnA-3, 2ZnB-1 and ZnAB-1.  Table A1. Structures with Zn at A (B) positions are represented by an empty circles (filled squares). The energy zero at each U(Zn d ) has been set at the most stable structure..   Figure C3. The energy zero at each U(Zn d ) has been set at the most stable structure. The line is a guide to the eye.