Numerical study of magnetic properties for Co-Fe-Nb nanofilms as promising materials for magnetoresistive memory

Spintronic devices are much more efficient than semiconductor electronics devices by using the spin of an electron rather than its charge to record information. They use much less energy, and the density of information recording is many times higher. Spintronics materials can be based on various mechanisms and phenomena, such as the effect of giant magnetoresistance[1]} or tunneling magnetoresistance. Currently, these effects are already used in tunneling magnetoresistive sensors in hard disk read heads. This has been responsible for a significant increase in the information storage density of hard disks[2]. In addition, more and more efforts are being made to develop different types of magnetoresistive random access memory (MRAM)[3-5].

One of the main tasks for spintronics is the search for new promising materials. An alternative name for spintronics is magnetoelectronics, as a result, materials for spintronic devices should have a number of magnetic properties. Therefore, ferromagnets in various heterostructured systems such as ferromagnet-paramagnet[6], ferromagnet-superconductor[7], and ferromagnet-antiferro-magnet[8] are considered as one of the most common materials.

At the moment, there is an active search for combinations of elements and materials for spintronic heterostructures. Antiferromagnets (AFM) in combination with ferromagnets, for example, give good results for writing information, although they do not have high efficiency in reading, and the ferromagnet-superconductor structure faces a number of limitations due to differences in operating temperatures and the Meissner effect. The purpose of this work is to study the formation processes, structure and properties of nanomaterials for spintronics devices by using theoretical methods and mathematical modeling. The work is a development of earlier studies of the formation processes for photovoltaic converters[9], multilayer superconductor-ferromagnet nanostructures[10-12], the growth of ordered arrays of luminescent matrices based on templates made of porous aluminum oxide[13], the technology of creating composite nanoparticles for special purposes and probabilistic analysis of the mechanisms of their formation and growth[14].

Magnetoresistive memory has proven to be one of the most efficient, least power-consuming, and quite compact memory technologies. Magnetic Random Access Memory (MRAM)[15-19] structure is based on a spin valve with giant magnetoresistance and consists of three layers: two ferromagnetic layers and an insulating nanofilm between them. In any of the presented varieties of memory it is possible to use the ferromagnet-superconductor structure as a spin valve or magnetic tunnel junction. The apparatus of mathematical modeling can be used to find a composition and structure of memory cells, in which their disadvantages can be reduced or minimized, and the characteristics of memory devices can be optimized.

Predictive modeling of the deposition process and magnetic behavior of atomic spins in spintronic nanofilms can significantly speed up the sample manufacturing process, find optimal conditions for the formation of nanomaterials, such as substrate temperature, deposition rate, conditions for minimal interdiffusion at nanofilm interfaces, and many other features.

The equations of motion of atoms and their spins can be obtained by applying Poisson brackets to the generalized Hamiltonian. This approach is described in detail in Ref.[20]. Taking into account relaxation processes leads to the model of stochastic combined (lattice and magnetic) molecular dynamics to the following form:

where \( i=1, 2, \ldots N \) and \( N \) is the number of atoms in the system; \( \{\ldots, \ldots\} \) is the Poisson bracket notation; \( \hat{H} \) is the general Hamiltonian; \( {\bf r}^{\,}_i \) and \( {\bf p}^{\,}_i \) are the radius vector and momentum vector of the \( i \)-th atom; \( m^{\,}_i \) is the mass of the \( i \)-th atom; \( U({\bf r}) \) is atomic interaction potential and \( {\bf r}=\big\{{\bf r}^{\,}_1, {\bf r}^{\,}_2,\ldots {\bf r}^{\,}_N\big\} \) is the generalized variable showing the dependence on the whole set of radius-vectors of atoms; \( J(r^{\,}_{ij})>0 \) is the exchange integral for a pair of spins located at a distance \( r^{\,}_{ij} = |{\bf r}^{\,}_{ij}|=|{\bf r}^{\,}_{i}-{\bf r}^{\,}_{j}| \); \( {\bf s}^{\,}_i \) and \( {\bf s}^{\,}_j \) are the unit vectors of the spins of the \( i \)-th and \( j \)-th atoms; \( {\boldsymbol \chi}(t) \) and \( {\boldsymbol \eta}(t) \) are vectors of random forces, the meaning of which will be explained further, \( \varkappa \) and \( \lambda \) are the viscous friction force parameter and damping spin coefficient, respectively; the meaning of the vector \( {\boldsymbol \omega}^{\,}_i \) will be explained further.

The change of atomic momenta in Eq(2) is due not only to the derivative of the interatomic potential, but also to the dynamics of magnetic forces depending on the spin configurations of the particles. In this case, only exchange interactions are considered.Consider a system consisting of \( N \) magnetic atoms, described by the generalized Hamiltonian in the following form[21-23]:

