2024 Volume 92 Issue 4 Pages 043024
A polarizable ion model was applied to the solid and molten alkaline-earth halide MX2, the parameters of which were determined by using first-principles calculations based on density functional theory, where M = Ca, Sr, Ba, and X = F, Cl, Br. The obtained parameters were used to evaluate the ionic conductivity, shear viscosity, and thermal conductivity in the molten and solid states by molecular dynamics simulations using the Green-Kubo relations. The calculated results were in good agreement with the experimental ionic conductivities and shear viscosities. The behaviors of all the calculated properties were well accounted for by ionic mass, number density, and packing fraction. Especially, the calculated thermal conductivities were well expressed by the empirical formula obtained for molten alkali halides. In addition, it was revealed that the reversal of cationic dependence in ionic conductivity of fluorides between solids and melts is due to the mass effects of carrier ions.
Alkaline-earth halides, not as much as alkali halides, are simple ionic substances with high melting points, and their melts are characterized by high thermal and electrochemical stability and low vapor pressure. Therefore, some of their melts are useful as solvents and electrolytes for many industrial applications. For example, CaF2 has been widely used as a flux to reduce the viscosity of slag melts. CaCl2 has been extensively investigated as one of the main components of electrolytes in molten salt electrolysis for novel production methods of minor metals such as Ti, Nb, Ta, and Dy. On the other hand, some of their crystals have unique properties that are unknown for alkali halides. Since CaF2 single crystal has good transmittance over a wide range of wavelength and unique dispersion characteristics of refractive index, it is in particular used as optical elements such as fluorite lenses and laser media (Yb:CaF2). Fluorite-type crystal structures such as CaF2 and SrCl2 are known to have superionic conductive phases as well.
For fundamental understanding of superionic conductive phases, extensive molecular dynamics (MD) simulations have been performed on CaF2 and SrCl2.1–12 In the early stages of classical simulation, Dixon and Gillan,3–9 and de Leeuw10–12 have reported results using various empirical potentials, since the pioneering work by Rahman.1,2 First, Rahman simulated CaF2 using the numerical values for the repulsive part of the potential obtained by Kim and Gordon.13,14 Dixon and Gillan then simulated CaF2 and SrCl2 using the rigid ion part of the shell model for alkaline-earth halide crystals developed by Catlaw et al.15,16 On the other hand, de Leeuw performed simulations using a modified CaF2 rigid ion model by Benson and Dempsey.17,18 Whereas Dixon and Gillan, and Rahman focused on studying the superionic phases, de Leeuw was more concerned with the structure and properties of the liquid state of SrCl2.
The shell model was used by Catlaw and Norgett15 to introduce polarization effects, but was not used by Dixon and Gillan as it was considered computationally expensive and insignificant. However, de Leeuw pointed out the necessity of introduction of polarization effects to improve the interionic potentials. Polarization effects in ionic condensed phases have been successfully formulated by Madden19 as a polarizable ion model (PIM) based on first-principles electrostatic calculations. Its importance has already been reputed for many ionic substances by MD calculations using PIM.20 For CaF2, Dent et al.21 performed MD calculations using PIM and compared the results with those of Gillan,5,6 and discussed why the empirical potential function used by Gillan can quantitatively reproduce the experimental results on superionic conduction without considering the polarization effect explicitly. As a result, it was concluded that, in the Gillan potentials, the unphysically large dispersion forces acting between fluoride ions substantially play the role of the polarization interactions.21
The purpose of this study is to systematically provide PIM force field parameters based on first-principles calculations of alkaline-earth halides, focusing on the evaluation of the transport properties mainly in the molten state. Although PIM parameters have been reported fragmentarily for these systems,21–23 there remains room to reexamine from systematic perspectives for a transferability in various molten salts, and for a further application in oxyhalide-based glass ceramics.24 What is meant by systematicity here is, for example, that the parameters among fluoride ions are common in a series of fluorides including alkali halides. An advantage of sharing parameters with like ions is that it facilitates future expansion to mixed systems. PIM has so far been applied to mixed systems such as LiF–BeF2,25 LiF–ZrF4,26 LiF–YF3,27 and the experimental results of the structure and dynamical properties of these melts have been successfully reproduced. One of the advantages of PIM is that it can be applied to different conditions such as high temperature and high pressure. In recent years, the phase diagrams of alkaline-earth halides at high pressure have been intensively studied from both standpoints of experiments and simulations.28–32 The results of this research are expected to be useful not only for application to these high-pressure phases, but also for applications of BaF2-LaF3 to fluoride ion batteries.
In the previous work,33 the force field for PIM has been successfully parameterized for a series of alkali halides. In this study, we systematically construct the parameters for MX2, where M = Ca, Sr, and Ba, and X = F, Cl, and Br, while maintaining consistency with those of alkali halides. These compositions cover all major salts with superionic conductive phases (CaF2, SrF2, BaF2, SrCl2, and SrBr2) in alkaline-earth halides. In order to verify the reliability of the obtained parameters, MD calculations were performed to investigate the shear viscosity, ionic conductivity, and thermal conductivity in the molten state. The calculated ionic conductivity and shear viscosity were compared with experimental results.34–36 Since no experimental result has been reported regarding thermal conductivity in the molten state, this estimated result would be useful as information comparable to experimental values. In addition, thermal conductivity and superionic conductivity in the solid state of some fluorides were also evaluated and compared with experimental and other calculated results.37–43 The behaviors of these calculated transport properties will be explained by their dependence on ionic mass, number density, and packing fraction derived from dimensional analyses of the properties.44,45
The PIM potential function consists of charge-charge, polarization, repulsive and two dispersion terms as
| \begin{align} \phi^{\text{tot}} & = \sum_{i} \sum_{j > i} \frac{q_{i}q_{j}}{r_{ij}} + \phi^{\text{pol}} + \sum_{i} \sum_{j > i} B_{ij}\exp(-\alpha_{ij}r_{ij})\\ &\quad - \sum_{i} \sum_{j > i} \left(\frac{C_{6}^{ij}}{r_{ij}^{6}}f_{6}^{ij}(r_{ij}) + \frac{C_{8}^{ij}}{r_{ij}^{8}}f_{8}^{ij}(r_{ij}) \right) \end{align} | (1) |
where qi is the charge of particle i: +2 e for cations and −1 e for anions. The polarization energy term ϕpol includes charge-dipole and dipole-dipole interactions as
| \begin{align} \phi^{\text{pol}} & = \sum_{i} \sum_{j > i} \left(\frac{q_{i}\boldsymbol{r}_{ij} \cdot {\boldsymbol{\mu}}_{j}}{r_{ij}^{3}}f_{4}^{ij}(r_{ij}) + \frac{{\boldsymbol{\mu}}_{i} \cdot q_{j}\boldsymbol{r}_{ji}}{r_{ij}^{3}}f_{4}^{ji}(r_{ij}) \right)\\ &\quad + \sum_{i} \sum_{j > i} \left(\frac{{\boldsymbol{\mu}}_{i} \cdot {\boldsymbol{\mu}}_{j}}{r_{ij}^{3}} + \frac{3({\boldsymbol{\mu}}_{i} \cdot \boldsymbol{r}_{ji})(\boldsymbol{r}_{ij} \cdot {\boldsymbol{\mu}}_{j})}{r_{ij}^{5}} \right) + \sum_{i} \frac{|{\boldsymbol{\mu}}_{i}|^{2}}{2\alpha_{i}} \end{align} | (2) |
where μi is the induced dipole moment of particle i. The fnij is damping function which corrects short-range errors in the charge-dipole and dispersion interactions,46
| \begin{equation} f_{n}^{ij}(r_{ij}) = 1 - c_{n}^{ij}\exp(- b_{n}^{ij}r_{ij}) \sum_{k = 0}^{n} \frac{(b_{n}^{ij}r_{ij})^{k}}{k!}. \end{equation} | (3) |
In the MD simulation, the induced dipoles are calculated at each time step by solving self-consistently the set of equations,
| \begin{equation} {\boldsymbol{\mu}}_{i} = \alpha \boldsymbol{E}_{i}(\{q_{j}\}_{j \neq i},\{\mu_{j}\}_{j \neq i}) \end{equation} | (4) |
where Ei is the electric field generated at ri by the entire set of charges and induced dipoles from the ions j ≠ i, and αi is the polarizability of ion i. These instantaneous dipole moments are practically determined by minimizing the total energy using the conjugate gradient method. The energies and forces on each ion from the dispersion, charge-charge, charge-dipole and dipole-dipole interactions are evaluated under the periodic boundary condition by using the Ewald method.47
In this work, the PIM parameters for a series of alkaline-earth halides were fitted by minimizing the deviations, χμ2 and χF2, from the dipole moments μiDFT and forces FiDFT respectively, obtained from a first-principles calculation based on density functional theory (DFT), on a series of typical liquid configurations and solid states having fast ion conductive phase which were produced using the empirical rigid ion models. The target functions are defined as
| \begin{equation} \chi_{\mu}^{2} = \frac{1}{N} \sum_{i = 1}^{N} \frac{|{\boldsymbol{\mu}}_{i}^{\text{DFT}} - {\boldsymbol{\mu}}_{i}^{\text{PIM}}|^{2}}{|{\boldsymbol{\mu}}_{i}^{\text{DFT}}|^{2}} \end{equation} | (5) |
| \begin{equation} \chi_{F}^{2} = \frac{1}{N} \sum_{i = 1}^{N} \frac{|\boldsymbol{F}_{i}^{\text{DFT}} - \boldsymbol{F}_{i}^{\text{PIM}}|^{2}}{|\boldsymbol{F}_{i}^{\text{DFT}}|^{2}} \end{equation} | (6) |
where μiPIM and FiPIM are the dipole moments and forces obtained from PIM, respectively, and N the number of particles treated in the DFT calculations. All the DFT calculations with plane-wave pseudopotential were carried out using the Car-Parrinello MD (CPMD) package.48 The Perdew–Burke–Ernzerhof generalized-gradient approximation was used for the exchange-correlation functional.49 The threshold in self-consistent procedure through a DFT calculation was set to 10−6 atomic unit (a.u.) in the maximum change of electronic gradient. The induced dipole moments μiDFT were obtained from maximally-localized Wannier functions.50,51
The polarizabilities and the other parameters obtained via the fitting procedure52 are collected in Tables 1 and 2, respectively. The polarizabilities of Ca2+, Sr2+, and Ba2+ ions were newly determined in this work, while those for halide ions were take from the previous work on alkali halides in order to hold the transferability between both the systems. Although, to the best of our knowledge, there are no experimental results on the polarizability of Sr2+ and Ba2+ ions, they are theoretically reported as 5.792 and 10.491 a.u., respectively.53 The polarizabilities obtained in this study are in good agreement with these values and experimental results presented in Table 1.53,54 In this work, the repulsive parameters between cations remain fixed regardless of cation species of M2+(Ca2+, Sr2+, Ba2+), since they are relatively insignificant in alkaline-earth halides. Then, as following the preceding manner of PIM force field,23 the repulsive interaction of large Ba2+ was set to be more emphasized than those of small Ca2+ and Sr2+, as shown in Ref. 23. Table 3 shows the final sets of the χμ2 and χF2 for a series of molten alkaline-earth halides.
| Ion pair | Bij | αij | b4ij = b4ji | c4ij | c4ji | C6ij | C8ij |
|---|---|---|---|---|---|---|---|
| F−-F− | 282.3 | 2.444 | 15 | 150 | |||
| F−-Ca2+ | 66.9 | 1.848 | 1.608 | 1.139 | 0.131 | 18 | 36 |
| F−-Sr2+ | 75.1 | 1.764 | 1.726 | 1.667 | −0.288 | 40 | 80 |
| F−-Ba2+ | 100.2 | 1.734 | 1.850 | 1.975 | −0.210 | 60 | 120 |
| Cl−-Cl− | 275.1 | 1.797 | 140 | 280 | |||
| Cl−-Ca2+ | 133.6 | 1.730 | 1.686 | 2.259 | 0.540 | 70 | 140 |
| Cl−-Sr2+ | 164.4 | 1.684 | 1.755 | 3.130 | 0.281 | 110 | 220 |
| Cl−-Ba2+ | 139.5 | 1.564 | 1.623 | 3.184 | −0.988 | 170 | 340 |
| Br−-Br− | 218.2 | 1.678 | 250 | 500 | |||
| Br−-Ca2+ | 88.9 | 1.565 | 1.849 | 5.318 | 0.982 | 40 | 80 |
| Br−-Sr2+ | 109.9 | 1.525 | 1.675 | 3.975 | −0.996 | 100 | 200 |
| Br−-Ba2+ | 156.8 | 1.521 | 1.809 | 5.751 | 1.196 | 200 | 400 |
| M2+-M2+ | 1.0 | 5.000 | 0 | 0 | |||
| Ba2+-Ba2+ | 3.0 | 5000.0 | 0 | 0 |
| Salt | χμ2 | χF2 | Salt | χμ2 | χF2 | Salt | χμ2 | χF2 |
|---|---|---|---|---|---|---|---|---|
| CaF2 | 0.037 | 0.051 | CaCl2 | 0.058 | 0.123 | CaBr2 | 0.060 | 0.116 |
| SrF2 | 0.029 | 0.031 | SrCl2 | 0.041 | 0.121 | SrBr2 | 0.064 | 0.112 |
| BaF2 | 0.031 | 0.086 | BaCl2 | 0.068 | 0.081 | BaBr2 | 0.060 | 0.102 |
In order to determine the dispersion parameters, C6ij and C8ij empirically, MD simulations were performed under NPT ensemble of ambient pressure so as to reproduce the experimental density.34 Note that these parameters for the anion pairs, e.g. F−–F−, are set to be the same as those for a series of alkali halides determined in the previous work33 to hold transferability between alkali and alkaline-earth halides. The optimized dispersion parameters are shown in Table 2. Figure 1 compares the densities of solids and liquids evaluated using the parameters together with experimental results.39 The deviations from the experimental values are at most 5 % for CaBr2 in solid states, and much better fitted in the other cases.

