Autonomous synthesis system integrating theoretical, informatics, and experimental approaches for large-magnetic-anisotropy materials

ABSTRACT We developed an autonomous and efficient system for synthesising ferromagnetic materials with large magnetocrystalline anisotropy by integrating theoretical, informatics, and experimental approaches. By combining the first-principles calculation of the magnetic anisotropy with Bayesian optimisation, we virtually screened candidate materials, comprising four elements and four-layer periods, from various magnetic multilayers. We employed the expected improvement as the acquisition function and Matern52 as the kernel function, to develop a robust machine learning model. We fabricated the top three predicted magnetic materials under laboratory conditions by monoatomic layer deposition and evaluated their magnetic anisotropy using a superconducting quantum interference device (SQUID). Ultimately, we demonstrated that [Fe/Co/Fe/Ni]13 is a novel ferromagnetic material whose magnetic anisotropy exceeds that of L10-FeNi- and L10-FeCo-type alloys. Furthermore, the origin of the perpendicular magnetic anisotropy was derived from the spin-conserving as well as the spin-flip terms. We determined that Bayesian optimisation offers promising configurability features in terms of the electronic structure that extend beyond the empirical knowledge and human intuition. Graphical abstract


Introduction
Novel magnetic materials can be synthesised using several element combinations, compositional ratios, multilayer structures, and growth temperatures. Additionally, these materials show complex interactions between their crystal structure and electronic state, which complicates the development of new materials with large magnetocrystalline anisotropy (MCA) for next-generation IoT devices. There is considerable scope for designing MCA materials by devising new synthetic strategies and discovering novel materials based on the literature. To date, materials science researchers have explored new material candidates via trial-and-error experiments through a series of calculations, synthetic methods, and material analyses. However, this approach consumes much time and requires high-skilled specialists. Recent developments in informatics and robotics can significantly aid in this effort. Data-driven materials design has been introduced in various fields, including the high-throughput synthesis of functional chemical materials and synthesis of photocatalytic materials by autonomous robots using Bayesian optimisation algorithms [1,2]. In this study, we introduced a state-of-the -art methodology combining theoretical, informatics, and experimental approaches for the autonomous synthesis of novel magnetic materials. In particular, we have developed an integrated materials synthesis system that predicts highly functional magnetic materials based on data-driven decisions and integrates experimental material synthesis.
The atomistic design approach is invaluable for the synthesis of novel magnetic materials. In particular, the MCA can vary greatly depending on the elemental composition/ratio and stacking order and is determined by the complex interactions between the crystal structure, magnetic moment, and electronic structure. Conventionally, materials scientists investigate novel materials by manually selecting target materials, performing theoretical calculations and evaluating the material functions, and individually adjusting the experimental conditions. However, in multi-element and multilayer systems, the number of candidate materials increases substantially, owing to which conventional material investigation methods present numerous limitations. Scientists in the field of magnetic multilayers have employed theoretical calculations and empirical experience to design the elements and stacking order on an atomic scale and have spent considerable effort in preparing evaluating individual materials. However, complex material systems result in combinatorial explosions, which significantly increase the effort requirement and hinder the discovery of novel materials. Currently, fabrication and characterisation of new materials are becoming more efficient through the introduction of data science. However, the theoretical calculation of magnetic anisotropy requires a highly accurate measurement of the electronic state of the configured structure, which is time-consuming and cost-ineffective even if a powerful computer is available. Similarly, the material fabrication process requires structural control at the atomic scale. The functional evaluation process is labour-intensive; individual magnetic hysteresis measurement is required for each specimen. Thus, there is great demand for the establishment of an innovative material synthesis technique to bypass such inefficient research processes. Spintronic devices, an example application of magnetic multilayers, are energy-saving and high-density memory devices that can contribute to the information society in the near future. In recent years, the development of magnetoresistive random access memory has progressed rapidly, and there is high demand for the development of novel ferromagnetic materials with large MCA [3]. To date, L1 0 -type ordered alloys with 2-element, 2-layer periods have been developed; L1 0 -FePt is a typical example, and recently, L1 0 -FeNi has been widely investigated theoretically and experimentally [4][5][6][7]. At present, L1 0type magnetic materials are based on two elements, and few studies have been conducted to extend the number of elements and their periodicity. Despite the ongoing investigation on ternary multilayer materials, the candidates explored thus far are few, leaving considerable space for further investigation [8][9][10].
In this study, we developed an autonomous and efficient strategy to synthesise magnetic materials with large MCA by integrating theoretical, informatics, and experimental approaches. The combination of materials science with machine learning is a data-driven approach for finding optimal solutions about big data and large-scale materials [11][12][13][14][15]. Recently, polymers, photocatalysts, and thermoelectric materials have been widely explored for application to various material fields [16][17][18]. Currently, there is great demand for an autonomous, machinelearning-robotics integrated material development method for predicting functional ferromagnetic materials [19]. In the case of magnetic materials, quantitative extraction of the factors governing the MCA behaviour is impractical because of the complex interactions between the lattice distortions, stacking order, number of multilayers, and electronic structure [20][21][22]. In addition, there are numerous possible combinations of constituent elements, rendering the candidate material selection process a bottleneck in the development of magnetic materials. Here, we employed the Bayesian optimisation method, which updates the learning model by sequentially predicting the output magnetic anisotropy from the input crystal structure. Bayesian optimisation is a prominent blackbox optimisation approach, which offers significant advantages in selecting ideal materials from a large sample space [23]. In this approach, Gaussian process regression predicts the shape of the function by calculating its variance and mean values from the training data. Acquisition function that controls the balance between exploitation and exploration is calculated to determine promising search points for the regression model. The new predicted data are then added to the training data to update the prediction model. By continually repeating this cycle, the global optimum can be obtained with a small number of iterations. This method is useful for research subjects whose cost of data collection is large or relationship between the explanatory and objective variables is unknown. Bayesian optimisation has been used to predict the microstructural heat transfer mechanism in thermal radiation multilayers and solar cells [24][25][26][27]. However, it has not yet been applied to predict the effects of complex interactions on MCA; hence, the application could present interesting research opportunities. By setting the MCA as the output of the Bayesian optimisation and the material configuration as the variable to be optimised, the target material configuration with a large MCA can be determined autonomously.