where \( \hat{H}^{\,}_{\rm lat} \) and \( \hat{H}^{\,}_{\rm mag} \) are the Hamiltonians of the lattice subsystem, which take into account the spatial motion of atoms, and the magnetic subsystem, respectively.

For the microcanonical ensemble (\( NVE \)), which is used for representing systems with fixed composition (\( N \)), volume (\( V \)), and total energy (\( E \)), the Hamiltonian can be written as follows:

Hamiltonian \( \hat{H}^{\,}_{\rm lat} \) defines the mechanical behavior of atoms.

A simple expression for the magnetic Hamiltonian of \( N \) interacting particles possessing spins is defined by the expression:

where \( \mu^{\,}_B \) is the Bohr magneton, \( \mu^{\,}_0 \) is the magnetic constant, \( g^{\,}_i \) is the Lande multiplier of the \( i \)-th spin, \( {\bf H}^{\,}_{\rm ext} \) is an external magnetic field strength, and \( {\bf B}^{\,}_{\rm ext}=\mu^{\,}_0{\bf H}^{\,}_{\rm ext} \) is the induction of an external magnetic field. This external field can be either constant or time-varying. The exchange Hamiltonian (the second term in Eq(6) provides a natural connection between the spatial and spin degrees of freedom through the exchange integral \( J(r^{\,}_{ij}) \), acting on the interatomic distance \( r^{\,}_{ij} \). The form of the magnetic Hamiltonian in Eq(6) is variational. In the general case, the magnetic Hamiltonian may include additional terms responsible for various spin effects.

In accordance with the approach[24, 25] proposed by Antropov and his co-authors, the system is considered as a set of atoms, each of which is associated with a spin vector \( {\bf s}^{\,}_i \). The connection between the ordinary spin vector with the normalized one is carried out by means of the relation:

where \( \hbar \) is the reduced Planck constant.

The relaxation processes in the magnetomechanical system were modeled using the Langevin equation[26, 27] and the stochastic Landau-Lifshitz-Hilbert equation[23, 28]. The joint solution of these equations is the basis of the spin dynamics model of particles used in this work.

The parameter \( {\boldsymbol \omega}^{\,}_i \) in Eq(3) is the gyromagnetic ratio multiplied by the local magnetic field in which the spin is located. The latter value can be calculated in accordance with the following expression[23]

It was shown[18, 29, 30] that thermal and magnetic fluctuations are well described in the context of the Langevin approach. In this approach the random forces \( {\boldsymbol \chi}(t) \) and \( {\boldsymbol \eta}(t) \) are given with their characteristic first and second moments:

where \( \alpha \) and \( \beta \) are vector components, \( t \) and \( t' \) are different time periods, \( \delta(t-t') \) is the Dirac delta function,

\[ \begin{eqnarray} \nonumber D^{\,}_{\rm lat} = 2\varkappa k^{\,}_B T^{\,}_{\rm lat} \quad \mbox{and} \quad D^{\,}_{\rm s} = \frac{2\pi\lambda k^{\,}_B T^{\,}_{\rm s}}{\hbar} \end{eqnarray} \]

are the amplitudes of coordinate and spin random fluctuations, \( k^{\,}_B \) is the Boltzmann constant, \( T^{\,}_{\rm lat} \) and \( T^{\,}_{s} \) are the lattice and spin temperatures, respectively.

The lattice and spin temperatures in the considered system provide averaged measures of the distribution of spatial and spin degrees of freedom of the atoms:

Various approaches for determining thermodynamic parameters are known from the literature[31-33]. The instantaneous lattice temperature in this work was determined as the average value of atom kinetic energy[34]. The expression proposed in Ref.[35] was used to calculate the instantaneous spin temperature.

In this work, we used the potential of the Modified Embedded Atom Method (MEAM) as one of the most promising and actively developing ones[36, 37]. The potential is based on the electron density functional theory and is a multi-particle potential. The essence of the MEAM potential is defined by the following expression:

where \( U^{\,}_i({\bf r}) \) is potential energy of individual atoms of the system, \( F^{\,}_i \) is embedding function of the \( i \)-th atom, \( \phi^{\,}_{ij}(r^{\,}_{ij}) \) is pair potential between two atoms under consideration. The embedding function depends on the background electron density \( \bar{\rho}^{\,}_i \), so it takes into account the locations of all neighboring atoms of the system.

The total background electron density \( \bar{\rho}^{\,}_i \) from Eq(11) takes into account the contribution of different types of electron orbitals \( s \), \( p \), \( d \), \( f \).