Simulated density of solids and melts of alkaline-earth halides compared with experimental results shown by solid lines: (a) fluoride, (b) chloride, and (c) bromide.
Using the force field parametrized as mentioned before, MD simulations were carried out for molten and superionic conductive states in order to evaluate transport properties of alkaline-earth halides. The molten states of alkaline-earth halides were prepared from initial configurations corresponding to their crystal structures. They were equilibrated in the NVT condition to reproduce molten states with the experimental density. The Nose-Hoover chain method suggested by Martyna et al. was employed as the thermostat with a coupling time of 20 ps.55 By contrast, for the superionic conductive phases of CaF2, SrF2, SrCl2, β-SrBr2, BaF2, and BaCl2, the fluorite crystal structure with experimental lattice parameters were relaxed under the ambient condition with the NPT ensemble. The isotropic barostat by Martyna et al. was also adopted here with a coupling time of 20 ps.55 After the careful equilibration for 2 ns, the temperature was gradually elevated up to the melting temperatures with each interval of 50 K for chlorides and bromides, and 50–200 K for fluorides. The orthorhombic or tetragonal crystals of CaCl2, CaBr2, α-SrBr2, and BaBr2 were also investigated in the same manner.
All the simulations to evaluate transport properties were performed under the NVT ensemble using the Nose-Hoover chain method. Total ion number N were set to 324 (216 anions and 108 cations) in the cell volume V determined by the experimental density at each temperature.34 The time step was set to 1.0 fs. After equilibration in each state, a series of 500 ps production runs were performed, until converged statistical averages could be obtained for the transport coefficients. As a result, 30 runs were performed for all the salts to evaluate the standard errors of transport coefficients, yielding a total production time of 15 ns for each of them.
The ionic conductivity σ, thermal conductivity λ and shear viscosity ηsv were evaluated by using the Green-Kubo formulae integrating the correlation functions,
| \begin{equation} \sigma = \frac{1}{T}L_{zz}, \end{equation} | (7) |
| \begin{equation} \lambda = \frac{1}{T^{2}}\left(L_{ee} - \frac{L_{ez}^{2}}{L_{zz}} \right), \end{equation} | (8) |
| \begin{equation} L_{\alpha \beta} = \frac{1}{kV} \int_{0}^{\infty} \langle j_{\alpha}^{a}(t)j_{\beta}^{a}(0) \rangle \mathrm{d}t,\quad (\alpha,\beta = e,z;\ a = x,y,z), \end{equation} | (9) |
| \begin{equation} \eta_{\text{sv}} = \frac{V}{kT} \int_{0}^{\infty} \langle \sigma^{ab}(t)\sigma^{ab}(0) \rangle \mathrm{d}t,\quad (a,b = x,y,z) \end{equation} | (10) |
where k is the Boltzmann constant, T the absolute temperature, jea and jza the a-component in the energy and charge currents, respectively, and σab the ab-component in the stress tensor. The microscopic expressions of the stress tensor and two currents are given elsewhere.52,56 Autocorrelation functions of the currents and the stress tensor were averaged over the three components (jαx, jαy and jαz) and the five components (σxy, σyz, σzx, σxx−yy and σ2zz−xx−yy), respectively, in order to hold good statistics. Since it is known that these transport properties obtained from MD calculations by using the Green-Kubo formulae depend only weakly on the cell size, in contrast to the diffusion coefficient,33,57–59 no correction was performed for the present results.
Figure 2 shows the temperature dependence of the transport coefficients using PIM at experimental densities corresponding to the temperatures in molten alkaline-earth halides. The behaviors between ionic conductivity and thermal conductivity were at first focused. Both relative magnitudes aligned in the same order of cations among the salts with the same anion such as CaF2, SrF2, and BaF2, and also in the same order of the halides such as fluorides, chlorides and bromides. These analogous behaviors between ionic conductivity and thermal conductivity can be explained by dimensional analysis of both properties, suggesting that both have the same dependence on ionic mass as mA−1/2 where mA denotes average mass of anion and cation, and have the same dependence on number density as (N/V)2/3, where N is the number of particles included in the system volume V.44,45 The temperature dependence of (N/V)2/3 is represented in Fig. 2d, in which the effect of ionic mass is originally excluded. These are the reasons why the analogous behaviors are observed in the cation series and the halides series.

