Abstract
At present, the midstory isolation (MSI) technology has great potential for application in historical buildings’ retrofitting and multifunction buildings. The coupling effect due to the variability of the location of the isolation layer may amplify the structural seismic response and is required for in-depth analysis. This paper aims to evaluate the magnitude of the coupling effect and delimitate the region of the coupling effect to be considered. Based on the complex mode superposition method, the explicit formulas for calculating the random response of the simplified model are deduced. The root-mean-square (RMS) ratio of the shear force coefficient of the upper isolation system is adopted as the performance indicator to evaluate the coupling amplification effect of the MSI system. Parameter analysis indicates that the coupling region is closely related to the mass ratio and frequency ratio of the upper and lower structures to the isolation layer. In general, the region of the coupling effect to be considered can be divided into two parts according to parameters of frequency ratios, depending on the thresholds of the performance indicator. As the mass ratio of the upper isolation system to the entire system increases, one of the coupling regions shrinks and eventually disappears, indicating that the coupling amplification effect in this region can be neglected under certain conditions. Finally, the time-domain analysis of three representative numerical cases of MSI buildings was performed to verify the reliability of the results obtained from the frequency-domain analysis. The research results can provide technical guidance for the preliminary design of the MSI buildings.
1. Introduction
The midstory isolation (MSI) technology is an extension of the base isolation technology and is proved to be effective in reducing seismic response [1–4]. The location of the isolation layer in MSI systems is not restricted. Collisions between the isolation layer and moat walls can be therefore be avoided. As the demand for maintenance of existing buildings increased, seismic isolation technology was applied to the reinforcement and renovation of existing buildings, which has contributed to the development of the MSI technology. Training Centers of Taisei Corporation, completed in 1997 [5], were reported as the first applications of retrofit with seismic midisolation method in Japan. Following this, the MSI technology began to be applied to high-rise multifunctional buildings; typical cases are 63 m Iidabashi First Building [6, 7], 120 m Shiodome Sumitomo Building [8, 9], 200 m Nakanoshima Festival Tower [10], and 230 m Roppongi Grand Tower [11]. By installing the upper seismic isolation system on the top of the existing buildings, the MSI technology is proved to be an effective strategy for adding the new function to existing buildings [12–15]. By inserting the isolation devices into the existing buildings, MSI technology has also been reported for the reinforcement of historic structures in accordance with the current design criteria [16]. Kim and Kang [17] installed magnetorheological (MR) dampers into the MSI system, which proved to be capable of reducing the displacement demand of the MSI system. Zhang et al. [18] incorporated the inerter system into the MSI system and verified its effectiveness in reducing the seismic response of the upper substructure and the deformation of the isolation layer. The relevant recent studies reported in [19–24] indicated that the incorporation of negative stiffness elements in isolation devices can also improve the seismic performance of isolated structural systems. In addition, a series of comprehensive studies on the optimal design of MSI systems have been conducted and reported in [25–29]. At present, due to its advantages in construction efficiency and adaptability to diverse needs, MSI technology has great potential for application in retrofitting and the reinforcement of historical buildings as well as in multifunctional buildings.
An MSI building consists of three segments: the lower structure, the isolation layer, and the upper structure. Due to the flexibility of the isolation layer, it was reported that the contribution of higher modes cannot be neglected, and the frequency of predominate modes may not be well separated, which in turn has an impact on the seismic response of the structure, defined as the modal coupling effect. Over the last few decades, a great deal of effort has been documented. Kobayashi et al. [30, 31] carried out a series of eigenvalue analyses and found that the response amplification on the upper structure of the isolated story was caused by modal coupling between modes of vibration on the upper and lower structures. Afterward, Wang et al. [32, 33] extended the validity of Kobayashi et al. [30, 31] methodology, adopting a 3DOF model to investigate the dynamic behavior of MSI structures. Results disclosed that seismic responses at the upper structure were enlarged due to the closeness of the two higher modal frequencies. A linear relationship between the frequency ratios of the upper and lower structures was put forward to predict the condition in which the modal coupling effect occurs. Based on the pole allocation method, Ikeda [34] proposed a closed-form expression to illustrate the tradeoff relationship between the damping ratios of the lower and upper substructures. Faiella and Mele [35, 36] utilized lumped mass multidegree-of-freedom models to identify the range of different behavioral modes and further verified that the modal coupling arises when the higher frequency of the upper structure is close to any frequency of the lower structure.
In the above studies related to the modal coupling effect, researchers mainly focused on the conditions in which modal coupling occurs and qualitatively illustrated the coupling amplification effect by comparing the peak transfer function or seismic response of the upper structure. Little has been documented for the magnitude of the modal coupling effect and the extent of its influence on the dynamic response of the structure, as well as the relationship between the modal coupling effect and the parameters of segments. Distinct from previous studies, this paper aims to evaluate the magnitude of the modal coupling effect with various design parameters and delimitate the region whether modal coupling should be considered based on the threshold values of performance indicator, which is essential for practical design. The explicit formulas for calculating the random response of the simplified model were derived by the residue theorem. At the same time, a series of dimensionless parameters were introduced for the parametric analysis of structural dynamic response. Then, the RMS ratio of the shear force coefficient of the upper isolation system was adopted as the performance indicator to evaluate the coupling amplification effect of the MSI system. Furthermore, for systems with different design parameters, the region of the modal coupling effect to be considered was delineated according to the performance indicator thresholds. Finally, the time-domain analysis of three representative numerical cases of MSI buildings was performed to verify the reliability of the results obtained from the frequency domain analysis.
2. Thend potential applications, the parameter ran MSI Systems
2.1. Model and Dynamic Equations
To gain insight into the model coupling effect of MSI systems, a simple linear three-degree-of-freedom (3-DOF) model is adopted for elementary analysis with reference to the 2-DOF model adopted in the base isolation system [32, 34, 37]. The MSI system consists of a lower structure, isolation layer, and upper structure, as shown in Figure 1(a). The lower structure can be simplified as an S-DOF system shown in Figure 1(b). With lower structure infinitely rigid, the 3-DOF model of the MSI system degenerates to the 2-DOF model of the base isolation system, as shown in Figure 1(c). The parameters of the models, including mass m, damping c, and stiffness k, are shown in Figure 1, and the subscripts ls, iso, and us denote the lower structure below the isolation layer, the isolation layer, and the upper structure above the isolation layer, respectively.