The advantages of using the modified embedded atom potential for molecular dynamics are severalfold. First, this approach allows for the influence of the environment on perturbations of the electron configurations of molecules, leading to more accurate calculations of interactions and energies under realistic environmental conditions. Second, by modifying the potential, it is possible to adaptively account for changes in the environment's configuration during molecule movement or reorganization, improving the predictive power for diffusion processes, reaction dynamics, and intermolecular interactions. Third, the embedded atom approach allows for the consideration of effects such as polarization and screening, which significantly affect binding energy and transient parameters, particularly in complex systems at phase boundaries or on surfaces. This enables more realistic simulations over a wide range of conditions and better predictions of properties that are critically dependent on the local environment.

The methodology employed involves modeling atomic magnetic spins associated with spatial displacements of atoms (lattice). Magnetic spins interact with each other and with the lattice via pairwise interactions. The general form of the Hamiltonian expression for describing the total energy of magnetic systems can be represented as several contributions, including Zeeman and exchange interactions, magnetic anisotropy, Dzyaloshinsky-Moriya interactions, magnetoelectric interactions, and dipole interactions[23, 38].

The value of the exchange integral \( J(r^{\,}_{ij}) \) in Eq(9) can be estimated using the Bethe-Slater dependence[23, 28]:

Depending on the crystal lattice of the material under study, different forms of magnetic anisotropy can occur. In Refs.[39, 40] approximations for modeling spin-orbit coupling were proposed. Thus, in particular, Ref.~[41] used the functions proposed by Neel to model bulk magnetostriction and surface anisotropy in cobalt. To describe the above phenomena, it is possible to use Neel's model of pairwise anisotropy between pairs of magnetic spins. To calculate the uniaxial magnetic anisotropy in the modeled systems, it is possible to use the considered Hamiltonian.

It is also well known that the combination of the exchange interaction and the spin-orbit interaction can lead to noncollinear spin states. The most common way to model such an effect is to relate the exchange interaction to another interaction, called the antisymmetric Dzyaloshinsky-Moria interaction[38]. In particular, it is known that the Dzyaloshinsky-Moria interaction is a key mechanism for the stabilization of magnetic skyrmions.

In magnetoelectric materials and multiferroics, there is a relationship between their magnetic and electrical properties. Therefore, magnetoelectric effects should be taken into account when considering these types of materials. According to Refs.[23, 42], the mechanisms of magnetoelectric effects can be accounted for through the antisymmetric spin-spin effective interaction as a special case of the Dzyaloszynski vector.

In magnetic systems below the paramagnetic limit, the intensity of the dipole interaction is usually much smaller than that of the other magnetic interactions considered. For this reason, this effect can be safely omitted.

In the numerical experiments related to spin dynamics, a simplified form of the magnetic Hamiltonian was used in this work, i.e., only the Zeeman and exchange interactions were considered.

The model and numerical algorithms described in this paper are implemented using the LAMMPS software package[43] and corresponding calculation scripts. This software package was created by a team of authors from Sandia National Laboratories and is distributed under the GPL license, i.e., it is available free of charge in the form of source codes. The additional package LAMMPS SPIN allows numerical studies of magnetic systems and calculation of spin dynamics of atoms[23, 41]. The scripts developed by the authors for the LAMMPS software package were used in the particular calculations and data analysis.

The computational experiments performed in this work represent three separate series of calculations. In the first case, the processes of formation and structuring in a multilayer cobalt-iron-niobium nanocomposite were considered. The modeling was aimed at investigating the mechanisms of attachment and interaction of atoms in three-component layered systems. The second example of numerical experiments analyzed the mutual self-ordering of directions and reorientation of spins in crystalline iron in the presence and absence of a magnetic field. In this study, the mechanism of mutual arrangement of atomic spins and the manifestation of spontaneous magnetization typical for ferromagnets, of interest. The third numerical experiment was focused on modeling the magnetic properties of a cobalt-iron bilayer film under conditions of a uniform external magnetic field. In this case, attention was focused on the study of the magnetic properties of the nanocomposite depending both on the mutual arrangement of atoms of different materials and on external magnetic factors. The parameters of the MEAM interaction potential for various atoms are given in Table 1.

Table 1

\[ \]

The MEAM potential parameters for niobium, cobalt and iron. Here \( B^{(k)}_0 \) are the empirical potential parameters for a certain type of atom (\( k=0, 1, 2, \) and \( 3 \)), \( t^{(k)}_0 \) are the weighting coefficients of the model, \( \rho^{0} \) is the background electron density of the initial structure, \( A \) is the empirical potential parameter, \( E^{0} \) is the sublimation energy.

(1)

The first computational experiment addresses the problem of studying the processes of ferromagnet-supercon-ductor structure formation, for a three-layer system based on cobalt, iron, and niobium. Cobalt and iron act as a substrate, while the niobium layer is formed by deposition on a thin iron nanofilm. Niobium under certain conditions exhibits the properties of a superconductor, which is why it was chosen as one of the simulated materials. The relevance of the study of such systems is shown, for example, in Ref.[7]. The problem was solved by the method of classical molecular dynamics; the spin behavior of atoms was not considered in this simulation.

