Experimental investigations of the perpendicular anisotropy in thin films that are interpreted using a macrospin model often require the inclusion of a higher-order anisotropy contribution. However, recent ferromagnetic resonance experiments on [Co/Ni]$N$ multilayers indicate that the macrospin model cannot explain the full angular dependence in this system. Using micromagnetic simulations of a system with lateral variations of the second-order perpendicular uniaxial anisotropy, we show that while the macrospin model is able to capture the average properties of the system along high symmetry orientations by including a higher-order anisotropy, the model fails to reproduce the full angle dependence. Our studies provide another indication that higher-order anisotropies may not be intrinsic to these systems but instead may be caused by inhomogeneities.

## I. INTRODUCTION

Spintronic devices require precise control of the magnetic anisotropy and the magnetization dynamics. The origin of the perpendicular anisotropy has been studied extensively.^{1–4} Perpendicular anisotropy in multilayers arises from the broken symmetry at the interfaces due to enhanced orbital and spin moments relative to the bulk system.^{5,6} From a phenomenological standpoint, perpendicular anisotropies can be classified by the power of the directional cosine along the symmetry axis. The first non-vanishing term is second-order, i.e., proportional to the square of the directional cosine, while the next higher-order term is of fourth-order. In some nomenclatures, these terms are also referred to as first-order and second-order anisotropies. Theoretically, higher-order anisotropy values are predicted to be significantly smaller than the second-order anisotropy values.^{7–10} However, unusually large higher-order anisotropy values have been reported experimentally in numerous publications.^{11–13} Despite this, the origin of the higher-order anisotropy is still not well understood. However, analytical models^{14,15} have shown that the lateral fluctuation of the second-order perpendicular anisotropy can give rise to higher-order anisotropies. Similarly, Sun^{16} has shown that an inhomogeneity along the film normal can give rise to an apparent higher-order anisotropy. Recently, using micromagnetic calculations, Mohammadi *et al.*^{17} have shown that lateral inhomogeneities of the second-order uniaxial anisotropy in thin film samples with perpendicular uniaxial anisotropy can give rise to significant higher-order anisotropy contributions.

In this article, we present the results of fully micromagnetic simulations using a periodic checkerboard structure with lateral variations of the second-order perpendicular uniaxial anisotropy with no higher-order anisotropy. We demonstrate that the polar angle dependence of these micromagnetic simulations can be described to some extent using a macrospin model with second- and fourth-order anisotropy terms. However, as we will show, this model does not always fully capture the micromagnetic results, thereby indicating limitations of the macrospin model, in particular, for cases when the lateral size of inhomogeneities of the second-order uniaxial anisotropy is larger than the exchange length of the system.

## II. MOTIVATION

Many experimental results of samples with a perpendicular anisotropy can be well explained using a macrospin model by including a second-order $K~2,eff$ and a fourth-order $K4,eff$ uniaxial anisotropy term,^{5,17}

where $\theta M$ is the angle of the magnetization measured with the magnetization easy axis (film normal). This macrospin model can be used to predict both quasistatic and dynamic properties of typical samples under investigation by fitting the second-order and fourth-order anisotropy constants. An example for this can be found in the supplementary material but also in numerous publications.^{18–26} However, careful inspection of the polar angle dependence of the ferromagnetic resonance data can reveal deviations of the experimentally observed dynamic response and the one predicted by the macrospin model; see, for example, Fig. 1(b) in the supplementary material. As the anisotropy constants are typically determined using a least-square fit, the remaining discrepancy between the macrospin model and experimental data may be an indication of a limitation of the model itself. However, typically the experimental data in itself offers not much insight why the model fails. Attempting to describe an experimental thin film system using full micromagnetic simulations would require numerous assumptions regarding the involved anisotropies, their distribution, and the length scale of their variation that are extremely difficult to determine independently. Therefore, to better understand this discrepancy, we carried out full micromagnetic simulations using our MATLAB based micromagnetic code M$3$ and analyzed the results using a macrospin model. This approach gives us complete control over the anisotropies in the micromagnetic model and thereby enables us to trace any discrepancy better than this can be done in an experiment.

## III. MICROMAGNETIC MODELLING