(a)

(b)

(c)
The dynamic equations of the 3-DOF system subjected to the ground excitation can be expressed as follows:where xls, xiso, and xus, denote the relative displacement of the three masses in 3-DOF. The dot on xls, xiso, and xus denote the differential operations with respect to the time t. For simplicity of derivation, the following parameters are defined:where r1 is the mass ratio of the upper isolation system to the whole system, reflecting the location of the isolation layer in the whole system. r2 is the mass ratio of the upper structure to the upper isolation system. ξls and ξus are the damping coefficient of the lower structure and the upper structure, respectively, and ξiso is the equivalent viscous damping ratio of the upper isolation system. εls and εus are frequency ratios between the different segments.
By introducing the above-defined parameters into equation (1), the dynamic equations for the isolation system can be rewritten as follows:where M, C, and K are the mass, damping, and stiffness matrices of the system, respectively, and E is the vector that each degree of freedom to the ground motion, and x is the inner-story relative displacement of different segments. The M, C, K, and E matrices related to the 3-DOF system are shown in the following:
Similarly, for the 2-DOF system and S-DOF system, the dynamic equations can also be written in the form of equation (3). The related matrices for the 2-DOF system are expressed as follows:and for the S-DOF system,
According to the above analysis, the relevant parameters of the structural dynamic equation (3) are the nondimensional parameters εls, εus, r1, and r2; the damping ratios ξls, ξiso, and ξus; and the frequency of upper isolation system ωiso. Accordingly, the dynamic characteristics and response of the system can be expressed as functions of the above parameters. The damping ratio of the upper and lower structure is generally a given constant that varies with the type of structure. To cover the range of practical applications [5–10, 36] and potential applications, the parameter ranges for related variables are defined as follows:
2.2. Dynamic Characteristics
Dynamic characteristics of 3-DOF MSI systems are discussed in this section. The frequencies of the 3-DOF systems are defined as ω1, ω2, and ω3, and the modal participating mass ratios are defined as L1, L2, and L3. According to the basic theory of dynamics, the i-th (i = 1, 2, 3) modal participating mass ratios can be calculated as follows:where , M and E can be obtained with equation (4), and is the i-th mode shape of the 3-DOF system.
Figures 2(a)–2(c) summarize the ratios of the fundamental frequencies to the upper isolation system and ratios of adjacent mode frequencies, denoted as ω1/ωiso, ω1/ω2, and ω2/ω3, respectively. From Figure 2(a), it can be seen that the frequency ratio ω1/ωiso increases with the increase of the frequency ratio εls, and εus. The first natural frequency of the 3-DOF model is close to the first natural frequency of the upper isolation system. As shown in Figure 2(b), with the increase of εls, the frequency ratio ω1/ω2 decreases after the increase in the first time and then tends to stabilize. The maximum value of ω1/ω2 is about 0.3, which validates the first and second frequencies are far apart. On contrary, as shown in Figure 2(c), it is apparent that values on the peak line of ω2/ω3 vary in a very narrow range with values close to 1, indicating that the modal coupling occurs between the second and third modes. Further numerical analysis shows that the frequency ratio corresponding to the peak line of ω2/ω3 satisfies the following relationship:which is proved to be consistent with the equation proposed in [32]. Interestingly, based on Kelly’s work [37], for the 2-DOF system shown in Figure 1(b), the approximation of second frequency is given by