At the initial time, the system contained a 2 nm thick cobalt substrate (8000 atoms, \( 20 \times 20 \times 5 \) crystalline layers) and a 0.9 nm thick iron nanofilm (3200 atoms) located on it. The substrate and the film have a crystalline structure. The bottom layer of substrate atoms is partially fixed, preventing movement in the vertical direction. Partial fixation means that the atoms of the lower layer cannot move along the \( z \) axis, but they can move freely in the horizontal direction. This is done in order to eliminate the displacement of the simulated system during the deposition of atoms on the substrate. Similar restrictions are not imposed on the movement of other atoms. In the horizontal directions (along the \( x \) and \( y \) axes), periodic boundary conditions were applied in the system, and reflection boundary conditions were applied at the top and bottom. Before deposition of niobium (3000 atoms), the cobalt substrate and the iron nanofilm were at rest, and there were no external forces in the system. The deposited niobium atoms were given an initial velocity of 0.05 nm/ps. Subsequently, the niobium atoms lost velocity when approaching the substrate after interaction with its atoms. The value of the initial deposition velocity of niobium atoms was tested for adequacy by a multiple decrease. For this purpose, a similar problem of niobium nanofilm formation was solved, but the atoms were deposited at a velocity of 0.025 nm/ps. The structure and characteristics of the niobium nanofilm were the same. The deposition flux density was regulated by the number of iterations after which a new deposited atom was added to the system. For the current problem, this parameter corresponded to a value of 29 iterations. The total simulation time was 200 ps (200 thousand iterations). The initial and thermostat temperatures were set at 300 K. A Nose-Hoover thermostat with a damping parameter of 0.1 ps (100 time steps) was used. The thermalization performed during the modeling process affected only the substrate, the iron nanofilm and the already deposited niobium atoms. The atoms in the gas phase were not subjected to temperature control. An integration step of 1 fs was used in the first computational experiment. The general scheme of the computational experiment and the spatial orientation of the coordinate axes are illustrated in Figure 1..

Figure 1.

Schematic diagram of the modeled process of niobium nanofilm deposition on a bilayer nanocomposite of cobalt and iron.

The result of modeling the deposition of niobium atoms on a layered cobalt and iron base is demonstrated in Figure 2.. The uneven surface topography of the deposited niobium nanofilm is clearly visible. The nanofilm is formed roughly, with height differences of several angstroms. Such an effect, according to the authors, can be explained by more intense interaction forces arising between niobium atoms compared to other types of atoms considered. Visual analysis of the system shows that the Co-Fe layers have not undergone significant rearrangement, and the structure of these films predominantly remained crystalline. The structure of the deposited niobium nanofilm is difficult to judge visually, but there is no clear crystalline structure in it. Also, no deep penetration of niobium atoms inside Co-Fe was observed during the deposition process. This result is important because in some cases during sputtering, mixing of materials occurs and diffuse interfaces between nanofilms are formed. Diffuse interfaces can lead to disruption of the magnetization mechanisms of the composite, and thus subsequently degrade its basic functional properties.

Figure 2.

Appearance and structure of cobalt-iron-niobium nanocomposite obtained by modeling of niobium nanofilm deposition by the molecular dynamics method.

Quantitative analysis of the formed nanocomposite composition was carried out layer by layer. For this purpose, horizontal thin layers of the system with a thickness of 0.2 nm were calculated in the vertical direction and the quantitative fractions of the studied types of atoms in each layer were calculated. The calculation started from the bottom layer and ended with the surface of the deposited niobium nanofilm. The fractions of deposited atoms (in percentage) by layer are shown in Figure 3..

Figure 3.

Height distribution of cobalt-iron-niobium nanocomposite composition, obtained after deposition of niobium nanofilm.

The composition study presented in Figure 3. indicates a fairly clear separation of the different material nanofilms and confirms the results of the visual analysis of the atomic structure shown in Figure 2.. Nevertheless, the formation of a more diffuse contact zone is observed between the iron and niobium nanofilms compared to the contact between the cobalt and iron layers. The deposited niobium atoms have high kinetic energy, resulting in a partial introduction into the surface layers of the iron substrate.

To evaluate the structure of the atomistic material, we used the lattice centrosymmetry parameter calculated according to the following expression:

where \( Z \) is the number of nearest neighbors for the atom under consideration; \( {\bf r}^{\,}_i \) and \( {\bf r}^{\,}_{i+Z/2} \) are radius vectors of the analyzed and one of the nearby atoms; \( a^{\,}_{\rm lat} \) is the lattice constant. In expression (14) the whole set of possible nearest neighbors of an atom is searched and the square of the distance between them is calculated in pairs. The obtained average value of the parameter characterizes the overall deviation of the nanomaterial structure from the ideal crystal structure.

