*Article* **Magnetic Transition State Searching: Beyond the Static Ion Approximation**

**Robert A. Lawrence , Scott J. Donaldson and Matt I. J. Probert \***

School of Physics, Engineering and Technology, University of York, Heslington, North Yorkshire YO10 5DD, UK

**\*** Correspondence: matt.probert@york.ac.uk

**Abstract:** The effect of structural relaxations on the magnetocrystalline anisotropy energy (MAE) was investigated by using density functional theory (DFT). The theory of the impact of magnetostructural coupling on the MAE was discussed, including the effects on attempt frequency. The MAE for ferromagnetic FePt (3.45 meV/formula unit) and antiferromagnetic PtMn (0.41 meV/formula unit) were calculated within the local density approximation (LDA). The effects of the structural relaxation were calculated and found to give a <0.5% reduction to the MAE for the ferromagnet and ∼20% for the antiferromagnet.

**Keywords:** magnetocrystalline anisotropy; FePt; PtMn; transition state search; magnetism; density functional theory

#### **1. Introduction**

The modern world has been transformed since the advent of the computer due to ever-improving data storage, enabling more complex devices and programmes to be manufactured and written. Huge advances have been made possible due to the development of and improvements in magnetic data storage from the primitive core-rope memory (∼4 Kb) of the Apollo era [1], to modern multi-Tb storage drives [2]. This progress has depended upon the construction of an ever-greater areal bit density (number of data storage "sites" per unit area) [2] and higher performance bits.

For high performance (high clock speed) memory, such as is found in modern random access memory (RAM), the storage must typically be electric rather than magnetic, because conventional magnetic storage techniques take too long to write a bit. This causes the memory to be volatile [3], and has a sizable effect on the power consumption of that memory, because constant refreshing of the bits is required. One potential way around this is antiferromagnetic magnetoelectric RAM (AFM-MERAM) [4], which makes use of the high-switching speeds of antiferromagnetic (AFM) materials to create memory that offers equivalent performance to standard charge-based RAM while still having the low energy requirements and nonvolatility of magnetic storage.

Now, considering the areal bit density, the greatest current limit is the requirement that the memory itself is thermally stable [5]; it is not much use if a warm day causes the data to be erased due to random bits flipping. Therefore, it is necessary to have a high barrier against bit-flipping, such that the probability of a random flip caused by thermal fluctuations is very small, even over long timescales (there must be a high activation energy for the bit flipping at *kBT* for typical operational temperatures).

In modern memory-storage devices that rely upon the magnetisation orientation of a magnetic material, this barrier to bit flipping is controlled by the energy required to reorient all of the spin moments within the bit. This barrier is therefore extrinsic, and increases with the size of the device. Any increase of the areal bit density, however, requires shrinking the size of a bit, and therefore to maintain a constant barrier to flipping needs a higher intrinsic barrier to flipping.

**Citation:** Lawrence, R.A.; Donaldson, S.J.; Probert, M.I.J. Magnetic Transition State Searching: Beyond the Static Ion Approximation. *Magnetochemistry* **2023**, *9*, 42. https://doi.org/10.3390/ magnetochemistry9020042

Academic Editors: C˘at˘alin-Daniel Constantinescu and Lucian Petrescu

Received: 30 November 2022 Revised: 16 January 2023 Accepted: 20 January 2023 Published: 27 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This intrinsic barrier to the flipping of individual bits is the magnetocrystalline anisotropy. This property, typically evaluated as a magnetocrystalline anisotropy energy (MAE), is a relativistic effect induced by spin-orbit coupling linking together orbital and spin angular momenta. One way of thinking about this effect is that, in the absence of relativity, spins exist in a spherical orbital. This may then be convolved with an electronic orbital to produce the standard "spin-orbitals" of, for example, Hartree–Fock theory, without affecting the spins. In the relativistic limit, however, *l* and *s* are no longer acceptable quantum numbers, and therefore we cannot neglect the effect of the convolution on the spin-like degrees of freedom. Because the original electronic orbitals are both directionally dependent and oriented with respect to the crystal lattice (due to the crystal field effect), this means that the effect of spin-orbit coupling is to make the spins experience an environment which is directionally dependent — magnetocrystalline anisotropy.

Much effort has gone into evaluating MAE over the years, both in ferromagnetic (FM) and antiferromagnetic (AFM) materials through both experimental [6,7] and theoretical [8–10] approaches based on *ab initio* calculation. For antiferromagnetic materials in particular, this is made especially difficult by their lack of apparent response to an applied field, and therefore the MAE of AFM materials tends to be inferred rather than directly measured [7].

Reliable simulation techniques have huge potential within the field of materials discovery because they enable a cheap initial search of likely candidate materials and novel structures [11,12], which may then be synthesised to confirm the findings of the simulations. For AFM materials in particular, in which the MAE is not directly accessible, this ability to directly access a magnetic property through a simulation is a potentially important avenue for materials discovery. Current simulation methods, however, typically significantly overestimate the anisotropy value compared with experiment [8,9].

One reason for this discrepancy between theory and experiment may be that the theoretical simulations do not calculate precisely the same observable that is measured experimentally, such as different timescales of experiment and theory. This sort of error cannot be corrected by using more sophisticated treatments of the electronic structure — such as additional treatments of strong correlations, or better exchange-correlation functionals — because these better methods for electronic structure would still be applied toward more accurately calculating the wrong number. Additionally, by finding a more experimentally relevant description it is possible to gain greater insight into the underlying physical processes, thereby both enhancing understanding and potentially opening up new routes for device optimisation and material parameters tuning.

A potential cause of the overestimation of MAE may be that these techniques typically rely on optimising the ground-state structure and then either sampling different constrained spin orientations [13] or calculating the spin-torque and using that to evaluate the anisotropy [14]. These methods, however, do not typically allow for a relaxation of the crystal structure in response to the changing spin orientation. In this paper, we propose a new strategy to evaluate the magnetic anisotropy which also includes the magnetostructural coupling and apply these to ferromagnetic FePt and antiferromagnetic PtMn.

#### **2. Theory**

#### *2.1. Magnetic Transition State Searching*

Rather than the narrow definition of MAE as the difference in energy between orientation of the spins along the magnetic easy and hard axes, it is possible to extend the definition to a higher dimensional case (i.e., also considering structural relaxations) by instead defining it as the maximum energy of the minimum energy pathway — the energy of the transition state — between the magnetic ground state and its symmetrically identical reversed state (for example turning all "up" spins into "down" spins and vice versa). When we only consider magnetic relaxations, this definition yields the same MAE as before. However, when simultaneous structural and magnetic relaxations are considered (as must be the case in a real experiment) the extra degrees of freedom make further, potentially lower, energy routes valid.

Transition state searching within chemistry, which considers only the movement of atoms, is a well-established field with many techniques for evaluating the energy barrier for transition between two local minima within a wider energy landscape (in chemical contexts, these are the reactant and product of a reaction step). Much work has been done on developing procedures for these, such as the nudged elastic band method [15,16]. However, an older, and more computationally efficient, method — known as linear synchronous transit/quadratic synchronous transit [17] (LST/QST) — may be sufficient for accurate magnetic transition state searching, and so it is used here.

In LST/QST, a linear pathway between the two minima is constructed by interpolation, and the maximum along this pathway is found through a bisection search as a first approximation to the transition state (this is the LST component). Because this state has no obligation to be a saddle point — which all transition states must be [18] — a further optimisation in directions orthogonal to the original search direction is performed until a minimum is found (the QST component). By identifying a maximum in the direction between minima which is a minimum in all other directions, the saddle point associated with the transition between local minima will have been identified.

For magnetic systems, we may consider the spin orientation of the system to be an orthogonal search direction to the structural configuration of the system (see Figure 1), because any spin configuration and any structural configuration may coexist (even if with a significant energy penalty). This suggests that performing an LST search in spin space, followed by a structural optimisation under the constraint of fixed spin orientations (in the manner of the QST component of the search) would lead to an improved estimation of the transition state. This QST style search also has benefits over standard structural-only LST/QST in that the QST relaxation is guaranteed as the constrained spin orientation cannot be altered by the change of atomic positions or lattice vectors. It is this QST-like optimisation of the crystal structure that we propose as a necessary minimum correction to evaluate the true transition state between two magnetic configurations (a more sophisticated evaluation of the transition state would require a more sophisticated technique).

**Figure 1.** A cartoon schematic of the energy landscape for inverting the magnetic ordering in a domain. Note that the lowest energy route passes through a saddle-point which has a different geometry to the ground state of the magnetic easy axis. Note that all magnetic degrees of freedom are compressed into a single axis, and all geometric degrees of freedom are compressed into another.

Furthermore, considering only the spin-like directions, we know that our initial and final states are identical by symmetry (a naïve rotation of the crystal by 180◦ maps one to the other — and this does not affect the internal energy of the system). Through Hammond's postulate [19], we find the result that for degenerate minima the transition state lies halfway along the reaction coordinate (which for our case ranges from *θ* = 0 to *θ* = 180◦ ) — thereby explaining why the easy and hard axes are perpendicular for a simple collinear magnet. This relationship is unaffected by structural relaxations due to the orthogonality of the spin and structural dimensions.

Finally, we note that the highest quality transition state search estimations of the barrier height also include a correction for the zero-point energy of the phonon mode [20,21] associated with the transition (see Figure 2). These zero-point energy corrections are typically on the order of meV, and may be safely neglected for chemical barriers (which are on the order of eV) but must be considered for these magnetic barriers, which are on the order of meV themselves.

**Figure 2.** The effect of zero-point energy and temperature on the MAE of an arbitrary barrier. The energy of the electronic state of the system is increased such that the bottom of the potential well is no longer accessible, leading to a decrease in the effective barrier, which change the MAE from being represented by arrow a to arrow b. Increasing the temperature will lead to arrow c.

#### *2.2. Quasiparticle Picture of Magnetic Anisotropy*

The standard picture of magnetic anisotropy — that of an external field "pushing" spins from one configuration to another — is somewhat troubling when one considers the linear magnetic field alignment in experimental measures of anisotropy [7], as this cannot produce an angular rotation of spins. This is because there is never a component of the applied field perpendicular to the spins. It may then be preferable to consider a quasiparticle picture of the anisotropy, in an analogous way to considering soft phonons to model phase transitions [22]. This does not change the underlying physics but simply uses different language to describe it, which makes the physical processes more apparent.

In this picture, any arbitrary magnetic structure can be written as a sum of the groundstate magnetic structure and a collection of symmetry-appropriate magnons. For the rotation of the ground state to its symmetric partner (exchanging up and down labels), this corresponds to a long-wavelength (*q* → 0) acoustic magnon.

This gives rise to the current strategy for evaluating MAE; calculate (or measure) the creation energy of this particular excitation to the spin structure with a static structure. However, if the minimum energy (transition) pathway (MEP) follows a route that requires a structural relaxation, then a quasiparticle responsible for this (a phonon) must be generated simultaneously to the magnon. This scenario is shown in Figure 3.

**Figure 3.** A schematic Feynmann diagram for the calculation of MAE. An incoming photon enables the creation of a magnon (and a phonon if magnetostuctural coupling provides the MEP), thereby driving the magnetic transition.

In terms of quasiparticles, the experimentally measured MAE is the energy loss due to the generation of either the magnon directly or the magnon–phonon pair (for the bound pair where they act as a single quasiparticle rather than as two coexisting quasiparticles, we use the term "magnetophonon"). Simulating the MAE then requires one to evaluate the creation energy of all of these pathways. The authors note that further expansions in terms of additional interactions are possible, such as multiple phonon and multiple magnon excitations; however the contribution of these higher order terms to the expansion are smaller, and therefore are neglected here.

One final subtlety of note is that magnons and phonons are Majorana particles — that is, they are their own antiparticles [23] — and accordingly, one can read the Feynman diagram in Figure 3 in multiple ways. The most obvious are the two solutions that come from considering the phonon as an emitted (see Equation (1)) or absorbed particle (Equation (3)). Additionally, one may consider the case where no phonon interactions are involved (Equation (2)). We have

$$E^{+} = E\_{magnon} + E\_{phonon} = E\_{magnetophonen} \tag{1}$$

$$E^0 = E\_{\text{magnon}} \tag{2}$$

$$E^{-} = E\_{magnon} - E\_{phonon}.\tag{3}$$

This makes the experimentally measured MAE an expectation value of the multiple pathways with appropriate weighting with the final experimentally measured barrier given as

$$
\langle MAE \rangle = w^+ E^+ + w^- E^- + w^0 E^0 \,\text{(4)}
$$

where *w* ± are the Arrhenius factors for each route and may be expressed as

$$w^i = A\_i e^{\frac{E^i}{K\_b T}},\tag{5}$$

where *A<sup>i</sup>* is the Arrhenius prefactor for route *i*. This gives an explicit temperature dependence for the MAE and means that when *kBT* < *E i* , the lowest energy route will dominate the kinetics. In order to evaluate the MAE at any arbitrary temperature (that is, below the critical temperature for the relevant magnetic phase), all that is required is to evaluate the energy of the individual routes and the Arrhenius prefactor. This prefactor may be thought of as the rate at which a system may attempt to overcome the transition barrier.

We also note that although a continuum of pathways for the spin flip exists, because the energy landscape is continuous, we only consider the routes that are lowest in energy for any given conditions. An experiment will average over all possible pathways, with those that are lower in energy dominating due to the exponential factor within Equation (5) (the Arrhenius equation).

#### *2.3. Attempt Frequency*

It is well known that a further important property when interpreting experimental data relating to the anisotropy is the attempt frequency [24,25]. This is typically stated to be a constant characteristic frequency for spin-flipping events (and therefore corresponds to the Arrhenius prefactor). It is notable that the attempt frequency for antiferromagnetic materials (10<sup>12</sup> Hz) [25] is three orders of magnitude higher than for ferromagnetic materials (10<sup>9</sup> Hz) [26] and also that these are of the same order of magnitude as the reciprocal phonon–phonon lifetimes *τph*−*ph* for acoustic and optical phonons, respectively [27].

A possible explanation for this is outlined in schematic form in Figure 4. In this scenario, the specific magnon mode commensurate with a spin-transition and its commensurate phonon (which must have the appropriate symmetry to couple strongly with this magnon mode) are abstracted from two general baths of magnons and phonons respectively, and specific transition rates between all four (and internal conversions for the bath) are considered. In the strong magnon–phonon coupling limit (*ke*,*ph* → ∞), these two modes form a single magnetophonon and may be considered as a single quasiparticle, which may be scattered by either magnons or phonons. By considering rate equations and requiring the system to be in a dynamic equilibrium, we may determine rates of conversion between and the populations of all levels.

**Figure 4.** Schematic of the possible interactions between the magnon associated with the selected spin transition, its commensurate phonon, and the baths of available phonon and magnon states.

First, we define an attempt as the creation of a magnon-like quasiparticle (either a magnon or magnetophonon) that is commensurate with changing the state of the electronic system to the barrier in spin-space between up and down configurations. Because these quasiparticles are bosonic, there is no limit to the amount of population of a single mode, so accordingly there are no density-of-states considerations in the transition rates, only requirements for conservation of energy and momentum which then yield selection rules.

By the principle of detailed balance, for a system in equilibrium (or near to equilibrium in the adiabatic limit) the rates of creation and annihilation of the magnon/magnetophonon mode must be equal (otherwise the population of the mode will change and the system is not in equilibrium). This suggests that the attempt frequency is the lifetime of the magnon/magnetophonon, which for the case of the magnetophonon may be calculated from the magnon and phonon scattering rates by using Matthiessen's rule,

$$\frac{1}{\tau\_{\text{magnetophono}}} = \frac{1}{\tau\_{\text{magnon}}} + \frac{1}{\tau\_{\text{phonon}}},\tag{6}$$

where *τquasiparticle* is the lifetime of the respective quasiparticle. The lifetime of the magnon seems unlikely to change significantly for different spin structures, as the magnon will remain acoustic (an optical mode corresponds to an FM-AFM transition). The magnon and phonon lifetimes will be dominated by scattering into their respective baths, and therefore the lifetime of the magnetophonon will be lower than the shorter of the phonon and magnon lifetimes. This magnetophonon picture (rather than the bare-magnon) explains the finding that AFM materials have a higher attempt frequency; they couple to optical phonons with a shorter lifetime, thereby increasing the attempt frequency, which may be represented as *τ* −1 *magnetophonon*.

#### 2.3.1. Symmetry Considerations

There are symmetry considerations with regard to the selection of which magnon and phonon modes couple together to form the magnetophonon associated with the MAE. For the magnon modes, the magnon drives the reorientation of the spins with respect to the crystal axes, and the lowest energy magnon mode is an acoustic mode in the long wavelength limit (*q* = 0). This mode rotates all spins across all space uniformly.

Any phonon mode coupled to this must affect all of space in the same way as the magnon (*q* = 0), but whether it is an acoustic or an optical mode depends on whether the system is ferromagnetic (in which case it will be acoustic) or antiferromagnetic (in which case the opposite responses of the different sublattices will lead to an optical phonon).

One result of this is that applying a magnetic field to a ferromagnetic system will also generate an acoustic phonon, leading to a volume change (magnetostriction). Conversely, applying a magnetic field to an antiferromagnetic system will generate an optical phonon (which conserves the volume) and so there will not be any magnetostriction. This behaviour is in agreement with previous studies on other magnetic systems [28].

It is also worth noting that the energy of a *q* → 0 acoustic phonon is vanishingly small, which suggests that the overall effect on the MAE of the magnon–phonon coupling will also be very small. However, this is not true for an optical phonon, and therefore the effect should be significantly larger in AFM materials than in FM materials.

#### 2.3.2. Conservation of Energy

The energy for the creation of the magnon (or magnetophonon) that drives the transition cannot come from the applied magnetic field. The total energy of a 1 T magnetic field, as typically used in experiments, in a volume of free space equal to the size of a crystallographic unit cell is ∼0.3 meV — and the energy per photon in this region will be even smaller. Hence, the external field only serves to couple different electronic states together, thereby enabling the stimulated rather than spontaneous generation of our commensurate magnon (magnetophonon). The energy to drive the transition, therefore, comes not from the external applied field but rather from internal sources of energy (and therefore, indirectly, the environment), namely the magnon heat bath and the phonon heat bath. Hence, this is a case of photon-assisted magnon–phonon conversion (in contrast to [29], which use phonons under illumination to generate magnons). The applied field enables the transition between different electronic states of the system to occur, whereas the energy comes from internal conversion due to scattering from phonon and magnon modes that have been thermally populated.

Considering that these two heat baths are connected [30,31] and the material is at a constant temperature (in thermal equilibrium), the magnon and phonon bath may be expected to have an equal distribution of energy, under the principle of equipartition. This suggests that there are two types of routes leading to the creation of the *q* → 0 acoustic magnon that corresponds to the MAE. The first of these routes is the "standard" *E* 0 route, which is purely magnonic, and involves no structural relaxations (except as a subsequent potential decay path for the commensurate magnon). The second type of route is the phononic routes identified in Section 2.2 (*E* <sup>+</sup> or *E* −). These routes involve scattering within

the phonon bath to populate a mode that obeys the same symmetries as the commensurate magnon (or in the magnetophonon case, is the same quasiparticle).

These two phononic routes differ in that the *E* <sup>+</sup> route treats the magnon and phonon as a single quasiparticle, whereas the *E* − route treats them as two independent quasiparticles. Which of these regimes dominates (and therefore controls the "fastest switching pathway" MAE for an AFM material) depends on the coupling strength between the magnon associated with the MAE and its associated phonon [31]. For *q* → 0 magnons and phonons, the low momentum leads to a large uncertainty in the position of both quasiparticles and therefore, within the relaxation time approximation, the coupling should be extremely large (→ ∞), leading to *E* <sup>+</sup> being the dominant route. In this strong coupling limit (*E* <sup>+</sup>), it does not make sense to consider magnons and phonons separately, but rather as a single magnetophonon [32].

#### *2.4. Recap*


#### **3. Methodology**

#### *3.1. Density Functional Theory Parameters*

The plane-wave DFT [33,34] code CASTEP [35] was used to simulate the MAE by evaluating the total energy of a conventional unit cell as the orientation of the spin moments on the non-Pt (Fe or Mn) ions — where the majority of the spin in the system is located — were rotated by applying spin constraints. Fully relativistic pseudopotentials from the CASTEP SOC19 pseudopotential library were used with vector spin-modelling and spinorbit coupling in order to capture this effect. Because both FePt and PtMn are metals, the Perdew–Zunger formulation [36] of the local density approximation (LDA) was used and may be expected to perform acceptably well.

Although more sophisticated treatments of exchange-correlation are available for the collinear limit, many do not yet have explicit noncollinear reformulations [37]. However, the methodological improvements suggested in this paper are independent of the choice of exchange-correlation functional, and therefore it is sufficient to demonstrate these benefits for the LDA, despite the fact that other choices might improve other aspects of the physical model.

Numerical convergence of the plane wave basis set was performed to ensure total energy convergence of 0.01 meV/atom (2100 eV plane wave cutoff), and the same tolerance was used to converge the k-space sampling (a 25 × 25 × 25 Monkhorst-Pack grid [38] for PtMn and 20 × 20 × 20 for FePt).

#### Applied Constraints

In order to sample the potential energy surface away from the global minimum, spin constraints were applied to ensure that the desired spin configuration was achieved through a method of Lagrangian minimization [39]. To this end, a constraint that vanishes at the desired spin-configuration but was positive (destabilising) elsewhere was applied. This ensured that the correct value of the potential energy landscape was sampled for the constraint configuration, but also that this was an artificial minimum in the energy landscape such that the minimisation of the total energy could identify this as the appropriate ground state of the system. This technique is conceptually similar to metadynamics [40], and ensures that both the true potential energy surface as well as non-ground-state spin configurations can be sampled.

#### *3.2. Geometry Optimisation Parameters*

The structures were relaxed by using BFGS minimisation until stresses were lower than <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> GPa and forces were less than <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> eV/Å, and a finite basis correction was used to ensure constant basis quality.

During the geometry optimisation, spatial symmetry operations were not enforced. This is because the magnetic ordering lowers the symmetry of the crystal structure (always for the AFM, and for the FM whenever the spins are not aligned along the c-axis), and it was necessary to remove any potential "false symmetries" that may be present which could lead to force cancellation by symmetry.

Additionally, a collinear LDA+U relaxation of the structures was performed, with an effective U correction of 4.5 eV applied to the Mn 3d orbitals and 3.5 eV to the Fe 3d orbitals. The MAE for this LDA+U relaxed structure was then evaluated by using the previous method with no Hubbard correction applied.

#### **4. Results and Discussion**

#### *4.1. Materials*

FePt is an exemplar high-anisotropy ferromagnetic material which has an L10 structure (see Figure 5a). PtMn is the AFM equivalent of FePt, with the only differences emerging from substituting Fe for Mn (see Figure 5b). This substitution yields significant differences to the magnetic structure (which changes from collinear FM to collinear C-type AFM) but not the crystal structure or the magnitude of the spin-orbit coupling introduced by the "heavy metal ion" Pt, which is constant between both.

**Figure 5.** (**a**) The ground-state magnetic and geometric structure of FePt. (**b**) The ground-state magnetic and geometric structure of PtMn. Brown atoms are Fe, purple atoms are Mn and silver are Pt. Note that the major difference between these structures is in the magnetic coupling between the Fe atoms and the Mn atoms.

#### *4.2. Magnetic-Only Transition State*

In order to perform the magnetic-only (LST-like) transition search — corresponding to the standard method with clamped ions — a series of "sweeps" were performed by varying the spin orientation from *θ* = 0 (+~*c*) to *θ* = 180◦ (−~*c*). These sweeps varied by their *φ* orientation (angle to ~*a*) in order to ensure a good coverage of the magnetic phase space. Due to the symmetry of the crystals, the *φ* variation is periodic with a period of 90 degrees. It is also a minor variation compared with the variation in *θ*, and accordingly in the interest of clarity only the *φ* sweeps at 0 and 45 degrees have been reported.

#### 4.2.1. FePt

For FePt, the sweeps are shown in Figure 6. We observe the expected *K* sin<sup>2</sup> (*θ*) behaviour (including higher order terms, as expected from an *ab initio* model). Additionally, the *φ* dependence is observed to be very weak, although with a very slightly softer *φ* = 45◦ axis, which is the minimum energy pathway through magnetic space only for the LDA unit cell. Nevertheless, this *φ*-variation is so weak that FePt may be said to exhibit uniaxial anisotropy with a hard plane.

**Figure 6.** Variation of total energy (in meV / f.u.) with angle from the easy axis, *θ*, for two values of the azimuthal angle, *φ* = 0 and *φ* = 45◦ for FePt. *θ* is expected to be the dominant term under the uniaxial approximation.