Workflow of material exploration by integrating theoretical, informatics, and experimental approaches
We combined the first-principles calculation of the MCA with Bayesian optimisation to determine the optimal candidate materials. Based on the predictions, we prepared nanoscale magnetic multilayers under experimental conditions and analysed their ferromagnetic properties to develop an actual material. We adopted four-layer periodic multilayers, each consisting of a non-rare-earth element (Fe, Co, Ni, and Cu), as a model case for experimental deposition. Furthermore, the applicability and performance of the proposed method with respect to determining perpendicular MCA materials was evaluated from theoretical and experimental viewpoints. Finally, we discussed the mechanisms underlying the improvements in the magnetic anisotropy of the fabricated materials from electronic-structure viewpoint. The workflow of this study is presented in Figure 1.
We initially generated the input data for the MCA energy calculation using the first-principles calculation method ( Figure 1). We then predicted the configuration of the multilayer with the highest improved MCA energy value using the Bayesian optimisation method ( Figure 1). The structural information (stacking order of elements) and MCA energy were the explanatory and objective variables, respectively. The information of the predicted multilayer composition was input again in step (a) to calculate the MCA energy. Steps (a) and (b) were repeated for each iteration to update the machine learning model and search for materials with large MCA energies. Furthermore, we examined several kernel functions that determine the similarities among the data to improve the search efficiency and optimise the machine learning model. To experimentally fabricate the predicted materials, it is important to validate the prediction performance of not only the first candidate but also the lower candidates. We then experimentally examined the feasibility of this method by fabricating and characterising the materials (Figure 1). The predicted magnetic multilayers were prepared by pulsed laser deposition (PLD). The deposition rate of the PLD system was regularly stabilised by automatically controlling the laser power. Accordingly, the stabilisation of the ablation process allowed for the deposition of a smooth thin film. A wide variety of multilayers can be fabricated by employing a programmable automatic target exchange and an automatic shutter open/close. The entire PLD system presented here is automated and programmatically controlled to achieve precise film growth while reducing human effort. Since two of the six PLD targets were used for substrates and caps, utmost four elements were set as the components of the multilayer. Although the area explored in this study is not vast, there has been no previous Bayesian exploration of the complex physical quantity 'MCA'; moreover, few studies have combined Bayesian optimisation, first-principles calculations, and experimentation. The development of fundamental principles and feasibility analyses in this study could be significant in exploring novel magnetic materials. Regarding the material properties, the surface structure was investigated by atomic force microscopy (AFM), whereas the crystal structure was analysed by X-ray diffraction (XRD). A superconducting quantum interference device (SQUID) magnetometer was used to evaluate the MCA energy. Thus, an autonomous magnetic material synthesis system combining theoretical, informatics, and experimental approaches was developed, and its feasibility was demonstrated. To comprehensively illustrate the mechanism underlying the improvement in the magnetic anisotropy, we analysed the electronic structure in detail. We evaluated the anisotropy of the orbital moments of each element to determine the contribution of the discovered samples to the MCA energy. We also incorporated the secondorder perturbation term to analyse the projected density of states (PDOS) near the Fermi level and related the origin of the spin-orbit interaction to the improvement in the MCA energy.