In general, the behavior of the lattice centrosymmetry parameter depends on many factors, including temperature, since thermal fluctuations affect the coordinate deviations of atoms from their ideal positions. Nevertheless, for solid crystalline materials the average value of this parameter is small, while for amorphous materials it has a large positive value. For the investigated cobalt-iron-niobium nanocomposite, the distribution of the lattice centrosymmetry parameter for different axial projections is illustrated in Figure 4..

Figure 4.

Distribution of the lattice centrosymmetry parameter for cobalt-iron-niobium composite atoms at the final time after deposition for axial projections: \( \bf a \) – the \( yoz \) plane (side view), \( \bf b \) – the \( yox \) plane (top view) for \( z=1.2\, \)nm (white dashed line in panel a).

Analysis of the formed nanocomposite atomic structure in Figure 4. shows that the smallest value of the parameter \( C^{\,}_{\rm sim} \) is found in surface niobium atoms. This effect is explained by the incomplete set of nearest neighbors of these atoms, so the lattice centrosymmetry parameter cannot be used to study the surface structure of nanomaterials. Cobalt and iron nanofilms have a small value of the lattice centrosymmetry parameter, which indicates that their structure is close to crystalline. The deposited niobium nanofilm is characterized by higher values of \( C^{\,}_{\rm sim} \). Its non-ideal structure is clearly visible in Figure 4.. The niobium layer undergoes significant rearrangement and subsequent slight compaction during deposition due to transformation processes occurring between atoms. Nevertheless, the final structure of niobium also exhibits lattice distortions. The dependence of the magnetic parameters of the sample on its atomic structure is shown, for example, in Ref.[44]. Defects in the structure and local arrangement of atoms arising during the deposition of nanofilms can cause deterioration of macroscopic magnetic parameters such as the magnetization modulus, magnetic permeability and its temperature coefficient, saturation induction, and others.

The averaged value of the lattice centrosymmetry parameter during the modeling process varied in the range from 0 to 15. Such a large value is primarily due to the influence of chaotically arranged deposited niobium atoms, which made a significant contribution to the calculated value. The niobium atoms appeared randomly in the deposition zone during the formation of the nanofilm, so their structure was far from the ideal crystalline structure. The value of the lattice centrosymmetry parameter of the deposited atoms subsequently decreased.

As part of the study of nanofilm deposition processes, computational experiments were performed in this work, in which the size of the systems was multiplied both in horizontal directions and with respect to the number of sputtered atoms. In all cases, similar results were obtained for the structure and composition of the formed nanocomposite. Thus, the influence of boundary effects on the properties of the formed three-layer sample was excluded. In case of increasing only the number of deposited niobium atoms, the thickness of the final nanofilm increased, but its structure and uneven surface structure remained.

Modeling of nanomaterials deposition processes in this work was considered for physical vapor deposition (PVD) technology. This technology includes a rather large group of methods for obtaining materials (thermal sputtering, molecular beam epitaxy, magnetron sputtering, laser beam or vacuum arc evaporation, focused ion beam heating) and is based on the transformation of the deposited material into the gas phase. The sputtering process by this technology, typically takes place in a vacuum environment. At this stage of research, the method of atom deposition by directed flow, characteristic, for example, of molecular beam epitaxy technology, was considered in this work. In the future, it is planned to develop and expand the study of PVD technology to deposition methods using electromagnetic fields, including consideration of the method of magnetron sputtering.

As noted earlier, the second computational experiment analyzed the mutual self-ordering of directions and reorientation of spins in crystalline iron. The modeling considered a nanosystem consisting of 20 elementary iron cells in each horizontal direction and 3 crystalline cells in the vertical direction, resulting in about 3200 atoms. The size of the system in nanometers was \( 5.7 \times 5.7 \times 0.9 \). The integration step of 0.1 fs was used in the simulation. The boundary conditions in all directions were periodic. The direction of the spatial axes coincided with the symmetry axes of the iron crystal lattice. The total simulation time was 10 psec or 100 thousand iterations. This stage of modeling is necessary to evaluate the correctness of the model and selected parameters of the computational experiment, as well as to investigate the influence of the external magnetic field on the considered system. In the second computational experiment, the temperature was maintained at 8 K using a Langevin thermostat. The thermostat damping parameter was 0.1 ps.

The second computational experiment included two parts: consideration of the crystalline iron nanosystem both in the absence and in the presence of an external magnetic field. As a result of the modeling, the spatial distribution of the iron atoms' spins during the entire time of the study was obtained. The direction of the spin vectors at the initial time was set random to minimize the probability of its influence on the final distribution of the spin vectors.

Figure 5.

Direction of crystalline iron spin vectors for the time instant of 10 ps in the absence of an external magnetic field (panels a and c) and in the presence of an external magnetic field oriented along the nanofilm surface (panels b and d).