The MAE from these results (3.45 meV/f.u., see Figure 6) were found to agree with certain literature results [41] but to be in relatively poorer agreement with others [8]. This was driven by the choice of lattice parameters — in this work we have used the relaxed LDA geometry with no Hubbard-U correction, but other works [8] have used values nearer the bulk experimental lattice parameters (the relaxed LSDA+U geometry). Using the same lattice parameters as these other works (|~*a*<sup>|</sup> <sup>=</sup> <sup>|</sup>~*b*<sup>|</sup> <sup>=</sup> 3.863 <sup>Å</sup> , <sup>|</sup>~*c*<sup>|</sup> <sup>=</sup> 3.783 <sup>Å</sup> ) yielded the same values as they reported (2.68 meV/f.u. using our method compared with 2.68 meV/f.u. for the calculation on the same cell [8]). Thus, the choice of geometry clearly has an effect on the evaluated MAE; the overbinding of LDA leads to shorter "bonds" and therefore increased effective spin-orbit interactions for the Fe orbitals (which have hybridised with Pt orbitals), and also a higher crystal field splitting of the orbitals — thereby increasing the calculated anisotropy. Nevertheless, in order to avoid mixing different approximations to the exchange-correlation functional, we have chosen to evaluate all of our MAEs at the bulk LDA lattice parameters with no Hubbard-U corrections.

In addition, we calculated the variation of the MAE with applied biaxial strain (keeping the~*<sup>c</sup>* lattice vector fixed and measuring the value of the MAE at several <sup>|</sup>~*a*<sup>|</sup> <sup>=</sup> <sup>|</sup>~*b*<sup>|</sup> values) as shown in Figure 7. This clearly demonstrates that the MAE is dependent upon the strain of the system. Although not perfectly linear, there is a clear linear trend, which may be indicative of decreasing effective spin-orbit coupling for electrons predominantly centred on Fe atoms as the overlap between Fe and Pt *d-*orbitals decreases.

**Figure 7.** Change of energy difference between~*c*- and~*a*-aligned spin structures with strain for FePt. This figure includes both the LDA (blue line) and LDA+U (yellow line) structures at their exact in-plane lattice parameter.

#### 4.2.2. PtMn

The angular sweeps for PtMn are shown in Figure 8. The predicted anisotropy per formula unit is lower than for FePt (as expected), and shows the same sin<sup>2</sup> (*θ*) behaviour with strongly uniaxial behaviour. In this case, the *φ* = 0 orientation provides a fractionally lower energy pathway; however, this is not resolvable in Figure 8. We also note that the energy per formula unit is not the same as the energy for the unit cell, although a conventional unit cell is required to model the magnetic structure the formula unit is explicitly 1 Pt and 1 Mn.

**Figure 8.** Variation of total energy (in meV) with angle from the easy axis, *θ*, for two values of the azimuthal angle, *φ* = 0 and *φ* = 45◦ .

In addition to the full angular sweep shown in Figure 8, the MAE was evaluated for various systems with biaxial strain (along~*a* and~*b* axes) applied in order to mimic the effects of growth seed layer choice upon the MAE. These results are shown in Figure 9. The total range of strains evaluated, of the order of ±3%, is large, but not unreasonable in the thin film limit typically used to grow samples. As for FePt, the MAE was evaluated for all systems (including the LDA+U geometry) without applied Hubbard-U. This does not, however, capture the direct effects of the applied-U to the MAE beyond the effects induced by the change in lattice parameters.

**Figure 9.** Change of energy difference between ~*c*- and~*a*-aligned spin structures with strain. Note that the MAE is the absolute value of this number, and the sign change indicates a strain-induced alteration of the easy and hard magnetic axes. The blue line represents the LDA lattice parameter, and the yellow line indicates the LDA+U lattice parameter.

The change of MAE with strain shown in Figure 9 indicates a strong dependence of the MAE on the underlying geometry. These are large changes — and so the underlying assumptions of the phonon-based model of small amplitude oscillations are not applicable — and accordingly a large change in anisotropy is seen. In addition to the obvious ability this gives to engineer thin-film anisotropy based upon a careful choice of seed layer, it should also be able to create electrical control of the MAE by coupling a PtMn layer with a piezoelectric layer.

We also note that these results are in qualitative agreement with other recent theoretical efforts which considered varying both the ~*a* and~*b* lattice parameters [42]. On the other hand, it seems likely that this change is driven by the changing crystal field splitting at the Mn sites driven by the change of interatomic distances rather than the |~*c*|/|~*a*| ratio *per se*.

#### *4.3. Magnetostructural Transition State*

Following the identification of the peak of the magnetic anisotropy (the hard axis), the structure was relaxed in order to perform the QST-like search throughout the orthogonal structure space (as opposed to the spin space). For a maximum, such as the magnetic-only transition state, there are formally no forces, because it is a stationary point. Accordingly, several small initial perturbations were applied to this structure before the relaxation was performed. This is to ensure that the symmetry of the system is broken to enable the unit cell to relax (as the forces are zero at a maximum).

#### 4.3.1. FePt

The relaxed geometry of L10 FePt was found for both the easy and hard magnetic configurations through geometry optimisation of the lattice parameters with the magnetic configuration constrained to the easy and hard axes respectively. The ions remained at high symmetry points. The lattice parameters given by these relaxation are given in Table 1.

**Table 1.** Lattice parameters for relaxed PtMn dependant upon the different spin configurations. LDA+U and experimental lattice parameters are also provided for comparison.


These relaxations show a small (0.02 Å 3 ) increase in total cell volume and a breaking of the symmetry of the unit cell for the hard magnetic configuration. This is in agreement with the previous arguments about magnetostriction in FM materials.

Following relaxation, the total energy of the system for all combinations of these spin and geometry configurations were evaluated (results shown in Table 2).

**Table 2.** Relative energies (in meV/formula unit) for the various configurations of spin and geometry in FePt. The easy geometry is the ground-state geometry of the system, and the hard geometry is the relaxed geometry of the system under the constraint that the system is aligned along the magnetic hard axis.


It is clear that the structural relaxation is a small effect (∼0.5% ) — in line with the ∼0.1% ~*a* lattice parameter change. However, this is per formula unit; for a macroscopic layer of FePt, this effect is multiplied up until these apparently small energy barriers become significant — with even the small difference between easy and hard geometries for the easy spin configuration becoming a barrier larger than *kBT* when 200 conventional unit cells are present.

Finally, it is worth noting that in addition to the thermodynamic arguments here, there are also dynamical considerations. As discussed in Sections 2.2 and 2.3, there are three possible routes from the magnetic easy to the magnetic hard axis: *E* 0 , corresponding to a fixed easy axis geometry; *E* <sup>+</sup> corresponding to a direct change from the easy geometry to the hard geometry; and *E* − corresponding to a fixed hard axis geometry. The attempt frequency for these different routes varies significantly and therefore, in a short time period such as in a data read or write process, the route with the highest attempt rate that is thermally accessible is likely to dominate. In a data-storage scenario, however, the long time period means that the lowest energy route is likely to be favoured.

For FePt these considerations are small. The difference in the barrier between *E* <sup>0</sup> and *E* − — the routes with the largest difference in energy — is less than 1%, making the choice of protocol a smaller effect than other errors associated with the choice of functional or geometry used. This is in accord with the very low energy of the *q* = 0 acoustic phonon mode. For AFM materials, however, the relatively large energy of the optical phonon mode must be considered.

#### 4.3.2. PtMn

The same procedure was carried out for antiferromagnetic L1<sup>0</sup> PtMn, and the structural changes are shown in Table 3. As anticipated, the observed magnetostriction was small (<2 mÅ) and may be attributed to remnant FM couplings between the component of the spins, which are canted by the Dzyaloshinskii–Moriya interaction (DMI). The magnitude of DMI is expected to be small because it leads to the presence of a spin-spiral ground state when it is strong [44]. The largest effect on the structure is a very small movement of the Pt ion in ±~*c*, which corresponds to the 22 meV transverse optical (TO) mode of the system. At first it may seem surprising that the Pt ions should move in response to the reorientation of the spins at the Mn sites, but this may be explained by considering that the movement of the Pt ions affects both the strength and symmetry of the crystal field experienced at the Mn sites. This perturbation changes the energy of individual bands, thereby leading to a reduction in energy of the system when the spins are aligned along the hard axis.

**Table 3.** Lattice parameters and atomic positions for relaxed PtMn dependant upon the different spin configurations.


The key change in the structure is the perturbation of the Pt atoms (unlike the FePt case in which they remained stationary at the high symmetry positions). This change, corresponding to a TO phonon mode of 22 meV, is in agreement with previously reported values of 18 meV (4.3 THz) [42]. DFT LDA calculations (as used here) overbind and hence give higher phonon frequencies, whereas PBE (as used in ref [42]) underbinds and hence gives slightly lower phonon frequencies.

With PtMn, unlike FePt, the energy of the system with spin aligned along the hard axis does not significantly change when the geometry is relaxed (in agreement with the fact that antiferromagnetic systems do not respond strongly to applied fields). However, when the hard-axis relaxed geometry is measured in the easy-axis configuration, a significant increase in the energy is observed, as shown in Table 4.

Now, the displacement that leads to this increase in energy is very small (≈2.9 mÅ). This is actually smaller than the zero-point amplitude of the phonon mode associated with (≈18.7 mÅ). The fact that this displacement is smaller than the zero-point motion which is an effective permanent population of the relevant phonon mode—means that this lower energy route, *E* −, is always achievable. Also unlike FePt, the effect of this small displacement is to significantly reduce the MAE from 0.41 meV/f.u. to 0.32 meV/f.u., which is an ∼20% reduction in the size of the theoretical barrier.

**Table 4.** Relative energies (in meV/f.u.) for the various configurations of spin and geometry in PtMn. The easy geometry is the ground-state geometry of the system, and the hard geometry is the geometry of the system having been relaxed under the constraint that the system is aligned along the magnetic hard axis.


Although the absolute magnitude of these effects will vary with the level of theory used in DFT for the system under investigation — for example, the use of a GGA or Hubbard-U correction might be expected to remove the overbinding of LDA and make the model of the system more accurate—the effect of the magnetostructural coupling on the MAE will remain. Although it is a negligible effect in ferromagnetic systems, this work shows that it is a significant effect for antiferromagnetic systems. Hence, it should be included in all calculations of MAE in AFMs, even with higher levels of theory, if the MAE is to be accurately modelled and not systematically overestimated.

#### **5. Conclusions**

The MAE was evaluated for ferromagnetic FePt and antiferromagnetic PtMn and found to be 3.45 meV/f.u. and 0.41 meV/f.u., respectively, for the LDA approximation in the static-ion limit. The choice of lattice parameters (for example LDA vs LDA+U geometries) was found to provide a significant change in the calculated MAE, up to an exchange of the easy and hard axes of the system. Relaxing the systems to find the magnetostructural rather than just the magnetic transition states yielded a ∼1% reduction in the barrier height of FePt (to 3.44 meV/f.u.), and a ∼20% reduction for PtMn (to 0.32 meV/f.u.) when the lowest energy pathway was considered. This difference is attributed to the difference in coupling between the magnon responsible for the transition and either an acoustic or optical phonon (for FePt and PtMn respectively). This shows the importance of considering magnetostructural effects when evaluating the MAE of systems in which the magnetic structure cannot be adequately expressed in the same primitive unit cell as the geometric structure, such as an antiferromagnet.

**Author Contributions:** Conceptualization and investigation, R.A.L.; writing—original draft preparation, R.A.L.; writing—review and editing, S.J.D.; visualization, R.A.L. and S.J.D.; supervision, resources, project administration and funding acquisition, M.I.J.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** Computational support was provided by the UK national high performance computing service, ARCHER2, for which access was obtained via the UKCP consortium and funded by EPSRC grant ref EP/P022561/1. The authors also acknowledge EPSRC grant EP/V047779/1 for financial support.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in the York Research Database at doi:10.15124/b1ec96de-4796-45a6-89f9-dc0c34d0f431.

**Acknowledgments:** The authors would like to thank Jerome Jackson, Kevin O'Grady and Gonzalo Vallejo Fernandez for their kind comments and valuable conversations in preparing this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Role of Disordered Precursor in L1<sup>0</sup> Phase Formation in FePt-Based Nanocomposite Magnet**

**Alina Daniela Crisan <sup>1</sup> , Ioan Dan <sup>2</sup> and Ovidiu Crisan 1,\***


**\*** Correspondence: ocrisan@infim.ro

**Abstract:** In order to prove the usefulness of having a structurally disordered precursor to the formation of FePt L1<sup>0</sup> phase and to facilitate the co-existence of exchange coupled hard and soft magnetic phases with optimized magnetic properties in various conditions of annealing, a Fe-Pt-Zr-B melt spun alloy has been synthesized and detailed structural and magnetic investigations have been undertaken to probe its phase evolution during annealing. The dynamics of formation of the hard magnetic L1<sup>0</sup> phase during the gradual disorder–order phase transformation has been monitored by using a complex combination of X-ray diffraction methods and <sup>57</sup>Fe Mössbauer spectroscopy methods, over a wide range of annealing temperatures. Multiple phases co-existing in the annealed sample microstructures, observed in XRD, have been reconfirmed by the Mössbauer spectra analysis and, moreover, accurate quantitative data have been acquired in what concerns the relative abundance of each of the observed crystalline phases in every stage of annealing. It is shown that the formation of the hard magnetic phase, emerging from the chemically disordered precursor, is gradual and occurs via complex mechanisms, involving the presence of a disordered Fe-Zr-B-rich intergranular region which contributes to an increase in the abundance of the L1<sup>0</sup> phase for higher annealing temperatures. Magnetic measurements have confirmed the good performances of these alloys in terms of coercivity and remanence. These results contribute to the development of these alloys as the next generation of rare earth, free permanent magnets.

**Keywords:** nanocomposite magnets; L1<sup>0</sup> phase; Mössbauer spectroscopy; magnetic properties

#### **1. Introduction**

Nanocomposite FePt-based magnets have been under scrutiny as a new class of permanent magnets, due to their high corrosion resistance and their high working temperature. The interest in using FePt alloys is based on the disorder–order structural phase transition [1] that occurs in the alloy. FePt usually has a disordered face-centeredcubic fcc A1 structure. Upon annealing, for quasi-equiatomic stoichiometry, the structure undergoes a structural phase transformation, towards the formation of an L1<sup>0</sup> ordered face-centered-tetragonal fct phase that has a high coercivity and a very large magnetocrystalline anisotropy. Non-equilibrium synthesis methods are appropriate for achieving the formation of an L1<sup>0</sup> phase in melt-spun FePt-based ribbons without the need for post-synthesis annealing [2]. Another way to create magnetic alloys with good magnetic properties is to achieve a hard–soft magnetic structure, where the hard magnetic FePt phase is exchange-coupled to soft magnetic ones, and to take full advantage from the high magnetization of the soft phase and the high coercivity of the hard magnetic phase, in order to maximize the energy product (BH)max. To achieve such a phase structure, the initial composition [3] must be modulated by addition of elements, which can provide the necessary additional soft magnetic phase such as boron. This a well-known glass-forming element which introduced in FePt-based metallic alloys allows formation of a disordered intermetallic precursor. Following appropriate annealing, this precursor may transform in

**Citation:** Crisan, A.D.; Dan, I.; Crisan, O. Role of Disordered Precursor in L1<sup>0</sup> Phase Formation in FePt-Based Nanocomposite Magnet. *Magnetochemistry* **2021**, *7*, 149. https://doi.org/10.3390/ magnetochemistry7110149

Academic Editor: Cătălin-Daniel Constantinescu

Received: 1 October 2021 Accepted: 9 November 2021 Published: 14 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

both hard and soft magnetic phases. The advantage of emerging form the same precursor is that the obtained nanocrystalline phases are crystallographically coherent and their exchange coupling is favored. In this way, suitable hard–soft arrangements of phases with well-dispersed nanometric-size grains upon annealing can be obtained. Arrangements of magnetically soft and hard phases for the improvement of hard magnetic properties has been already employed by other researchers. For instance, Chrobak et al. [4] obtained ultra-high coercivity of (Fe86−*x*Nb*x*B14)0.88Tb0.12 bulk nanocrystalline magnets where such arrangements of different magnetic phases were achieved.

Recently, there have been reports on the preparation and characterization of such ternary Fe-Pt-B alloys, using a plethora of non-equilibrium synthesis methods, to obtain both layered structures, such as thin films and multilayers, and bulk alloys, such as melt spun ribbons [5–7]. Magnetic composites of the multilayer type, made of FeB pre-alloy co-doped within FePt deposited films, were studied by Tsai et al. [8] Such exchangecoupled hard–soft system have shown a perpendicular anisotropy with a single switching field. In sandwiched systems made of soft (CoFeB) and hard (FePtB) layers, it has been shown [9] that oxidation during annealing induces stresses through the capping layer, which facilitates greatly the formation of the L1<sup>0</sup> phase in the FePtB hard layer. FePtB melt-spun alloys of various compositions have also been studied in [10]. In this case, the precursor alloy is synthesized in an as-cast state and is then annealed in order to obtain the crystalline state. In this case, hard, tetragonally ordered L1<sup>0</sup> FePt and soft Fe2B phases were formed, but the composition modulation and Pt deficit strongly alters the magnetic performances. Co-sputtered FePt-B multilayers have also been studied [11] and it was found that rapid annealing induces ordering of L1<sup>0</sup> FePt via intermixing with B at the interfaces as well as by the introduction of 0.1 at% Ag. Other synthesis procedures, such as pulsed laser deposition, have also been used to construct exchange-coupled FePtB-based composites [12] and it was found that best magnetic performance has been found for very thick (1.7 microns) FePtB where both hard FePt and soft iron boride phases co-exist. Even lower ordering temperatures with higher thicknesses (3 microns) of co-sputtered FePtB are shown to also produce good magnetic results [13]. Optimizing the microstructure is a factor that strongly influences the overall magnetic behavior. For this, the performed annealing is of great importance. The effect of annealing time on the obtained magnetic properties in FePtB melt-spun alloys was investigated in [14]. The compositional effect of introducing boron on the exchange coupling effect in FePtB has also been investigated [15]. It was shown [16] that from application point of view, such exchange-coupled FePt-FePtB composites show good potential as magnetic recording media.

In our previous work [17], we have shown that in Fe-rich Fe-Pt-Nb-B alloys the microstructure is highly sensitive to the stoichiometry of the as-cast alloy [18]. Boron content of about 8–9% is proven to be not enough to ensure an amorphous-like as-cast state and the samples are obtained in their nanocrystalline A1 structure, while for B content of 18–20%, a disordered structure is obtained in the as-cast state, with large Bragg lines of the A1 structure, as seen in the X-ray diffraction. For samples with a lower Fe content, (Fe0.65Pt0.35)78+*x*Nb2B20−*x*, we have proven that, upon appropriate annealing, very good exchange spring properties are obtained with an energy product (BH)max of 70 kJ/m<sup>3</sup> [19].

In the case of systems exhibiting an L1<sup>0</sup> phase, the high ordering temperature hinders their potential as future nanocomposite RE-free exchange-coupled magnets. It is thought that the reduction of this value can be achieved by adding other elements to the composition, such as Au, Ag, Nb or Zr, with the aim of promoting earlier ordering by segregation to the FePt grain boundaries [20–24] or a glass-forming element (B) to allow the formation of a chemically disordered phase as a precursor for both hard and soft magnetic phases, exchange-coupled to optimize in this way the magnetic properties [25–35]. We have previously shown [2] that the Fe-Pt-Ag-B melt-spun alloy shows a direct formation of the L1<sup>0</sup> phase from the as-cast state without subsequent annealing. It seems quite clear that the different magnetic features in L10-based nanocomposite magnets are strongly influenced by their multiphase nanostructure features.

The present work is dedicated to the detailed structural and morphological study of the microstructure modifications and phase evolution during annealing, as well as the magnetic properties of a Fe-Pt-Zr-B melt-spun alloys, co-existence of exchange coupled hard and soft magnetic phases and the magnetic properties in various conditions of annealing.

#### **2. Experimental**

An alloy with nominal composition Fe65Pt15Zr3B<sup>17</sup> has been synthesized by rapid solidification of the melt. The alloy was synthesized starting with elemental powders and flakes of high purity. They were melted together in an induction furnace with a controlled melting temperature (set at 1300 ◦C). The primary alloy re-melted 3 times to prevent element segregation and to improve its chemical homogeneity. A total amount of 5 g has been used for each sample. The rapid solidification of the melt is performed on a Buhler Melt Spinner SC with protective Ar atmosphere. The obtained melt is purged onto the surface of a Cu wheel. The wheel has 40 cm in diameter and rotates with 2000 rot/min. The melt is flown away from the quartz tube through a circular nozzle of 0.5 mm using Ar pressure of 40 kPa. The size of the nozzle dictates the width of the obtained ribbons. Since the total mass of the experiment was only 5 g, it was necessary to use such a small nozzle. The melt solidifies with a cooling rate of about 10<sup>6</sup> K/min. Away from the wheel, continuous and homogeneous ribbons are obtained. They are about 30 microns thick, 2–3 mm wide and several decimeters long.

In order to obtain the full crystallization and formation of the hard magnetic phase, the as-obtained ribbons were subjected to isothermal annealing. The procedure has been performed at various temperatures, and were chosen at 100 ◦C intervals. The chosen temperatures of annealing were: 500 ◦C, 600 ◦C, 700 ◦C and 800 ◦C. The heating rate was established at 5 K/min and annealing time (time spent at maximum temperature) was set to 30 min. The annealing was performed under a high vacuum (10−<sup>4</sup> mbar) in order to avoid oxidation.

The detailed structural and morphological study of the microstructure modifications, phase evolution during annealing, and magnetic properties of the annealed alloys, was performed by X-ray diffraction (XRD), Mössbauer spectrometry (MS) and magnetic characterizations. For XRD, a Bruker D8 Advance (Bruker AXS GmbH, Karlsruhe, Germany) with Cu K<sup>α</sup> radiation wavelength of 1.54 Å was used. The <sup>57</sup>Fe Mössbauer spectrometry is performed at 300 K and 77 K using a conventional set-up in transmission geometry with a <sup>57</sup>Co source in a Rh matrix. Magnetic characterization has been performed with a SQUID (Superconducting QUantum Interference Device) unit of an MPMS (Magnetic Properties Measurement System) from Quantum Design (Quantum Design Europe GmbH, Darmstadt, Germany), under an applied field of up to 5.5 Tesla, in parallel and perpendicular geometry and temperatures ranging from 4.2 K to 300 K.

#### **3. Results and Discussion**

#### *3.1. As-Cast State*

The ribbon as-cast state has been investigated using XRD and Mössbauer spectrometry. The X-ray diffractogram for the as-cast Fe65Pt15Zr3B<sup>17</sup> sample exhibits very broad Bragg lines. They have been indexed as belonging to the fcc A1 phase. As proven by the lines wide broadening, there is chemical disorder within the A1 fcc structure. The full-profile analysis of the XRD patterns has been analyzed using MAUD (Materials Analysis Using Diffraction) software. The structure of the as-cast sample has disordered fcc A1 symmetry, as resulted from the fittings. Figure 1 shows the X-ray diffractogram fitted with MAUD. It is probable that this structure could be a precursor for the formation of FePt and iron boride phases after appropriate crystallization treatment. We have calculated, based on the integral breadth approach, the crystallite size for all the Bragg peaks indexed in the fcc A1 system. The mean grain size is found to be 3.8 ± 1.4 nm.

based on the integral breadth approach, the crystallite size for all the Bragg peaks indexed

in the fcc A1 system. The mean grain size is found to be 3.8 ± 1.4 nm.

**Figure 1.** X-ray diffractogram of as-cast Fe65Pt15Zr3B17 sample fitted with the MAUD software. The indexed lines belong to fcc A1 FePt. Bottom graph is the difference between fitted and experimental **Figure 1.** X-ray diffractogram of as-cast Fe65Pt15Zr3B<sup>17</sup> sample fitted with the MAUD software. The indexed lines belong to fcc A1 FePt. Bottom graph is the difference between fitted and experimental data.