Magnetic anisotropy energy calculation by first-principles calculations
The electronic structures of the candidate materials were determined via first-principles calculations based on density functional theory (DFT) using Quantum Espresso (QE), a first-principles calculation package [28]. The MCA energy calculations were performed after structural relaxation to determine the mechanically stable structure for a given film configuration ( Figure 1). Each unit cell had a tetragonal structure with the period of four layers. The results obtained herein for three known materials were compared with those of previous studies, and the Perdew-Becke-Ernzerhof (PBA) generalised gradient approximation (GGA) was applied as the exchangecorrelation function for the structural relaxation, with a k-point of 12 × 12 × 8 [29]. Table 1 shows that the structural relaxation with the GGA for FeNi and CoNi is in good agreement with that for WIEN2K [30]. For FeCo, the B2 structure, rather than the L1 0 type, is the stable structure because it is the equilibrium state of the FeCo system (c/a = 0.707). The MCA energy was defined as the difference between the total energy when the magnetisation was oriented in-plane [100] and out-of-plane [001], and it was considered positive when the magnetisation easy axis was perpendicular to the plane. The results of appropriate exchange-correlation functions (GGA and local density function (LDA)) for MCA calculations are listed in Table 2 [31]. In Table 2, results from previous studies using all-electron calculation, SPR-KKR, WIEN2K, FP-LMTO, and two types of calculations using QE are presented [30,32]. GGA tends to underestimate the MCA energies, whereas calculations using QE LDA agree well with those in previous studies. In this study, we adopted LDA as the exchange-correlation function for MCA calculations with a k-point of 24 × 24 × 18. The tetrahedron method was used to efficiently perform the Brillouin zone integration [33]. The projector augmented wave (PAW) method was used to describe the core and valence electron potentials in these calculations [34]. The cut-off energy of the plane-wave basis was set to 75 Ry.