Micromagnetic calculations are carried out using our MATLAB based finite differences code M$3$, which uses a Fast Fourier Transform (FFT) method to calculate the magnetostatic interaction field, utilizing Newell’s formulation of the demagnetizing tensor at short distances and a point-dipole approximation for the far field. This code also offers different implementations of the exchange interaction, including 6, 12, and 26 neighbor methods. We can assign arbitrary anisotropies (unidirectional, uniaxial, cubic, etc.) to the simulation volume.

We created lateral variations of the second-order perpendicular uniaxial anisotropy in the form of a periodic checkerboard structure where the anisotropy is assumed to be constant across the film thickness. The gyromagnetic ratio was set to $\gamma =2.21\xd7105$ mA$\u22121$s$\u22121$, the damping constant $\alpha =0.007$, the saturation magnetization $Ms=1000$ kA/m, and the exchange constant $A=1\xd710\u221211$ J/m. Simulations were carried out with a second-order perpendicular anisotropy $K2,A=\u22121.5\xd7106$ J/m$3$ for regions A, which is sufficient to overcome the demagnetizing field of an infinite film resulting in an easy axis along the film normal. The second-order perpendicular anisotropy for regions B was $K2,B=0,\u22120.15\xd7106,\u2026,\u22121.5\xd7106$ J/m$3$ for different simulations. For all micromagnetic simulations, no fourth-order aniosotropy was included; i.e., $K4=0$ for regions A and B. However, the difference between the second-order anisotropies of the two regions causes frustration, which leads to the emergence of a fourth-order anisotropy contribution when one tries to describe the system using a macrospin model. The wavelength $Lx,y$ of the checkerboard pattern is the same along both $x$ and $y$ directions. The cellsize throughout the simulation volume is $1\xd71\xd71nm3$. These dimensions are smaller than the exchange length,^{17,27} calculated as $\lambda A=A/|K2,A|\u22483$ nm for a system with second-order perpendicular uniaxial anisotropy in regions A $(K2,A=\u22121.5\xd7106$ J/m$3)$.

## IV. MAGNETIZATION DYNAMICS

The magnetization dynamics in ferromagnetic systems is described by the Landau–Lifshitz–Gilbert equation of motion,^{28–30}

where $M\u2192$ is the magnetization vector; $Ms$ is the saturation magnetization; $H\u2192eff$ is the effective magnetic field, which includes external, exchange, anisotropy, and demagnetization fields; $\gamma =g\mu B\u210f$ is the gyromagnetic ratio; and $\alpha $ is the damping constant. In the above equation, the first term represents the precession, while the second term represents the dissipation or damping.

The dynamic response of the system in our simulation is probed by applying a $2$ T static magnetic field in the $yz$ plane. Prior to recording the dynamic response, the system is relaxed in the presence of $6$ mT field along the $x$-direction, i.e., always perpendicular to the static magnetic field as shown in Fig. 1. Next, this small field is removed at $t=0$, and the time evolution of the magnetization, described by the Landau–Lifshitz–Gilbert (LLG) equation of motion, is recorded. Fundamental resonance frequencies $(fres)$ at various polar angles $(\theta H)$ are extracted using a Fast Fourier Transform (FFT) of the time evolution of the magnetization. The details of this method have been discussed in a previous publication.^{17} Fundamental resonance frequencies at various external static field angles from the checkerboard structure surface normal at a constant static field are shown in Fig. 2. For these simulations, we always used the same second-order uniaxial anisotropy in regions A with $K2,A=\u22121.5\xd7106$ J/m$3$ but used three different values of the second-order uniaxial anisotropy constant for regions B. First, we assigned the same value of the second-order uniaxial anisotropy to regions B of a checkerboard structure: $K2,B=\u22121.5\xd7106$ J/m$3$; i.e., no inhomogeneity is introduced in the system, which is represented by green circles. Second, as shown by red circles, we assigned half of the second-order uniaxial anisotropy of regions A to regions B: $K2,B=\u22120.75\xd7106$ J/m$3$; i.e., a moderate inhomogeneity is introduced in the system. Last, we assigned no second-order uniaxial anisotropy to regions B: $K2,B=0$; this is the most inhomogeneous system studied and represented by black circles.