data. The 300 K and 77 K Mössbauer spectra of the as-cast sample (Figure 2) exhibit broad magnetic sextets, typical of distributed Fe environments encountered in Fe-rich amorphous ribbons. The shape of the spectra confirms the XRD results where a chemically disordered A1 phase has been identified. The hyperfine field (HF) distributions derived from the fitting of the Mössbauer spectra (Figure 3) show a bimodal-type large Gaussian profile distribution, which is characteristic of a disordered Fe environment, with two main environments for Fe. From the numerical fitting of the size distribution (Figure 3) with two Gaussian profiles, we have obtained for the low field mode an average HF of 16 T and 18 T for 300 K and 77 K respectively and for the high field mode average HF of 26 T and 29 T for 300 K and 77 K respectively. The relative proportion of the high-field to low-field relative contributions to the HF distributions is about 3:1. From the HF values we can presume that the low field mode corresponds to a disordered precursor that would give The 300 K and 77 K Mössbauer spectra of the as-cast sample (Figure 2) exhibit broad magnetic sextets, typical of distributed Fe environments encountered in Fe-rich amorphous ribbons. The shape of the spectra confirms the XRD results where a chemically disordered A1 phase has been identified. The hyperfine field (HF) distributions derived from the fitting of the Mössbauer spectra (Figure 3) show a bimodal-type large Gaussian profile distribution, which is characteristic of a disordered Fe environment, with two main environments for Fe. From the numerical fitting of the size distribution (Figure 3) with two Gaussian profiles, we have obtained for the low field mode an average HF of 16 T and 18 T for 300 K and 77 K respectively and for the high field mode average HF of 26 T and 29 T for 300 K and 77 K respectively. The relative proportion of the high-field to low-field relative contributions to the HF distributions is about 3:1. From the HF values we can presume that the low field mode corresponds to a disordered precursor that would give rise upon annealing to a boride (possibly Fe2B) phase, while the high field mode would give rise to an fcc FePt phase.

#### rise upon annealing to a boride (possibly Fe2B) phase, while the high field mode would give rise to an fcc FePt phase. *3.2. Isothermal Annealing*

#### 3.2.1. XRD Analysis

*3.2. Isothermal Annealing*  3.2.1. XRD Analysis In order to obtain the full crystallization and formation of the hard magnetic phase, isothermal annealing has been performed at various temperatures, chosen at 100 °C intervals, for observing better the degree of formation of the crystalline phases as well as to monitor the evolution of the disorder–order phase transformation. The chosen temperatures of annealing were: 500° C, 600° C, 700° C and 800° C. The heating rate was estab-In order to obtain the full crystallization and formation of the hard magnetic phase, isothermal annealing has been performed at various temperatures, chosen at 100 ◦C intervals, for observing better the degree of formation of the crystalline phases as well as to monitor the evolution of the disorder–order phase transformation. The chosen temperatures of annealing were: 500 ◦C, 600 ◦C, 700 ◦C and 800 ◦C. The heating rate was established at 5 K/min and annealing time (time spent at maximum temperature) was set to 30 min. The annealing was performed under a high vacuum (10−<sup>4</sup> mbar) in order to avoid oxidation.

lished at 5 K/min and annealing time (time spent at maximum temperature) was set to 30 min. The annealing was performed under a high vacuum (10−4 mbar) in order to avoid oxidation. The structure of the resulting samples was characterized by XRD, Mössbauer spectrometry and their magnetic properties by vibrating sample magnetometry and SQUID. The X-ray diagrams of the annealed samples are plotted in Figure 4. For a better comparison, the as-cast X-ray diffractogram is also added. In the as-cast state, the broad Bragg peaks of the disordered A1 FePt-rich solid solution have been identified, as previously discussed. The broad line centered at around 25◦ is due to the sample holder and it can be seen that this broad feature is preserved at all annealing temperatures, until 800 ◦C,

where the intensity of the main Bragg peak is very high and the broad feature appears strongly diminished due to scaling. At 500 ◦C, the A1 fcc pattern is better formed, with a smaller linewidth than in the as-cast state. Besides these lines, the main superlattice peaks of the tetragonal L1<sup>0</sup> phase appear already formed. At 600 ◦C the main superlattice peaks of the tetragonal L1<sup>0</sup> phase are more pronounced and the process of peak narrowing is furthermore continued, as expected. Starting with 600 ◦C a small Bragg peak, attributed to the main Bragg reflection of the Fe2B tetragonal phase, is observed. *Magnetochemistry* **2021**, *7*, x FOR PEER REVIEW 5 of 13 *Magnetochemistry* **2021**, *7*, x FOR PEER REVIEW 5 of 13


**Figure 2.** 300K and 77K Mössbauer spectra of as-cast Fe65Pt15Zr3B17 sample. **Figure 2.** 300K and 77K Mössbauer spectra of as-cast Fe65Pt15Zr3B<sup>17</sup> sample. **Figure 2.** 300K and 77K Mössbauer spectra of as-cast Fe65Pt15Zr3B17 sample.

**Figure 3.** 300K and 77K Mössbauer HF distribution of as-cast Fe65Pt15Zr3B17 sample. **Figure 3.** 300K and 77K Mössbauer HF distribution of as-cast Fe65Pt15Zr3B17 sample. **Figure 3.** 300K and 77K Mössbauer HF distribution of as-cast Fe65Pt15Zr3B<sup>17</sup> sample.

The structure of the resulting samples was characterized by XRD, Mössbauer spectrometry and their magnetic properties by vibrating sample magnetometry and SQUID.

The structure of the resulting samples was characterized by XRD, Mössbauer spectrometry and their magnetic properties by vibrating sample magnetometry and SQUID. The X-ray diagrams of the annealed samples are plotted in Figure 4. For a better compar-

ison, the as-cast X-ray diffractogram is also added. In the as-cast state, the broad Bragg peaks of the disordered A1 FePt-rich solid solution have been identified, as previously discussed. The broad line centered at around 25° is due to the sample holder and it can be seen that this broad feature is preserved at all annealing temperatures, until 800 °C, where the intensity of the main Bragg peak is very high and the broad feature appears strongly diminished due to scaling. At 500 °C, the A1 fcc pattern is better formed, with a smaller linewidth than in the as-cast state. Besides these lines, the main superlattice peaks of the tetragonal L10 phase appear already formed. At 600 °C the main superlattice peaks of the

peaks of the disordered A1 FePt-rich solid solution have been identified, as previously discussed. The broad line centered at around 25° is due to the sample holder and it can be seen that this broad feature is preserved at all annealing temperatures, until 800 °C, where the intensity of the main Bragg peak is very high and the broad feature appears strongly diminished due to scaling. At 500 °C, the A1 fcc pattern is better formed, with a smaller linewidth than in the as-cast state. Besides these lines, the main superlattice peaks of the tetragonal L10 phase appear already formed. At 600 °C the main superlattice peaks of the

phase.

tetragonal L10 phase are more pronounced and the process of peak narrowing is further-

**Figure 4.** X-ray diffractograms of the Fe65Pt15Zr3B17 sample, in its as-cast state as well as annealed at **Figure 4.** X-ray diffractograms of the Fe65Pt15Zr3B<sup>17</sup> sample, in its as-cast state as well as annealed at various temperatures (500 ◦C, 600 ◦C and 800 ◦C).

various temperatures (500 °C, 600 °C and 800 °C). Then at 800 °C the lines of the crystalline phases are even narrower, and the tetragonal L10 phase is the most predominant in the sample. The XRD diagram at 800 °C shows a completely crystallized pattern and besides the L10 phase, the Fe2B is also visible in the Then at 800 ◦C the lines of the crystalline phases are even narrower, and the tetragonal L1<sup>0</sup> phase is the most predominant in the sample. The XRD diagram at 800 ◦C shows a completely crystallized pattern and besides the L1<sup>0</sup> phase, the Fe2B is also visible in the diagram.

diagram. It has to be mentioned that the two main superlattice peaks of tetragonal L10, the (0 0 1) and (1 1 0) reflections occurs for the Cu Kα radiation we used, at about 24° and 32° respectively. From the observation of the Bragg line linewidths, it can be seen that the main Bragg lines observed in the fcc FePt solid solution in the as-cast state undergo a refinement process as the annealing temperature is increased, from bottom to top. The main 5 Bragg lines of the fcc A1 FePt solid solution in the as-cast state, from lower to higher angles, are, in order, attributed to: (1 1 1), (2 0 0), (2 2 0), (3 1 1) and (2 2 2). Upon annealing, these main reflections become overlapped with the main Bragg reflections of the L10 FePt tetragonal phase. It is very interesting to observe that already, from the first annealed sample at 500 °C, the (2 0 0), (2 2 0) and (3 1 1) peaks split into two. It seems that It has to be mentioned that the two main superlattice peaks of tetragonal L10, the (0 0 1) and (1 1 0) reflections occurs for the Cu Kα radiation we used, at about 24◦ and 32◦ respectively. From the observation of the Bragg line linewidths, it can be seen that the main Bragg lines observed in the fcc FePt solid solution in the as-cast state undergo a refinement process as the annealing temperature is increased, from bottom to top. The main 5 Bragg lines of the fcc A1 FePt solid solution in the as-cast state, from lower to higher angles, are, in order, attributed to: (1 1 1), (2 0 0), (2 2 0), (3 1 1) and (2 2 2). Upon annealing, these main reflections become overlapped with the main Bragg reflections of the L1<sup>0</sup> FePt tetragonal phase. It is very interesting to observe that already, from the first annealed sample at 500 ◦C, the (2 0 0), (2 2 0) and (3 1 1) peaks split into two. It seems that the peaks from the fcc phase are now spectrally separated from those of the tetragonal phase.

the peaks from the fcc phase are now spectrally separated from those of the tetragonal For the other two annealing treatments, at 700 °C (Figure 5) and 800 °C, the separation continues and the peaks become more and more narrow and better separated between one and another (see for instance the fitting from Figure 5 for the sample annealed at 700 For the other two annealing treatments, at 700 ◦C (Figure 5) and 800 ◦C, the separation continues and the peaks become more and more narrow and better separated between one and another (see for instance the fitting from Figure 5 for the sample annealed at 700 ◦C). The specific signatures of the tetragonal L1<sup>0</sup> phase, the superlattice peaks (0 0 1) and (1 1 0), are clearly observed for all the annealed samples at around 24◦ and 33◦ (in 2-theta), respectively.

°C). The specific signatures of the tetragonal L10 phase, the superlattice peaks (0 0 1) and (1 1 0), are clearly observed for all the annealed samples at around 24° and 33° (in 2-theta), respectively. The good separation of the spectral lines allowed us an accurate determination by numerical fitting of the peak positions and linewidths, a calculation of the lattice parameters, the tetragonality or ordering parameter c/a and of the mean grain sizes for each of the observed crystalline phases. The obtained results are schematized in Table 1.

It can be seen that the lattice parameters, as revealed from the fit of the experimental spectra, do not change drastically with the increase of the annealing temperature; however, a steady increase of the ordering parameter c/a, associated with the tetragonality or the degree of ordering of the L1<sup>0</sup> phase, is observed upon increasing the annealing temperature.

Fe65Pt15Zr3B17

**Figure 5:** MAUD fitting of the XRD diagram of the sample annealed at 700 °C for 30 min. **Figure 5.** MAUD fitting of the XRD diagram of the sample annealed at 700 ◦C for 30 min.

observed crystalline phases. The obtained results are schematized in Table 1.

The good separation of the spectral lines allowed us an accurate determination by numerical fitting of the peak positions and linewidths, a calculation of the lattice parameters, the tetragonality or ordering parameter *c/a* and of the mean grain sizes for each of the



**(nm) (nm) a (Å) c (Å) c/a a (Å)**  As-cast - - - 3.8411 ± 0.02416 - 3.8 ± 1.4 500 °C 3.8542± 0.0043 3.7127± 0.0067 0.9632 3.8315 ± 0.0308 15 ± 2 17 ± 3 600 °C 3.8537 ± 0.0027 3.7134 ± 0.0042 0.9636 3.8268 ± 0.0284 30 ± 3 36 ± 4 700 °C 3.8524 ± 0.0049 3.7147 ± 0.0071 0.9642 3.8379 ± 0.0247 46 ± 5 49 ± 5 800 °C 3.8516± 0.0052 3.7166± 0.0051 0.9649 3.8292 ± 0.0303 57 ± 3 54 ± 3 The grain size for both L10 and A1 FePt phases in the annealed samples increases continuously from the as-cast values 15 nm (17 nm respectively) as the annealing temper-The grain size for both L1<sup>0</sup> and A1 FePt phases in the annealed samples increases continuously from the as-cast values 15 nm (17 nm respectively) as the annealing temperature increases. At 800 ◦C annealing, the grain size in the L1<sup>0</sup> FePt phase is about 54 nm, while for the same annealing in the A1 FePt phase the average grain size is estimated to be around 57 nm. In conclusion, in the annealed samples, using structural data, we have unambiguously proven the coexistence of the hard and soft magnetic phases with well-refined grain sizes. These observations may be made quantitative with the help of the Mössbauer spectroscopy, where the relative fraction of each of the phases present in the samples can be obtained.

#### ature increases. At 800 °C annealing, the grain size in the L10 FePt phase is about 54 nm, 3.2.2. Mössbauer Analysis

while for the same annealing in the A1 FePt phase the average grain size is estimated to be around 57 nm. In conclusion, in the annealed samples, using structural data, we have unambiguously proven the coexistence of the hard and soft magnetic phases with wellrefined grain sizes. These observations may be made quantitative with the help of the Mössbauer spectroscopy, where the relative fraction of each of the phases present in the samples can be obtained. <sup>57</sup>Fe Mössbauer analysis on all the samples has been performed using a conventional setup in transmission geometry. The Mössbauer source we used was a <sup>57</sup>Co source embedded in a Rh matrix. Measurements were carried out at an ambient temperature. The annealed Fe65Pt15Zr3B<sup>17</sup> samples have been investigated and their Mössbauer spectra, recorded at −10 +10 mm/s velocity range, have been fitted, with a procedure involving several subspectra, magnetic sextets. All the hyperfine parameters, i.e, the hyperfine field (HF), quadrupole splitting (QS) and the isomer shift (IS), were kept as free parameters during the fitting. In Figure 6 the Mössbauer spectra of three of the annealed Fe65Pt15Zr3B<sup>17</sup> samples are presented together with the subspectra, as identified from the fitting. All the hyperfine field parameters obtained from the fit are depicted in Table 2. The components used for the fitting are of the same color as their hyperfine parameters, for easier identification.

3.2.2. Mössbauer Analysis

cation.

**Figure 6:** Mössbauer spectra of the Fe65Pt15Zr3B17 samples, annealed at 500 °C, 600 °C and 700 °C together with the fitted subspectra. **Figure 6.** Mössbauer spectra of the Fe65Pt15Zr3B<sup>17</sup> samples, annealed at 500 ◦C, 600 ◦C and 700 ◦C together with the fitted subspectra.

57Fe Mössbauer analysis on all the samples has been performed using a conventional setup in transmission geometry. The Mössbauer source we used was a 57Co source embedded in a Rh matrix. Measurements were carried out at an ambient temperature. The annealed Fe65Pt15Zr3B17 samples have been investigated and their Mössbauer spectra, recorded at −10 +10 mm/s velocity range, have been fitted, with a procedure involving several subspectra, magnetic sextets. All the hyperfine parameters, i.e, the hyperfine field (HF), quadrupole splitting (QS) and the isomer shift (IS), were kept as free parameters during the fitting. In Figure 6 the Mössbauer spectra of three of the annealed Fe65Pt15Zr3B17 samples are presented together with the subspectra, as identified from the fitting. All the hyperfine field parameters obtained from the fit are depicted in Table 2. The components used for the fitting are of the same color as their hyperfine parameters, for easier identifi-

**Table 2.** The obtained R-factors, as resulted from the full-profile refinement XRD patterns of annealed samples. Profile Rp, weighted profile Rwp and goodness-of-fit χ2 are shown. **Table 2.** The obtained R-factors, as resulted from the full-profile refinement XRD patterns of annealed samples. Profile Rp, weighted profile Rwp and goodness-of-fit χ <sup>2</sup> are shown.


800 °C 3.54 2.05 1.29

Mössbauer spectra of the annealed samples show different and more complex hyperfine features than the as-cast ones. Indeed, the Mössbauer patterns consist of sextets with much narrower lines, which indicate that the samples are crystallized. The fitting model chosen to adjust the spectra of annealed samples takes into account the XRD results where a microstructure consisting of multiple phases was proven. The spectra have been fitted Mössbauer spectra of the annealed samples show different and more complex hyperfine features than the as-cast ones. Indeed, the Mössbauer patterns consist of sextets with much narrower lines, which indicate that the samples are crystallized. The fitting model chosen to adjust the spectra of annealed samples takes into account the XRD results where a microstructure consisting of multiple phases was proven. The spectra have been fitted with a number of subspectra that were indexed and assigned to the phases identified in the XRD.

The Mössbauer spectra of the sample annealed at 500 ◦C has been fitted with 5 subspectra, a fitting model that corresponds to the situation observed in the XRD data. The two subspectra designed by the red line in Figure 6 have been assigned to the A1 fcc FePt phase. The hyperfine parameters are depicted in Table 3. The quadrupole shift value suggests a cubic symmetry, while the hyperfine field of 30 T is typical for A1 FePt phase which is consistent with the results on the as-cast sample. The relative proportion of the fcc phase, as resulted from the fit, is 18%. The HF parameters obtained are consistent with the results reported in the literature [36] and with those obtained from the investigation on the as-cast sample. The subspectrum designed by the blue line in Figure 6 has been assigned to the L1<sup>0</sup> FePt phase. The quadrupole splitting of 0.27 mm/s suggests a tetragonal symmetry and the hyperfine field of 27.6 T corresponds to the value found for the L1<sup>0</sup> phase in the literature. For the annealing at 500 ◦C, this phase is the most abundant in the sample. Its

relative proportion, as revealed from the fit, is 28%. This indicates that the disorder–order phase transition has started already at 500 ◦C and an important fraction of the initial fcc FePt solid solution is already transformed into tetragonal L1<sup>0</sup> FePt. The pink lines are the components of the boride phases. Compositionally, there is an Fe excess which does not crystallize in the FePt phases; therefore, in principle, during annealing there are conditions that Fe and B form borides, as there is a large affinity for these two elements to alloy. Fe2B has been already observed in the sample annealed at 500 ◦C and Fe3B is a metastable phase that may be formed in intermediate annealing temperatures. Their cumulate relative intensity amounts to 14%. The central part of the spectrum was fitted with an HF distribution between 3.5 and 18.5 T. These contributions are represented to be convoluted as one contribution (green) in Figure 6.


**Table 3.** Hyperfine parameters as resulted from fitting of the Mössbauer spectra (Figure 6).

The estimated errors are: ±0.02 mm/s for IS and QS/2ε, ±0.1 T for Bhf and ±1 for the relative proportion.

The average hyperfine field of these low-field contributions is 9.9 T and its relative intensity amounts to 40%. From the HF value we can assume that this low-field contribution to the Mössbauer spectrum may be attributed to the disordered Fe-Zr-B-rich phase. The quadrupole shift close to 0 sustains the argument of a cubic symmetry of this phase. The low hyperfine field of this contribution suggests that it is a Fe-poor phase. From all of the arguments we consider that the green contribution may be a disordered Fe-Zr-B phase with few Pt atoms interstitially inserted or randomly occupying the Fe sites. Most of this phase is probably located in interfacial regions between the nanocrystalline grains and the in-grain boundaries. We have to emphasize that the disordered phase, as fitted in the Mössbauer spectra, accounts for all Fe-distributed ionic chemical environments. For the sake of fitting consistency, every contribution that does not belong to the magnetic subspectra (sextets) is included in one disordered sublattice. In Mössbauer terms this is a hyperfine field distribution. Usually, a hyperfine field distribution is used to fit the Mössbauer spectra of Fe-disordered components or Fe-amorphous alloys. Therefore, the disordered part, illustrated by the green component, represents in fact the contribution of all of the chemically distributed Fe sites, each of them with different probabilities, as shown, for instance, in Figure 3. We have named it disordered Fe-Zr-B, to account in fact for all of the disordered Fe chemical environments. This is in agreement with other interpretations of Mössbauer spectra in Fe-Zr-B alloys [37].The Mössbauer spectra of the

sample annealed at 600 ◦C and 700 ◦C have the same complex pattern as the previous one, but with better resolved lines. The spectra were fitted using the same fitting model. The parameters were consistent throughout the series of measurements. From Table 3 one can observe the phase evolution with the annealing conditions, as well as the relative abundance of the A1 and L1<sup>0</sup> FePt phases. The tetragonal phase increases slightly up to 44%, while Fe3B decomposes. It can be seen that Fe2B relative abundance does not increase, therefore we may assume that from the decomposition of Fe3B, the obtained Fe helps the formation of L1<sup>0</sup> FePt and B is incorporated into the Fe-Zr-B-rich residual phase. With an increasing annealing temperature, all the phases that decompose help the formation of the L1<sup>0</sup> phase, which is the most ordered one, in the ribbon microstructure.

#### 3.2.3. Magnetic Properties

The hysteresis loops for the Fe65Pt15Zr3B<sup>17</sup> samples, as-cast and annealed at 600 ◦C and 700 ◦C, have been recorded at 300 K with a magnetic field applied parallel to the ribbon plane, and are plotted in Figure 7. The hysteresis loop for the as-cast ribbons is typical for soft ferromagnets. The magnetization saturates almost immediately after applying a small magnetic field. The saturation reaches 1.3 T, while the loop shows a very low hysteresis (around 30 mT). The annealed samples show, however, a completely different behavior. The magnetization has a slower approach to saturation than in the as-cast state. This may be explained by the presence of a significant fraction of hard L1<sup>0</sup> grains. The value of the maximum magnetization is about 15% lower than that of the as-cast sample, reaching about 0.9 T. Upon decreasing the field, a very large value of remanent magnetization is noticed, as is the case for hard magnetic materials. The sample annealed at 600 ◦C shows a large openness of the loop. *Magnetochemistry* **2021**, *7*, x FOR PEER REVIEW 11 of 13 highest (BH)max values of 64 kJ/m3 are obtained for the sample annealed at 700 °C. All the magnetic parameters obtained for these samples at 300 K in the parallel field are depicted in Table 4.

**Figure 7:** RT hysteresis loops for the as-cast and annealed samples. **Figure 7.** RT hysteresis loops for the as-cast and annealed samples.

**Table 4.** Magnetic parameters obtained for annealed samples at 300 K in parallel applied field. **Annealing Conditions Hc (kA/m) μ0MS (T) μ0Mr (T) Mr/Ms (BH)max (kJ/m3)**  As-cast 2.94 1.26 1.2 0.95 - 500 °C 194 1.1 0.85 0.77 61 600 °C 847 0.9 0.63 0.70 64 700 °C 877 0.9 0.63 0.70 64 The remanent magnetization is in this case 0.63 T, thus giving a remanence-to-saturation ratio of about 0.7. The shape of the hysteresis loop presents two vaguely observed inflection points (at around −190 and −840 kA/m applied field), typical for systems with both hard and soft magnetic phases in co-existence; however, here the hard magnetic phase is predominant and the exchange coupling between the hard and the soft magnetic phases is not complete. As a consequence, a quite strong coercive field is observed. The coercivity reaches in this case 847 kA/m, which is more than 3 times increased, compared with the annealing at 500 ◦C. The switching field of the soft magnetic phase is smaller than the coercive field, which indicates a better exchange coupling between the grains. The value

the non-equilibrium melt spinning method. The stoichiometry was chosen in such a way as to obtain a chemically disordered as-cast state. From this state, by appropriate annealing, we have been able to obtain a mixture of hard L10-FePt and soft A1 FePt and boride magnetic phases. The dynamics of formation of the hard magnetic L10 phase during the gradual disorder–order phase transformation has been monitored by using a complex combination of X-ray diffraction methods and 57Fe Mössbauer spectroscopy methods over a large range of annealing temperatures. Multiple phases co-existing in the annealed sample microstructures, observed in XRD, have been reconfirmed by the Mössbauer spectra analysis and moreover, accurate quantitative data have been acquired in what concerns the relative abundance of each of the observed crystalline phases in every stage of annealing. It has been shown that the formation of the hard magnetic phase, emerging from the chemically disordered as-cast precursor, is gradual and occurs via complex mechanisms involving the presence of a disordered Fe-Zr-B-rich intergranular region, as well as the decomposition of iron borides, providing excess Fe, which contributes furthermore to the increase of the hard magnetic phase abundance, for higher annealing temperatures. Magnetic measurements have confirmed the good performances of these alloys in terms of

**4. Conclusions** 

of the switching field is 155 kA/m for the soft magnetic grains, increasing substantially compared with the annealing at 500 ◦C. The sample annealed at 700 ◦C presents an almost similar hysteresis loop as the sample annealed at 600 ◦C. While the coercivity is slightly higher than before (877 kA/m), the remanent and the maximum magnetization are almost identical to the previous annealing, with a remanent-to-saturation ratio of about 0.7. The energy product for all the samples was calculated using the B-H loops. The highest (BH)max values of 64 kJ/m<sup>3</sup> are obtained for the sample annealed at 700 ◦C. All the magnetic parameters obtained for these samples at 300 K in the parallel field are depicted in Table 4.