(a)

(b)

(c)
On the other hand, for the S-DOF system shown in Figure 1(c), we have
Comparing equations (9)–(11), it can be concluded that the model coupling occurs when the second frequency of the upper isolation system is close to the first frequency of the lower structure, which is similar to the tuning effect of TMD.
Substituting equation (9) into equation (3), the second frequency of the 3-DOF system can be obtained as , and the modal shape vectors can be expressed as [1 0–1/r2]T. Correspondingly, the second modal participating mass ratio defined in equation (8) can be calculated as follows:
Accordingly, the modal participating mass ratios L1 and L3 can be approximated as follows:
As indicated in equation (12), it is apparent that L2 increases smoothly with increasing mass ratio r2 and decreases with increasing mass ratio r1. Moreover, the above analysis clearly shows that when modal coupling occurs, the vibration of the structure is characterized by the first and the second modal vibration, where the second modal participating mass ratio affects the magnitude of the coupling effect. When coupling occurs, the larger the L2, the more significant the coupling effect. As shown above, equations (12) and (13) are calculated with the condition of equation (9); it is worthwhile to investigate the influence of related parameters on the modal participating mass ratios under general conditions.
The modal analysis indicates that the first mode of the 3-DOF is characterized as the first-order vibration of the upper isolation system, and the second and third orders are characterized as the second mode of the upper isolation system or the first mode of the lower structure, depending on the value of the frequency ratio. The variation of modal participating mass ratios also reflects the transition of the dominant modes. Figure 3 presents the comparison of the modal participating mass ratios L1, L2, and L3 corresponding to the variation of the mass ratios r1 and r2 and the frequency ratio εls, when the frequency ratio εus = 4, for example. As depicted in Figure 3(a), it can be seen that the first modal participating mass ratio L1 may be significantly affected by the mass ratio r1. The increase in r1 resulted in the increase of L1, and L1 gradually decreases and approaches r1 with the increase of the frequency ratio εls. As for L2 and L3, it is apparent that at εls = 8.95, the effective mode switches from the second to the third mode. Moreover, L2 + L3 increases slightly and approaches to (1 – r1) as εls increases. It is reasonable to assume that the mass ratio r1 plays an important role in the dynamic response of the MSI system.

(a)