Temperature dependence of (a) ionic conductivity, (b) thermal conductivity, (c) shear viscosity of molten alkaline-earth halides calculated with PIM. The error bars are provided with the standard error over the results of 30 runs, and then, the standard errors of ionic conductivity well converge inside the filled symbols. The temperature dependence of number density, which is common to all the transport properties, is shown in (d), and the temperature dependence of number density and packing fraction, which is specific to ionic conductivity, is shown only for fluorides in (e).
On the other hand, the temperature dependence indicated in Figs. 2a and 2b is totally opposite between both properties. Our previous study on the thermal conductivity of molten alkali halides shows that it is practically independent of temperature.60 The apparently negative temperature dependence shown in Fig. 2b is attributed to the number density dependence, i.e. the thermal expansion as shown in Fig. 2d. By contrast, according to our preliminary evaluation, the ionic conductivity is proportional to temperature. Furthermore, a strong negative packing-fraction dependence of ionic conductivity contributes positively to the temperature dependence, as η−2 where η is the packing fraction calculated using Shannon’s ionic radii,61 although the number density dependence of ionic conductivity contributes negatively to the temperature dependence: the former overwhelms the latter as totally shown in Fig. 2e. These are the reasons why they have apparently opposite temperature dependence.
The behavior of shear viscosity shown in Fig. 2c seems to be a little complicated. However, this also can be accounted for by dimensional analysis: shear viscosity depends on mA1/2 and (N/V)2/3, while ionic conductivity and thermal conductivity depend on mA−1/2 and (N/V)2/3. That is, whereas shear viscosity has the same number-density dependence as ionic conductivity and thermal conductivity, but it has a mass dependence inverse to those. These are the main reasons why the shear viscosity shows apparently complicated behavior.
Figure 3 compares the calculated ionic conductivity and shear viscosity including those for alkali halides33 with the experimental results.34–36 The deviations of the calculated values from the experimental values are depicted as correlations between ionic conductivity and shear viscosity. These indicate that the present PIM parameters for alkaline-earth halides have almost the same quality as those for alkali halides obtained by our previous study.33 While the Fumi-Tosi model62,63 gives the results with systematic and significant deviations from the experimental results in alkali halides,33 the PIM gives those with only small deviations scattered uniformly even in alkaline-earth halides as well as alkali halides,33 as shown in Fig. 3c.