**Table 4.** Magnetic parameters obtained for annealed samples at 300 K in parallel applied field.

#### **4. Conclusions**

We have synthesized a nanocomposite magnet based on the FePt-Zr-B alloy, by using the non-equilibrium melt spinning method. The stoichiometry was chosen in such a way as to obtain a chemically disordered as-cast state. From this state, by appropriate annealing, we have been able to obtain a mixture of hard L10-FePt and soft A1 FePt and boride magnetic phases. The dynamics of formation of the hard magnetic L1<sup>0</sup> phase during the gradual disorder–order phase transformation has been monitored by using a complex combination of X-ray diffraction methods and <sup>57</sup>Fe Mössbauer spectroscopy methods over a large range of annealing temperatures. Multiple phases co-existing in the annealed sample microstructures, observed in XRD, have been reconfirmed by the Mössbauer spectra analysis and moreover, accurate quantitative data have been acquired in what concerns the relative abundance of each of the observed crystalline phases in every stage of annealing. It has been shown that the formation of the hard magnetic phase, emerging from the chemically disordered as-cast precursor, is gradual and occurs via complex mechanisms involving the presence of a disordered Fe-Zr-B-rich intergranular region, as well as the decomposition of iron borides, providing excess Fe, which contributes furthermore to the increase of the hard magnetic phase abundance, for higher annealing temperatures. Magnetic measurements have confirmed the good performances of these alloys in terms of coercivity and remanence, a maximum energy product of 64 kJ/m<sup>3</sup> being obtained for annealing at 700 ◦C. These findings open good perspectives for a future class of permanent magnets made from these alloys.

**Author Contributions:** Conceptualization, O.C.; methodology, I.D.; writing—original draft preparation, A.D.C.; writing—review and editing, O.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded through POC Project: "Materiale multifunct,ionale inteligente pentru aplicat,ii de înaltă tehnologie (MATI2IT)" (code MySMIS 105726) funded by the Romanian Ministry of European Funds.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are not publicly available due to IPR protection measures.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Electric Field Control of Magnetic Properties by Means of Li<sup>+</sup> Migration in FeRh Thin Film**

**Gengfei Li 1,2, Yali Xie 1,2,\* , Baomin Wang 1,2 , Huali Yang 1,2 and Run-Wei Li 1,2,3**


**Abstract:** Recently, the electric control of magnetism by means of ion migration has been proven to be effective with nonvolatility and low energy consumption. In this work, we investigated the control of the magnetic properties of FeRh films by means of Li<sup>+</sup> migration in FeRh/MgO heterostructures. We found that the migration of Li<sup>+</sup> could reduce the phase transition temperature by 2 K with an applied voltage of 1 V. Meanwhile, the voltage-dependent saturated magnetization exhibited a repetitive switching behavior from high to low magnetization values while the voltage was switched from 4 to <sup>−</sup>4 V, indicating that the migration of Li<sup>+</sup> in the FeRh film can be reversible. This provides a means to control the magnetic properties of FeRh films.

**Keywords:** electric control of magnetic phase transition; FeRh; Li<sup>+</sup> migration

#### **1. Introduction**

Fallot et al. discovered that the CsCl-type FeRh alloy undergoes a first-order phase transition from an antiferromagnetic (AF) to ferromagnetic (FM) phase at around room temperature accompanied by an approximately 1% volume expansion of the crystal lattice [1]. Due to this magnetic transition, FeRh has attracted extensive attention for its potential applications in heat-assisted magnetic recording, AF memory, and magnetic refrigeration [2–7]. In order to realize its applications under different conditions, many studies have been conducted to tune the magnetic properties of the FeRh film by means of the strain, magnetic field, electric field, etc. [8]. However, the manipulation of FeRh properties by means such as the strain and magnetic field is difficult to implement in spintronics devices, limiting its practical applications [9,10]. In recent years, the electric control of magnetism has been proven to be effective with low energy consumption and repeatability [11–16]. For example, Cherifi et al. found that a moderate electric field can produce a large magnetization variation in FeRh/BaTiO<sup>3</sup> [14]. Lee et al. reported a large resistivity modulation by applying an electric field to FeRh/PMN-PT heterostructures [16]. However, these methods usually require a high-quality interface or very thin film. Recently, the electric control of magnetism by means of ion migration has been proven to be effective with nonvolatility and low energy consumption [17–19]. For instance, Dhanapal et al. demonstrated reversibly controlled magnetic domains of Co via electric-field driven oxygen migration at the nanoscale [18]. In this work, we investigated the electric control of magnetic properties of FeRh films by means of Li<sup>+</sup> migration in FeRh/MgO heterostructures. We found that this method can not only manipulate the phase transition temperature, but also change the saturated magnetization of FeRh film.

**Citation:** Li, G.; Xie, Y.; Wang, B.; Yang, H.; Li, R.-W. Electric Field Control of Magnetic Properties by Means of Li<sup>+</sup> Migration in FeRh Thin Film. *Magnetochemistry* **2021**, *7*, 45. https://doi.org/10.3390/ magnetochemistry7040045

Academic Editors: Cătălin-Daniel Constantinescu and Lucian Petrescu

Received: 1 March 2021 Accepted: 25 March 2021 Published: 26 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **2. Experiment**

First, 30 nm thick FeRh films were deposited by sputtering from an FeRh alloy target using a DC power of 100 W and an argon pressure of 5 mTorr. The (100)-oriented MgO substrates were preheated to 530 ◦C for 1 h in a vacuum chamber and then held at this temperature during deposition. After growth, the samples were annealed at 700 ◦C for about 1 h at a base pressure below 1.0 <sup>×</sup> <sup>10</sup>−<sup>5</sup> Pa in order to obtain the chemically ordered CsCl-type FeRh. The film thicknesses were calibrated by X-ray reflectivity (XRR) and the rate of growth was 5 nm/min. The surface topographies of FeRh films were measured using a Bruker Icon atomic force microscope. The crystal structure was characterized using a Bruker Discover X-ray diffractometer. The temperature-dependent magnetizations of the FeRh films were determined using a superconducting quantum interference device (SQUID) magnetometer. The electric field was applied using a commercial three-electrode experimental setup [20]. In this work, FeRh/MgO, a standard calomel electrode, and a platinum sheet were used as the working electrode, reference electrode, and counter electrode, respectively. In addition, 0.1 M of LiClO<sup>4</sup> dissolved in propylene carbonate was used as the electrolytic solution. When voltage was applied, the Li<sup>+</sup> could migrate into or out of the FeRh films. When we performed magnetic measurements, the FeRh/MgO sample was taken out of the experimental setup and the electrode setup was also removed.

#### **3. Results and Discussion**

An AFM image of the FeRh films is displayed in Figure 1a, which indicates that the roughness of the film was very small (~0.3 nm). Figure 1b shows the XRD *θ*–2*θ* pattern of the FeRh film. As seen in Figure 1b, only the reflections at 29.90◦ and 62.11◦ could be observed, which demonstrates that the obtained FeRh film had the required CsCl-type structure. The X-ray ϕ-scan of FeRh/MgO in Figure 1c indicates that the FeRh films were epitaxially grown with a 45◦ in-plane lattice rotation with respect to the MgO(001) substrates. The temperature-dependent magnetization of the FeRh film measured with an in-plane magnetic field of 2 kOe is plotted in Figure 1d. When warming up, the magnetization steeply increased from 100 to 1.1 <sup>×</sup> <sup>10</sup><sup>6</sup> A/m at around 350 K, exhibiting a typical characteristic of the transition from the AF to the FM phase in FeRh. By contrast, the magnetization of the FeRh film gradually decreased to the original value of the AF state in the cooling process. A temperature hysteresis in the magnetization of about 25 K could be observed, indicating the first-order phase transition of FeRh.

**Figure 1.** (**a**) The AFM image of the FeRh films. (**b**) *θ*-2*θ* scan and (**c**) ϕ-scans of FeRh/MgO. (**d**) The temperature-dependent magnetization of the FeRh film measured with an in-plane magnetic field of 2 kOe.

Figure 2a shows the temperature-dependent magnetization of the FeRh film at 2 kOe without and with applied voltages of 1 V through the electrochemical method. It can be seen that with an applied voltage of 1 V, the behavior of the temperature-dependent magnetization was almost unchanged compared with that for 0 V. Meanwhile, the phase transition shifted to the lower temperature after the applied voltage of 1 V. The voltagedependent critical phase transition temperatures of FeRh, which were defined as *T*AFM-FM on heating and *T*FM-AFM on cooling processes, are plotted in Figure 2b. *T*AFM-FM exhibited a decrease of about 2 K with an applied voltage of 1 V, indicating that the Li<sup>+</sup> migrated into the FeRh film. However, no further change was observed after increasing the applied voltages. *T*FM-AFM behaved similarly to that of *T*FM-AFM. This kind of behavior may be ascribed to the small size of Li<sup>+</sup> . When voltage was applied, the Li<sup>+</sup> could migrate into the FeRh film and the lattice volume could expand, which have been confirmed in our previous work [21]. The magnetization and phase transition temperature of FeRh are strongly dependent on the lattice parameter [22–24]. Consequently, the *T*AFM-FM and *T*FM-AFM change due to the migration of Li<sup>+</sup> . However, the lattice volume cannot change by very much, due to the small size of Li<sup>+</sup> , leading to a modulation of about 2 K in *T*AFM-FM and *T*FM-AFM when a voltage of 1 V was applied, as well as an easy saturation.

**Figure 2.** (**a**) Temperature-dependent magnetization without and with an applied voltage of 1 V for the FeRh thin film at 2 kOe. (**b**) Voltage-dependent *T*AFM-FM and *T*FM-AFM of FeRh on heating and cooling processes.

Figure 3a shows the reversibility of the temperature-dependent magnetization of the FeRh film measured with applied voltages of 4 and −4 V. The applied magnetic field was 2 kOe. It can be seen that when 4 V voltages were applied, the saturated magnetization of the first temperature-dependent magnetization (1) at 400 K could reach about 0.8 <sup>×</sup> <sup>10</sup><sup>6</sup> A/m. Meanwhile, when voltages of <sup>−</sup>4 V were applied, the saturated magnetization of the first temperature-dependent magnetization (2) at 400 K decreased to about 0.7 <sup>×</sup> <sup>10</sup><sup>6</sup> A/m, indicating that applying an electric field of <sup>−</sup>4 V can lead to the decrease in the saturated magnetization compared with that for a applied voltages of 4 V. When applied voltages of 4 and −4 V were applied again, the change in the saturated magnetization (3 and 4) was repeated. This kind of behavior can be seen more clearly in the inset of Figure 3a. It can also be ascribed to the migration of Li<sup>+</sup> . When 4 V voltages were applied, the Li<sup>+</sup> can migrate into the FeRh film. However, when voltages of −4 V were applied, the Li<sup>+</sup> migrated out of the FeRh film, leading to a decrease in the saturated magnetization. The saturated magnetizations under different voltages are summarized in Figure 3b. It can be seen that the voltage-dependent saturated magnetization exhibited a repetitive switching behavior from high to low magnetization values when switching the voltage from 4 to <sup>−</sup>4 V repetitively, indicating that the migration of Li<sup>+</sup> in FeRh film can be reversible.

**Figure 3.** (**a**) Temperature-dependent magnetization with applied voltages of 4 and −4 V for FeRh thin films. The inset shows the enlarged figure of the temperature-dependent magnetization. (**b**) Voltage dependence of saturated magnetization measured with a magnetic field of 2 kOe in FeRh thin films.

#### **4. Conclusions**

In summary, the magnetic properties of FeRh films can be controlled by means of electrochemistry in FeRh/MgO heterostructures. Both the phase transition temperature and saturated magnetization of FeRh films can be modulated by varying the voltage. This kind of modulated behavior can be ascribed to the migration of Li<sup>+</sup> ions by applying voltage. This work provides a new method to control the magnetic properties of FeRh films.

**Author Contributions:** Conceptualization, Y.X. and B.W.; Methodology, Y.X. and B.W.; Investigation, G.L., Y.X., H.Y.; Writing—Original Draft Preparation, G.L.; Writing—Review and Editing, Y.X. and B.W.; Supervision, B.W.; Funding Acquisition, Y.X., B.W. and R.-W.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a National Key Technologies R&D Program of China (2016YFA0201102), the National Natural Science Foundation of China (51871233, 51871232, 51931011, and 51525103), the Youth Innovation Promotion Association of the Chinese Academy of Sciences (2019299), the Ningbo Science and Technology Bureau (2018B10060), and the Ningbo Natural Science Foundation (2019A610054, 2019A610051).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data generated or analyzed during this study are included in this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Feasibility Study of Cooling a Bulk Acoustic Wave Resonator by Nanoparticle Enhanced Phase Change Material**

**Mohammad Yaghoub Abdollahzadeh Jamalabadi**

Chabahar Maritime University, Chabahar 99717-56499, Iran; my.abdollahzadeh@cmu.ac.ir

**Abstract:** In the current study, the coupling of a cooling problem with the electromagnetic resonance of a bulk acoustic wave (BAW) material is investigated. As well, a new cooling method by the addition of nanoparticles to a phase change material surrounding the BAW resonator is presented. To solve the governing equations of piezoelectric charge and momentum balance, thermal balance, and fluid flow a code with the method of finite element is introduced. After validation of various features of the code with melting profile, heat generation, charge curve, and dispersion curve with benchmarks, the eigenfrequency analysis of the system is done. The thermal behavior of the system at first mode and various boundary conditions are studied. As well, the effect of nanoparticles in fastening the cooling of the BAW resonator is demonstrated.

**Keywords:** bulk acoustic waves; filter design; thermal management; MEMS; phase change material; nanofluid; electronic packaging

**Citation:** Abdollahzadeh Jamalabadi, M.Y. Feasibility Study of Cooling a Bulk Acoustic Wave Resonator by Nanoparticle Enhanced Phase Change Material. *Magnetochemistry* **2021**, *7*, 144. https://doi.org/ 10.3390/magnetochemistry7110144

Academic Editor: C˘at˘alin-Daniel Constantinescu and Lucian Petrescu

Received: 6 September 2021 Accepted: 15 October 2021 Published: 28 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Nanotechnology could be a technique that searches about the behavior of materials on a nanometric scale to improve their abilities [1–5]. Nanomaterials have many industrial and engineering applications. One of the current challenges related to nanotechnology lies in how nanomaterials are produced [2]. In quickly creating a universe of IoT, the RFFE of shrewd gadgets should deal with higher information rates and access the full data transmission of 4G/5G remote innovation. The higher bandwidth and data rate of 5G requires smarter technology. The explanation behind this is the developing requests of pervasive low dormancy information (data should transfer fast) at higher working frequencies needed to oblige improved information transmission abilities and quickly develop quantities of clients [3]. Usual electro-acoustic resonators utilized for RF channel applications depend on two strategies for the excitation and proficient energy catching of an acoustic wave (1200 to 4000 m/s) proliferating in a piezoelectric material. The first strategy and technique depends on SAW, where the surface acoustic wave is created on a piezoelectric substrate by metal IDTs by all accounts. The second strategy and technique depends on BAW, where the mass acoustic wave in the bulk of the material is energized by the use of an electric field through terminals above and under a dainty piezoelectric plate [4]. Key execution markers and quality index for resonator configuration are frequency of resonance, quality factor, coupling coefficient, and impedance for RF channel applications [5].

A bulk acoustic wave (BAW) device is usually composed of a piezoelectric material (i.e., zinc oxide or aluminum nitride) between two electrodes [6–10]. The BAW devices are currently used in Wi-fi systems, wireless positioning systems, cell phones, and satellite navigation. The range of working frequency (10<sup>8</sup> Hz to <sup>2</sup> <sup>×</sup> <sup>10</sup><sup>10</sup> Hz) benefits the 6G technologies in small sizes (10−<sup>6</sup> m) [7]. For acoustic isolation of the BAW device from the surrounding or substrate, two methods usually are used: air cavity (free-standing) and (solidly mounted). The free-standing method is cheaper than using Bragg reflector [8]. Such acoustic mirrors usually used a series of coupled low and high acoustic impedance [9]. The thickness of the Bragg reflector is equal to the quarter wavelength for maximum acoustic reflectivity [10].

Many numerical methods are used in the layout design of the BAW devices to reduce the production costs and increase the accuracy of the device for requested resonance modes [11]. The network of BAW resonators makes a frequency filter that protects the electronic device from unwanted frequencies. When electrostatic discharge can occur in the system, in case of high power requests, lower insertion loss requests, high-frequency requests, (*<sup>f</sup>* <sup>&</sup>gt; 1.5 2.5 <sup>×</sup> <sup>10</sup><sup>9</sup> Hz), or in case of steep stopband attenuation requests, the SAW device is replaced by a BAW device [12]. The design of a BAW device requires a scattering profile, impedance, stage, and Q-factor [13]. The scattering profile distinguishes all the potential modes, energy-related with every mode, the get over between the primary method of interest and other plate modes when the BAW resonator is energized. The impedance and stage contour measure the non-idealities in the BAW resonator plan and the general effect on the Q-values [14]. BAW resonator-based circuits require energy control in the piezoelectric layer dependent on the ideal choice of Bragg reflector layer, piezoelectric film, and terminal thickness.

Figure 1 shows the BAW device geometry which characterizes the ideal resilience on the piezoelectric film and electrodes. The thickness of cathode and anode could be different [15]. That thicknesses will affect the frequency resonances of the system. Five significant decisions in the plan and advancement of high-quality acoustic wave resonator are numbered in the list as follows:


**Figure 1.** Schematic of piezoelectric and electrodes in a BAW configuration.

The solution of governing equations is not straightforward. This is because of the intricacy of the piezoelectric conditions that depict the gadget and its stacking conditions, delivering scientific arrangements of the conditions for non-trivial 2D and 3D calculations [16]. The first challenge in the plan of thin-film resonators is the concealment of side resonances that can be energized around the recurrence of the ideal mode shapes. The BC applied by IDTs and electrodes can add many extra modes. Some tough layers are added at the end of fingers to control such effects. In any case, the anode measurements directed by electrical coordinating regularly are with the end goal that precise outcomes do not warrant 1D or 2D display. By 1D or 2D analysis, the effect of fingers and spurious modes or fake resonances cannot be detected [17]. Along these lines, one should examine the fake resonances through full 3D FEM calculations. Be that as it may, the solid advantage of FEM is in its natural ability to oblige convoluted calculations, various materials, piezoelectricity, and full precious stone anisotropy. Subsequently, just a few disentangling approximations are required, bringing about an exact and adaptable reproduction strategy [18].

Various nanoparticles such as Ag, TiO2, CuO, ZnO, Fe3O4, Indium, SiO2, SWCNTs, and MWCNTs, can be mixed by a base fluid throughout sonication to make some improvements in thermo-physical properties [9]. Usual PCMs in the industry are paraffin, wax, and poly-a-olefin [19]. In the process of sonication which takes from 30 min up to 24 h, some surfactants (such as ethanol and sodium dodecyl sulfate) are added for better mixing. In Figure 2, the geometry of the BAW resonator assisted by PCM box is presented.

**Figure 2.** Geometry of the BAW resonator assisted by PCM box in 2D simulation.

As the thermal management and temperature compensation are important, the thermal design of BAW resonators is considered in its design [20]. Since the purpose of the current work is the use of NEPCM for cooling of BAW resonators in its design has significance in BAW development [21]. The main aim of the work is the investigation of nano-added PCM in thermal management of BAW which has the principal conclusions of durability and stability in electronic packaging application [22]. A novel design for PCM is a topic in BAW cooling [23]. Phase change of NEPCM could be used in BAW cooling [24].

Rational design calculations of surface acoustic wave devices were performed in many types of research [25]. However, the use of latent heat of phase change in a thermal energy storage systems especially PCM including nanoparticles for cooling on a smallscale is a new idea. In some batteries for temperature control purposes, the PCM and heat pipes are used throughout the discharge–charge cycle [26]. Property (temperaturedependent) variations in thermal energy storage and Brownian motion can affect freezing and melting problems in PCM. The big issues of thermal stability and compatibility within the context of materials are critical in latent heat storage systems. Before using for practical cases life-cycle assessment of phase change material is required. In such systems which are usually assisted by finned heat pipes (with highly conductive metal foams), hightemperature is a common problem [27]. Nowadays phase change materials are used as thermal energy storage for buildings. Thermal performance enhancement of NePCMs for solar energy heat exchangers and accumulators were tested. Since they can work with multiple phase change materials [28]. In this paper, a finite element method (FEM) code is developed for the mathematical arrangement of the electro–elastic conditions that oversee the directly constrained piezoelectric vibrations. Both methodologies, that of taking care of the field issue (consonant investigation) and that of taking care of the relating eigenvalue issue (modular examination), are portrayed. A FEM programming bundle has been made without any preparation. Significant angles fundamental to the proficient execution of FEM are clarified, for example, memory of the executives and tackling the summed up piezoelectric eigenvalue issue. Calculations for diminishing the necessary PC memory through enhancement of the framework profile, just as matrix calculation for the arrangement of the eigenfrequency issue are connected into the product from outer mathematical estimations. Current FEM programming is applied to to solve the mathematical governing equations of a thin-film mass BAW. Specifically, 3D reproductions are utilized to explore the impact of the top terminal shape. The legitimacy of the displayed strategy is shown by contrasting the recreated and estimated removal profiles at a few frequencies. The outcomes show that valuable data on the exhibition of the flimsy BAW acquired by generally large cross-sections and, subsequently, some assets. The novelty of the current study is that the NEPCM has not been used previously for cooling thin-film bulk acoustic wave resonance. The performance of that method of cooling and its side effects on electrical and mechanical parts are investigated in this paper.

#### **2. Materials and Methods**

In this section governing equations and boundary conditions of the system are presented.

#### *2.1. Cooling Part*

The governing equations of fluid, the non-slip boundary conditions of fluid on walls, and the boundary condition for thermal energy are presented in Tables 1–3 respectively. Enhancement of solidification rate of latent heat thermal energy storage using corrugated fins (shell-and-tube) was proved but here just a simple PCM cube is used. In some systems the PCM slurries' density changes at phase change temperature, but here a Newtonian model is used. The nanoparticles dispersed in a phase change material can also improve melting characteristics which is neglected here.

**Table 1.** Governing equations of fluid [2].


**Table 2.** Non-slip boundary conditions of fluid on walls [4].


**Table 3.** Thermal boundary conditions of fluid [8].


Figure 1 shows the resonator configuration. The geometry of the BAW cell is equal to one of that circuit lines. Figure 2 illustrate the resonator and cooling part. The effective characteristics are shown in Table 4 while the material thermophysical properties are presented in Table 5. The effect of the addition of nanoparticle on the thermal behavior of PCM is obtainable from Table 4. Many fluids and solids such as organic materials, carb-oxy-methyl cellulose carbon foams, paraffin, and expanded graphite are used as phase change material previously for a thermal storage system. However, here TH29 was selected as [2,29].

Graphs of heat flow (in joule per second) as a function of temperature which is called in literature the differential scanning calorimeter (DSC), usually provided by PCM producers. TH29 (Calcium-chloride hexahydrate) performance in DSC form is presented in [29] for the pure material and its mixing with cellulose. As well, the TH29 shows good performance in aged cases which exposes to ambient air for two weeks. From the mass measurements, the

hydrated salt-based TH29 absorbed moisture at about 50% of their weight and lost their ability in eight days.

As well, various nanoparticles were used such as Ag, carbon nanofiber, Al, aluminum, and carbon nanotubes. However, here copper nanoparticle was used. The governing equations and thermal boundary conditions for heat transfer in solid part are presented in Table 6.



**Table 5.** Material thermo-physical properties [2].


Usually, the thermal delay method is used to determine the effective characteristics of micro-encapsulated phase change material. Table 5 mentions TH29, a commercially available PCM.

**Table 6.** Governing equations and thermal boundary conditions for heat transfer in solid part [6].


#### *2.2. Piezoelectric Induced Wave*

The piezoelectric constitutive conditions that couple the mechanical and electrical amounts in the piezoelectric material are communicated utilizing the coupling equations. The conservation of momentum for the bulk acoustic wave in an electrical part is presented in Table 7.



The electrical boundary condition for the piezoelectric wave is depicted in Figure 3. The wave motion (*ρ ∂* 2*ui ∂t* <sup>2</sup> = *∂Tij ∂x<sup>j</sup>* ) which comes from conservation of momentum and presented in in Table 7. The coupled constitutive equations in strain-charge form [6] are