(b)
As shown in Figure 3(b), the first modal participating mass ratio L1 remains stable as r2 varies. From the comparison of L2 and L3 in Figure 3(b), it can be seen that r2 has a great influence on the values corresponding to the switching points. The frequency ratios εls for the effective modal switching point are 5.65, 7.30, and 12.65, corresponding to r2 of 0.5, 0.7, and 0.9, respectively. In summary, it should be noted that both mass ratios r1 and r2 have an evident influence on modal participating mass ratios: r1 affects the stable value of the first mode, while r2 affects the position of the switching point between the second and third modes. The application of the MSI technique causes the fundamental change in the basic dynamic characteristics of the system.
3. Modal Coupling Effect of MSI Systems
3.1. The Absolute Acceleration Response of n-DOF (n = 1, 2, and 3) Systems
The MSI technique protects the upper structure by mitigating the vibration induced by the lower structure. Meanwhile, the vibration of the lower structure decreases with the building mass damper effect induced by the upper isolation system. The vibration mitigating effect can be visually reflected in the acceleration of different segments. Therefore, the absolute acceleration is of importance and is therefore analyzed in this study.
Substituting the base acceleration in equation (3), the absolute acceleration transfer function of all masses in the n-DOF (n = 1, 2, 3) systems can be derived from equation (3). Moreover, according to the integral formula for random vibration analysis [38], the variance of acceleration response σ2 of the n-DOF system under white-noise random excitation (with the constant power spectral density S0), can be described as follows:where
H(ω) is the frequency response function of the absolute acceleration for different segments, which can be calculated according to Laplace transformation of equation (3). ε is given by ε = ω/ωiso, which represents the frequency ratios of excitation to the upper isolation system. B0-Bm-1 and A0-Am are time-invariant real constants. Taking absolute acceleration of the upper mass in the 3-DOF model as an example, B0-Bm-1 and A0-Am are given as follows:
Substituting equations (16) and (17) in the resultant equations in Appendix, RMS of the absolute acceleration σ for the upper mass in the 3-DOF model can be obtained as an explicit function of the design parameters r1, r2, εls, εus, ξls, ξiso, and ξus. The calculations required for random vibration analysis can be significantly simplified with equation (14) and the equations in Appendix. Similar expressions can be obtained for the lower mass and isolation mass in the 3-DOF model and the corresponding mass in 2-DOF and S-DOF models.
3.2. Performance Indicator for the Modal Coupling Effect
Based on the effectiveness of the n-DOF (n = 1, 2, 3) models in capturing the acceleration response for MSI systems, base isolation systems, and fixed-base systems, the coupling amplification effect can be conducted in the frequency domain. In this section, comparisons in the response of the MSI systems and the evaluation of the coupling effect are conducted.
To analyze the absolute acceleration amplification of the lower mass, isolation mass, and upper mass in the 3-DOF model, ratios are defined as follows:where σ is the RMS value of the absolute acceleration response, which can be calculated with equation (14). The subscripts III, II, and I represent 3-DOF, 2-DOF, and S-DOF systems, respectively.
Figures 4(a)–4(c) illustrate the RMS ratios Rals, Raiso, and Raus of the absolute acceleration with respect to different frequency ratios εls and εus (r1 = 0.8, r2 = 0.7, and ξiso = 0.20). It can be observed that Rals for the lower mass varies between 0.3 and 0.8, and Raiso and Raus for the isolation mass and upper mass vary between 1.0 and 1.6. This indicates that the response of the lower structure can be significantly reduced compared to the corresponding fixed-base structure with the installation of the upper isolation system. On the contrary, due to the flexibility of the lower structure, the response of the upper isolation system is amplified compared with the base isolation system. Furthermore, it is worth noting that there are apparent peak lines in Figures 4(a)–4(c), and it is verified that the relationship of frequency ratios corresponding to the peak lines is approximately consistent with equation (9).

(a)

(b)

(c)
This study focuses on the modal coupling amplification effect of the MSI system caused by the flexible isolation layer. The above analysis indicates that the absolute acceleration amplification effect occurs in the upper mass and isolation mass. Therefore, the RMS ratio Raui is adopted as the performance indicator to analyze the coupling amplification effect, defined as follows:where aui is the shear force coefficient of the upper isolation system, defined as follows:where σ is the RMS of the shear force coefficient, which can be calculated with equations (14) and (20). The subscripts III and II represent 3-DOF and 2-DOF systems, respectively. The design conditions are defined as follows:where [RA] is the design threshold, which can be determined according to the seismic design performance requirements. To examine the influence of the upper and lower structures on the dynamic performance of the 3-DOF MSI systems, the boundary curve with [RA] = 1.2 is adopted for illustration. The contours of variation of Raui with frequency ratios εls and εus are plotted in Figures 5 and 6 for r2 = 0.7 and ξiso = 0.20 and r2 = 0.9 and ξiso = 0.20, respectively.

(a)

(b)

(c)

(d)

(a)