Comparison of (a) ionic conductivity and (b) shear viscosity of molten alkaline-earth halides calculated with PIM together with the experimental values. Solid line is a criterion of which slope is unity. Open circles show the results for a series of molten alkali halides. The deviations of calculated values from experimental ones are shown in (c) as a correlation between both properties.
Contrary to these properties, there has been no available and reliable experimental data for the thermal conductivity of molten alkaline-earth halides. Thus, the present calculated thermal conductivities would be expected to be the most reliable data for them in light of the accuracy shown above for both the ionic conductivity and shear viscosity. Figure 4 compares the calculated thermal conductivity with the estimated values using two different empirical equations such as (a) the Gheribi equation64 and (b) our equation.60 Although the Gheribi equation was based on a finding that the thermal conductivity of alkali halides has practically no temperature dependence which was disclosed by our previous work,60 there seem to be significant differences between both the results using these equations. The slope of the Gheribi equation is about less than half of that for the calculated thermal conductivity of alkaline-earth halides. Since the Gheribi equation is given as
| \begin{equation} \lambda(T) = \lambda_{\text{fus}} + b(T - T_{\text{fus}}), \end{equation} | (11) |
where λfus and Tfus are the thermal conductivity and temperature at melting point, respectively, the coefficient b which is related with the thermal expansivity and Gruneisen parameter at melting point seems to be underestimated. Although our equation is simpler than Eq. 11, it expresses the present results for alkaline-earth halides within the errors and dispersions which are almost the same uncertainty for all the experimental, calculated, and estimated thermal conductivity for alkali halides.60,65–68