$$
\begin{bmatrix}
\begin{matrix} S\_{11} \\ S\_{12} \\ S\_{13} \\ S\_{22} \\ S\_{23} \\ S\_{33} \end{matrix} \\
\begin{bmatrix}
\mathbf{s}\_{11}^{E} & \mathbf{s}\_{12}^{E} & \mathbf{s}\_{13}^{E} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{s}\_{22}^{E} & \mathbf{s}\_{23}^{E} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{s}\_{33}^{E} & \mathbf{s}\_{33}^{E} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{s}\_{44}^{E} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{s}\_{55}^{E}
\end{bmatrix}
\begin{bmatrix}
T\_{11} \\ T\_{12} \\ T\_{13} \\ T\_{22} \\ T\_{23} \\ T\_{23}
\end{bmatrix}
\begin{bmatrix}
0 & 0 & d\_{31} \\ 0 & 0 & d\_{32} \\ 0 & 0 & d\_{33} \\ 0 & 0 & d\_{33} \\ 0 & 0 & d\_{33} \\ 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\frac{\partial \phi}{\partial t} \\ \frac{\partial \phi}{\partial t} \\ \frac{\partial \phi}{\partial t} \\ \frac{\partial \phi}{\partial t}
\end{bmatrix}
\end{bmatrix}
$$

and

$$
\begin{bmatrix} D\_1 \\ D\_2 \\ D\_3 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 & d\_{15} & 0 \\ 0 & 0 & 0 & d\_{25} & 0 & 0 \\ d\_{31} & d\_{32} & d\_{33} & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} T\_{11} \\ T\_{12} \\ T\_{13} \\ T\_{22} \\ T\_{23} \\ T\_{33} \end{bmatrix} - \\ \tag{29}
$$

$$
\begin{bmatrix} \varepsilon\_{11}^T & 0 & 0 \\ 0 & \varepsilon\_{22}^T & 0 \\ 0 & 0 & \varepsilon\_{33}^T \end{bmatrix} \begin{bmatrix} \frac{\partial \phi}{\partial x} \\ \frac{\partial \phi}{\partial y} \\ \frac{\partial \phi}{\partial z} \end{bmatrix}
$$

Consider next the case where the slab is electrically excited by applying a harmonic voltage, v = *V*0*e jt*, across the electrodes. In this case, it is convenient to consider a slab of length L with ends at x = L/2. The stresses at the ends are zero. The fields of interest are expressed by eigenfrequency analysis (−*ρω*2*<sup>u</sup>* <sup>=</sup> *∂T*<sup>11</sup> *∂x* ). The dielectric permittivity, elastic and piezoelectric constants of ZnO are considered in this study.

As the velocity of the electromagnetic wave is faster than the acoustic wave, the electrostatic equation of Laplace is applied for electric potential (∇2*<sup>φ</sup>* <sup>=</sup> 0). As well, the electrical displacement is the function of strain and the electric field (*D<sup>i</sup>* = *diklSkl* + *ε T ikE<sup>k</sup>* ).

**Figure 3.** Schematic of BAW piezoelectric material (ZnO) and periodic condition (for real and imaginary k), terminal, and ground boundary condition at aluminum surfaces.

#### **3. Results**

Here, predominantly, the instance of voltage stacking is considered, i.e., one terminal is electrically grounded, and a sinusoidal voltage with endorsed consistent plentifulness is associated between the grounded and "hot" anode. Arrangement of the eigenvalue issue at that point brings about arrangement reverberation frequencies (short-circuited terminals) at which the electrical impedance expects low qualities. The equal reverberation frequencies with enormous estimations of impedance can be acquired from the open-circuited case. For the open-circuit limit conditions, the potential is set consistent on the "hot" anode, yet not fixed (drifting potential).

In this section which describes the results, first, the confrontation of the obtained results with other research are presented through validations. Such a suitable comparison significantly raises the meaning of the presented paper. Then the effect of adding PCM to BAW is discussed.

#### *3.1. Validation of Used Code*

In this section, various parts of the code are validated against well-known benchmarks and commercial software.

#### 3.1.1. Validation of Fluid Part

For the purpose of the validation of the fluid part, the previous work of Jamalabadi and Park [2] is used. The vertical axis in Figure 4 is solid content that has no unit. As revealed in Figure 4 the solid percent of PCM decreased by the increase of time (melting process). As seen in Figure 4 both methods are in good agreement.

**Figure 4.** Validation of fluid part by comparing the solid percent of PCM versus time in a melting process.

#### 3.1.2. Validation of Heat Generation

In this section, the validation of heat generation caused by resonance with the analytical solution is performed. The loss mechanisms in the thin film BAW resonators include leakage to the substrate, laterally escaping waves, ohmic losses, viscous losses, dielectric losses, and wave scattering. Additionally mentioned in the literature are eddy current losses and losses due to the formation of conducting channels in the high resistivity substrate interface. The analytical solution of the piezoelectric governing equation for heat generation is given in Appendix A. To calculate the heat generation of the electric source, both mechanical part (elastic energy) and electric part (electric displacement energy) should be considered. The Poynting vector (W/m<sup>2</sup> ) which is the directional energy flux (the energy transfer per unit area per unit time) usually used as a symbol of conservation of electromagnetic energy (*S* = *E* × *H*). The basic idea is that Poynting energy losses some parts as passed through the piezoelectric material and electrodes. Here the power loss distribution when the slab is electrically excited by an applied voltage is considered. The electrical excitation with the amplitude of 10 V is applied across the terminals. Material and geometrical parameters for the case of heat generation are presented in Figure 5 and Table 8.

**Table 8.** Material and geometrical parameters for the case of heat generation.


**Figure 5.** Geometrical parameters for the case of heat generation in a slab.

The comparison of the current code with the analytical solution which is given in Appendix B is plotted in Figure 6. As seen in Figure 6 both numerical and analytical results are in good agreement. The validation results of Figure 6 has the error of simulation near zero as it designed to respond in 1D and the effect of transverse direction is not too much as the electrode extends all over the top and bottom surface. In engineering applications, the electrodes cover the surface partially and some materials are added to help the system to have these ideal mechanical support connections.

**Figure 6.** Heat generation in axial direction calculated by FEM and analytical solution.

#### 3.1.3. Validation of Charge Curve

In this section, the 3D problem is considered while in another part of the paper you can consider the 2D model as the third direction has not affected the results. Simulation of a conventional BAW resonator is performed in the appendix code presented at the end of the paper. After a variable definition on substrate, electro, piezo, and resonator parameters. Element types of solid226 and solid95 are used to model the piezoelectric parts. Materials parameters of aluminum nitride for various directions as anisotropic material in the form of tabular data are given. Some damping parameters are also considered for mechanical and electrical losses. Silicon and molybdenum data are entered as well. For building geometry, the bottom substrate, bottom electrode, piezo layers, and top electrode simple block are used. 3D regular mesh is created to apply the mechanical boundary conditions (two symmetry faces, right surface from top view symmetry, and up to surface from top view symmetry, down surface from top view hard boundary as a fixed plate) and electrical boundary conditions (electrode bottom, electrode top). The results of harmonic analysis in a range of 1100 and 1250 Mhz are used to validate the charge collected on the electrode surface. Material parameters definition and geometry of the compared case are presented in Tables 9 and 10, respectively. As seen the piezoelectricity and full anisotropy of the materials are taken into account. The comparison of the current code with the ANSYS-APDL code results which the code is provided in Appendix B is plotted in Figure 7. As seen in the figure, both results are in good agreement. The validation results of Figure 7 have the error of simulation near zero. As both used the same method (FEM) these results were anticipated.


**Table 9.** Material parameters definition.


**Figure 7.** Accumulated charge on electrode as a function of frequency estimated by FEM and ANSYS.

#### 3.1.4. Validation of Dispersion Curve

The behavior of the side resonances in a multilayer structure is a function of Lamb modes curves traveling in a multilayer structure is investigated in this section. The comparison of a dispersion curves presented by Makkonen et al. [3] is given in Figure 8. As depicted in the Figure 8, both methods are in good agreement.

#### *3.2. Eigenfrequency Analysis*

The first six modes of the BAW problem are presented in this part. The real part of the first six modes of BAW device is 221.4200 222.02, 222.97, 224.16, 225.51, and 226.96 MHz. The imaginary parts are three orders of magnitude less than the real parts with the values of 83.593, 84.982, 84.029, 90.701, 90.162, and 90.142 kHz. Detailed parameters are given in Table 11. To have a better understanding of given parameters they have been plotted in Figure 9. As shown, the coefficients have a fluctuation from one frequency to another. In participation factors of x and z, even frequencies (2, 4, 6, . . .) have higher values and y participation factors have maximum in odd frequencies (1, 3, 5, . . .).


**Table 11.** Eigenfrequency study of six modes of BAW resonator.

**Figure 9.** Eigenfrequency, participation factors, and effective modal masses of the first six modes of BAW piezoelectric in eigenfrequency evaluation.

Figures 10–15 show horizontal and vertical displacement of mode 1 to mode 3. Each of propagating modes and evanescent modes can also be visualized by its displacement field. As seen in Figure 10, the horizontal displacement of mode 1 has two peaks at the middle of terminals. Figure 11 presents four peaks in the vertical displacement of mode 1 where one of them is in opposite direction. Horizontal displacement of mode 2 as expected has four (2 × *N*) peaks which are symmetric around the middle of piezoelectric which is visible in Figure 12. As illustrated in Figure 13, vertical displacement of mode 2 has four maximum and four minimum which are aligned in three parallel lines. The presence of propagating and evanescent modes is clear. The exhibition of six (2 × *N*) simple maximum of mode 3 of horizontal displacement of mode 3 is demonstrated in Figure 14. As shown the eigenfunctions are aligned in three parallel lines and are symmetric conditions. Finally, the eigenfunctions of mode 3 of vertical displacement with six visible maximum are depicted in Figure 15. A detailed BAW design, general discussion on BAW resonators, and electric potential modes are discussed in [1,6]. In surface plotted figures including Figures 10–15 the color is automatically proportional to surface height.

**Figure 10.** Horizontal displacement of mode 1.

**Figure 11.** Vertical displacement of mode 1.

**Figure 12.** Horizontal displacement of mode 2.

**Figure 13.** Vertical displacement of mode 2.

**Figure 14.** Horizontal displacement of mode 3.

**Figure 15.** Vertical displacement of mode 3.

The damping ration, quality factor, and *S*<sup>11</sup> are plotted in Figures 16–18 respectively. Damping ratio (ratio of imaginary of the frequency to its absolute value) and quality factor (ratio of absolute value to two times of the imaginary of the frequency) as a function of frequency are plotted in Figure 16. The quality factor for frequency as a function of frequency is depicted in Figure 17. The real part and the imaginary part of the scattering parameter as a function of frequency are illustrated in Figure 18 for the layer structure of a BAW resonator. Here the quality plot is against all selected intervals of frequency. The dispersion curve of imaginary (evanescent mode) and real (propagating mode) parts of the wavenumbers are figured in Figure 19 for a thin-film BAW composite resonator. The dispersion schematic includes both imaginary and real parts of the wavenumber. The plots actually are the contours of admittance |*Y*11| for real part of k, and contours of admittance |*Y*22| for imaginary the part of k. Figure 19 shows the calculated dispersion curve. The results for the the real and the imaginary parts of the wavenumber are combined by taking the positive *<sup>k</sup>*-axis with an optimum near (<sup>6</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> , 2 <sup>×</sup> <sup>10</sup><sup>8</sup> ) for the real part and the negative *<sup>k</sup>*-axis with a maximum near (−<sup>2</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> , 2 <sup>×</sup> <sup>10</sup><sup>8</sup> ) for the imaginary part. The horizontal axis of dispersion curve is *k*, wave number, with the unit of <sup>1</sup> *<sup>m</sup>* while the vertical axis is *f* , frequency, with aunit of <sup>1</sup> *s* . They correspond to the lower branch for the *TS*<sup>2</sup> mode (thickness-shear), the upper branch for the *TE*<sup>1</sup> mode (thickness-extension), and the left branch for the evanescent modes.

**Figure 16.** Damping ratio (blue dashed-line) and quality factor (red line) as a function of frequency.

**Figure 17.** Quality factor for frequency as a function of frequency.

**Figure 18.** Real part (blue line) and imaginary part (red dashed-line) of scattering parameter as a function of frequency.

**Figure 19.** Dispersion curve for imaginary wave numbers (evanescent waves) and for real wave numbers (propagating waves).

Total power dissipation as a function of frequency is plotted in Figure 20 for the specific device structure thin-film BAW composite resonator considered here. As shown, the maximum value of heat generation happened in the first mode while is in accordance with maximum admittance (*Y*11). This graph will be used in the next part as the input for heat generation.

**Figure 20.** Total power dissipation as a function of frequency.

#### *3.3. Fluid and Thermal Analysis*

Profile development and contours of velocity magnitude distribution in the liquid part of the PCM cell are plotted in Figure 21. As shown, the liquid profile is developed fast near the initial condition as the absorbed heat in the PCM is converted to phase change in the liquid–solid boundary. Figure 21 does not have a color legend and is focused on showing the melting profile development. In the red zone, the velocity magnitude is maximum while in blue zones (near walls) the velocity magnitude is near zero. As time goes on, the melting process developed in the PCM box. As shown, the velocity is minimum near walls and is maximum at the middle of two flow channels created horizontally at the PCM cell. Detailed velocity distribution in PCM cell with cooling at the end of melting process is plotted in Figure 22. As shown in velocity vectors and magnitude in the PCM field, the maximum velocity happens at the middle of two flow channels with the horizontal direction. At the end of the flow stream, it turns vertically to turn inside by natural convection.

**Figure 21.** Profile development and contours of velocity magnitude in liquid part of PCM cell (**a**) t = 10, (**b**) t = 100, (**c**) t = 200, (**d**) t = 500, (**e**) t = 1000, (**f**) t = 2000, (**g**) t =5000, (**h**) t = 10,000.

**Figure 22.** Velocity vectors and magnitude in PCM field.

The power distribution inside BAW slab for first and second mode are plotted in Figure 23. The distribution is axially and matches the analytic solution given in Appendix A.

The temperature distribution inside BAW slab for first and second mode are plotted in Figure 24. The distribution is axially and for the case of constant temperature BC at both sides.

**Figure 23.** Power distribution inside BAW slab for first and second mode.

**Figure 24.** Temperature distribution inside BAW slab for first and second mode.

Temperature difference (∆*T* = *T* − *T*∞) contours are plotted in Figure 25. Figure 25 does not have a color legend and is focused on showing the melting profile development. In the red zone, the velocity magnitude is maximum while in blue zones (near walls) the velocity magnitude is near zero. The figure shows the temperature profiles through the time for various BC. The initial condition of PCM cell is temperature equals the environmental temperature (*T* = *T*∞). As time goes on, the cold zone is shortened through the volume, and the heated zone expands. If both sides are kept at environmental temperature the temperature profile inside the slab is parabolic. The case of constant temperature BC at both sides happens when the periodic structure of the BAW + PCM cell forced the melting temperature at both side until the end of phase change. The case of adiabatic BC at the left happens when each BAW + PCM cell is considered as an isolated domain from the other. Such conditions force the melting temperature on the right side until the end of phase change while adiabatic condition at the left of the cell. As shown, the periodic structure causes less temperature increase.

Temperature distribution inside BAW + PCM cell is illustrated in Figure 26. Time dependency of the temperature in the BAW + PCM cell shows that, by considering PCM, the thermal inertia happens in the system.

**Figure 25.** Temperature profiles inside liquid part of PCM cell (**a**) constant temperature BC at the both sides, (**b**) maximum temperature while constant temperature BC at the both sides, (**c**) adiabatic BC at the left, (**d**) maximum temperature while adiabatic BC at the left.

**Figure 26.** Temperature distribution inside BAW + PCM cell (**a**) profile plot (**b**) surface plot.

The effect of nanoparticle concentration on time dependency of maximum temperature is illustrated in Figure 27. By increasing the nanoparticle the rate of interface convection increased and faster phase change happens.

**Figure 27.** Effect of nanoparticle concentration on maximum temperature.

#### **4. Conclusions**

In this paper, a FEM code is introduced for the mathematical arrangement of the electro–elastic conditions that oversee the directly constrained vibrations of piezoelectric media. A consonant time reliance is accepted. Both of the methodologies, that of taking care of the field issue (consonant investigation) and that of taking care of the relating eigenvalue issue (modular examination), are portrayed. A FEM programming bundle has been made without any preparation. Significant angles fundamental to the proficient execution of FEM are clarified, for example, the memory of the executives and tackling the summed up piezoelectric eigenvalue issue. Calculations for diminishing the necessary PC memory through enhancement of the framework profile, just as Lanczos calculation for the arrangement of the eigenvalue issue are connected into the product from outer mathematical libraries. Current FEM programming is applied to nitty-gritty mathematical demonstrating of thin-film bulk acoustic wave (BAW) composite resonators. Correlation of results from 2D and full 3D recreations of a resonator is introduced. Specifically, 3D reproductions are utilized to explore the impact of the top terminal shape on the resonator electrical reaction. The legitimacy of the displaying strategy is shown by contrasting the recreated and estimated removal profiles at a few frequencies. The outcomes show that valuable data on the exhibition of the flimsy film resonators can be acquired even with generally coarse cross-sections and, subsequently, moderate computational assets. To highlight the advantages and disadvantages of using NEPCM for cooling application of BAW resonator it is good to notice that when comparing the current devised solution with other solutions from the scientific literature, the efficiency of the current design is more clear. While in common circuits the heat is removed from the semiconductor by conduction to support or with convection or radiative heat transfer modes to the environment, here the phase change plays the role of the heat sink with more thermal inertia. As well, the use of nanoparticles facilitates the heat removal, and also decreases the delay time and thermal inertia of the cooling system. As it would be useful to add information on further research of the authors related to the continuation of this research topic, future research could be conducted to consider the 3D design of a thin-film bulk acoustic wave composite resonators in conjunction with another electrical part to see the effectiveness of current study in an engineering case.

**Funding:** This research received no external funding.

**Acknowledgments:** The author want to acknowledge Humboldt officials for Research Fellowship Program for Experienced Researchers and Ing. Ulrich Nieken, head of Institut für Chemische Verfahrenstechnik at Universität Stuttgart, for hosting during writing of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Nomenclature**


)


#### **Appendix A. Analytical Solution in One Dimension**

The current model as presented in Figure 5 simulates a uniform layer structure extending horizontally to infinity. The computed dispersion diagram of such an infinite plate can be used to optimize device designs. For this tutorial, one can take a thin slice from the center of the complete device and use the periodic boundary condition to extend it laterally to infinity. If we apply *φ*(*z* = *h*) = *φme <sup>j</sup>ω<sup>t</sup>* at the top electrode while bottom electrode is ground (*φ*(*z* = 0) = 0) then from Laplace equation in electrostatic approximation

*∂*

$$\frac{\partial^2 \phi}{\partial z^2} = 0 \tag{A1}$$

we have

$$
\phi = \phi\_m e^{j\omega t} \frac{z}{h} \tag{A2}
$$

thus the electric field in z-direction is

$$E\_3 = -\frac{\phi\_m e^{j\omega t}}{h} \tag{A3}$$

if we differentiate from the first equation of piezoelectric

$$S\_{11} = s\_{11}T\_{11} - d\_{31} \frac{\partial \phi}{\partial z} \tag{A4}$$

respect to *x* direction, then

$$\frac{\partial S\_{11}}{\partial \mathbf{x}} = s\_{11} \frac{\partial T\_{11}}{\partial \mathbf{x}} \tag{A5}$$

Considering the definition of strain from displacement field (*S*<sup>11</sup> = *<sup>∂</sup><sup>u</sup> ∂x* ) the following equation for the differential of the stress respect to *x* direction is obtained

$$\frac{\partial T\_{11}}{\partial \mathbf{x}} = \frac{1}{s\_{11}} \frac{\partial^2 u}{\partial \mathbf{x}^2} \tag{A6}$$

substituting in right hand side of the dynamic equation

$$-\rho\_s \omega^2 u = \frac{\partial T\_{11}}{\partial x} \tag{A7}$$

one can obtain

$$-\rho\_s \omega^2 \mu = \frac{1}{s\_{11}} \frac{\partial^2 \mu}{\partial \mathbf{x}^2} \tag{A8}$$

The above equation can solve by the free stress boundaries at the ends (*T*11(*x* = ±*L*/2) = 0)

$$\mathcal{S}\_{11}(\mathbf{x} = \pm \mathcal{L}/2) = \frac{\partial u(\mathbf{x} = \pm \mathcal{L}/2)}{\partial \mathbf{x}} = -d\_{31} \frac{\partial \phi}{\partial z} \tag{A9}$$

gives the displacement field as

$$u(\mathbf{x},t) = -\frac{d\_{31}\phi\_m e^{j\omega t}}{h\sqrt{\rho\_s s\_{11}}\omega} \frac{\sin(\sqrt{\rho\_s s\_{11}}\omega \mathbf{x})}{\cos(\sqrt{\rho\_s s\_{11}}\omega L/2)}\tag{A10}$$

Considering the definition of strain from the displacement field

$$S\_{11}(\mathbf{x},t) = -\frac{d\_{31}\phi\_m e^{j\omega t}}{h} \frac{\cos\left(\sqrt{\rho\_s s\_{11}}\omega \mathbf{x}\right)}{\cos\left(\sqrt{\rho\_s s\_{11}}\omega L/2\right)}\tag{A11}$$

Considering the definition of stress from the first piezoelectric equation

$$T\_{11}(\mathbf{x},t) = -\frac{d\_{31}\phi\_m e^{j\omega t}}{h s\_{11}} (\frac{\cos(\sqrt{\rho\_s s\_{11}}\omega \mathbf{x})}{\cos(\sqrt{\rho\_s s\_{11}}\omega L/2)} - 1) \tag{A12}$$

Considering the definition of electric displacement from the second piezoelectric equation

$$D\_3 = d\_{31} T\_{11} - \varepsilon\_{33}^T \frac{\partial \phi}{\partial z} \tag{A13}$$

The analytical solution of electric displacement is obtained

$$D\_3 = -\frac{d\_{31}^2 \phi\_m e^{j\omega t}}{h s\_{11}} \left(\frac{\cos(\sqrt{\rho\_s s\_{11}} \omega x)}{\cos(\sqrt{\rho\_s s\_{11}} \omega L/2)} - 1\right) - \varepsilon\_{33}^T \frac{\phi\_m e^{j\omega t}}{h} \tag{A14}$$

or

$$D\_3 = -\frac{\phi\_m e^{j\omega t}}{h} \left[ \frac{d\_{31}^2}{s\_{11}} \left( \frac{\cos(\sqrt{\rho\_s s\_{11}} \omega \chi)}{\cos(\sqrt{\rho\_s s\_{11}} \omega L/2)} - 1 \right) + \varepsilon\_{33}^T \right] \tag{A15}$$

If piezoelectric material coefficients are real when there will be no losses. Internal losses are attributed to the imaginary part of them, as complex piezoelectric coefficients are:

$$s\_{11} = \bar{s}\_{11} + j\dot{s}\_{11}'\tag{A16}$$

$$d\_{31} = \bar{d}\_{31} + jd\_{31}'\tag{A17}$$