In the process of modeling, as shown in Figure 6., the reorientation of atoms' spins is observed both in the absence of an external magnetic field and in its presence. Figure 5.a and Figure 5.c show the formation of vortex structures (skyrmions), \( i.e. \) regions of spontaneous homogeneous magnetization, in which the magnetic moments of Fe atoms are co-directed, although no magnetic field was applied to the structure. Skyrmions are vortex magnetic textures that can be viewed as topological objects[45]. They are quasiparticles that can move in a magnetic system, similar to electrons in a conductor. Skyrmions are found both in thin magnetic nanofilms[46] and in superconductors[47]. Materials in which skyrmions arise are promising for spintronics devices, including as an element of track memory[48] — the next generation of magnetic memory.

Figure 6.

Time variations of the absolute value of magnetization (panel a) and the magnetic energy of the crystalline iron system (panel b) during the simulation in the absence and in the presence of the external magnetic field of 0.1 T. In the LAMMPS software package, the unit for magnetization is \( e/({\mathring{A}\cdot ps}) \), which is equal to \( 1.60\cdot 10^3 \) A/m (here \( e=1.60\cdot 10^{-19} \) C is the elementary charge)

The coordinated reorientation of atomic spin vectors indicates the occurrence of spontaneous magnetization in iron. Such behavior of spins means that the material has ferromagnetic properties. The manifestation of ferromagnetic properties of iron is a known fact, which confirms the adequacy of the considered model, as well as the used parameters of potentials and magnetic interactions. Figure 5.b and Figure 5.d show the distribution of magnetic moments in the system in an external magnetic field. The direction of the magnetic induction vector for this system coincided with the direction of the ox axis. The magnitude of the magnetic induction vector in the simulation was 0.1 T. Figure 5. shows that the application of an external magnetic field resulted in the displacement of the magnetic vortex regions. In the system with an applied external magnetic field, as well as in the system in the absence of a magnetic field, the formation of skyrmions is observed.

The overall magnetic behavior of the material can be analyzed using the system magnetization vector and magnetic energy. The variation of these parameters during the simulation is shown in Figure 6.. The magnetization norm in the Euclidean space is the modulus of the magnetization vector: \( |M|=\sqrt{M_x^2+M_y^2+M^2_z} \).

As can be seen from Figure 6.b, the total magnetic energy of the system during the computational experiment stabilizes around 280 eV both in the presence of an external magnetic field and without it. The abrupt change in the magnetic energy at the initial stage of the calculations is due to the nonequilibrium initial state of the system under consideration, caused by the random distribution of spin vectors at the initial moment of time.

A more detailed behavior of the magnetization vector of the iron nanofilm is illustrated in Figure 7., where the dynamics of the individual components of this vector is presented. The behavior of the magnetization vector projections differed significantly between the case with an applied magnetic field and the case without it.

Figure 7.

Components of the magnetization vector along the \( x \)-axis (panel a), the \( y \)-axis (panel b), and the \( z \)-axis (panel c) in the absence/presence an external magnetic field

The graphs in Figure 7. show that the magnetic moments of atoms are reoriented under the influence of an external magnetic field. Since the magnitude of the magnetization vector components along the \( y \)- and \( z \)-axes approaches zero (Figure 7.b and Figure 7.c), we can conclude that the spin vectors of atoms are reversed along the \( y \)-axis. The revealed formation of skyrmions and their behavior under the influence of a magnetic field allows us to discuss the possibilities of promising use of crystalline iron nanofilms for memory devices functioning on the basis of the principle of controlled displacement of vortex magnetic regions under the influence of an external magnetic field.

Despite the fact that the anisotropy energy is not introduced as a separate contribution in the mathematical model, this model allows describing some asymmetric magnetic effects. Thus, in the second problem, in the absence of an external magnetic field, the general magnetization vector of the iron nanofilm has the direction \( {\bf M}\Big|_{{\bf B}_{\rm ext}=0} = (M^{\,}_x, M^{\,}_y, M^{\,}_z) = (13, 4, -5.5) \). The asymmetric location of the magnetization vector of the sample relative to the coordinate axes is explained by the fact that the paper considers a thin iron film, in which the formation of magnetic vortices of the spin vectors of atoms is observed. When considering iron nanofilms of different thicknesses, the direction and norm of the magnetization vector can change. In computational experiments studying bulk iron samples, magnetic vortices are not observed, the spin vectors of atoms are oriented predominantly in one direction, and the magnetization vector reaches a certain saturation. As can be seen from Figure 7., under conditions of an external magnetic field \( {\bf B}_{\rm ext}=0.1\, \textrm{T} \), the general magnetization vector of the iron nanofilm changes its direction to \( {\bf M}\Big|_{{\bf B}_{\rm ext}\neq 0} = (6.5, 1.5, 1) \). In this case, the magnetization vector is predominantly directed along the ox axis, that is, in the direction of the external magnetic field. This was not observed in computational experiments without an external magnetic field. The decrease in the magnetization norm of the sample under conditions of an external magnetic field \( {\bf B}^{\,}_{\rm ext}=0.1\, \textrm{T} \), compared to the magnetization norm without an external field, is associated with a decrease in the absolute value of the components of the magnetization vector \( M^{\,}_x \) and \( M^{\,}_y \). These are the directions that lie not along the vector of the external magnetic field, but are perpendicular to it.