## V. DATA ANALYSIS

In the following, we will attempt to describe the full micromagnetic simulations using the macrospin model introduced in Sec. II. This model has been used extensively to describe experimental data of the magnetization dynamics of systems with perpendicular anisotropy^{5,31,32} including polar angle dependent measurements.^{18–26,33,34} As noted earlier, in order to obtain a reasonable description of the angular dependence of the magnetization dynamics in the systems with perpendicular anisotropy, a second-order $K~2,eff$ and fourth-order $K4,eff$ anisotropy terms are required in the energy density. $K~2,eff$ contains the contribution from shape anisotropy and the uniaxial perpendicular anisotropy, which can be expressed as $K~2,eff=\mu oMs22+K2,eff$ if one assumes that the film is infinite and homogeneously magnetized. $K4,eff$ is the fourth-order uniaxial anisotropy term. In this notation, the negative values of the anisotropy constants indicate that the easy axis is along the film normal.

Using the energy density from Eq. (1) for the macrospin model, a least-square fitting code calculates the resonance frequency at a constant external static magnetic field using the Smit–Beljers relation^{30,35}

where $\omega =2\pi f$ is the angular frequency, $\theta $ is the polar angle of magnetization with respect to film normal, $\varphi $ is the azimuthal angle of magnetization, and $E$ is the magnetic free energy density. The derivatives are taken at the points of equilibrium $(\theta M,\varphi M)$ of the magnetization for which

To illustrate the importance of the information contained in the full angular dependence, we use two different approaches to calculate the second-order $K~2,eff$ and fourth-order $K4,eff$ anisotropy constants of the macrospin model. In the first approach, we only use two data points along two high symmetry directions of the system, namely, $0\xb0$ and $90\xb0$. In this case, one obtains analytical in-plane (ip) and out-of-plane (oop) Kittel equations with second- and fourth- order anisotropy that enable a straightforward calculation of the anisotropy values,^{5,17,36}

where $\gamma \u2032=\gamma 2\pi $ is the reduced gyromagnetic ratio, and the effective magnetizations are

Using Eqs. (7) and (8) in (5) and (6), we find $K~2,eff$ and $K4,eff$ using only the in-plane ($90\xb0$) and out-of-plane ($0\xb0$) frequency data points. With the effective anisotropy values determined this way, one can simulate the overall polar angle dependence, shown as dashed lines in Fig. 2. As these calculated anisotropy values only depend upon the high symmetry directions, it is not surprising that they do not capture the full angular dependence of the micromagnetic data.

Hence, we use a second approach where the anisotropy constants of the macrospin model are determined by fitting the entire polar angle dependence of the micromagnetic simulations, and this is shown as solid lines in Fig. 2. As this approach takes into consideration the entire angular dependence to determine the second-order $K~2,eff$ and the fourth-order $K4,eff$ uniaxial anisotropy constants, one expects a better agreement of the macrospin model with the micromagnetic simulations.

As shown in Fig. 2, we considered three different cases of inhomogeneity to illustrate the adequacy of the macrospin fitting models in explaining micromagnetic simulations. In the case of no inhomogeneity in the system, i.e., both regions A and B having same second-order anisotropy constants, both approaches capture the micromagnetic simulation exactly as illustrated by the green dashed and solid lines. As expected, if there is no inhomogeneity in the system, a macrospin model is sufficient to explain these micromagnetic simulations. In the second case, where we assigned a moderate inhomogeneity to the system with half the value of second-order anisotropy in regions B as compared to A, the macrospin model using anisotropy values based on high symmetry directions (red dashed line) perfectly explains the in-plane ($90\xb0$) and out-of-plane ($0\xb0$) micromagnetic simulation data as one expects. However, one notes that this approach does not capture the remaining micromagnetic simulation data points adequately. The solid red line, which represents the macrospin model fit using the entire angular dependence, also cannot fully describe the micromagnetic simulations. In the third case, where we assigned full strength of a second-order anisotropy value for regions A while no second-order anisotropy for regions B, i.e., the most inhomogeneity in the system, the deviations of the macrospin model become even more pronounced. These observations indicate that while the macrospin model is able to reasonably describe systems with small variations of the second-order uniaxial anisotropy by introducing a fourth-order anisotropy, this approach fails for systems with significantly different second-order anisotropies in different regions of the sample.