(a) Comparison of the calculated thermal conductivity of molten alkaline-earth halides by using PIM together with estimated values. Solid line is a criterion of which slope is unity. Open circles show the results for a series of molten alkali halides. (b) A plot of the calculated thermal conductivity vs. our scaling equation using mass, number density and packing fraction for alkali halides. Solid line is a guide to the eye.
Figures 5a and 5b illustrate the ionic conductivity and thermal conductivity calculated for fast ion conductive phases with some experimental values,35,37–42 respectively. The quantitative agreement with the experimental results for the ionic conductivity is fairly good, comparing the differences among the experimental results.37–39 The calculated ionic conductivity shows the opposite order in cation series between solid and liquid phases. Considering that the charge carriers in the fast ion conductive phase are fluoride ions, there should be no mass dependence in the phase. Furthermore, according to our preliminary evaluation, the dependence of ionic conductivity on packing fraction η is strongly negative: η−2. Figure 2e shows the temperature dependence of (N/V)2/3η−2. That is, in the fast ion conductive phase, the strong negative packing-fraction dependence dominates the cationic order in the magnitude of ionic conductivity of fluorides, whereas, in the molten states, the ionic mass dependence caused by cations inverts the order of ionic conductivity. The calculated thermal conductivity was in good agreement with the results of first-principles calculation,43 implying that the PIM reflects the DFT calculations well. On the other hand, it seems to be a little underestimated compared with the experimental results at ambient temperature.40–42