(b)
For simplicity thereafter, the regions delineated in Figures 5 and 6 are summarized in Figure 7. According to the given [RA], the relevant parameter range can be divided into four regions A–D, as shown in Figure 7. Specifically, there is a peak line in region C, which is consistent with Figure 5. Regions A and C with amplification factor Raui > [RA] are defined as the coupling regions, while regions B and D with Raui < [RA] can be defined as design regions. As depicted in Figures 5 and 6, the amplification patterns of regions A and C are completely different. Combined with the dynamic characteristics analysis in Section 2.2, it can be inferred that for region A, the amplification effect is caused by modal coupling between first and second modes, while for region C, it is caused by modal coupling between second and third modes.

As shown in Figure 5, in the four cases, the scope of region A remains stable and can be expressed as εls < 3 approximately. The value of Raui in region A tends to decrease rapidly with the increase of εls and remains essentially constant with εus. The scope of coupling region C gradually decreases with the mass ratio r1 varies from 0.2 to 0.7. In particular, when r1 is 0.8, region C disappears, and regions B and D are combined. This indicates that the coupling amplification effect caused by the coupling between second and third modes is gradually weakening with the increase of mass ratio r1. In addition, the peak line in region C is approximately parallel to the boundary line of the region and is independent of the frequency of the upper isolation system ωiso. The results show that the seismic response of the upper isolation system can be reduced by slightly adjusting the frequency ratios εls and εus according to the characteristics of different regions.
The comparison of Figures 5(b) and 6(a) shows that as r2 increases, the scope of region C shrinks rapidly, and the values on the peak line decrease accordingly. As for the case r1 = 0.7 and r2 = 0.9, region C disappears as shown in Figure 6(b), indicating that the coupling effect in region C can be ignored with the increase of both mass ratios r1 and r2. In addition, the comparison of the peak line in region C in Figures 5 and 6 reveals that the slope of the peak line decreases with the increase of r2. Further numerical analysis indicates that the peak line is parallel to the line corresponding to equation (9), which can be expressed as follows:where d is the deviation between results of the peak line obtained from response analysis of equations (14) and (19) and the peak line of ω2/ω3 deduced from eigenvalue analysis. For different cases, the values of d are shown in Table 1. It can be observed that d varies between 0.18 and 0.26 and increases with the increase of r1 and r2, indicating that the peak line of Raui deviates from the line obtained from the nominal frequency relation, while equation (9) can only approximately predict the peak line.
To investigate the influence of damping ratio ξiso on the response of the upper isolation system, the variation surface of Raui as a function of damping ratio of the upper isolation system ξiso and the frequency ratio εls for εus = 4 is plotted in Figure 8. It is seen that when εls is greater than 5, the ratio Raui increases with the increase of damping ratio, which means the decrease of the isolation effect.

Figure 9 shows the contour of the system with ξiso = 0.30 versus the frequency ratios εls and εus. Compared with the system with ξiso = 0.20 shown in Figure 6(b), the slope of the peak line remains the same while the position is shifted upward. Further analysis shows that as the damping ratio ξiso increases from 0.2 to 0.3, the value of Raui decreases in region A and increases in region C. Correspondingly, the scope of region C increases.

4. Application to Case Buildings
For application, two actual cases, the Iidabashi First Building [9, 36] (denoted as Case 1) and the Shiodome Sumitomo Building [9, 36] (denoted as Case 2), and one calculation case (denoted as Case 3) are selected as cases to investigate the coupling effect of the MSI systems with frequency ratios in different regions. Case 3 is a numerical calculation model of a ten-story building designed so that the structural parameters are located on the peak line of Raui. The damping ratios of the upper structure and lower structure ξls and ξus for the three cases are 0.02. The 3-DOF model for MSI structure is utilized, and the related original and calculated parameters for the three cases are shown in Table 2.
4.1. Frequency Response Analysis
The contours of the performance indicator Raui derived from the given mass ratios r1 and r2 and damping ratio ξiso, as well as the location of the design parameters adopted in practical applications (denoted as Els and Eus, marked by the red circle), are shown in Figures 10–12 for Cases 1–3, respectively. It is seen that the design parameters of Cases 1–3 correspond to different regions of coupling effect. According to the delineation of the regions in Figure 7, it can be seen that Case 1 is located in the coupling region A, Case 2 in the intersection of regions A and C, and Case 3 on the peak line of the coupling region C.