To further investigate the limitations of the macrospin model, we have studied the length scale dependence of the inhomogeneities. For this study, we are considering the case with the largest inhomogeneity in the system by using the full value of second-order perpendicular anisotropy for regions A, i.e., $K2,A=\u22121.5\xd7106$ J/m$3$ and no second-order perpendicular anisotropy for regions B, i.e., $K2,B=0$, where for $Lx=Ly=50$ nm, the deviations of the macrospin model from the micromagnetic results were well pronounced, as discussed previously. As shown in Fig. 3(a), when the anisotropy variation occurs on a small length scale $Lx=Ly=2$ nm, the macrospin model can describe the micromagnetic simulations very well. When the length scale of the anisotropy inhomogeneity is increased to $Lx=Ly=20$ nm, the deviation of the macrospin model description from the micromagnetic simulations becomes noticeable. This deviation further increases when the length scale is increased to $Lx=Ly=100$ nm. With increasing length scale of the anisotropy inhomogeneity, the magnitude of the fourth-order anisotropy required in the macrospin model also increases. The micromagnetic simulations further show that with increasing inhomogeneity length scale, the resonance becomes increasingly localized, and as expected, the relevant length scale for this localization is determined by the exchange length of the system. For the smallest length scale of the anisotropy variation studied here, the contributions of both regions to the power spectral density of the fundamental resonance are equal and do not change with the angle of the applied field [black data points in Fig. 3(b)]. In this case, the system to a good approximation can be described with a macrospin model. However, with increasing length scales $Lx,y$, both regions no longer contribute equally to the power spectral density of the fundamental resonance of the system for the out-of-plane and in-plane configurations. Furthermore, as shown in Fig. 3(b), the power spectral density contribution now strongly depends on the angle of the applied field; that is, the spatial localization of the resonance is now angle dependent. Given the highly nonlinear nature of this localization, it is not surprising that a macrospin model is unable to reproduce the angle dependence of the ferromagnetic resonance if one limits the effective anisotropies to the lowest two contributions.

## VI. CONCLUSIONS

We have discussed the limitations of the macrospin model with a second- and fourth-order anisotropy term to explain the magnetization dynamics of a complex inhomogeneous system using a periodic checkerboard structure with lateral variations of the second-order uniaxial anisotropy but no fourth-order anisotropy. The variation of the second-order uniaxial anisotropy in regions A and B causes frustration in the system giving rise to a fourth-order anisotropy term when describing the system using the macrospin model. With increasing frustration, the deviation between the full micromagnetic simulations and the macrospin model description of the data increases. With increasing length scale of the inhomogeneities in the system, these deviations become increasingly apparent in the polar angle dependence of the resonance. While the macrospin model is able to reproduce the resonance condition for the in-plane and out-of-plane orientation, it is not able to reproduce the full polar angle dependence using second- and fourth-order anisotropies. Analysis of the micromagnetic data shows that in these cases, the power spectral density contribution of the different sample regions shows a complex polar angle dependence. Our results suggest that the inability to describe experimental polar angle dependent ferromagnetic resonance data using a macrospin model may be another indication that fourth-order anisotropies observed in these systems are caused by inhomogeneities rather than being an intrinsic property of the materials. This also implies that full polar angle dependent measurements are crucial to determine the origin and obtain a better understanding of the higher-order anisotropy.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the experimental ferromagnetic resonance results of [Co/Ni]$N$ multilayers.

## ACKNOWLEDGMENTS

Part of this research was supported by the Defense Advanced Research Projects Agency (DARPA) program on Topological Excitations in Electronics (TEE) under Grant No. D18AP00011. A.S. and C.M. also acknowledge support through National Science Foundation (NSF)-CAREER under Grant No. 1452670. Furthermore, we would also like to acknowledge the support of the University of Alabama High Performance Computing facilities and staff.

## DATA AVAILABILITY

The data that supports the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Handbook of Nanomagnetism*(Pan Stanford Publishing, 2015), pp. 71–96.