Temperature dependence of the calculated (a) ionic conductivity and (b) thermal conductivity in solid states of alkaline-earth fluorides.
The polarizable ion model has been parameterized for alkaline-earth halides, MX2, where M = Ca, Sr, and Ba, and X = F, Cl, and Br. The parameters were determined to reproduce the forces and dipole moments calculated using DFT. The calculated results of the ionic conductivity and shear viscosity were in good agreement with the experimental ones. The PIM enables us to predict the transport properties of alkaline-earth halides in the solid state as well as in the molten states. The calculated transport properties for both alkali and alkaline-earth halides could be well expressed by a scaling equation using ionic mass, number density, and packing fraction.
This study was partially supported by a Grant-in-Aid for Scientific Research (C) (Grant No. JP22K05033), a Grant-in-Aid for Scientific Research (B) (Grant No. JP22H01788), and JST ACT-X (Grant No. JPMJAX23D3). MD calculations were performed using the supercomputers of Research Center for Computational Science at Okazaki (22-IMS-C067 and 23-IMS-C059), Kyushu University (ITO), and Hokkaido University (Grand Chariot) partially through the HPCI System Research Project (project IDs: hp220068 and hp230083).
Yoshiki Ishii: Data curation (Lead), Formal analysis (Lead), Investigation (Lead), Methodology (Lead), Visualization (Lead), Writing – review & editing (Lead)
Sataro Kiko: Data curation (Supporting), Formal analysis (Supporting)
Norikazu Ohtori: Conceptualization (Lead), Supervision (Lead), Writing – original draft (Lead)
The authors declare no conflict of interest in the manuscript.
Japan Society for the Promotion of Science: JP22K05033
Japan Society for the Promotion of Science: JP22H01788
JST ACT-X: Grant No. JPMJAX23D3
N. Ohtori: Author to whom all correspondence should be addressed.
Y. Ishii and N. Ohtori: ECSJ Active Members