Prediction of large MCA energy materials by Bayesian optimisation method
Bayesian optimisation is a sequential strategy for global optimisation. In this study, the function model f of the MCA energy y, and the material configuration x is constructed. Bayesian optimisation approaches have two important steps: the definition of a prior probability and acquisition function. Regarding the first step, from a Gaussian process, the prior distribution of f(x) is described by the following parameters [35].
where the Gaussian process has two parameters: the mean value μ x ð Þ and the kernel function k x; x 0 ð Þ. The kernel function k x; x 0 ð Þ characterises the smoothness and other properties of the functional f. The prior distribution is updated to a posterior distribution using the available data. A relationship between these probability terms is given by Bayes' Theorem in Equation (2).
where P f jx n ð Þ is the posterior distribution, P x n jf ð Þ is the likelihood function, and P f ð Þ is the prior distribution. In this study, the elements of each layer were defined as simple descriptors: Fe={1}, Co={2}, Ni={3}, and Cu={4}. Such material descriptors have already found various applications, from atomic-scale crystal structures to multilayers [24][25][26]. The number of elements in the machine learning vector is defined by the periodicity of the specimen and is independent of the number of elements. This reduces computational costs and is suitable for material exploration in small  laboratories. Further, using continuous variables signifies the value itself; however, the elements are indexed by atomic number, which functions as a reasonable descriptor in terms of the number of electrons. Because the periodicity of the thin film was set to four layers, the explanatory variables indicating the film composition are represented by a fourdimensional array. Four initial training data were randomly selected, and the variance, σ 2 , and mean, μ, required for estimating the functional form were calculated using Equation (3): where k X 1:t ; x ð Þ and k x; X 1:t ð Þ are the vectors of the kernel functions for some value of x and the existing data, and k X 1:t ; X 1:t ð Þ contains all previously evaluated combinations of values. The estimation of the functional form depends on the kernel function that determines the similarity between the data. An efficient material search cannot be performed without selecting an appropriate kernel function. Therefore, in this study, four kernel functions were used to compare the search efficiency: exponential, Matern32, Matern52, and the radial basis function (RBF) [35].
The second step is to define the acquisition function. The calculated expectation and variance were used to compute an acquisition function that determines the degree of improvement from the existing data set. The acquisition function used here was the expected improvement (EI), shown in Eq. (4), which considers not only the probability of improvement but also the amount of improvement, allowing for an automatic balance between exploitation and exploration [36].
where ϕ, Φ, and Z represent the probability density function, cumulative distribution function, and random variables of the normal distribution, respectively. By calculating this function, we can predict the optimal material configurations. Four initial training data points were selected by random sampling (Figure 1). Bayesian optimisation was applied to this data set to determine the next search point, candidate materials. Using the predicted structure, the structural relaxation and MCA energy were calculated by first-principles calculations, and Bayesian updating was performed to provide new training data. Finally, we created an environment in which the above-mentioned processes were carried out automatically until all material candidates were proposed. The python implementation of the Bayesian optimisation method uses the GPyOpt package [37]. The total computation time for structural relaxation and MCA calculations, for all candidates, was approximately 28 d. Structural relaxation was carried out multiple times, with different initial structures, to determine the stable structure. The calculations were performed using 30 threads of parallel computing (Processor: AMD® Ryzen Threadripper 1950 × 16-core processor × 32, Memory size: 64 GB).