Therefore, in order to compare the coupling effect in different regions, the frequency ratio εls is selected as the variable for Case 1 and the frequency ratio εus for Cases 2 and 3, as shown in Table 3.
The RMS ratios Ra for the upper mass, the isolation mass, and the lower mass of Case 3 are shown in Figure 13. As depicted, for different segments of the structure, the peak response is located at εus = Eus. It should be noted that in general, the RMS ratios Ra for the upper mass and isolation mass are greater than 1, while Ra for the lower mass are less than 1, and all of them gradually approach 1 as the frequency ratio εls increases, indicating that the response of the MSI structure is consistent with that of the base-isolated structure when the lower structure is sufficiently rigid.

4.2. Time-History Analysis
To further confirm the validity of the obtained analytical results, a series of time history analyses subjected to earthquake excitations were carried out for the three cases. Three real earthquake records were selected as input ground motions: the EW component of El Centro 1940, the NS component of Taft 1952, and the NS component of Hachinohe 1968. The magnification factor M, defined as ratios of maximum acceleration responses to peak ground acceleration (PGA), is shown in Figures 14–16.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)
As shown in Figures 14–16, for Case 1, located in the coupling region A, the magnification factor M decreases for the upper and isolated masses and increases for the lower mass as the frequency ratio εls increases. For Case 2, located in the intersection of regions A and C, the increase of εus above the peak line in region C leads to an increase of the magnification factor M of the upper and isolated masses. For Case 3, located in the coupling region C, the maximum magnification factor M for all masses is obtained when εls = Eus, which indicates that either increasing or decreasing the stiffness of the lower structure here can reduce the structural response. These results are consistent with those obtained from the frequency domain analysis.
5. Conclusion
Based on random vibration analysis, this study proposed a method to evaluate the coupling effect of the MSI systems. Explicit dynamic response formulas are derived for the MSI systems (represented by the F model), the base isolation systems (represented by the 2DOF model), and the fixed-base systems (represented by the S-DOF model) under white-noise excitation. The RMS ratio of the shear force coefficient of the upper isolation system is adopted as the performance indicator to evaluate the coupling amplification effect of the MSI system. In addition, the coupling effect region corresponding to a certain threshold can be obtained with different structural parameters. Finally, the time-domain analysis was performed to verify the reliability of the results obtained from the frequency domain analysis. Some conclusions can be drawn as follows:(1)The frequency ratios εls and εus and mass ratios r1 and r2 have a significant effect on the dynamic characteristics and dynamic response of the MSI systems. The mass ratio r1 mainly affects the distribution of contribution between the first and higher modes, while the mass ratio r2 affects the distribution of contribution between the second and the third modes.(2)The coupling effect between adjacent modes amplifies the seismic response of the upper isolation system compared to the corresponding base isolation system. In contrast, the coupling effect effectively reduces the seismic response of the lower structure compared to that of the fixed-base system.(3)In general, the region where the coupling effect should be considered can be divided into two parts, regions A and C, according to the thresholds of the performance indicator. As mass ratios r1 increases, region C shrinks and eventually disappears, indicating that the coupling amplification effect in this region is negligible under certain conditions.(4)In region A, where the amplification effect is caused by modal coupling between the first and second modes, the effect of frequency ratio εls is more significant than frequency ratio εus. The performance indicator Raui decreases rapidly with increasing εls and remains essentially constant with increasing εus. In region C, where the amplification effect is caused by modal coupling between second and third modes, there is a peak line on which the maximum of indicator Raui can be obtained. The peak line in region C is approximately parallel to the boundary line of the region. The results show that the seismic response of the upper isolation system can be reduced by slightly adjusting the frequency ratios εls and εus according to the characteristics of different regions.(5)The effect of the damping ratio on the amplification effect varies in different regions and needs to be determined according to the frequency ratios εls and εus and mass ratios r1 and r2.
The proposed approach does not require detailed dynamic simulations of the system components and can be applied in the preliminary design of the MSI structures.
Appendix
where
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The research presented herein was sponsored by the Japan Society for the Promotion of Science (JSPS; Grant no. P19777), the Japan Society for the Promotion of Science (Kakenhi; no. 18K04438), and the Tohoku Institute of Technology Research Grant. The financial support is gratefully acknowledged.