Reproducibility of the obtained results of computational experiments is due to the fact that at any initial distributions of velocities and directions of atomic spins the system comes to a single equilibrium physical state. To confirm this fact, additional computational experiments with alternative distributions of initial velocities and spin vectors of atoms were carried out.

The third computational experiment focused on investigating the magnetic properties of a layered nanocomposite of cobalt (8000 atoms) and iron (3200 atoms) under conditions of a constant magnetic field. The total size of the system in nanometers was \( 5.7 \times 5.7 \times 2.9 \), where the height of the cobalt nanofilm was 2.0 nm, and the height of the iron nanofilm was 0.9 nm. The integration step of 0.1 fs was used in the simulation. The direction of the spatial axes coincided with the symmetry axes of the crystal lattice of the simulated materials. The total simulation time was 10 ps or 100 thousand iterations. Such a system was considered in the first calculation as a basis for the deposition of the top niobium nanofilm.

Figure 8.

Distribution of atomic magnetic moments in thin film Co-Fe nanocomposite for a simulation time of 10 ps.

The magnetic induction vector \( {\bf B}^{\,}_{\rm ext}=0.1\, \textrm{T} \) was directed along the ox axis as shown in Figure 8.. The materials had a crystalline structure and the MEAM potential was used for the interaction between atoms. Along the horizontal directions, the computational cell had periodic boundary conditions, and in the oz axis direction, reflection boundary conditions were in effect. The initial velocity and spin vectors of atoms were set randomly in accordance with the initial lattice and spin temperatures \( T^{\,}_{\rm lat}=T^{\,}_s=8\, \textrm{K} \), which coincided with the temperature values from the previous numerical experiment. To maintain the lattice and spin temperatures, a Langevin thermostat was used, the operating principle of which is described in Eqs(15)-(11). The damping parameters of the thermostats were the same at a level of 0.1 ps.

During the simulation, the lattice and spin temperatures were kept at the same level by using the Langevin dynamics. This value of temperatures was chosen specifically below the superconducting transition temperature of niobium (9.25 K). It is known that in superconductors, including niobium, the phenomenon of complete or partial expulsion of the magnetic field from the bulk due to the Meissner effect occurs during the transition to the superconducting state[49]. For this reason, the niobium layer was not considered in the calculations performed.

Figure 9.

Formation of magnetic domains in a cobalt nanofilm at a simulation time of 10 ps.

The result of the spin dynamics of the two-component Co-Fe systems at a simulation time of 10 ps is shown in Figure 8.. At the initial time, the magnetic moments of atoms in both the cobalt and iron nanofilms were differently oriented. Subsequently, a joint change in the orientation of the atoms' spins was observed. The Fe layer was more prone to forming vortex magnetic regions and skyrmions, which were previously obtained by modeling a single layer of crystalline iron. The magnetic behavior of the cobalt nanofilm differed from the pattern of mutual orientation of spins in iron. Sufficiently well-defined magnetic domain regions were obtained, which are easily identified in Figure 9..

The magnetic domain zones shown in Figure 9. and selected by different geometrical shapes have different spatial orientation of spin vectors from each other. At the same time, within the selected domains, a consistent uniform orientation of atomic spins is observed. The shape of the domains differs. At the junction of magnetic domains, the spin vectors of atoms are rotated. However, the overall magnetization of the system is low due to the absence of a distinct priority direction of magnetic moments. The zone structure of cobalt nanofilm magnetism can be used in the creation of new film structures and magnetic nanoobjects in superdense recording and information storage devices.

Figure 10.

Time evolution of the spin and lattice temperatures for Co-Fe nanocomposite.

The average dynamics of atomic motion and changes in their spins can be estimated by calculating the temperatures of the system. The change of lattice (\( T^{\,}_{\rm lat} \)) and spin (\( T^{\,}_s \)) temperatures of the investigated cobalt-iron composite in the process of modeling is illustrated in Figure 10.. Analysis of the temperature graph shows that significant changes and jumps of these parameters are observed at the initial moments of time. Such an effect is explained by stochastic initial distributions of velocities and magnetic moments of atoms. Subsequently, the velocities and directions of spins, which were unstable in the initial state, are rearranged, and fluctuations of the values become moderate. The temperature dynamics reaches stationary regimes corresponding to thermostat values of 8 K. Minor fluctuations of temperatures near the target value indicate that the composite is in an energy stable state, and the lattice and spin thermostats function adequately in the system.