Material search results using Bayesian optimisation
A robust material search method is used to discover promising materials efficiently, independent of the type of initial training data. This is because the initial training data for machine learning are selected by either empirical artificial sampling or random sampling. In this study, empirical methods are inappropriate for sampling the initial training data because manual prediction of the parameters of the fourdimensional variables is difficult. Here, we use random sampling to input four initial training data and perform a material search using Bayesian optimisation. Because the appropriate kernel function for the MCA data set of ordered alloys is unknown, we optimised the kernel function by comparing the prediction accuracy of Bayesian estimation for four typical kernel functions (exponential, Matern32, Matern52, and RBF). We also compared the results of randomly selected candidates with the results of the Bayesian optimisation. The search region had a maximum of four elements and four layers, so the total number of candidates was 256 (4 4 ); however, by considering the c-axis inversion and translation symmetries, duplicates were eliminated and the number of candidates was reduced to 55. We prepared 100 sets of initial training data and performed Bayesian optimisation predictions for each of them ( Figure 2). The horizontal axis shows the number of iterations of the Bayesian optimisation, and the vertical axis shows the maximum MCA energy at that time. All candidates in Figure 2(a) were selected by random sampling, whereas Figure 2(b-e) show the search results after Bayesian optimisation, wherein only the grey areas were selected as the initial training data by random sampling. Compared with random sampling, Bayesian optimisation tended to find the maximum MCA energy in less than half the number of iterations for any kernel, regardless of the type of initial data. Region (I) in Figure 3 illustrates that the maximum MCA energy was discovered at an earlier stage. In region (II), the average MCA value was relatively high because the search was focussed on similar structures as those of the material found in region (I). The average MCA value decreased significantly in region (III), whereas the average value tended to increase in region (IV). The reason for this behaviour was ascribed to the discovery of zero anisotropy, rather than to the prediction of large anisotropy. This trend with the increasing and decreasing average MCA energy value in Figure 3 was confirmed for all candidate materials and kernels (Figure 3)). This search model was demonstrated to find materials with large MCA values in the first half, and materials with negative or nearzero anisotropic values in the second half. These results are expected to vary significantly depending on the selection of the acquisition function, which determines the balance between exploitation and exploration. To illustrate the search process, we focussed on the number of findings regarding the top two materials in terms of anisotropy (Figure 3)). The projected material candidates with large MCA values were localised in region (I) in Figure 3 for all kernel functions, indicating that the finding of the largest MCA material once will lead to a series of materials with similar composition in any kernels. Assuming that the best material, as inferred from the result, cannot be fabricated experimentally, it is useful to develop a model that suggests multiple material candidates, including second and third candidates. Table 3 lists the average number of iterations required to find the best and second-best candidates in terms of anisotropy. As presented in Table 3, compared to random sampling techniques, Bayesian optimisation can find materials with large anisotropy even in notably earlier iterations. This indicates that Bayesian optimisation can successfully establish the relationship between the film structure and the MCA of the ordered alloy, even if this relationship is a black box. In particular, here, Bayesian optimisation using Matern kernels could find a material composition with a large MCA energy in less than 20% of the total number of iterations. Promising materials can be discovered in approximately 6 d with a total calculation time of 28 d. This represents a fivefold time-saving for material search in small laboratories with limited resources. Moreover, the difference in the results between the Matern32 and Matern52 functions was negligible because both functions use discrete values to define the film features. By using discrete instead of continuous values as features of the film composition, the efficiency of the material search becomes independent of the kernel function; consequently, a highly efficient and stable material search can be performed. By combining these methods, we successfully developed a robust magnetic material exploration model.
Our predictions using theoretical calculations and machine learning showed that the compositions of the materials with the largest and second-largest MCA values were (Fe/Cu/Fe/Cu) and (Fe/Cu/Co/Cu), respectively. The non-magnetic element, Cu, constituted half of the composition of both materials. Because of the lack of magnetic moment, intuitively predicting the contribution of non-magnetic elements to the improvement of anisotropy is difficult. In addition, the non-magnetic elements were sandwiched between the ferromagnetic elements in the structure. The material with the third-largest MCA energy was (Fe/Co/Fe/Ni), which has a thermodynamically stable structure with large MCA energy and saturation magnetisation values, as per theory [38]. However, the top two materials ((Fe/Cu/Fe/Cu) and (Fe/Cu/Co/Cu)) have not been previously reported either experimentally or theoretically. Our proposed model could effectively search for new materials without being restricted by empirical rules. Furthermore, MCA strongly depends on lattice distortion [22,[39][40][41]. The relationship between lattice distortion, c/a, and MCA for all candidates is shown in Figure 4. From Figure 4, the face-centred tetragonal (fct) structure (c/a ≧ .8) promotes the appearance of perpendicular MCA above 1.00 MJ/m 3 . In particular, FeCo binary ordered alloys have a stable structure with a compressed body-centred tetragonal (bct) structure (c/a < 2.0) in the c-axis direction, regardless of the film configuration; this was disadvantageous for perpendicular MCA. Screening candidates by eliminating bct-structured materials from the structural relaxation results would improve the efficiency of the machine learning process.  Label descriptors assist candidate selection, consequently minimising calculation times and reducing memory consumption. One-hot encoding would facilitate a detailed causal analysis of MCA and component elements, and elemental contributions could be weighted to identify factors improving MCA. In future research, descriptors should be designed according to the purpose of the analysis. These results provide a guideline for extending the search space and the exploration of novel, large MCA materials.

Theoretical analysis of the origin of magnetocrystalline anisotropy
For the next step, we performed theoretical analyses on the perpendicular MCA for the materials with the largest MCA energies (as predicted by the Bayesian optimisation and DFT calculations). In Table 4, we present the MCA energy and orbital moment anisotropy (OMA, Δm orb ) of the materials with the top three MCAs, namely, (Fe/Cu/Fe/Cu), (Fe/Cu/Co/ Cu), and (Fe/Co/Fe/Ni) [42][43][44]. The Δm orb was defined by the difference in the orbital moments between the magnetisation directions [001] (m orb 001 ½ �) and [100](m orb 100 ½ �), that is, Δm orb ;m orb 001 ½ � À m orb 100 ½ �, and was proportional to the MCA energy under certain conditions, referred to as the Bruno formularisation [45]. Since the total OMA values of the materials in Table 4 were positive, the Δm orb values were consistent with the MCA results. However, the order of material dependence on Δm orb was (Fe/Cu/Fe/Cu)>(Fe/Co/Fe/Ni)> (Fe/ Cu/Co/Cu), which was different from the corresponding order of material dependence on the MCA energy. This implied that an additional factor plays an important role in the perpendicular MCA. Regarding the atomic OMAs, the Δm orb of Fe in (Fe/Cu/Fe/Cu), Δm orb of Co in (Fe/Cu/Co/Cu), and Δm orb of Fe and Co in (Fe/Co/Fe/Ni) showed larger values than the other constituent elements, indicating their higher contributions to the perpendicular MCA.
Subsequently, we performed second-order perturbation analyses of the spin-orbit interaction to further understand the MCA mechanism. The theoretical details of the perturbation calculations can be found in Ref. [4]. Figure 5 indicates the contributions of the perturbation processes to the MCA between the initial and final spin-state,σ, through the intermediate spin states, σ 0 , for each atom. The MCA contribution of Δm orb in Table 4 corresponds to E ## MCA term and does In terms of MCA, the fct structure with extended structure in the c-axis direction was determined to be favourable. not include other terms (E "" MCA , E "# MCA , E #" MCA ). Figure 5 shows that in the case of (Fe/Cu/Fe/Cu), which had the largest MCA energy, the contribution of E ## MCA term of Fe atom to the MCA was significant, while that of E "# MCA (spin-flip) term surpassed that of E ## MCA term. Thus, both the spin-conserving term (E ## MCA ) and the spin-flip term (E "# MCA ) contributed positively to the MCA energy, which was the origin of the large perpendicular MCA energy of (Fe/Cu/Fe/ Cu). Figure 5 shows that, in the case of (Fe/Cu/Co/ Cu), which had the second-largest MCA energy, the spin-conserving term (E ## MCA Þ of Co primarily contributed to the perpendicular MCA, as opposed to that of Fe, which was consistent with the Δm orb result in Table 4. However, the spin-flip term (E "# MCA ) of Co was also positive, leading to the larger perpendicular MCA of (Fe/Cu/Co/Cu). In the case of (Fe/Co/Fe/Ni), which showed the third-largest MCA energy, the spinconserving terms of Fe and Co predominantly contributed to the perpendicular MCA. On the other hand, the spin-flip term of Ni had a large negative contribution, resulting in the reduction of the perpendicular MCA in (Fe/Co/Fe/Ni). In Table 4, the material dependence of the total Δm orb was not consistent with that of the MCA energy, wherein the total Δm orb of (Fe/Co/Fe/Ni) was almost comparable to that of (Fe/Cu/Fe/Cu). In parallell, the MCA of (Fe/Co/Fe/ Ni) was approximately half that of (Fe/Cu/Fe/Cu). Figure 5 clarifies that this was due to the large negative spin-flip term (E "# MCA Þ of Ni. In Figure 6, we show the PDOS results on thed orbitals of Fe in (Fe/Cu/Fe/Cu) and Co in (Fe/Cu/Co/Cu). According to the formulation of the second-order perturbation analysis shown in Equation (1) of Ref. [4], the electronic states around the Fermi level play an important role in the MCA. In Figure 6, the minority-spin states of Fe d(3z 2 -r 2 ) in (Fe/Cu/Fe/Cu) showed a peak around the Fermi level. From the analysis of the matrix elements in the second-order perturbation, we found out that the spin-flip terms between the majority-spin d(yz), and the minority-spin d(3z 2 -r 2 ) significantly contributed to the perpendicular MCA of (Fe/Cu/Fe/Cu). On the other hand, there were four d orbitals around the Fermi level ( Figure 6), namely, d(3z 2 -r 2 ),d(zx), d(yz), and d(xy), that provided a positive spin-conserving term (E ## MCA ) and two spin-flip terms (E "# MCA ,E #" MCA ), which contributed to the perpendicular MCA of (Fe/  Cu/Co/Cu). In general, to enhance the MCA of transition metal alloys, it is important to design electronic structures with peaks of PDOS around the Fermi level. However, the intuitive human-driven design of such electronic structures with strong perpendicular MCA is difficult because the PDOS around the Fermi level can contribute to both positive and negative MCA energies. In the present study, our machine learning using the Bayesian optimisation offered ferromagnetic materials with perpendicular MCA energies, wherein both the spin-conserving term (E ## MCA ) and the spinflip term (E "# MCA ,E #" MCA ) contributed to the positive MCA energy. It is usually difficult to increase the MCA by simultaneously considering these two terms. This signifies that our machine-learning using the Bayesian optimisation was successfully carried out, and its effectiveness was demonstrated from the electronic-structure viewpoint.

Experiments
A material search using Bayesian optimisation assumes a perfectly ordered crystal structure and uses theoretical data from first-principles calculations. To demonstrate this method experimentally, we fabricated the first-, second-, and third-best material candidates in terms of MCA energy under laboratory conditions and evaluated their magnetic anisotropy. The samples were fabricated by PLD using a Nd:YAG laser (wavelength: 266 nm, pulse width: 6 ns, repetition frequency: 10 Hz), and the monoatomic layers were stacked alternately. To prevent droplet generation and to ablate the target sufficiently, the laser power was controlled to 7.0 ± 0.07 mJ (LabVIEW TM , National Instruments). First, a Cu layer (25 nm) was deposited on a singlecrystal MgO (100) substrate as the base layer and heated at 300°C for 30 min. In total, 52 [Fe/Cu/Fe/Cu] 13 , [Fe/Cu/Co/Cu] 13 , and [Fe/Co/Fe/Ni] 13 monoatomic layers were deposited at room temperature (22.1±1 °C), followed by a layer of Ru (10 nm) as an antioxidant layer.
The surface roughness and morphology of the samples were observed by AFM (Dynamic mode, NaioAFM, nanosurf). Reflection high-energy electron diffraction (RHEED) and XRD (SmartLab, Rigaku) were used to confirm the epitaxial growth. The magnetic properties were measured using a SQUID (DC mode, MPMS3, Quantum Design). The maximum magnetic field applied in the in-plane and perpendicular directions was ±30 kOe. The MCA (K u ) was calculated from the difference in the area enclosed by the magnetisation curves in the in-plane and perpendicular directions and was corrected for shape anisotropy using Eq. (5): where H, M, and M s are the applied magnetic field, magnetisation, and saturation magnetisation, respectively.

Magnetic property analysis
XRD spectra of the crystal structures showed only fundamental peaks along the orientation of the Cu underlayer for all samples. RHEED patterns had sharp streaks, indicating a smooth surface condition; AFM images also indicated that a sufficiently smooth surface was obtained. The crystal structure characterisation by XRD and surface morphology analysis by RHEED and AFM of the fabricated specimens are presented in Appendix A. Single-crystal thin films with magneticresearch-grade quality were obtained, that is, the quality of the obtained multilayer was equivalent to that in previous studies demonstrating MCA improvements [6]. Figure 7 shows the magnetisation curves of the fabricated thin films. The spontaneous easy axis of magnetisation was along the in-plane direction in both samples. The magnetic moment (M s ) and K u values of the three thin films studied here were determined as follows: [Fe/Cu/Fe/Cu] 13 (Figure 7). These were the largest magnetic properties reported by our group thus far ( Figure 7) [6,7]. The improvement in K u was attributed to the low content of non-magnetic Cu and the increase in the demagnetisation energy term (2πM s 2 ) in Equation (5). In [Fe/Co/Fe/Ni] 13 , out-of-plane magnetisation curves showed a sudden rise near zero field, which was followed by a gradual, discontinuous increase. Further, the saturation of the in-plane magnetisation of this sample was more difficult than that of other samples. This result could suggest that the sample begins to partially acquire spontaneous perpendicular magnetisation, reflecting an improvement in MCA. In general, such a discontinuous magnetisation curve cannot be explained by a simultaneous magnetisation rotation based on the Stoner-Wohlfarth model, suggesting that the inhomogeneities in the magnetic domain structure and grain boundaries could play a role. It is hypothesised that regions with large perpendicular components are locally generated in the sample. Microscopic analysis using photoemission electron microscopy or magnetic force microscopy could elucidate this [21,46]; this investigation would be conducted in the future.
For applications, we utilise the effective magnetic anisotropy energy, K u eff , considering the effect of shape anisotropy (Table 5). K u eff was obtained using Equation (5), òH in dM À òH out dM. The spontaneous easy axis of magnetisation was in the in-plane direction; hence, K u eff was negative for all samples. The negative trend in K u eff was first due to the square increase in M s in shape magnetic anisotropy and second due to the influence of order-disorder structural changes on MCA, which were examined by simulation. We considered disordered cases by using the akai-KKR-CPA code (cpa2002v010) [47,48], wherein the conditions for calculating MCA in the ordered cases are the same as those when using Quantum ESPRESSO. In the KKR-CPA calculations, LDA was used for the exchange-correlation function, and the coherent potential approximation (CPA) was used to account for substitution-type disorder effects [49,50]. In general, MCA increases with the degree of order; hence, two types of calculations were used to investigate the contribution of the degree of order. Table 6 shows the calculated MCA for ordered and disordered systems, based on theoretical and experimental lattice constants. For all candidates, the MCA is nearly zero in the disordered phase (i.e. the degree of order is zero), indicating that a large MCA is achieved   by structural ordering. The same behaviour was confirmed by the theoretical and experimental lattice constants, and these results agree with those of previous studies on FeNi and FeCo [51,52]. Therefore, it is suggested that the MCA can be further improved by increasing the degree of order through process optimisation, such as growth temperature optimisation. In view of the above, this study successfully demonstrates the effectiveness of our novel ferromagnetic material exploration system that integrates Bayesian optimisation and first-principles calculations. We anticipate our system to offer an optimal platform for developing novel ferromagnetic materials with high MCA. We performed material prediction using computation and informatics in virtual space and integrated this with experimental materials fabrication in real space. We have determined the fundamental principles for discovering novel magnetic materials and have experimentally demonstrated the proof-ofconcept of the proposed method based on real materials. Expanding the material search space by increasing the number of elements and layers and by design additional features in the temperature space would help improve the structural ordering. This development would facilitate further enhancement of MCA and offer a unique approach to the discovery of unexplored magnetic materials.

Conclusions
In this study, we developed a novel synthetic method for fabricating ferromagnetic materials with large magnetic anisotropy, by combining first-principles calculations and Bayesian optimisation. We investigated the optimum machine learning parameters for achieving an efficient candidate material search, which is independent of the initial training data. By using the EI as the acquisition function and Matern52 as the kernel function, we could find a material with a large MCA energy within less than 20% of the total number of iterations. Following our strategy, we experimentally fabricated the [Fe/Cu/Fe/Cu] 13 , [Fe/Cu/Co/Cu] 13 , and [Fe/Co/Fe/Ni] 13 ferromagnetic alloys and demonstrated that [Fe/Cu/Co/Cu] 13 and [Fe/Co/Fe/Ni] 13 exhibit particularly large K u values that exceed those of L1 0 -FeNi-and L1 0 -FeCo-type alloys. In addition, the second-order perturbation analysis of the discovered candidates ((Fe/Cu/Fe/Cu) and (Fe/Cu/Co/Cu)) inferred that both the spin-conserving and spin-flip terms were the origin of the perpendicular magnetic anisotropy. This result suggests that Bayesian optimisation can discover promising materials in terms of electronic structure, which extend beyond empirical experience and human intuition. This technique is applicable to advanced magnetic materials with highly complicated electronic correlations, such as Heusler alloys and spin-thermoelectric materials. Finally, the integration of our strategy with the field of robotic materials, coupled with advanced synchrotron radiation analysis, could further promote the development of autonomous material exploration systems.