$$
\varepsilon\_{33}^T = \mathfrak{E}\_{33}^T + j\varepsilon\_{33}^{T'} \tag{A18}
$$

where *s*¯11, ¯*d*31, *ε*¯ *T* <sup>33</sup> denote real and *s* 0 <sup>11</sup>, *d* 0 <sup>31</sup>, *ε T*0 <sup>33</sup> denote imaginary parts, respectively. An expression for the power dissipation density in the slab is obtained by considering hysteresis loops. Integrating to determine the areas under the loops and averaging over time yields:

$$Q = \frac{-1}{T} \int\_{0}^{T} D\_{\overline{3}} \frac{\partial E\_{3}}{\partial t} dt + \frac{1}{T} \int\_{0}^{T} T\_{11} \frac{\partial S\_{11}}{\partial t} dt \tag{A19}$$

where *Q<sup>E</sup>* is

$$\begin{split} Q\_{E} &= \frac{-1}{T} \int\_{0}^{T} D\_{3} \frac{\partial E\_{3}}{\partial t} dt \\ = \frac{-1}{T} \int\_{0}^{T} -\frac{\phi\_{m} e^{j\omega t}}{h} \left[ \frac{d\_{31}^{2}}{s\_{11}} \left( \frac{\cos(\sqrt{\rho\_{5} s\_{11}} \omega \mathbf{x})}{\cos(\sqrt{\rho\_{5} s\_{11}} \omega \mathbf{L}/2)} - 1 \right) + \varepsilon\_{33}^{T} \right] \\ &\qquad \frac{\partial}{\partial t} (-\frac{\phi\_{m} e^{j\omega t}}{h}) dt \\ = \frac{-\phi\_{m}^{2}}{Th^{2}} \int\_{0}^{T} j\omega e^{2j\omega t} \left[ \frac{d\_{31}^{2}}{s\_{11}} (\frac{\cos(\sqrt{\rho\_{5} s\_{11}} \omega \mathbf{x})}{\cos(\sqrt{\rho\_{5} s\_{11}} \omega \mathbf{L}/2)} - 1) + \varepsilon\_{33}^{T} \right] dt \end{split} \tag{A20}$$

and *Q<sup>M</sup>* is

$$Q\_{M} = \frac{1}{T} \int\_{0}^{T} T\_{11} \frac{\partial S\_{11}}{\partial t} dt$$

$$= \frac{1}{T} \int\_{0}^{T} -\frac{d\_{31} \phi\_{m} e^{j\omega t}}{l s\_{11}} (\frac{\cos(\sqrt{\rho\_{s} s\_{11}} \omega x)}{\cos(\sqrt{\rho\_{s} s\_{11}} \omega L/2)} - 1)$$

$$\frac{\partial}{\partial t} (-\frac{d\_{31} \phi\_{m} e^{j\omega t}}{h} \frac{\cos(\sqrt{\rho\_{s} s\_{11}} \omega x)}{\cos(\sqrt{\rho\_{s} s\_{11}} \omega L/2)}) dt \tag{A21}$$

$$= \frac{\phi\_{m}^{2}}{T h^{2}} \int\_{0}^{T} e^{2j\omega t} \frac{d\_{31}^{2}}{s\_{11}} (\frac{\cos(\sqrt{\rho\_{s} s\_{11}} \omega x)}{\cos(\sqrt{\rho\_{s} s\_{11}} \omega L/2)} - 1)$$

$$j\omega \frac{\cos(\sqrt{\rho\_{s} s\_{11}} \omega x)}{\cos(\sqrt{\rho\_{s} s\_{11}} \omega L/2)} dt$$

if one replaces the complex piezoelectric coefficients in *Q<sup>E</sup>* the following is obtained

$$\begin{aligned} Q\_E &= \frac{-\phi\_m^2}{Th^2} \int\_0^T j\omega e^{2j\omega t} \left[ \frac{(\bar{d}\_{31} + jd\_{31}')^2}{\bar{s}\_{11} + j\dot{s}\_{11}'} \right] \\ & \left( \frac{\cos(\sqrt{\rho\_s(\bar{s}\_{11} + j\dot{s}\_{11}')} \omega x)}{\cos(\sqrt{\rho\_s(\bar{s}\_{11} + j\dot{s}\_{11}')} \omega L/2)} - 1 \right) + \bar{\varepsilon}\_{33}^T + j\dot{\varepsilon}\_{33}' \frac{T'}{33} \right) dt \\ &= \frac{-\phi\_m^2}{Th^2} \int\_0^T j\omega \left( \cos(2\omega t) + j\sin(2\omega t) \right) \left( \frac{(\bar{d}\_{31} + jd\_{31}')^2}{\bar{s}\_{11} + j\dot{s}\_{11}'} \right) \\ & \left( \frac{\cos(\sqrt{\rho\_s(\bar{s}\_{11} + j\dot{s}\_{11}')} \omega x)}{\cos(\sqrt{\rho\_s(\bar{s}\_{11} + j\dot{s}\_{11}')} \omega L/2)} - 1 \right) + \bar{\varepsilon}\_{33}^T + j\dot{\varepsilon}\_{33}' \frac{T'}{33} \Big] dt \end{aligned} \tag{A22}$$

As

$$\begin{split} \frac{1}{T} \int\_{0}^{T} \cos(2\omega t)dt &= 0\\ \frac{1}{T} \int\_{0}^{T} \sin(2\omega t)dt &= \frac{1}{4\pi} \end{split} \tag{A23}$$

The *Q<sup>E</sup>* is

