I. INTRODUCTION
Thin films based on magnetic materials are of considerable scientific and technological interest. Particular attention has been devoted to materials that exhibit a preferred magnetization direction, or magnetic anisotropy. For instance, a number of materials whose magnetization direction is perpendicular to the substrate have found multiple applications in the field of magneto-optical and perpendicular magnetic recording. Following the global downsizing trend, the road toward high-density data storage requires further reduction of the grain size. While for longitudinal magnetic recording, thermal fluctuation is identified as the main limit to further size reduction,Reference Weller and Moser1 for perpendicular recording the necessity to reduce the flying height of the read/write head with respect to the media surface results in various tribological issues such as friction heat, nanometer scratches leading to plastic deformations, etc. The mechanical constraints produced by such events can potentially alter the magnetization of the magnetic layer. For instance, lattice distortions are known to affect the magnetocrystalline anisotropy of Co-based films.Reference Lee, He, Lamelas, Vavra, Uher and Clarke2–Reference Andersson and Hjörvarsson6
The observation that the magnetic properties decrease due to the scratch damage of the CoCrPt magnetic layers in a magnetic memory storage device inquires an emergency question; what is the origin of the changes in the magnetic anisotropy energy (MAE) under the deformation for the magnetic layer in the devices? Experimental works have reported that under normal strain the MAE trends of Co-based films and hcp Co are similar.Reference Shimatsu, Okazaki, Sato, Muraoka, Aoi, Sakurai, Okamoto, Kitakami, Tanii and Sakuma7,Reference Wang, Sakurai, Oikawa, Ishida, Kikuchi, Okamoto, Sato, Shimatsu and Kitakami8 This is most likely due to the fact that the hexagonal close-packed (hcp) Co sustains the largest MAE among 3d transition metals, resulting in a large MAE for the Co-based magnetic film. As a result, to improve the efficiency and reliability of magnetic memory devices, it is important to understand the influence of strain deformations on the MAE of the hcp Co.
During the last two decades, many theoretical works have pointed out the difficulty to reproduce the experimental value of the MAE in bulk hcp Co. Originally, the experimental work of Paige et al.Reference Paige, Szpunar and Tanner9 reported a MAE of 65 μeV/atom for hcp Co at low temperature (T = 4.2 K), with the magnetization easy axis oriented along the [0001] direction and the hard axis along the [
$10\bar 10$] direction. Several theoretical works attempted to calculate the MAE of hcp Co at equilibrium and under strain. Daalderop et al. used the linear muffin-tin orbital method (LMTO) in the atomic-sphere approximation (ASA) within the local density approximation (LDA) and obtained a MAE magnitude of 29 μeV/atom including the spd orbitals, which is in the order of the experiment. However, the easy axis of magnetization in hcp Co was incorrectly predicted.Reference Daalderop, Kelly and Schuurmans10 In an attempt to find the origin of the disagreement with the experiment, Daalderop et al. estimated the MAE for different crystal volumes and noted a dependency of the MAE with strain. The correct easy axis of magnetization was later predicted by Trygg et al.Reference Trygg, Johansson, Eriksson and Wills11 by means of an LDA approach in combination with an orbital polarization (OP) correction method.Reference Brooks12,Reference Eriksson, Brooks and Johansson13 The addition of the OP correction to the LDA resulted in significant improvements of the MAE and the metal orbital moments. Furthermore, they showed the importance of the k-point sampling for accurate determination of total energy differences. Beiden et al.Reference Beiden, Temmerman, Szotek, Gehring, Stocks, Wang, Nicholson, Shelton and Ebert14 used a fully relativistic and spin-polarized extension of the locally self-consistent multiple-scattering method. Stiles et al.Reference Stiles, Halilov, Hyman and Zangwill15 used a spin-other-orbit interaction method within LDA. Both studies did not succeed to find the correct easy axis of magnetization. Recently, a calculated MAE of 65 μeV/atom for hcp Co has been reported by using the strong electronic correlations LDA+U approach (U = 1 eV) within a tight-binding linear muffin-tin orbital method in the atomic sphere approximation (TB-LMTO-ASA).Reference Xie and Blackman16 The change in the MAE for the 3d transition metals (Fe, Co, Ni) under deformation was reported by Burkert et al.Reference Burkert, Eriksson, James, Simak, Johansson and Nordström17 with a full-potential LMTO method and the force theorem. It was found that for the tetragonal Co, a reorientation of the magnetization easy axis appears in the range of c/a < 1/
$\sqrt 2$ and c/a > 1. However, little is understood about the change in the magnetic anisotropy of the hcp Co under deformation and requires further investigations with respect to the increasing interest in applications of magnetic films based on Co or Co alloys.
In this study, with the aim to understand the effect of normal and shear strain on the MAE of hcp Co, we perform ab initio density functional theory (DFT) calculations to obtain the landscape of the MAE and the magnetic moments (MMs) as the variation in the lattice parameters a and c/a ratio around equilibrium. This paper is organized as follows. Section II explains the computational details of the ab initio electronic structure calculations. Results and discussion are given in Sec. III. Section IV concludes the paper.
II. METHODOLOGY AND SIMULATION PROCEDURE
The pseudopotential plane-wave (PP-PW) method performs electronic structure calculations, where only the valence electrons are explicitly treated with the frozen-core approximation. This method enables efficient calculation of force exerted on atoms and stress; i.e., it is suitable for evaluation of properties under deformation. We used the Vienna Ab-initio Simulation Package (VASP)Reference Kresse and Hafner19,Reference Kresse and Furthmüller20 for the pseudopotential calculations based on the projector augmented wave method (PAW).Reference Blöchl21 We adopted the generalized gradient approximation (GGA), and two parameterizations by Perdew and Wang (PW91)Reference Perdew and Wang22 and by Perdew–Burke–Ernzerhof (PBE)Reference Perdew, Burke and Ernzerhof23 were used. For comparison, calculations using the LDA were also performed. Since it was pointed out that the number of valence electrons can strongly influence the MAE, the exchange-correlation function DFT-pv and DFT-sv were also considered in test calculations. In the DFT-pv pseudopotentials, the semi-core p states are treated as valence, whereas the semi-core s states are treated as valence for the DFT-sv.
The bulk hcp Co crystal contains two atoms per unit cell, and the MAE is obtained as the difference in total energy between the [
$10\bar 10$] and [0001] magnetization directions (ΔE = 
$E_{10\bar 10} - E_{{\rm{0001}}}$). Test calculations for the hcp Co have shown that a difference of 2 μeV in MAE was obtained between 700 and 500 eV cut-off energies suggesting that a cut-off energy of 500 eV is a reasonable compromise between accuracy and efficiency. All the calculations are performed in two steps. First, relaxation with the optimization of atom positions and cell dimensions is carried out to determine the geometry and magnetic ground state at equilibrium. The self-consistent electronic steps were conducted until the energy change between subsequent iterations becomes less than 1 × 10−7 eV. A minimum of 50 self-consistent electronic steps were carried out to ensure a good convergence of the total energies. The atomic relaxation was performed by means of a conjugate gradient scheme until the maximum force acting on each atom became less than 0.01 eV/Å, simultaneously the cell dimensions were optimized until all the stress components come within ±12.5 MPa. Second, the structure obtained in the first step is used to perform the total energy calculations including spin–orbit coupling, with the magnetic axis oriented along the [0001] and [
$10\bar 10$] directions. All the calculations are performed by two Monkhorst–PackReference Monkhorst and Pack24 k-point meshes of 16 × 16 × 10 and 20 × 20 × 14, corresponding to 2560 and 5600 k-points, respectively. The self-consistent electronic calculations and geometry optimizations were performed without symmetry constraints. Therefore, the numbers of k-points were not reduced.
The DFT+U approach including the effect of the strong electronic correlations was also used for comparison. We adopt the DFT+U approach formulated by Dudarev et al.Reference Dudarev, Botton, Savrasov, Humphreys and Sutton25 as implemented in VASP. Within this method, the difference between the Coulomb and the exchange parameters is expressed by a single effective Hubbard term U eff = (U − J). To clarify the effect of the +U term on the equilibrium lattice parameters, we calculated a and c/a as a function of (U − J).
After choosing the exchange-correlation function for a proper description of the MAE, we investigated the effect of strain on the magnetic properties. First, the effect of normal strain on the MAE and MM with varying lattice parameters, a and c/a (2.47 < a < 2.51 Å, 1.595 <c/a < 1.625), were calculated. We then applied shear strain, γxy and γzx, to examine the effect of shear on the magnetic properties. Here we adopt the notation of engineering shear strain. Increasing shear strain, γxy or γzx, was applied to the cell keeping the volume of the cell constant, and atom positions were optimized at each strain. A range of shear strain from γ = 0.0 to γ = 0.08 was investigated.
III. RESULTS AND DISCUSSION
A. Comparison of the MAE for various exchange-correlation functions and initial parameters
Since the MAE for the hcp Co is of the order of several μeV per atom, ab initio calculations have to be performed with very high accuracy. Various technical options are available and should be carefully chosen to ensure that the initial setting parameters will properly describe the material properties. Such options are, e.g., the exchange-correlation functions, the cut-off energies, the k-point mesh, the values of U for the Hubbard +U method, etc.
As it is important to choose the proper functionals to examine the MAE of hcp Co, our first approach of the problem was to investigate the influence of the exchange-correlation function. Figure 1 presents several test calculations of MAE/atom for choosing the proper exchange-correlation function. The figure compares the MAE (
$E_{10\bar 10} - E_{{\rm{0001}}}$) at equilibrium state for various exchange-correlation functions with 16 × 16 × 10, 20 × 20 × 14, and 22 × 22 × 16 k-point mesh corresponding to 2560, 5600, and 7744 k-points in the irreducible wedge of the Brillouin zone, respectively. Calculation of the MAE using the experimental lattice parameters are given for comparison. As seen in the figure, the GGA functionals (PW91, PBE, PBE-pv, and PBE-sv) succeed in predicting the correct easy axis of magnetization in Co, namely [0001] at ground states. However, the LDA, LDA-sv, and LDA-pv functionals predict the [
$10\bar 10$] direction as the easy axis of magnetization at ground states. The calculated MAEs with GGA yield values between 16 to 26 μeV/atom, which are fairly low compared to the experiment (65 μeV/atom) but in good agreement with other theoretical works.Reference Daalderop, Kelly and Schuurmans10 It is interesting to note that the LDA functionals predict the correct easy axis with the experimental lattice parameter; however, the calculated MAE is strongly underestimated with approximately 5 μeV/atom. It was found that the residual stresses in the case of LDA functionals are up to 20 GPa, indicating that the LDA is strongly affecting the structure of the crystal. As shown in Table I, the lattice parameters a for the relaxed structures are underestimated with GGAs by only 0.7%, whereas a 3% underestimation is obtained with LDAs, suggesting that a correct description of the crystal structure is essential for calculating the MAE of hcp Co. Hereafter in this paper, two parameterizations are chosen to perform the calculations, GGA-PW91 and GGA-PBE-pv. The calculation with the GGA-PBE-pv gives the highest MAE (26 μeV/atom) with 5600 irreducible k-points. The lattice parameters a and c/a are very close to the experiment with an underestimation of only 0.5%. Although the MAE and lattice parameters with GGA-PW91 are somewhat smaller than with GGA-PBE-pv, the GGA-PW91 functional serves as a point of reference in our study since it is a widely used approximation for ab initio calculations, which has proven to be mature enough to yield reliable results. Furthermore, as Perdew et al.Reference Perdew, Burke and Ernzerhof23 mentioned, the PBE is a simplification of PW91 and, in general cases, gives equivalent results for energies, lattice parameters, bulk modulus, etc. However, Mattsson et al.Reference Mattsson, Armiento, Schultz and Mattsson26 reported that in some complex cases, for instance system involving surfaces, the PBE and PW91 give nonequivalent results. It is therefore important to keep both parameterizations since MAE calculations are nontrivial.
Comparison of MAE = (
$E_{10\bar 10} - E_{0001}$)/atom for different exchange-correlation functions and k-point meshes. Green triangles, blue solid circles, and red solid squares refer to k-point meshes of 16 × 16 × 10, 20 × 20 × 14, and 22 × 22 × 16, respectively; 2560, 5600, and 7744 correspond to the number of irreducible k-points. Calculation of the MAE with the experimental lattice parameter and 5600 irreducible k-points is represented by black squares. Dashed lines are only a guide to the eye.
Lattice parameters for hcp Co with different exchange correlation functions compared with the experimental results by Ono (1988).Reference Ono and Maeta18
a Extrapolated values.
Regarding the number of k-points, the calculated MAE using GGA-PW91 with 2560, 5600, and 7744 irreducible k-points produced an MAE of 16.5, 16.0, and 14.0 μeV/atom, respectively. For the GGA-PBE-pv, the obtained MAE with 2560, 5600, and 7744 irreducible k-points are 7, 26, and 24.5 μeV/atom, respectively. To keep the computational cost relatively low and a reasonable accuracy, 2560 irreducible k-points with the GGA-PW91 and 5600 irreducible k-points with the GGA-PBE-pv were found to be sufficient for predicting the trend of the MAE under deformation.
The calculations with GGA-PW91 and LDA with 2560 k-points were carried out by taking into account the effect of the strong electronic correlations Hubbard U and Hund’s J in an attempt to improve the quantitative estimation of MAE. Figure 2 presents the MAE as well as the lattice parameters a and c/a for the GGA-PW91 and the LDA as a function of U eff. As shown in Fig. 2, the calculations with GGA-PW91 + U (U eff = 3.2 eV) and LDA+U (U eff = 4.2 eV) lead to an agreement with the experimental value (65 μeV/atom). However, strong deviations of the corresponding lattice parameters occur. In a previous investigation of fcc Ni using a DFT+U approach, Yang et al.Reference Yang, Savrasov and Kotliar27 found that the MAE value for fcc Ni is varying very rapidly as a function of U. We obtained a similar finding from our study for hcp Co with LDA+U and GGA-PW91 + U.
The MAE of hcp Co versus (U − J) within (a) GGA-PW91 and (b) LDA; the corresponding lattice parameters a and c/a versus (U − J) within (c) GGA-PW91 and (d) LDA.
Moreover, although the DFT+U method is known to cure the underestimation of the MAE that occurs in the conventional DFT, the choice of the value of the +U term is usually not obvious. Stiles et al.Reference Stiles, Halilov, Hyman and Zangwill15 noted that the U chosen for Ni by Yang et al.Reference Yang, Savrasov and Kotliar27 is not compatible with the value of U obtained from a realistic approach.Reference Wu and Freeman28 Xie et al.Reference Xie and Blackman16 used a value of U = 1.77 eV and J = 0.75 eV for Co, and Steiner suggested the value of U = 2.40 eV and J = 0.50 eV.Reference Steiner, Albers and Sham29 In this range of U, our calculations using the GGA-PW91+U or the LDA+U method led to the wrong easy axis. In this paper, we decided to use the conventional DFT-GGA method to evaluate the MAE and its response to strain because (i) the crystal structure is extremely sensitive to the U value and a “reasonable” crystal structure is needed to apply mechanical deformations, and (ii) the realistic U value suggested by Steiner gives a good crystal structure and a stable MAE, however the wrong easy axis is obtained. As our main interest is in finding, the response of the MAE to the variation of the lattice parameter around equilibrium, the choice of the GGA (not +U) is the most reasonable option.
As for previous studies of MAE using different functionals, Wu and Freeman found that the difference between the LDA and the GGA is not very large for the MAE calculation of a Co monolayer arranged in a square lattice. However, it became significant in the case of Ni and Fe.Reference Wu and Freeman28 They related this discrepancy to the differences in the atomic structure which may affect the magnetic properties of the material. Such differences are finding their origin in the severe underestimation of the lattice parameters with LDA for the 3d transition metals. GGA is known to improve the calculations.Reference Perdew, Chevary, Vosko, Jackson, Pederson, Singh and Fiolhais30 A similar structural-magnetic property relationship was observed in our study for hcp Co.
More generally, our results indicate that the magnitude of the MAE is underestimated regardless of the choice of exchange-correlation function. As discussed by Jansen,Reference Jansen31 this discrepancy is attributed to the neglect of orbital correlations in the functionals. In most cases, the OPReference Brooks12,Reference Eriksson, Brooks and Johansson13 correction has been a method of choice for improving the orbital moments and MAE. Although the OP method appears to be quite successful for calculating the MAE value of bcc Fe and fcc Co, it produced limited success for hcp Co and fcc Ni.Reference Trygg, Johansson, Eriksson and Wills11
For investigating the effect of mechanical deformation, it is of primary importance to obtain a good description of the crystal structure at equilibrium state. Therefore in the following, we chose not to include the +U method to investigate the MAE landscape and the effect of shear strain on the magnetic properties of the bulk Co. We use the GGA-PW91 with 2560 k-points, as well as the GGA-PBE-pv with 5600 k-points for comparison. Even though the calculated value of the MAE for hcp Co is lower than the experimental one, the qualitative trend for the variation of the MAE as a function of strain is assumed to be correctly described.
B. Magnetic properties at equilibrium and under normal strain
In Fig. 3, the MAE landscape, the MM, and the total energy as a function of the lattice parameters, a and c/a, are presented for the GGA-PW91 and the GGA-PBE-pv functionals. For the GGA-PW91 and 2560 k-points, the minimal total energy is obtained for a = 2.4936 Å and c/a = 1.6150 with a MAE of ∼16 μeV/atom. For the GGA-PBE-pv and 5600 k-points, the minimal total energy is obtained for a = 2.4938 Å and c/a = 1.6158 with a MAE of ∼28 μeV/atom. The trends observed here for both functionals are similar to each other, proving again the equivalence of the PBE and PW91 parameterizations. The perturbations observed on the MAE landscape for GGA-PW91 is a consequence of the difficulty in estimating the MAE since we are dealing with numbers of the order of a few μeV, and therefore, slight deviations in the MAE can be obtained locally.
(Top panels) Calculated MAE versus a and c/a for the (a) PW91 and (b) PBE-pv functionals; (Middle panels) the MMs versus a and c/a for the (c) PW91 and (d) PBE-pv functionals; (Bottom panels) the energy landscape of a and c/a for the (e) PW91 and (f) PBE-pv functionals. A dotted line indicates zero stress σz.
The reduction of the MMs with decreasing crystal volume, as seen in Figs. 3(c) and 3(d), is in good agreement with other theoretical works obtained for various magnetic materials.Reference Wang, Albina and Umeno32–Reference Bagayoko and Callaway38
Our investigation suggests that the MAE of a Co film can be increased by expanding the lateral lattice parameter a and reducing the c/a ratio following the line of zero stress in the z direction (σz = 0) as plotted in Figs. 3(e)–3(f). These findings are consistent with recent experimental observations.Reference Shimatsu, Okazaki, Sato, Muraoka, Aoi, Sakurai, Okamoto, Kitakami, Tanii and Sakuma7,Reference Wang, Sakurai, Oikawa, Ishida, Kikuchi, Okamoto, Sato, Shimatsu and Kitakami8 They found that the MAE of epitaxially grown Co and Co-based alloy thin films strongly depends on the axial ratio c/a of the hcp crystal lattice. It was observed that the MAE in hcp Co is enhanced by the decrease of c/a. The agreement between our results and the experiments supports the validity of the calculated response of the MAE with respect to normal strains at least qualitatively.
C. Magnetic and mechanical properties under shear
1. Effect of shear strain on (1000) along [
$10\bar 10$]
Figure 4 presents the change in MAE of hcp Co under the shear strain, γxy, with the GGA-PW91 (2560 k-points) and with the GGA-PBE-pv (5600 k-points). Two energy differences (ΔE) between the magnetic directions, [
$10\bar 10$] and [0001], [1000] and [0001], are presented. As shown in Fig. 4(a), the magnetization direction along [1000] should follow the atomic displacement under strain. However, since we are dealing with very small shear strain amplitudes, the orientation of the magnetization direction under shear is very close to [1000]. All the ΔE values decline with increasing shear strain. The reduction of ΔE is less than 25 μeV/atom with the PW91, whereas for the PBE-pv, the reduction of ΔE is twice as large with approximately 50 μeV/atom.
(a) Shear strain γxy in the hexagonal structure; magnetic anisotropy energy as a function of γxy with (b) GGA-PW91 and 2560 k-points, and (c) GGA-PBE-pv and 5600 k-points; stresses as a function of γxy with (d) GGA-PW91 and 2560 k-points, and (e) GGA-PBE-pv and 5600 k-points.
Above a shear strain of γxy = 0.05, the easy axis reorients itself to [1000]. The critical shear strain of 0.05 corresponds to a shear stress of τxy = 6.65 GPa [see Figs. 4(d) and 4(e)]. Note, however, that large normal stresses are superimposed due to the method of shear with constant volume, which may also slightly affect MAE.
2. Effect of shear strain on (0001) along [1000]
Figure 5(a) illustrates the crystal structure under the shear strain, γzx. Similar to the case of the shear on the (1000) plane, the MAE decreases with increasing shear strain. While the GGA-PW91 predicts a reorientation of the magnetization easy axis occurring around a shear strain γzx = 0.07, the critical strain estimated by the GGA-PBE-pv is slightly higher (>0.08).
(a) Shear strain γzx in the hexagonal structure; magnetic anisotropy energy as a function of γzx with (b) GGA-PW91 and 2560 k-points, and (c) GGA-PBE-pv and 5600 k-points; stresses as a function of γzx with (d) GGA-PW91 and 2560 k-points, and (e) GGA-PBE-pv and 5600 k-points.
From Figs. 5(d)–5(e) one can estimate that the critical shear strain (γzx) for the switching of the magnetization easy axis is ∼6 GPa. Note that normal stress components are small even though the volume of the cell was kept constant during shear (i.e., normal strains are not adjusted).
Magnetic disk devices utilize the magnetic anisotropy of magnetic metal for data storage. Co-based alloys such as CoCrPt are typical materials for the magnetic layer in the devices, where the c axis is oriented perpendicular to the disk plane for the perpendicular magnetic recording, while the magnetic layer is in a polycrystal form. Obviously, our findings of Co single crystal cannot directly be extrapolated to alloys, but the characteristic of Co single crystal may account for the feature of the magnetic layer because Cr and Pt tend to be segregated at grain boundaries and Co composes a substantial fraction of bulk.Reference Honda, Funamono, Hasegawa, Kawasaki, Kugiya, Koizumi, Yoshida and Tonomura39,Reference Bernards, Schrauwen and Luitjens40 In addition, Wang et al. reported a similar trend of MAE versus c/a was found in hcp Co, Co-Pt, and Co-Pd alloys.Reference Shimatsu, Okazaki, Sato, Muraoka, Aoi, Sakurai, Okamoto, Kitakami, Tanii and Sakuma7,Reference Wang, Sakurai, Oikawa, Ishida, Kikuchi, Okamoto, Sato, Shimatsu and Kitakami8 Therefore, our DFT calculations for MAE of hcp cobalt single crystal and its change under strain might provide useful insights to the feature of the device.
It is known that the device performance may deteriorate during operation, and it is suspected, empirically, that mechanical impact to the magnetic layer may be a culprit for the degradation.Reference Xu, Furukawa, Nakamura and Honda41 In particular, shear deformation is likely to occur because the disk in fast rotation is hit by a magnetic head above. Our simulation results suggest that the magnitude of MAE decrease or even disappear due to shear strain seems consistent with the empirical speculation that the performance of the magnetic layer is affected by mechanical damage. However, theoretical calculations of Co-based alloys are needed for more direct evidence of the mechanism of degradation.
IV. CONCLUSION
We carried out ab initio electronic calculations of the MAE (
$E_{10\bar 10} - E_{0001}$) of bulk hcp Co at equilibrium and under strains (normal and shear). In search of the ideal exchange correlation function for describing the MAE, it was found that the GGA could predict the correct easy axis but underestimate the MAE with respect to the experimental value. Although the experimental value of the MAE could be obtained by taking into account the strong electronic correlations Hubbard U, the results of the DFT+U calculations have proven to be unreliable for predicting the behavior of MAE under deformation due to strong lattice parameter distortions. Our calculations with two GGA parameterizations suggest that a pure isotropic tensile or compressive strain will lead to a reduction of the MAE. However, the simultaneous reduction of the c/a ratio and increase of the lateral lattice parameter a will strongly enhance the MAE of the material. This is consistent with recent experimental observations of Co thin films, supporting the validity of our calculations.
Finally, two shear strains, γxy and γzx, were applied on the bulk hcp Co with GGA (PW91 and PBE-pv). We found substantial effect of shear strain on the MAE; shear strain leads to the reduction of MAE or even switching of the magnetic easy axis. Although the physical complexity of magnetic memory devices based on Co alloys is far greater than that of the simple bulk material, our results seem to support the speculation that the mechanical damage in the magnetic layer can deteriorate the performance of the device.
ACKNOWLEDGMENTS
The present work was partly supported by the GMSI (Global Center of Excellence for Mechanical Systems Innovation) program at the University of Tokyo: A GCOE (Global Center of Excellence) program by Ministry of Education, Culture, Sport, Science and Technology (MEXT). J.W. gratefully acknowledges the MEXT scholarship program.