Figure 11.

Variation of the magnetization norm in layers of cobalt, iron of different thicknesses and nanocomposite as a whole.

For the layered cobalt-iron nanocomposite, the change in the magnetization norm was calculated (Figure 11.). The magnetization rate was determined both for the composite as a whole [labeled Co+Fe (\( d \)=1 nm)] and separately for cobalt [labeled Co] and iron [labeled Fe (\( d \)=1 nm)]. Additionally, the magnetic behavior of the same system but with a two-fold increased thickness of the iron nanofilm [labeled Fe (\( d \)=2 nm)] was investigated. The magnetic moments of the domain structure arising in cobalt are multidirectional, which causes the magnetization norm value of this nanofilm to be close to zero. As can be seen from Figure 11., the domains are built rather quickly (within the first picosecond of modeling), and further the magnetization norm of cobalt changes insignificantly.

The dynamics of the magnetization norm of iron nanofilms has a more variable character. This is due to the fact that the vortex orientation of the atomic spins requires a longer time to occur. In addition, already after the formation of skyrmions in iron, some displacement can occur, which also affects the change of the magnetization norm. The skyrmions in iron have a well-defined magnetic moment, which leads to larger values of the magnetization rate of iron, compared to the cobalt nanofilm. This effect is clearly visible in Figure 11.. In the variant of the computational experiment with increased Fe thickness, vortex structures appear not only in the plane of the nanofilm, but also unfolded in the volume, which causes reorientation of the magnetic moments of atoms and, as a consequence, an increase in the magnetization norm of the material.

The use of skyrmions, obtained from the modeling in the second and third problems, and Josephson junctions[50] is a very effective direction in the creation of fast and energy-efficient memory devices, as well as in the development of superconducting qubits and quantum circuits focused on next-generation of quantum CPUs[51]. Such nanostructures can also be used as tunable kinetic inductors, which are designed to implement and control artificial neural networks with magnetic data representation[52]. However, the creation of nanoscale thin-film multilayer materials and precise control of their magnetic states requires a thorough elaboration of their fabrication technologies, functioning processes, as well as expanding the understanding of the fundamental properties of nanoobjects.

Numerical study of niobium deposition processes on a bilayer composite of cobalt and iron has shown that cobalt and iron nanofilms have small values of the lattice centrosymmetry parameter. This fact indicates that their structure is close to the crystalline structure. The deposited niobium nanofilm is characterized by higher values of the lattice centrosymmetry parameter, non-ideal structure and uneven surface. A layer-by-layer study of the composition shows a fairly clear separation of nanofilms of different materials, but between the layers of iron and niobium the formation of a more blurred contact zone is observed. The deposited niobium atoms have high kinetic energy, which results in their partial introduction into the surface of the iron, on which the sputtering is carried out.

In the computational experiment investigating the magnetic behavior of iron crystalline nanofilm atoms under conditions of external field and its absence, it was obtained that in both considered variants the reorientation of spins was observed, the formation of spontaneous homogeneous magnetization regions and vortex currents occurred. The coordinated reorientation of spin vectors indicates the occurrence of spontaneous magnetization in iron and the manifestation of ferromagnetic properties. The presence of an external magnetic field led to the displacement of magnetic vortex regions, the reversal of atomic spins, and a decrease in the norm of the total magnetization of the nanomaterial.

As a result of modeling the magnetic properties of a layered nanocomposite of cobalt and iron under conditions of a constant magnetic field, the formation of skyrmions in the iron nanofilm and domain regions in cobalt was observed. The shape of the domains differed. The magnetic moments of the domain structure appearing in cobalt were multidirectional, which caused the magnetization norm value of this nanofilm to be close to zero. The iron layer in the nanocomposite had larger values of the magnetization norm compared to cobalt, since the formed skyrmions had a distinct magnetic moment. In a variant of the computational experiment with increased iron thickness, vortex structures were obtained not only in the plane of the nanofilm, but also deployed in the volume.

The work was carried out within the framework of the state assignment of the Ministry of Science and Higher Education of the Russian Federation (topic No. FUUE-2025-0001 Modeling and experimental studies of the structure and magnetic parameters of a Josephson contact made of a multilayer nanostructure, state registration No. 125020501412-9) (methodology for modeling magnetic properties) and with the financial support of the Ministry of Science and Higher Education, project No. 075-15-2025-010 (research at ultra-low temperatures) and the project No. 020201 of the Moldova State Program Nanostructures and advanced materials for implementation in spintronics, thermoelectricity and optoelectronics (samples preparation and investigation of Nb/Co interfaces).

Corresponding author: Aleksey Fedotov

orcid.org/0000-0002-0463-3089

alezfed@gmail.com