$$Q\_E = \frac{\phi\_m^2 j\omega}{h^2} \left[ \frac{(\bar{d}\_{31} + jd\_{31}')^2}{\bar{s}\_{11} + j\dot{s}\_{11}'} \left( \frac{\cos(\sqrt{\rho\_s(\bar{s}\_{11} + j\dot{s}\_{11}')} \omega x)}{\cos(\sqrt{\rho\_s(\bar{s}\_{11} + j\dot{s}\_{11}')} \omega L/2)} - 1 \right) \tag{A24}$$
 
$$+ \mathbb{E}\_{33}^T + j\varepsilon \frac{\boldsymbol{r}'}{r\_{33}}]$$

where the real part is

$$\begin{split} \bar{Q}\_{E} &= \frac{\rho\_{m}^{2}\omega}{h^{2}} \left[ \varepsilon^{T'}\_{\;33} + \\ \arg \left( \frac{(d\_{31} + jd\_{31}')^{2}}{\bar{s}\_{11} + js\_{11}'} \frac{\cos(\sqrt{\rho\_{s}(\bar{s}\_{11} + js\_{11}')} \omega \mathbf{x})}{\cos(\sqrt{\rho\_{s}(\bar{s}\_{11} + js\_{11}')} \omega \mathbf{L}/2)} \right) \right] \end{split} \tag{A25}$$

#### **Appendix B. ANSYS APDL Input Files**

! ==== 3D bulk acoustic wave resonator model ======= !Parameters xsubs=2e-4 ysubs=2e-4 zsubs=5e-7 xelectrobot=2e-4 yelectrobot=2e-4 zelectrobot=26e-8 xelectrotop=15e-5 yelectrotop=15e-5 zelectrotop=26e-8 xpiezo=2e-4 ypiezo=2e-4 zpiezo=1e-6 !Element type /PREP7 ET,1,solid226,1001 ET,2,solid95 ET,3,solid95 MP,DENS,1,3.26e3 MP,PERX,1,8 MP,PERY,1,8 MP,PERZ,1,10 TB,PIEZ,1 TBDATA,1,0,0,-0.58 TBDATA,4,0,0,-0.58 TBDATA,7,0,0,1.53 TBDATA,10,0,0,0 TBDATA,13,0,-0.48,0 TBDATA,16,-0.48,0,0 TB,ANEL,1 TBDATA,1,3.45e7,1.25e7,1.2e7 TBDATA,7,3.45e7,1.2e7 TBDATA,12,3.95e7 TBDATA,16,1.1e7 TBDATA,19,1.18e7 TBDATA,21,1.18e7 DMPRAT, 0.0003

mp,damp,1,1e-8 MP,DENS,2,2181.5 MP,EX,2,17652e7 MP,EY,2,17652e7 MP,EZ,2,17652e7 MP,PRXY,2,0.22 MP,PRYZ,2,0.22 MP,PRXZ,2,0.22 MP,GXY,2,72345e6 MP,GYZ,2,72345e6 MP,GXZ,2,72345e6 MP,DENS,3,10096 MP,EX,3,23236e7 MP,EY,3,23236e7 MP,EZ,3,23236e7 MP,PRXY,3,0.362 MP,PRYZ,3,0.362 MP,PRXZ,3,0.362 MP,GXY,3,8530e7 MP,GYZ,3,8530e7 MP,GXZ,3,8530e7 ! Mesh BLOCK,0,xelectrobot/2,0,yelectrobot/2,0,zsubs z0=zsubs+zelectrobot BLOCK,0,xelectrobot/2,0,yelectrobot/2,zsubs,z0 layers = 1 \*DO,i,1,layers BLOCK,0,xpiezo/2,0,ypiezo/2,z0+(i-1)\*zpiezo/layers, z0+(i)\*zpiezo/layers \*ENDDO z1=z0+zpiezo h1=(ypiezo-yelectrotop)/2 h2=(xpiezo-xelectrotop)/2 BLOCK,h2,h2+xelectrotop/2,h1,h1+yelectrotop/2, z1,z1+zelectrotop VGLUE,all ALLSEL VSEL,S,LOC,Z,0,zsubs VATT,2, ,2,0 ALLSEL VSEL,S,LOC,Z,zsubs,z0 VATT,3, ,3,0 ALLSEL VSEL,S,LOC,Z,zsubs+zelectrobot,zsubs+zelectrobot+zpiezo VATT,1, ,1,0 ALLSEL VSEL,S,LOC,Z,z0+zpiezo,z0+zpiezo+zelectrotop VATT,3, ,3,0 size=1e-2 ALLSEL ESIZE,size mshkey,0 MSHAPE,1,3D VMESH,all ! BC NSEL,s,LOC,x,xpiezo/2 NSEL,R,LOC,y,ypiezo/2 D,ALL,Ux,0 D,ALL,Uy,0

NSEL,s,LOC,x,xpiezo/2 D,ALL,Ux,0 NSEL,sR,LOC,y,ypiezo/2 D,ALL,Uy,0 ALLSEL NSEL,R,LOC,x,0, D,ALL,Ux,0 D,ALL,Uy,0 D,ALL,Uz,0 NSEL,s,LOC,y,0, D,ALL,Ux,0 D,ALL,Uy,0 D,ALL,Uz,0 NSEL,S,LOC,z,zsubs+zelectrobot CM,electromasasup,NODE CP,1,VOLT,ALL \*GET,electremasasupnumber,NODE"NUM,MIN NSEL,s,LOC,z,zsubs+zelectrobot+zpiezo NSEL,R,LOC,x,h2,h2+xelectrotop/2 NSEL,R,LOC,y,h1,h1+yelectrotop/2 CM,electrodetop,node CP,2,VOLT,ALL \*GET,electrodetopnumer,NODE"NUM,MIN ALLSEL D,electromasasup,VOLT,0.0 D,electrodetop,VOLT,1.0 !solution /SOLU ANTYPE,HARMONIC HARFRQ,1100E6,1250E6 NSUB,151 OUTRES,ALL,ALL KBC,1 SOLVE ! plot /post26 rforce,2,ELECTRODETOPNUMER,chrg"Ladug PI=3.141592 cfact,1,0,0,2\*PI,0,0 prod,3,2,1"current quot,4,1,1"volt quot,5,4,3"impedance quot,6,3,4"admittance quot,7,4,5"power plvar,freq,power

#### **References**


**Iliana Naumova Apostolova <sup>1</sup> , Angel Todorov Apostolov <sup>2</sup> and Julia Mihailova Wesselinowa 3,\***


**\*** Correspondence: julia@phys.uni-sofia.bg

**Abstract:** The magnetic, optical, and phonon properties of ion-doped CuAlO<sup>2</sup> nanoparticles on the Cu or Al site are theoretically investigated. The room temperature ferromagnetism in CuAlO<sup>2</sup> nanoparticles can be due to the surface, size, and doping effects. The magnetization increases with the decreasing nanoparticle size. The different radii of the transition metal ion and the host Cu ion lead to compressive strain, to the enhancment of the exchange interaction constants, and to increased magnetization *M<sup>s</sup>* and Curie temperature *TC*. By substitution with Mn or Cr on the Al site, tensile strain, a decrease in *Ms*, and an increase in dopants are observed. The size and ion-doping influence on the band-gap energy is also discussed. The phonon energy *ω* decreases, whereas the phonon damping *γ* increases with increasing temperature and decreasing NP size. They show a kink around *T<sup>C</sup>* ∼ 400 K. The behavior of *ω* and *γ* for different ion dopings is observed.

**Keywords:** CuAlO2; ion doping; magnetization; band gap; phonon energy; microscopic model

**Citation:** Apostolova, I.N.; Apostolov, A.T.; Wesselinowa, J.M. Size and Ion-Doping Effects on Magnetic, Optical, and Phonon Properties of CuAlO<sup>2</sup> Nanoparticles. *Magnetochemistry* **2022**, *8*, 169. https://doi.org/10.3390/ magnetochemistry8120169

Academic Editors: C˘at˘alin-Daniel Constantinescu and Lucian Petrescu

Received: 17 October 2022 Accepted: 22 November 2022 Published: 25 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Diluted magnetic semiconductors (DMS) play an important role in interdisciplinary materials science and future spintronics. Ferromagnetic DMS have been studied from first principles within mean field approximation [1,2]. Ferromagnetism in 3d transition metal (TM)-doped II–VI and III–V based DMS is predicted [3]. Kizaki et al. [4] have investigated the magnetism in a new DMS, namely CuAlO<sup>2</sup> (CAO), by the Korringa–Kohn–Rostoker method. Using density functional theory (DFT) calculations, Iordanidou et al. [5] have examined the hole doping effect on the magnetic and electronic properties of CAO. It is known that pure CAO is a p-type width band gap semiconductor with *E<sup>g</sup>* of about 3.5 eV [6]. It has the hexagonal delafossite structure and the space group of R3¯m. CAO can be applied as a DMS when doped with TM ions. Kizaki et al. [4,7] have investigated the electronic and magnetic properties under carrier doping treatment in TM-doped CAO on the Cu site (TM = Fe, Co, Mn, and Ni). The Curie temperature *T<sup>C</sup>* of the doped examples increases with increasing dopant concentration. Moreover, Fe or Co doping on the Al site of CAO bulk and thin films increases the spontaneous magnetization [8–12].

The origin of ferromagnetism in DMS is not yet clear, even if various methods have been proposed. First-principal calculations showed that (Cu,Fe)AlO<sup>2</sup> can be a candidate of ferromagnetic DMS [13]. From the calculated density of states, it seems that the double exchange interaction is the dominant exchange mechanism in (Cu,Fe)AlO<sup>2</sup> [13]. It was recently found that the weak magnetism of these materials may derive from the polarised unpaired electrons around impurities [14]. The magnetic properties of Mn-doped CAO, Cu(Al,Mn)O<sup>2</sup> have been reported by Zhang et al. [15]. The magnetization decreases with increasing Mn concentration. The influence of Cd impurity at Cu and Al sites on the electronic properties of CAO from first-principles calculations is discussed in Ref. [16].

The TM doping on the Al site also influences the optical properties of CAO thin films and nanofibers [17–21]. Raman spectra of pure and ion-doped CAO bulk and NPs are studied in [18,22–26]. So, carrier concentration and ion doping together significantly affect the oxide-based DMS.

The aim of the present paper is to investigate the magnetic, optical, and phonon properties of ion-doped CAO nanoparticles (NPs) for the first time using a microscopic model and Greens's function theory in order to clarify the origin of room temperature ferromagnetism in these systems on microscopic level.

It should be noted that the most theoretical papers consider the magnetic properties of CAO NPs using DFT. The DFT is a very powerful tool in investigation of many body problems. However, DFT is mostly concerned with ground-state properties at zero temperature. In our approach, we are able to cover the whole temperature regime. It is a finite temperature analysis including the entire excitation spectrum. In particular, the method allows us to study the total phase diagram, which is based on the different excitation energies realized in the system. The disadvantage of our approach consists of the consideration of collective properties from the beginning. Our basic quantities are not the "naked" electrons but effective spins of the underlying quasi-particles. Whereas within DFT all of the parameters of the system can be—at least in principle—calculated, we are forced to use additional models to find out those parameters. We are convinced that both approaches, DFT and Green's function method, are appropriate and, to a certain extent, are an alternative in describing many body systems.

#### **2. The Model and Green's Functions**

The room temperature ferromagnetism in bulk CAO can be caused by ion doping and can be described by the Heisenberg hamiltonian *H<sup>d</sup>* :

$$H\_d = \sum\_{i,j} \mathbf{x} I\_{\rm{dij}}(\mathbf{S}\_j \cdot \mathbf{S}\_j) - \sum\_i D\_i(S\_i^z)^2 - \mathbf{g}\mu\_B h \sum\_i S\_i^z. \tag{1}$$

**S***i* is the Heisenberg spin-operator of the TM at site *i*. *J<sup>d</sup>* is the exchange interaction between the TM ions. Kizaki et al. [4] found that *J<sup>d</sup>* between the magnetic ions in the same Cu-plane are ferromagnetic, whereas those between the Cu-planes can be neglected. *D* is the single-ion anisotropy, *h* is an external magnetic field, and *x* is the doping concentration.

In the case of a CAO NP, there also appears a ferromagnetism from surface and size effects due to the uncompensated Al spins on the surface described by *H<sup>s</sup>* as well as due to ion-doping effects (see *H<sup>d</sup>* ):

$$H\_{\mathbf{s}} = -\sum\_{i,j} (1-\mathbf{x}) f\_{sij}(\mathbf{S}\_j \cdot \mathbf{S}\_j). \tag{2}$$

From the spin Green's function defined as *Gij* = *S* + *i* ; *S* − *<sup>j</sup>* the magnetization *M* for arbitrary spin *S* is calculated:

$$M = \frac{1}{N^2} \sum\_{k} \left[ (S + 0.5) \coth[(S + 0.5)\beta E\_m(\mathbf{k})] - 0.5 \coth(0.5 \beta E\_m(\mathbf{k})) \right],\tag{3}$$

where *β* = 1/*kBT*, *Em*(**k**) is the spin-wave energy.

The spin–phonon and phonon–phonon interactions are described by:

$$H\_{sp-ph} = \frac{1}{2} \sum\_{i,j,k} F(i,j,k) Q\_i S\_j^z S\_k^z - \frac{1}{4} \sum\_{i,j,r,s} R(i,j,r,s) Q\_i Q\_j S\_r^z S\_s^z + h.c.\tag{4}$$

The normal coordinate *Q<sup>i</sup>* can be expressed in terms of phonon creation *a* <sup>+</sup> and annihilation *a* operators, *Q<sup>i</sup>* = (2*ω*0*i*) <sup>−</sup>1/2(*a<sup>i</sup>* + *a* + *i* ). *F* and *R* are the spin–phonon coupling constants in the first and second order, respectively.

$$H\_{ph-ph} = -\frac{1}{2!} \sum\_{i} \omega\_{0i} a\_i a\_i^+ + \frac{1}{3!} \sum\_{i,j,r} B(i,j,r) Q\_i Q\_j Q\_r$$

$$+ \quad \frac{1}{4!} \sum\_{i,j,r,s} A(i,j,r,s) Q\_i Q\_j Q\_r Q\_{s\cdot r} \tag{5}$$

where *ω*0*<sup>i</sup>* is the frequency of the lattice mode.

From the poles of the phonon Green's function *G*˜ *ij*(*t*) = hh*ai*(*t*); *a* + *j* ii the phonon energies are observed:

$$\begin{aligned} \omega\_{\dot{ij}}^2 &= -\omega\_0^2 - 2\omega\_0 \Big( M\_i M\_{\dot{j}} R\_{\dot{ij}} \delta\_{\dot{ij}} \\ &- \frac{1}{2N'} \sum\_r A\_{\dot{ij}r} (2\bar{N}\_r + 1) - B\_{\dot{ij}} \langle Q\_{\dot{ij}} \rangle \delta\_{\dot{ij}} \Big), \end{aligned} \tag{6}$$

$$\langle Q\_{\rm ij} \rangle = \frac{M\_i M\_{\bar{j}} F\_{\bar{i}\bar{j}} \delta\_{\bar{i}\bar{j}} - \frac{1}{N'} \sum\_{r} B\_{\bar{i}jr} (2\bar{N}\_r + 1)}{\omega\_0 - M\_i M\_{\bar{j}} R\_{\bar{i}\bar{j}} \delta\_{\bar{i}\bar{j}} + \frac{1}{N'} \sum\_{r} A\_{\bar{i}jr} (2\bar{N}\_r + 1)}. \tag{7}$$

The phonon correlation function *N*¯ *<sup>r</sup>* = h*a* + *r ar*i is obtained via the spectral theorem. The phonon damping *<sup>γ</sup>* = *<sup>γ</sup>sp*−*ph* + *<sup>γ</sup>ph*−*ph* is also calculated taking into account anharmonic spin–phonon and phonon–phonon interactions.

The band-gap energy *E<sup>g</sup>* of a CAO NP is defined by the difference between the valence and conduction bands:

$$E\_{\mathcal{S}} = \omega^{+}(\mathbf{k} = 0) - \omega^{-}(\mathbf{k} = \mathbf{k}\_{\sigma}).\tag{8}$$

The electronic energies

$$
\omega^{\pm}(k) = \epsilon\_k - \frac{\sigma}{2}I \langle \mathbb{S}^z \rangle \tag{9}
$$

are observed from the Green's function *g*(*k*, *σ*) = *ck*,*σ*; *c* + *<sup>k</sup><sup>σ</sup>* , *σ* = ±1, *c* + *iσ* and *ci<sup>σ</sup>* are Fermi operators, and *I*-s-d is the interaction constant [27].

#### **3. Numerical Results and Discussion**

A CAO NP with a cubo-octahedral shape is defined by fixing the origin at a certain Al spin in the center of the particle and including all other spins within the particle into shells, which are numbered by *n* = 1, . . . , *N* from the central to the surface shell. For the numerical calculations, we use the following model parameters: *J<sup>d</sup>* (Fe-Fe) = 195 meV, *Jd* (Co-Co) = 117 meV, *J<sup>d</sup>* (Mn-Mn) = 78 meV, *J<sup>d</sup>* (Ni-Ni) = 97.5 meV [4], *J<sup>s</sup>* = 20 meV, *<sup>D</sup>* <sup>=</sup> <sup>−</sup>0.068 meV, *<sup>F</sup>* = 23 cm−<sup>1</sup> , *<sup>R</sup>* <sup>=</sup> <sup>−</sup>18 cm−<sup>1</sup> , *A* = 6.61 cm−<sup>1</sup> , and *<sup>B</sup>* <sup>=</sup> <sup>−</sup>2.94 cm−<sup>1</sup> . The exchange interaction *Jij* = *J*(*r<sup>i</sup>* − *rj*) depends on the distance between the spins, i.e., on the lattice parameters. It is inverse proportional to the lattice parameters. Therefore, the exchange interaction constant on the surface *J<sup>s</sup>* and of the doped states *J<sup>d</sup>* can be changed or can increase or decrease in dependence on the strain in the lattice. So, we take into account the direct connection between the microstructure and magnetic behavior.

Firstly, we will consider the spontaneous magnetization *M<sup>s</sup>* in a pure CAO NP. It must be noted that bulk CAO is a p-type wide band gap semiconductor, and it is a nonmagnetic compound. However, due to the uncompensated surface spins of the Al3<sup>+</sup> ions on the surface, a spontaneous magnetization *M<sup>s</sup>* appears in the CAO NP, contrary to the bulk case, where *M<sup>s</sup>* = 0. We obtain that *J<sup>s</sup>* increases and the spontaneous magnetization *M<sup>s</sup>* of a CAO NP increases with decreasing NP size. The result is presented in Figure 1. A finite magnetization *M<sup>s</sup>* is also observed experimentally in pure CAO NPs [12,28] and in pure CAO thin films [29], due to the existence of point defects. It must be noted that we obtain

similar behavior for the size dependence of the phase-transition temperature *TC*(*N*), *TC*, which increases with decreasing NP size.

**Figure 1.** Size dependence of the spontaneous magnetization of a CAO NP for *T* = 300 K.

The origin of ferromagnetism in DMS is not yet clear. In order to investigate the observed room temperature ferromagnetism (RTFM) in ion-doped CAO-bulk and NPs, we have replaced some Cu-ions with Fe-ions, i.e., we consider Cu1−*x*Fe*x*AlO<sup>2</sup> NPs for *x* = 0–0.25. There are experimental data that show that in this concentration interval the samples are single-phase, whereas for larger *x* values, the second phase appears [4]. It must be noted that by the doping of TM on the Al site, the limit value of the ion doping concentration is smaller; it is around 0.05–0.1 [8,10,15,28,30], which is different for the different doping ions. The exchange interaction constant *J* depends on the distance between the spins, i.e., on the lattice parameter, on the different strain, on the lattice symmetry, and on the number of next neighbors. *J* is an inverse square function of the distance between two neighboring spins, i.e., *J* ∝ 1/(*r<sup>i</sup>* − *rj*) 2 . From the structural analysis, Cu1<sup>+</sup> is replaced by Fe3<sup>+</sup> or Fe2+, which leads to cation vacancies increasing because of charge variation [4,7,9]. Through the different radii of the dopants, lattice defects can be introduced or intrinsic host-lattice defects can be activated when different ions, such as Fe, occupied the Cu sites. Chen et al. [9] have shown that Fe element exists in the ion form, with the valence of +2 and/or +3. Since the ionic radius of Fe2<sup>+</sup> (0.92 *A*˙) or Fe3<sup>+</sup> (0.65 *A*˙) is always smaller than Cu1<sup>+</sup> (0.95 *A*˙), the lattice parameters decrease with increasing the Fe dopants as supported by Chen et al. [9]. The reduction of the lattice constants in Fe-doped CAO nanostructures leads to a compressive strain, i.e., in our microscopic model, it leads to an increase in the exchange interaction constant of the doped states *J<sup>d</sup>* with an increase in the Fe ion concentration because it is inverse proportional to the lattice parameters . So, we obtain an increase in the spontaneous magnetization *M<sup>s</sup>* and the Curie temperature *TC*. The ferromagnetic coupling between the doping Fe-ions contributes additively to this increase. Fe3+-ion (Fe2+-ion) has 5*d* (6*d*) electrons with the total spin of *S* = 5/2 (*S* = 2), whereas the Cu1+-ion has *S* = 0. Thus, Fe3<sup>+</sup> or Fe2<sup>+</sup> substitution into Cu1<sup>+</sup> induces an extra magnetic moment. The spontaneous magnetization *M<sup>s</sup>* and the Curie temperature *T<sup>C</sup>* in dependence on the Fe3<sup>+</sup> doping concentration are calculated. We observe that both increase with increasing Fe doping concentration *x*. The result for *TC*(*x*) is demonstrated in Figure 2, curve 1. We would obtain a similar increase in the phase-transition temperature *TC*, for example, for Ni-, Mn-, or Co-doped CAO, (Cu,TM)AlO<sup>2</sup> (Figure 2, curves 2–4), where the doping ions have a smaller ionic radius in comparison with that of the Cu ion. This behavior of the doping dependence of *TC*(*x*) is in good qualitative agreement with the experimental data [8–12].

It must be noted that for example, in Mn- or Cr-doped CAO on the Al site, where the ionic radius of Mn3+(0.72 *A*˙) or Cr3+(0.63 *A*˙) is larger than that of Al3+(0.51 *A*˙), the lattice parameters increase in agreement with the experimental data of Zhang et al. [15]. A tensile strain appears, i.e., *J<sup>d</sup>* decreases, which leads to a decrease in *M<sup>s</sup>* with increasing Mn or Cr ion doping concentration (see Figures 3 and 4, respectively). Similar behavior is reported for Cu(Al,Mn)O<sup>2</sup> by Zhang et al. [15]. Unfortunately, there are not experimental data for *M*(*x*) of Cr-doped CAO. A decrease in the magnetization is observed additively with increasing Mn content *x* due to the antiferromagnetic coupling among the Mn spins (contrary to the case of Mn doping at the Cu site). By the substitution of the Al ion with a rare earth ion (RE) whose radius is larger than that of Al, tensile strain and a reduction of the spontaneous magnetization *M<sup>s</sup>* with an increase in the RE ion doping concentration appears again. Unfortunately, there are not experimental data for RE doped CAO.

**Figure 2.** (Color online) The Curie temperature *T<sup>C</sup>* as a function of the ion doping concentration *x* for different ions substituted the Cu ion: (1) Fe; (2) Co; (3) Mn; and (4) Ni.

**Figure 3.** (Color online) Magnetic field dependence of the magnetization *M<sup>s</sup>* of Cu(Al,Mn)O<sup>2</sup> for different Mn-ion doping concentrations *x*: (1) 0.01; (2) 0.03; and (3) 0.05.

**Figure 4.** Dependence of the magnetization *M<sup>s</sup>* of Cu(Al,Cr)O<sup>2</sup> on the doping Cr ion concentration.

From Equation (8), we have calculated the band-gap energy *E<sup>g</sup>* for pure and iondoped CAO NPs (see Figure 5). It must be noted that there is a broadening of the band gap *E<sup>g</sup>* of pure CAO NPs with a decrease in the NP size. We now consider the case of a Cr3+-doped CAO NP, CuAl1−*x*Cr*x*O2. The lattice parameters increase with the increase in the Cr-doping concentration [17] because the ionic radius of the doping Cr ion (0.63 *A*˙) is larger than that of the host Al ion (0.51 *A*˙). A tensile strain appears, i.e., *J<sup>d</sup>* decreases with increasing *x*, which leads to a decrease in the magnetization (see Figure 4) and an increase in *E<sup>g</sup>* (curve 1). The observed behavior is in agreement with the experimental data of Jiang et al. [17]. We would also obtain a similar enhanced *E<sup>g</sup>* within our model by doping with Sb3<sup>+</sup> (0.76 *A*˙) ions or Y3<sup>+</sup> (1.04 *A*˙) ions, which also causes a tensile strain, i.e., an increase in the lattice parameters as well as a decrease in the exchange interaction constant *J<sup>d</sup>* and therefore reduced magnetization *M<sup>s</sup>* and enhanced *Eg*, reported in [18,31]. Otherwise, by doping with the Mg ion, we observe the contrary result: a decrease in the band-gap energy (curve 2), in agreement with [28,30,32]. The optical absorption of the CAO increases due to the reduction in *E<sup>g</sup>* by Mg doping. A smaller *E<sup>g</sup>* is also reported by Bi and Eu doping at the Al site in CAO [21,33].

Finally, we will discuss the temperature, size, and ion-doping dependence of the phonon energy *ω* and damping *γ* in CAO NPs. The temperature dependence for the A1*<sup>g</sup>* mode *ω*<sup>0</sup> = 767 cm−<sup>1</sup> of the phonon energy *ω* and phonon damping *γ* is shown in Figure 6. We take into account the anharmonic spin–phonon interaction *R* < 0, as well as the three-phonon *A* > 0 and four-phonon *B* < 0 anharmonic interactions. Let us emphasize that in order to obtain the experimentally observed softening of *ω*, we have to chose *R* < 0. For *R* > 0, we would obtain a hardening, an increase in *ω* with *T* [34]. The damping *γ* is proportional to *R* <sup>2</sup> and so independent of the sign of *R*. *γ* corresponds to the full width at half-maximum of the Raman lines and is indirectly proportional to the phonon life time. The damping is a direct consequence of the anharmonic coupling. The phonon energy decreases, whereas the phonon damping increases with increasing temperature, which is in agreement with [22,26]. It can be seen from Figure 6 that around *T<sup>c</sup>* ∼ 400 K, we observe a kink due to the spin–phonon interaction *R*. For *R* = 0, the kink disappears. Moreover, the kink is size-dependent; it shifts to smaller *T* values with increasing NP size. In Figure 1, we have shown that with decreasing NP size spontaneous magnetization *M<sup>s</sup>* with the corresponding *T<sup>C</sup>* value appears. A similar anomaly in the phonon energy and damping of CAO thin films is observed experimentally between 423 and 573 K by Singh et al. [22] and around *T* ∼ 300 K by Pellicer-Porres et al. [26]. The authors have proposed that these changes may be related to a negative thermal expansion (NTE) in the delafossite structure

along the O-Cu-O linkage. A similar NTE is also reported in other AMO<sup>2</sup> (A = Cu or Ag; M = Al, Sc, In, or La) [35–37]. NTE means a volume contraction, i.e., there is a change in the structure, for example, by transitions between ferro- and paramagnetic phases.

**Figure 6.** (Color online) Temperature dependence of the phonon energy *ω* (1) and phonon damping *γ* (2) of a CAO NP, *N* = 10.

The size dependence of the phonon energy *ω* and the phonon damping *γ* in a CAO NP is presented in Figure 7. It can be seen that the phonon energy decreases whereas the damping increases with decreasing nanoparticle size, which is in agreement with the experimental data of Yassin et al. [23].

Number of shells N

**Figure 7.** (Color online) Size dependence of the phonon energy *ω* and phonon damping *γ* of a CAO NP, *T* = 300 K.

The ion-doping also influences the phonon properties of CAO NPs. It may cause lattice distortions and unrelaxed strains. To show this, we will consider the Fe3<sup>+</sup> doping at the Al3<sup>+</sup> site. The ionic radius of the Fe ion (0.65 *A*˙) is larger than that of the Al ion (0.51 *A*˙), i.e., we observe tensile strain, an increase in the lattice parameters that leads to a decrease in *J<sup>d</sup>* and *R<sup>d</sup>* and to a decrease in the phonon energy with increasing Fe ions, see Figure 8, curve 1. This corresponds to a red-shift in the peak positions. The contribution to the damping of the ion doping is in addition to those of the bulk, surface, spin–phonon, and phonon–phonon interactions, so that *γ* increases with increasing Fe concentration (see Figure 8, curve 2). This means that the full-weight of the half maximum increases, i.e., the Raman lines are broadener. This is in agreement with the experimental results of Fe-doped CAO [38]. We obtain a similar decrease in *ω* and an increase in *γ* for Eu-doped CAO NP as reported by Liu et al. [21].

**Figure 8.** (Color online) Fe ion-doping dependence of the phonon energy *ω* and phonon damping *γ* of a Cu(Fe,Al)O<sup>2</sup> NP, *N* = 10, and *T* = 300 K.

#### **4. Conclusions**

To conclude, we have studied the size and ion-doping effects on the magnetic, optical, and phonon properties of transition metal (TM = Fe, Co, Mn, and Ni) ion-doped CAO NPs on the Cu and Al site using a microscopic model and Green's function theory. The RTFM in bulk CAO can be due to ion doping, whereas CAO NPs make an additive contribution in terms of surface and size effects . Due to surface defects, a magnetization is obtained that increases with decreasing NP size. The difference in the radii of the TM3<sup>+</sup> and the host Cu<sup>+</sup> ions leads to a compressive strain, an increase in the exchange interaction constants, and increased spontaneous magnetization *M<sup>s</sup>* and Curie temperature *TC*. The doping TM3<sup>+</sup> ion has an additional spin value compared to the zero spin value of the Cu<sup>+</sup> ion, which also enhances *M<sup>s</sup>* . On the other hand, a tensile strain is observed by substitution with Mn3<sup>+</sup> or Cr3<sup>+</sup> on the Al3<sup>+</sup> site, where the radius of the Mn or Cr ions is larger than that of the Al ion, which leads to a decrease in *M<sup>s</sup>* with an increase in the ion-doping concentration. The band-gap energy increases with decreasing NP size and by doping with Cr or Sb ions and decreases by Mg ion doping. The phonon energy *ω* decreases, whereas the phonon damping *γ* increases, with increasing temperature. They show a kink around the phasetransition temperature *T<sup>C</sup>* ∼ 400 K, caused by the anharmonic spin–phonon interaction. Due to surface and size effects, *ω* decreases, whereas *γ* increases, with decreasing NP size. The phonon energy decreases, whereas the damping increases, by Fe and Eu doping on the Al site, i.e., a red shift and broadening of the Raman peaks is obtained. The observed results are in good qualitative agreement with the existing experimental data.

Let us emphasize that in order to take into account the carrier doping treatment in transition-metal-doped CAO on the Cu site, which plays an additive role in the effect of the different ionic radii, we have to expand our model and consider the s-d model by also taking into account the electron–phonon interaction, which is important for CAO [39]. This could be done in the next paper.

**Author Contributions:** All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The raw data that support the findings of this study are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Mesostructure and Magnetic Properties of SiO2-Co Granular Film on Silicon Substrate**

**Natalia A. Grigoryeva1,**<sup>∗</sup> **, Victor Ukleev <sup>2</sup> , Alexey A. Vorobiev <sup>3</sup> , Alexander I. Stognij <sup>4</sup> , Nikolay N. Novitskii <sup>4</sup> , Leonid V. Lutsev <sup>5</sup> and Sergey V. Grigoriev <sup>6</sup>**


**Abstract:** Granular films SiO2(Co) exhibit unusual magnetic and magnetotransport properties which are strongly dependent on the composition of the film and material of a substrate. For example, the injection magnetoresistance (IMR) coefficient reaches a giant (GIMR) value of 105% at room temperature in SiO2(Co) films on an n-GaAs substrate. However, the IMR effect is negligible in the case of a similar granular film deposited on the n-Si substrate. In this report, the structural and magnetic properties of granular film SiO2(Co) on Si substrate are studied with the aim to understand the cause of the difference in IMR coefficients for SiO2(Co) thin film deposited on n-GaAs and on n-Si substrates. Investigations were carried out using complementary methods of Polarized Neutron Reflectometry, Grazing Incidence Small-Angle X-ray Scattering, X-ray Reflectometry, Scanning Electron Microscope, and SQUID magnetometry. It is shown that the interface layer between the granular film and Si substrate exhibits metallic rather than magnetic properties and eliminates the GIMR effect. This interface layer is associated with the Si diffusion to Co nanoparticles and the formation of the metallic cobalt silicides.

**Keywords:** granular film; nanoparticles; grazing-incidence small-angle X-ray scattering; polarized neutron reflectometry

**PACS:** 75.70.Cn; 75.75.–c; 75.50.Lk

#### **1. Introduction**

Granular films (GFs) are nanocomposites that consist of metallic magnetic nanoparticles incorporated in an insulating matrix. GFs with a very well-controlled nanostructure have promising technological applications in microwave devices, spintronics, medicine and biology [1–7]. Structural investigations of GFs play a crucial role and often facilitates the interpretation of their electrical and magnetic properties as well as their interconnection. A complex interplay of intrinsic properties can be caught through their structural features also accounting for finite-size effects, size distribution, surface effects, and nanoparticle interactions.

The effect of giant injection magnetoresistance (GIMR) has been observed in the GFs SiO2(x at.% Co) deposited on GaAs substrate [8,9]. The GIMR phenomenon is a magneticfield-induced suppression of the spin-polarized electron current from the GF into the semiconductor substrate. The IMR coefficient has maximum values for structures with Co concentrations in the range 54 ÷ 75 at.%. In the case of the opposite current direction (electrons drift from the semiconductor into the granular film) the magnetoresistance effect

**Citation:** Grigoryeva, N.A.; Ukleev, V.; Vorobiev, A.A.; Stognij, A.I.; Novitskii, N.N.; Lutsev, L.V.; Grigoriev, S.V. Mesostructure and Magnetic Properties of SiO2-Co Granular Film on Silicon Substrate. *Magnetochemistry* **2022**, *8*, 167. https://doi.org/10.3390/ magnetochemistry8120167

Academic Editors: C˘at˘alin-Daniel Constantinescu and Lucian Petrescu

Received: 25 September 2022 Accepted: 15 November 2022 Published: 24 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

becomes less pronounced. The GIMR coefficient reaches 10<sup>5</sup> % in a narrow temperature range 240 K < T < 300 K at the magnetic field of H = 10 kOe applied in the sample plane and at the voltage U = 70 V applied perpendicular to the sample plane.

The extreme value of this effect at room temperature and the possibility of controlling the resistivity of heterostructure by varying the magnetic field could be used in the development of spintronic devices such as spin injectors, spin diodes, and transistors. Unlike crystalline silicon, the gallium arsenide substrate has a number of disadvantages in practical use: it is mechanically fragile and thermally unstable above 600 ◦C due to As evaporation, which imposes restrictions on the synthesis technology. The use of silicon substrates would eliminate these problems and provide additional advantages, such as its intrinsic carrier concentration being n*<sup>i</sup>* = 1.02 <sup>×</sup> <sup>10</sup><sup>10</sup> cm−<sup>3</sup> and the energy gap being E*<sup>g</sup>* = 1.12 eV (compared to GaAs <sup>−</sup> <sup>n</sup>*<sup>i</sup>* = 2.1 <sup>×</sup> <sup>10</sup><sup>6</sup> cm−<sup>3</sup> and E*<sup>g</sup>* = 1.43 eV, respectively). However, while the GIMR effect reaches 10<sup>5</sup> % in the case of the GF deposited on GaAs substrate, the GIMR is negligible in the case of a similar granular film deposited on n-Si substrate [8], and the intrinsic magnetoresistance of SiO2(x at.% Co)/Si films has even negative values. In the present work we studied Au/SiO2(Co x at.%)/Si heterostructures with cobalt content of x = 54, 60, 70, 75, 82 at.%. The aim of the study is to investigate the structural features and the associated magnetic properties of the heterostructures Au/GF/Si and to reveal the role of Si substrate in the negligible GIMR effect. For this purpose, we used the Polarized Neutron Reflectometry (PNR) technique, which we combined with complementary techniques of Grazing Incidence Small Angle X-ray Scattering (GISAXS), X-ray Reflectometry (XRR), Scanning Electron Microscope (SEM), and Superconducting Quantum Interference Device Magnetometry (SQUID).

The paper is organized in the following way. Section 2 gives details of the sample preparation and experimental results obtained via different methods to reveal structural and magnetic properties. Section 3 provides a discussion and interpretation of the results combined. A combination of different methods makes it possible to build a model explaining the negligible value of the GIMR effect, providing that almost all structural and magnetic properties of GF are practically the same for the films disposed on GaAs and Si substrates. It is shown that the interface layer between the GF and Si substrate shows metallic rather than magnetic properties and eliminates the GIMR effect. This interface layer is associated with the Si diffusion to Co nanoparticles and the formation of metallic cobalt silicides. Section 4 is the conclusion.

#### **2. Experiment**

#### *2.1. Sample Preparation*

The granular films SiO2(x at.% Co) were synthesized by ion beam co-sputtering [10,11] of the composite cobalt-quartz target on the commercial n-Si substrate heated to 200 ◦C and misoriented by 2 degrees from (100) plane to (110) plane and with the thickness of 0.4 mm. The resistivity of the Si substrate was measured by the dc four-probe method at room temperature and was equal to 3.7 Ω·cm. The concentration of cobalt nanoparticles in SiO<sup>2</sup> was set by the ratio of cobalt and quartz target areas. The film was capped by 2 nm of the gold layer as a conductive contact for magnetoresistance measurements. The sputtering method and choice of components are described in ref. [12]. A detailed description of the characteristics of granular films is given in ref. [8].

The granular films SiO2(x at.% Co) can be differentiated by their conductivity into three classes depending on the Co concentration.


particles. The concentration of Co x ≈ 40 at.% is a percolation threshold, below which only the tunneling process contributes to the conductivity of the granular metal.

(3) Class of GFs (SiO2(x at.% Co), at x > 85 at.%) can be attributed to metals as the volume fraction of the metal x is large, metal grains merge and form a metallic continuum with dielectric inclusions.

For this study, the granular films SiO<sup>2</sup> (x at.% Co) with Co concentrations of x = 54, 60, 70, 75, 82 at.% deposited on Si substrate were chosen. A total of 5 batches of SiO<sup>2</sup> (x at.% Co)/Si heterostructures with the same concentrations were synthesized and studied to ensure the repeatability of synthesis results. Figure 1 demonstrates SEM images for SiO<sup>2</sup> (54 at.% Co)/Si (a), and for SiO<sup>2</sup> (82 at.% Co)/Si (b). Field Emission SEM (FE-SEM) Zeiss Leo 1530 capable of 2.5 nanometer resolution equipped with an In-Lens secondary electron detector has been used. The average size of cobalt nanoparticles has been obtained from a set of linear measurements using ImageJ software and equals about 75 ± 15 Å. The microscope achieves a resolution of approx. 1 nm.

**Figure 1.** SiO2(54 at.% Co)/Si (**a**) and SiO2(82 at.% Co)/Si (**b**) images in the secondary electron mode with the InLens detector, mounted above the films. Light contrast corresponds to the SiO<sup>2</sup> matrix and dark areas correspond to Co inclusions.

#### *2.2. Grazing Incidence Small Angle X-ray Scattering and X-ray Reflectometry*

The Grazing Incidence Small Angle X-ray Scattering (GISAXS) and X-ray Reflectometry (XRR) experiments were carried out at ID10 beamline of European Synchrotron Radiation Facility (ESRF, Grenoble, France) [13]. The geometry of the GISAXS (Figure 2) is described by the grazing incidence angle *α<sup>i</sup>* and two scattering angles *α<sup>f</sup>* (out of sample plane) and *ϕ* (in-plane). Information on the electron density distribution in the z direction is obtained by measuring the scattering intensity as a function of the momentum transfer component Q*<sup>z</sup>* perpendicular to the sample plane. Meanwhile, the lateral structure of the granular film, its surface and its interface with the silicon substrate can be studied by measuring intensity as a function of the component Q*<sup>y</sup>* lying in the sample plane. A detailed description of the theory and methodology of the GISAXS technique can be found in the review [14–16]. The GISAXS experiment with the GF SiO2(x at.% Co) on GaAs substrate can be found in [17]. The collimated X-ray beam of the size 1 <sup>×</sup> 0.1 mm<sup>2</sup> and wavelength *<sup>λ</sup>* = 0.95 Å, impinging the sample surface at the grazing angle *α<sup>i</sup>* = 0.32◦ , was used for the GISAXS experiment. The scattering patterns were recorded with the two-dimensional position-sensitive detector MARCCD-133. The central part of the detector was protected from high-intensity direct and specular beams by a lead absorber.

**Figure 2.** Geometry of the GISAXS experiment: **k<sup>i</sup>** and **k<sup>f</sup>** are the wave vectors of the incident and scattered beams, respectively; *α<sup>i</sup>* , *α<sup>f</sup>* and *ϕ* angles determine the momentum transfer components Q*x*, Q*y*, Q*z*.

The two-dimensional GISAXS intensity maps from the samples Au/GF/Si are shown in Figure 3a–c. Several features can be distinguished. Firstly, the presence of a broad elliptic diffuse halo whose intensity is increased at *α<sup>f</sup>* = *α<sup>c</sup>* = 0.25◦ (where *α<sup>c</sup>* is the critical angle of the Total External Reflection) due to the Yoneda anomalous scattering and has zero intensity at *α<sup>f</sup>* = 0, i.e., below the sample horizon. This broad elliptic diffuse halo is a typical manifestation of scattering from an isotropic three-dimensional system of scatterers (cobalt nanoparticles) located at a distance *l*<sup>1</sup> = 2*π*/*Q*<sup>1</sup> from each other, where Q<sup>1</sup> expressed in momentum transfer units is the radius of the diffuse halo.

The second feature of the scattering pattern is the presence of two symmetrical peaks at angles *ϕ* = ±0.29◦ . Such peaks correspond to the scattering from the nanoparticles spatially arranged within the monolayer at the interface with the substrate [14–16,18]. By taking the peaks coordinate Q<sup>2</sup> it is easy to calculate the average interparticle distance at the interface: *l*<sup>2</sup> = 2*π*/Q2. To determine both *l*<sup>1</sup> and *l*<sup>2</sup> distances the cuts along *ϕ* direction of the two-dimensional intensity maps at *α<sup>f</sup>* = 0.25◦ were taken (Figure 3d).

**Figure 3.** (**a**) The two-dimensional GISAXS intensity maps from the samples SiO2(38 at.% Co)/Si (**a**), SiO2(54 at.% Co)/Si (**b**), and SiO2(82 at.% Co)/Si (**c**). The dashed horizontal lines indicate the integration regions for the one-dimensional profiles shown on (**d**) panel. (**d**) The scattering intensity profiles along *φ* direction at *α<sup>f</sup>* = 0.25◦ for SiO2(38 at.% Co)/Si—curve 1, SiO2(54 at.% Co)/Si—curve 2, and SiO2(82 at.% Co)/Si—curve 3. The gray area corresponds to the beamstop shadow.

In XRR experiments, the intensity of X-rays scattered in the direction *α<sup>i</sup>* = *α<sup>f</sup>* and *ϕ* = 0 was measured as a function of the incidence angle *α<sup>i</sup>* . The measurements were performed with a 1 <sup>×</sup> 0.1 mm<sup>2</sup> beam with a wavelength of 1.56 Å. To detect the reflected beam, we used an mBraun-50M linear position-sensitive detector, which made it possible not only to measure the desired signal but also to simultaneously determine the background, which was then subtracted. The specular reflection curve R(Q*z*) thus obtained is related to the one-dimensional in-depth distribution of the electron density of the sample *ρ*(*z*) using the Parratt method [19]. The results of the XRR experiment are presented in Figure 4, where the intensity of the reflected beam is shown as a function of the momentum transfer for Au/SiO2(60 at.% Co)/Si and Au/SiO2(82 at.% Co)/Si.

**Figure 4.** Synchrotron X-ray reflectometry from Au/SiO<sup>2</sup> (60 at.% Co)/Si (**a**) and Au/SiO<sup>2</sup> (82 at.% Co)/Si (**b**) heterostructures. The gray points are experimental data, the solid line is the theoretical curve corresponding to the model of the electron-density distribution shown in the inset with the indication of electron densities plotted for the pure components of the sample.

The experimental reflectivity data shows clear maxima with oscillating intensity (Kiessig fringes) caused by the finite thickness of the film. The experimental curves (Figure 4) are sufficiently well reproduced by the model containing two layers on a substrate Si. The parameters of the model include the depth profile of the electron scattering length density (SLDX-rays), layer thickness (d) and interface (or surface) roughness (*σ*). The first layer corresponds to the technical layer of gold, the second layer corresponds to the granular SiO2(x at.%Co) film with a thickness of 450 Å. The fitting parameters for each layer adjusted to minimize the value of reduced *χ* <sup>2</sup>—a weighted measure of goodness of fit are presented in Table 1. The absence of fast oscillations corresponding to the total thickness of the bilayer Au/GF cannot be explained by the insufficient resolution of the instrument. Measurements were taken at a step of 0.02 degrees (50 steps/1 degree). This is enough to resolve a film with a thickness of about 450 Å (the oscillation period from 450 Å is 0.1 degrees). The reduction in the oscillation amplitude can be attributed to the interface roughness or, rather, to smeared SLDX-rays contrast (a contrast of the electron density) between a granular film and a substrate.


**Table 1.** Parameters of the model used to approximate the experimental X-ray reflectometry curves.

#### *2.3. Polarized Neutron Reflectometry*

A detailed description of the polarized neutron reflectometry technique can be found elsewhere [20]. The experiment was carried out at the SuperADAM reflectometer at Institut Laue-Langevin (Grenoble, France) [21–23]. Typical geometry of the PNR experiment is shown in Figure 5.

**Figure 5.** Typical geometry of the polarized neutron reflectometry experiment.

We exploited a neutron beam with the wavelength *λ* = 4.41 Å and the initial polarization P<sup>0</sup> = 0.98. The lateral components of neutron momentum transfer Q*<sup>x</sup>* and Q*<sup>y</sup>* were equal to zero, and the vertical component Q*<sup>z</sup>* was equal to 4*π*sin(*α<sup>f</sup>* )/*λ*, where *α<sup>f</sup>* is the angle of the scattered beam with wave vector **k<sup>f</sup>** . Intensities of the reflected beams with initial polarization +**P<sup>0</sup>** and −**P0**, i.e., along and against the direction of the external magnetic field **B**, were measured one by one. The magnetic field from 0 to 320 mT was applied parallel to the film surface and perpendicular to the incident neutron beam with wave vector **k<sup>i</sup>** . The direction of **P<sup>0</sup>** with respect to **B** was switched by a spin Flipper. Specularly reflected neutrons were detected by the <sup>3</sup>He position-sensitive detector. Measurements were performed at temperatures, where the GIMR effect is not detected (T = 120 K), reaches its maximum (T = 300 K) and again is not detected (T = 420 K) [8]. The data for every temperature were taken after the demagnetization process. No difference between reflectivity curves R(Q*z*, +**P0**) and R(Q*z*, −**P0**) was observed at B = 0 mT, while the application of the external field resulted in a pronounced splitting of two components (Figure 6).

**Figure 6.** The experimental (symbols) and fitted model (lines) reflectivity curves for the samples Au/SiO2(60 at.% Co)/Si (**a**–**c**) and Au/SiO2(70 at.% Co)/Si (**d**–**f**) at fields B = 0 and B = 240 mT for T = 120 K (**a**,**d**), 300 K (**b**,**e**), and 420 K (**c**,**f**). Measured and fitted model curves at B = 240 mT are multiplied by 10 for convenience.

For fitting the Polarized Neutron Reflectometry experimental data the Parratt algorithm was used [19]. The following parameters were included in the fitting routine: the beam width, the sample length, the resolution, the background and the detector offset. The scattering length density (SLD) profile (Figure 7) represents the distribution of the scattering potential of the sample into the depth of the film (in **Z** direction Figure 5). The results for all studied samples are shown to coincide qualitatively with insignificant numerical variations associated with the concentration of cobalt in the granular film (Figure 7). All PNR curves at different temperatures and **B** were fitted with the same structural parameters: thickness (d), roughness (*σ*), and nuclear scattering length density (N*n*). The magnetic scattering length density (N*m*), which depends on the applied magnetic field, varied in the corresponding datasets (Tables 2 and 3).

From Figure 7, it is clear that a good quantitative agreement between model calculations and experimental data has been achieved. The difference in the granular film thickness, d, defined by PNR and XRR methods is related to different batches of samples.

**Figure 7.** The SLD profiles *ρ* <sup>+</sup> (dot line) and *ρ* − (dashed line) of the samples Au/SiO<sup>2</sup> (60 at.% Co)/Si (**a**–**c**) and Au/SiO<sup>2</sup> (70 at.% Co)/Si (**d**–**f**) at B = 240 mT for T = 420 K (**a**,**c**), 300 K (**b**,**e**), and 120 K (**c**,**f**). The solid line corresponds to the SLD profile at H = 0 T for every polarization of neutrons.




**Table 2.** *Cont.*

**Table 3.** Parameters of the model used to approximate the PNR curves for Au/SiO<sup>2</sup> (70 at.% Co)/Si.


#### *2.4. SQUID Magnetometry*

Static magnetic properties of the samples were characterized by SQUID magnetometry (Quantum Design MPMS-5S) at the Institute of Condensed Matter Physics (Braunschweig, Germany). The magnetization curves M(H) were recorded up to 5 T. The magnetic field was applied parallel to the plane of the film. As an example, M(H) of the sample SiO<sup>2</sup> (75 at.% Co)/Si taken at T = 200 K and 300 K are shown in Figure 8 after subtraction of the diamagnetic contribution of the Si substrate and the SiO<sup>2</sup> matrix. A shift of the magnetic hysteresis loops along the field axis and a change of the coercive field as a function of temperature are observed. (insert (a) in Figure 8). In the fields region of 0.8 ÷ 1.8 T there are specific hysteresis loops—"pockets" (inserts (b,c) in Figure 8).

The actual magnetic saturation value of the films is quite difficult to detect using SQUID techniques because the overall high-field signal includes the dominant diamagnetic contribution of the substrate, appearing due to its relatively large volume. If the granular film itself shows a small linear response in the high field regions, this is convoluted with an overwhelming diamagnetic signal and it is challenging to reliably subtract the latter contribution in order to separate the magnetic properties of the nanolayer itself. The saturation magnetization value of the SiO2(Co) film obtained by SQUID magnetometry (Figure 8) M*<sup>s</sup>* <sup>≈</sup> 60 emu/cm<sup>3</sup> differs significantly from the saturation magnetization of bulk cobalt M*Co* = 1422 emu/cm<sup>3</sup> .

**Figure 8.** Magnetic field dependence of the SiO<sup>2</sup> (75 at.% Co)/Si sample magnetization at magnetic field from −5 T to 5 T. Inserts: (**a**)—region of the coercive fields, (**b**,**c**)—regions of the specific hysteresis loops.

#### **3. Discussion**

The GISAXS data from the samples Au/GF/Si has clearly shown two nanoparticle subsystems in GFs with different average interparticle distances *l*<sup>1</sup> and *l*<sup>2</sup> (Figure 3). The interparticle distance *l*<sup>1</sup> for the first subsystem increases with the increasing cobalt concentration: 64 ± 2 Å for SiO<sup>2</sup> (38 at.% Co)/Si, 70 ± 2 Å for SiO2(54 at.% Co)/Si, and 88 ± 2 Å for SiO<sup>2</sup> (82 at.% Co)/Si (Figure 3d). It is clear that the determined interparticle distance *l*<sup>1</sup> for these concentrations of Co should be equal to the size of nanoparticles because of the percolation limit. As one can see, the value of *l*<sup>1</sup> obtained by GISAXS is in good agreement with the particle size observed by SEM on the top surface of the GF, where percolation is clearly visible (Figure 1).

In contrast, the second subsystem with a characteristic distance *l*<sup>2</sup> = 300 nm is not revealed by SEM measurements. This is due to the fact that it is most likely located at the burred GF/Si interface. A similar structure containing two subsystems of nanoparticles in the granular films SiO<sup>2</sup> (x at% Co) deposited on GaAs substrate was observed in [17,24]. There it was shown that the specific interface layer of cobalt nanoparticles is formed between the GF and GaAs substrate with an in-plane interparticle distance ∼ 320 Å and the thickness ∼70 Å [17,24].

The formation of such a specific layer on the GF/Semiconductor interface may be caused by a number of reasons. The deposition of cobalt on a substrate at the initial stage leads to the clustering of a thin cobalt film due to the dynamics of Co on Si, or GaAs [25–28]. The activation energy for clustering was found to be 0.3 ± 0.2 eV [25], indicating that a Co surface diffusion process dominates the kinetics. The surface energy of cobalt is about 2500–2900 erg/cm<sup>2</sup> , while the surface energy of silicon corresponds to 1230 erg/cm<sup>2</sup> for the (111) plane, 1510 erg/cm<sup>2</sup> for the (110) plane, and 2130 erg/cm<sup>2</sup> for the plane (100). The surface energy of GaAs is 860 erg/cm<sup>2</sup> [29,30]. Another reason for the formation of a specific layer at the interface is presumably related to the process of heteroepitaxial growth of the GF by ion beam co-sputtering directly onto a semiconductor substrate, while the next layers are grown in the homoepitaxial regime [31].

However, the X-ray and Polarized neutron reflectometry data (Figures 4, 6 and 7) do not confirm the coexistence of two nanoparticle subsystems in GFs, or they are not sensitive enough to the existence of a specific layer at the GF/Semiconductor interface in SiO<sup>2</sup> (x at.% Co)/Si samples. In contrast to the GISAXS results, the film model which fits the XRR and PNR data in the best way consists of only two layers on the Si substrate: a granular film layer and the capping Au layer, i.e., there is no specific layer on the GF/Semiconductor interface. Why do the PNR method, which is highly sensitive to changes in magnetic properties, and the XRR method, which is highly sensitive to changes in electron density contrast, fail to detect an interface layer?

The next question deals with the study and interpretation of the magnetic properties of the granular films. The magnetic scattering length density N*<sup>m</sup>* has the opposite sign for neutrons with initial polarization +**P**<sup>0</sup> and −**P**<sup>0</sup> (along and against **B**) and its value is directly proportional to the magnetization component **M** parallel to **B**:

$$M = \left(\frac{m\_n \mu\_n}{2\pi\hbar^2}\right)^{-1} \frac{N\_m}{4\pi} \,\mathrm{\,} \tag{1}$$

where m*<sup>n</sup>* and *µ<sup>n</sup>* are the mass and magnetic moment of the neutron. Comparing SLD profiles *ρ* <sup>+</sup> = N*<sup>n</sup>* + N*<sup>m</sup>* and *ρ* <sup>−</sup> = N*<sup>n</sup>* - N*<sup>m</sup>* corresponding to R(Q*z*, +**P**0) and R(Q*z*, −**P**0), respectively, one can obtain a value of the magnetization in every layer forming the sample film.

From the non-zero values of N*<sup>m</sup>* at the SLD profiles (Figure 7) one can conclude that only GF is magnetized under the applied field and tends to saturate the magnetic moment at B ≥ 300 mT. Magnetizations of the GF in SiO2(60 at.% Co)/Si and SiO<sup>2</sup> (70 at.% Co)/Si samples calculated from N*<sup>m</sup>* according to Equation (1) are presented in Figure 9 as a function of the applied magnetic field for T = 120, 300 and 420 K. In Figure 9, the magnetic SLD N*<sup>m</sup>* is normalized to magnetic SLD of bulk cobalt *<sup>N</sup>m*(Co) = 4.24 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>Å</sup>−<sup>2</sup> . From Figure 9, it follows that the GF is in the ferromagnetic state at all three temperatures. Its magnetization decreases with increasing temperature: at the maximum applied field B = 320 mT we found that N*m*/N*m*(Co) = 0.21 at T = 120 K while decreasing down to 0.17 at T = 300 K and to 0.14 at T = 420 K for SiO<sup>2</sup> (60 at.% Co) and decreasing down to 0.14 at T = 300 K and to 0.09 at T = 420 K for SiO<sup>2</sup> (70 at.% Co). Why is the magnetization, determined by the PNR method so much smaller than the magnetization of bulk Co? It is smaller by a factor of 5 to 10 at different temperatures.

Let us analyze the nuclear scattering length densities of neutrons and X-rays for all samples in order to answer these questions. The fit parameters in PNR and XRR models show excellent agreement between the measured and expected nuclear SLD values of Si (2.073 <sup>×</sup> <sup>10</sup>−6Å−<sup>2</sup> and 2.024 <sup>×</sup> <sup>10</sup>−5Å−<sup>2</sup> , respectively). The fits of the nuclear SLD values of Au in the PNR and the XRR experiments are slightly lower than the expected bulk values (4.5 <sup>×</sup> <sup>10</sup>−6Å−<sup>2</sup> and 12.47 <sup>×</sup> <sup>10</sup>−5Å−<sup>2</sup> , respectively). The reason could be the discontinuous layer of gold since its thickness is below 50 Å. A number of studies

demonstrate changes in the physical and chemical properties of thin gold films associated with their island structure [32–34].

**Figure 9.** Magnetic field dependence of fitted magnetic SLD value N*m* normalized to magnetic SLD of bulk cobalt N*m*(Co) for the ferromagnetic films SiO<sup>2</sup> (70 at.% Co) (**a**) and SiO<sup>2</sup> (60 at.% Co) (**b**).

The measured nuclear SLD of the GF layer is substantially higher for neutrons and lower for X-rays as compared to expected values (Tables 1–4). Table 4 shows the calculated N*<sup>n</sup>* and SLDX-ray values for a given ratio of Co/SiO<sup>2</sup> in GF under the assumption that cobalt can interact with oxygen from the silicon oxide dielectric matrix or oxidize in air. Unfortunately, there is no reliable way to determine the thickness of the oxide layer on the surface of cobalt nanoparticles in the present study. It is obvious that the cobalt grains must be partially oxidized for the nuclear SLD to be higher for neutrons0 or lower for X-rays, compared to the values for a given ratio of Co/SiO<sup>2</sup> in GF. All cobalt oxides are antiferromagnets, in contrast to pure cobalt (Table 5). Magnetization measurements (Figure 8) reveal that for SiO2(x at.%Co)/Si heterostructures exchange biasing of the magnetic hysteresis loop along the field axis is observed. The exchange bias describes a phenomenon associated with interfacial coupling between ferromagnetic and antiferromagnetic materials [35]. In our investigation, the exchange biasing of the magnetic hysteresis loop confirms the oxidation of cobalt nanoparticles. The exchange biasing field, like the coercive field, changes as a function of temperature and is about 10 percent of the coercive field (insert (a) in Figure 8).

Another small but well-distinguished feature of the magnetization curve is the small hysteretic pockets in the field range from 0.8 to 1.6 T. They are well-detected both in the positive and negative halves of the magnetization curve. They were also detected for GF on the GaAs substrate [36]. An interesting feature of this hysteresis loop is counterintuitive

large values of the magnetization when the field increases and smaller magnetization values when the field decreases. This feature can be referred to as the same biasing phenomenon of the Co (ferromagnetic) particles with the oxidized (antiferromagnetic) surface, to the so-called "core-shell" model. In this case, the magnetization of the Co particle is responsible for the large ferromagnetic response of the GF system, while the antiferromagnetic shell may show a spin flop in a certain field range, which may be energetically more favorable than the collinear spin arrangement. Considering its really small (4% only) contribution imposed on the top of the magnetization curve, we point out the biased (ferro-antiferro) magnetic nature of this phenomenon but do not go into any further details.

**Table 4.** The calculated N*<sup>n</sup>* and SLDX-rays values for a given ratio of Co/SiO<sup>2</sup> in GF under the assumption that cobalt can be oxidized.


**Table 5.** The calculated N*<sup>n</sup>* and SLDX-rays values are given with Density (D), Molecular weight (MW) and Number of molecules per Å<sup>3</sup> (n) for compounds (Comp) that may be formed in the granular film, or at the interfaces.


Figure 10 demonstrates the neutron spin asymmetry calculated by:

$$SA = \frac{R(Q\_{z\nu} + \mathbf{P}\_0) - R(Q\_{z\nu} - \mathbf{P}\_0)}{R(Q\_{z\nu} + \mathbf{P}\_0) + R(Q\_{z\nu} - \mathbf{P}\_0)},\tag{2}$$

As the neutron spin asymmetry is a very reliable method to study the intrinsic nanomagnetic features [37,38], it provides direct evidence that these magnetic properties originate from the composite thin film region and the substrate and surface play no detectable role. The spin asymmetry is purely related to the local magnetic moment in the sample and characterizes the in-depth distribution of the magnetization in GF directly. The observed peak SA occurs at Q*<sup>z</sup>* = 0.0136 <sup>±</sup> 0.004 Å−<sup>1</sup> and exceeds 80% (Figure 10a), which competes with that of standard magnetic alloys [39] and is large compared to the spin asymmetries typical of dilute magnetic systems [40]. The magnetic moment of the sample proportional to the SA peak area increases with the applied external field and saturates at the experimental field of 240 ÷ 320 mT (insert in Figure 10a). In contrast to SQUID

magnetometry, PNR is unaffected by complications of diamagnetic or weak paramagnetic responses because they produce very low magnetizations per volume. Instead, the PNR signal is dominated by scattering effects from large local net magnetization, typically caused by a well-aligned collinear spin structure. Consequently, the small value of the neutron spin asymmetry at low fields (insert in Figure 10a) confirms the traits noted in the hysteresis from the SQUID data (Figure 8)—a low remanence and a low coercive field (B*<sup>c</sup>* < 1.4 mT) which are attributed to either superparamagnetic behavior or a wandering anisotropy axis. This shows that, without a sufficient applied field, the magnetic moments in GF do not stay well aligned in any specific direction. The possible reason for this is that there are inhomogeneous magnetically ordered patches in the GF layer that experience a distribution of local spin anisotropies and consequently have a distribution of magnetization angles that are only co-aligned in a magnetic field and relax under zero field. In addition, the magnetic moment in GF changes moderately with temperature (insert in Figure 10b), indicating a strong magnetic exchange coupling.

**Figure 10.** The neutron spin asymmetry, SA, for Au/SiO<sup>2</sup> (70 at.% Co)/Si fitted (gray line) to the chemical and magnetic model previously described in the text (Figure 7): (**a**)—under applied magnetic fields up to 320 mT, (**b**)—for different temperatures at B = 240 mT. Insert in panel (**a**)—the field-dependence of the integral intensity of the peak spin asymmetry for different temperatures. Insert in panel (**b**)—the temperature dependence of the integral intensity of the peak spin asymmetry at B = 240 mT. For clarity, data for (**a**) at B = 40, 100, and 320 mT are omitted.

As it is shown in Figure 10, the SA changes its sign upon the increase in the momentum transfer Q*z*. We interpret this observation as the presence of the "magnetically dead" layer in the GF at its interface with the Si substrate. The thickness of this layer is difficult to determine exactly, but it should be of the order of (1/10) of the thickness of the granular film itself. A similar effect of the magnetically dead layers in the PNR was observed in the multilayered structures (see, for example, [41,42].)

Based on the above, it should be assumed that the oxidation of cobalt nanoparticles is the reason for such a low magnetization value of the granular film with cobalt concentrations from 54 at.% to 82 at.% compared to pure cobalt. The same difference of the magnetic SLD value normalized to magnetic SLD of bulk cobalt was observed for GF deposited on GaAs substrate in [43].

The answer to the question of why the highly sensitive PNR and XRR methods fail to detect an interface layer between GF and Si substrate will also allow us to answer the main question of this article: what is the reason for the dramatic difference in IMR coefficients for identical SiO2(x at.% Co) granular films deposited on GaAs, or Si substrates.

It was noted above that Co due to the higher surface free energy as compared to Si would tend to grow as islands at the beginning of the deposition. The phase composition and magnetic properties of cobalt films deposited on the silicon substrate were studied in [36]. It was shown that as the film thickness increases, an interface cobalt silicide and an island film of a solid solution of silicon in cobalt are successively formed. The growth of metallic cobalt begins only after the deposition of a Co layer with a thickness of more than 7 Å. The calculated N*<sup>n</sup>* values are given in Table 5 for chemical compounds that may be formed in the granular film, or at the interfaces due to boundary interlayer diffusion leading to the formation of an interface of variable binary composition—cobalt silicides. Table 5 shows that cobalt silicides and cobalt oxides have similar SLDX-rays values, which does not allow the GF and an additional interface layer to be distinguished in the Z direction (Figure 2) by XRR experiments. Wherein in the Y direction (Figure 2), an additional interface layer with large periodicity *l*<sup>2</sup> = 300 ± 10 Å is detectable by GISAXS due to the high contrast of SLDX-rays between the multiphase cobalt granules and the SiO2. Thus, both the XRR and PNR experimental results are well described by the model of a uniform granular cobalt film with a noticeably oxidized surface of the granules and a high surface roughness at the interface with the silicon substrate.

Increasing the silicon concentration to 30% in the Co-Si solid solution makes the large cobalt granules practically non-magnetic due to the charge transfer effect. Moreover, due to the formation of cobalt silicide at the GF/Si interface, the conductivity of the granular film will be determined by the chains "granule (Co)-metal (Co-Si)-semiconductor (Si)", and not by electron tunneling between ferromagnetic granules through the SiO<sup>2</sup> dielectric matrix into the Si semiconductor. This leads to a strong decrease in the magnetoresistive coefficient to a few percent in SiO<sup>2</sup> (x at.% Co) heterostructures on a Si substrate compared to the really giant IMR value in SiO<sup>2</sup> (71 at.% Co) on a GaAs substrate [8].

#### **4. Conclusions**

In the present paper, we studied the structure and magnetic properties of the maze structure with interconnected cobalt particles SiO<sup>2</sup> (x at.% Co) deposited on Si substrate at 54 at.% ≤ x ≤ 82 at.%). Investigations were carried out using Polarized Neutron Reflectometry, Grazing Incidence Small Angle X-ray Scattering, X-ray Reflectometry, Scanning Electron Microscopy, and Superconducting Quantum Interference Device Magnetometry.

The reason for the dramatic difference in IMR coefficients for SiO<sup>2</sup> (x at.% Co) heterostructures deposited on GaAs or Si substrates was established. Despite the fact that both granular films have very similar values for thickness, interparticle distances, and magnetizations, as well as the additional layer of Co at the interface granular film/substrate, there is a difference in the interface morphology for SiO<sup>2</sup> (x at.% Co)/Si. On the interface granular film/Si substrate there is boundary interlayer diffusion of Si atoms into Co, leading to the formation of an interface of variable binary composition - cobalt silicides. The formation of cobalt silicide at the GF/Si interface leads to the conductivity of the granular film determined by the electron movement via chains "ferromagnetic granule (Co)-metal (Co-Si)-semiconductor (Si)", and not by the electron tunneling between ferromagnetic granules through the SiO<sup>2</sup> dielectric matrix into the Si semiconductor. A high value of the GIMR effect in SiO<sup>2</sup> (x at.% Co) on GaAs [8,44] substrate is probably related to the spin-dependent potential barrier formed in the accumulation electron layer in the

semiconductor near the interface. The action of the spin-dependent potential barrier is amplified by the avalanche process and by the electron accumulation in the quantum well in the semiconductor interface region induced by the backscattering process of injected electrons on exchange-split levels [44]. In the granular film SiO<sup>2</sup> (x at.% Co) on the Si substrate, the cobalt silicide does not allow the formation of the region at the GF/Si interface saturated with electrons capable of polarizing in the direction of the applied magnetic field. Therefore, further efforts to form a quantum well with exchange-split levels on the well top at the surface of Si are necessary to gain a high positive magnetoresistance in SiO<sup>2</sup> (x at.% Co)/Si heterostructure. For example, this can be achieved by introducing a buffer layer between the granular film and substrate [45,46].

**Author Contributions:** Conceptualization, N.A.G.; methodology, N.A.G.; validation, N.A.G. and S.V.G.; investigation, N.A.G., V.U., A.A.V.; resources, A.I.S., N.N.N. and L.V.L.; data curation, V.U. and A.A.V.; writing—original draft preparation, N.A.G. and V.U.; writing—review and editing, S.V.G.; visualization, A.A.V.; supervision, S.V.G.; project administration, L.V.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is financially supported by the Ministry of Science and Higher Education of the Russian Federation in the framework of Agreement No. 075-15-2022-830 (Prolongation of Agreement No. 075-15-2021-1358 from 12 October 2021).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** Authors would like to thank Institut Laue-Langevin, European Synchrotron Radiation Facility and Institute of Condensed Matter Physics of TU Braunschweig for the provided experimental opportunities. We would like to acknowledge M. Wolf, A.Mistonov and D.Menzel for the assistance with the neutron and magnetometry measurements. The authors are grateful to the staff of the Interdisciplinary Resource Center for Nanotechnology and the Center of X-ray diffraction studies at the Research park of Saint Petersburg State University for *preliminary research* of Au/SiO2(Co)/Si heterostructures, as well as Saint Petersburg State University for financial support (Activity 6—grant for academic mobility 2019, ID:41160111).

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

#### **References**

