Thermal conductivity minimum of graded superlattices due to phonon localization

Heat conduction engineering in semiconductor microstructures and nanostructures has been widely investigated in the past decades based on the particle picture of phonons.^1–5 As phonons are essentially wave-packets with finite width, the wave picture provides another important degree to engineer heat conduction when the characteristic size of nanostructures becomes comparable to the coherence length of phonons.^6–9 In recent years, there are mainly two categories of nanostructures proposed for the heat conduction tuning via the wave nature of phonons: nano-phononic meta-materials (pillared thin films,^10,11 nanowire cages,¹² etc.) and nano-phononic crystals [superlattices (SLs),^13–17 periodic porous nanostructures,^18–20 etc.] based on resonance and band-folding, respectively.

Localization is a well-known wave phenomenon that significantly impedes transport since the pioneering work of Anderson.²¹ It denotes the absence of wave diffusion in a disordered medium due to the coherent backscattering and destructive interference²² and has been widely studied with electrons,²³ electromagnetic waves,²⁴ and elastic waves.²⁵ The localization of thermal phonons, induced by the disorder of bond strength and/or atomic mass, actually widely exists in amorphous materials²⁶ and alloys.^27,28 However, the localization effect on heat transport is invisible or unknown in these systems.

The phonon transmission and localization in disordered SL systems were studied based on elastic continuum models in the low-frequency limit.^29–31 The influence of phonon localization on thermal transport properties has been more extensively explored by atomistic simulation methods via either the ballistic non-equilibrium Green’s function (NEGF)^32–35 or molecular dynamics (MD).^36–40 The localization effect was shown to be not observable in the thermal conductivity of isotope-disordered nanotubes due to the small contribution of high-frequency localized phonons at temperatures low enough to preserve phonon coherence from anharmonic scattering.^32,33 Important progress was made to explicitly demonstrate the phonon localization by the thermal conductivity maximum of periodic SLs with randomly embedded nanodots,³⁴ which has been experimentally observed recently.⁴¹ The decrease in thermal conductivity with total length is attributed to the exponential decay of phonon transmission to zero,^32–34 which is more directly evidenced in a recent modal-level NEGF study of randomly aperiodic SLs.³⁵ The phonon localization in aperiodic SLs in the presence of anharmonicity and interfacial mixing is thoroughly studied by MD simulation,³⁸ where the thermal conductivity maximum appears at a sufficient level of periodicity disorder. There are also some MD studies of phonon localization in random multilayer systems,^36,40 disordered porous nano-phononic crystals,³⁷ and heterogeneous interfaces with atomic inter-diffusion.³⁹ However, the thermal conductivity maximum is only obtained in the coherent part,^36,37 while any measurable evidence of phonon localization still remains ambiguous in these systems.

To sum up, the previous works on phonon localization in heat conduction are mainly featured by (i) random disorder, (ii) exponential decay of phonon transmission to zero with increasing system length, and (iii) a thermal conductivity maximum with length. Feature (ii) is a signature of strong localization⁴² caused by feature (i) and finally leads to feature (iii).³⁴ The graded systems with deterministic disorder have been considered for the manipulation of heat transport in some recent works.^43–45 In this work, we demonstrate phonon localization in graded SLs through a ballistic NEGF simulation. Although the configuration of atoms in graded SLs is deterministic as it follows a specific mathematical rule, the gradual variation of the period length represents similar to a moderate disorder in effect. Therefore, we obtain a phonon localization with very special features: (i) short-range order and long-range disorder, (ii) exponential decay of phonon transmission with system length to a non-zero constant, and (iii) a thermal conductivity minimum with length. This kind of localization is not strong as the level of disorder is not sufficiently high. As a result, the phonon transmission does not fully decay to zero with increasing system length. The decay of transmission to a non-zero constant was actually shown in few systems recently.^37,38 In this work, we aim to draw attention to such partial localization³⁸ of broadband phonons, which was seldom considered in heat conduction.

The remainder of the article is organized as follows. The physical model of graded Si/Ge SLs and the corresponding NEGF methodology will be introduced in Sec. II. The results and discussions of phonon localization in the graded SL system will be given in Sec. III. Finally, the concluding remarks will be made in Sec. IV.

In this section, first, we explain the detailed configurations of graded Si/Ge SLs in Sec. II A. In Sec. II B, we introduce the ballistic phonon NEGF formalism, together with the calculation of several important variables from the Green’s functions.

We design a series of graded Si/Ge SLs as exemplified in Fig. 1(a), with reference periodic Si/Ge SLs shown in Fig. 1(b). In both the graded and periodic SLs, the interfaces are atomically smooth. The graded SLs are made up of several different SL units, with each SL unit containing the same number (Np) of periods of each period length. In the example shown in Fig. 1(a), there are three SL units with the period lengths p₁ = 1 uc, p₂ = 2 uc, and p₃ = 3 uc and the number of periods Np = 2 in each SL unit. Here, 1 uc denotes the length of one conventional unit cell of Si. Other configurations with more SL units are given in Table I. The reference periodic SLs in Fig. 1(b) have a period length of p = 1 uc. The present graded SLs have both short-range order and long-range disorder. With the increasing number of SL units, the level of disorder elevates. For graded SLs with a specific number of SL units, the total system length L increases as Np increases. For simplicity, we assume the lattice constant and force constant of Ge to be the same as those of Si, with only the atomic mass difference taken into account.^46,47 All the considered period lengths lie within a range of about 1–4 nm, which are comparable to the dominant phonon coherence length in Si/Ge SLs^38,48 and ensure the coherent condition for phonon localization.

As phonon localization in heat conduction is coherent in nature, we adopt the ballistic NEGF formalism,^49–51 with the definition of real-time retarded Green’s function as

where “i” and 2πℏ are the imaginary unit and the Planck constant, respectively, θ(τ) is a step function of time, and $unατ=rnατMn$ with $rnατ$ being the atomic displacement of atom n along the α direction at time τ and M_n being the atomic mass. The brackets “[]” and “⟨⟩” in Eq. (1) denote the boson commutator and the ensemble average, respectively. For the steady-state heat transport in the device region with transverse periodicity in this work, the retarded Green’s function is computed in matrix notation by

with I being the unity matrix, $ω;q⊥$ denoting the frequency and wave vector dependences along the transport and the periodic transverse direction, respectively, and $Φ̃(q⊥)$ being Fourier’s representation of the harmonic force constant matrix.^47,51 The retarded self-energy matrix only includes the contribution from two contacts: $ΣRω;q⊥=Σ1Rω;q⊥+Σ2Rω;q⊥$⁠, which are related to the surface Green’s function of contacts computed by the decimation technique.⁵² The greater/lesser Green’s function is computed by^47,53

where the advanced Green’s function G^A is the Hermitian conjugate of G^R, and the greater/lesser self-energy matrix also only includes the contribution from the two contacts: $Σ>,<ω;q⊥=Σ1>,<ω;q⊥+Σ2>,<ω;q⊥$⁠, which are related to the retarded contact self-energy matrices and Bose–Einstein equilibrium distributions.⁴⁷ Note that in ballistic NEGF, it is usually not necessary to calculate G^>,<, which are required instead in the present work to compute the local phonon number density as introduced later.

Once the Green’s functions are resolved, the transmission through the SLs is obtained based on the Caroli formula,^49–51 

where N is the number of transverse wave vectors, “Tr” denotes the trace of a square matrix, and the broadening matrix is defined as $Γ1(2)=iΣ1(2)R−Σ1(2)A$⁠, with the advanced contact self-energy matrix Σ^A being the Hermitian conjugate of Σ^R. The thermal conductance of the SLs is calculated based on Landauer’s formula as^49–51 

where s is the area of transverse cross section and f_BE is the Bose–Einstein distribution. The thermal conductivity of the SLs is related to its thermal conductance as κ = σL. We also introduce the normalized thermal conductivity accumulation function vs frequency,

which provides a direct visualization of the frequency range where phonon localization takes place.

The local phonon number density is related to the diagonal blocks of G^< as^47,54

where R_n is the position of atomic site n and ɛ ≡ ω² is the eigenvalue of the harmonic force constant matrix, with the relation $ρω,Rn=ρε,Rnω/π$⁠. The local phonon number density provides an intuitive real-space representation of phonon localization: The phonon number density of localized states is concentrated around some atomic sites, as to be shown later. The local density of states (LDOS) is defined as^47,54

with the spectral function matrix computed by

In addition, we have the relation for LDOS as $LDOSω,Rn=LDOSε,Rnω/π$⁠. The overall DOS of the system is calculated as

with N_a being the total number of atoms in the system. Once we obtain the LDOS throughout the system, the spectral participation ratio (PR) is introduced to characterize the level of phonon localization,^55,56

To account for the influence of the number of states at a specific frequency, we also introduce a weighted PR defined as the product of PR and DOS.⁵⁶ 

In terms of numerical implementation, we adopt a recursive algorithm to solve the retarded and greater/lesser Green’s functions in Eqs. (2) and (3) within a massively parallelized framework.^47,57 The harmonic force constants of Si are obtained from the first-principles calculation, with the details given in our previous work.⁴⁷ The optimized size of one unit cell of Si is 1 uc = 5.4018 Å. A mesh of N_m = 200 for the frequency points and a mesh of N = 20 × 20 for the transverse wave vector points are adopted after careful independence check.

In this section, first, we demonstrate the appearance of a thermal conductivity minimum with the length of the graded Si/Ge SLs due to phonon localization in Sec. III A. A detailed analysis of the length-dependent transmission and participation ratio is presented as the evidence of phonon localization in Sec. III B. Finally, the real-space localization pattern is shown via the phonon number density distribution.

First, the results of thermal conductance and conductivity of periodic Si/Ge SLs with p = 1 uc at 200 K are given in Fig. 2 as a reference for later discussion. The reason of considering the temperature of 200 K will be given later. With increasing number of periods, the thermal conductance decreases and rapidly converges to a constant value after about 4–5 periods, as shown in Fig. 2(a). This is due to the band-folding process in SLs where several periods are required for the reflections and interferences of lattice waves to form the SL eigenmodes.⁴⁶ Once the folded modes are well formed, they transport ballistically across the SLs as in a homogeneous medium with a conductance independent of the system length.¹⁴ The thermal conductivity of periodic SLs increases almost linearly with the system length, as shown in Fig. 2(b), which also indicates ballistic heat transport. The present trend of thermal conductance and conductivity for the periodic SLs with smooth interfaces agrees quite well with that in previous NEGF studies.^46,58

Then, we discuss the results of thermal conductance and conductivity of graded Si/Ge SLs at 200 K in Fig. 3. For all the graded SLs with different numbers of SL units, generally the thermal conductance rapidly decreases with total length and gradually converges at Np ≥ 4–5, as shown in Fig. 3(a). The decrease in thermal conductance is more significant compared to the periodic SL case. With the increasing number of SL units, the thermal conductance converges to lower values. This is caused by the elevating level of long-range disorder in the graded SLs, which induces elastic backscattering and destructive interference of phonon wave-packets. As the number of SL units increases from 3 to 5, the elastic backscattering becomes stronger to approach the quantum diffusion limit (σ ∼ 1/L)^21,22 at the initial stage (Np ≤ 4–5), as indicated in the double-log-scale inset of Fig. 3(a). As a result, the thermal conductivity of the graded SLs transits from an increasing trend with length (the case of three and four SL units) to a length-independent regime (at Np ≤ 4–5 in the case of five SL units), as it is demonstrated in Fig. 3(b). As the number of SL units further increases, the elastic backscattering is so strong that the decrease in thermal conductance with length at the initial stage is faster than the quantum diffusion limit, i.e., σ ∼ 1/L^α with α > 1. Thus, the thermal conductivity even decreases with length (at Np ≤ 4–5 in the case of six and seven SL units), as shown in Fig. 3(b), which is a signature of phonon localization. For all the cases, the thermal conductivity starts to increase with length at Np ≥ 4–5 since the thermal conductance gradually converges. Therefore, a thermal conductivity minimum with length appears at Np = 4–5 for the graded SLs with six and seven SL units. The minimum shall come from a competition between the phonon localization and ballistic transport. As the level of disorder in the graded SLs is moderate compared to the random disorder (e.g., in aperiodic SLs^35,38), a strong localization⁴² cannot be reached and the ballistic transport wins out at Np ≥ 4–5. As a comparison, the thermal conductivity maximum reported in the strong localization regime^34,35,38,41 results from the ballistic to localized transition as the system length increases to be comparable to and larger than the phonon localization length. The present minimum is also different from the previous thermal conductivity minimum^13,15–17 due to the coherent-to-incoherent transition when increasing the period length of SLs under a fixed total length. The heat transport in graded SLs is somehow in a partial localization regime, which was introduced for the high-frequency phonons in aperiodic SLs due to decoherence as the period length is comparable to or even larger than the coherence length of those phonons.³⁸ We will provide deeper analysis and more direct evidence of localization via the phonon PR, transmission, and number density distribution in Subsections III B and III C.

The temperature dependence of the thermal conductance and conductivity of the graded Si/Ge SLs with six SL units is plotted in Figs. 4(a) and 4(b), respectively. As shown in Fig. 4(a), the thermal conductance decreases with decreasing temperature. The dip of thermal conductivity minimum with length becomes smaller as the temperature is lower, and fully vanishes at 10 K, as shown in Fig. 4(b). This indicates that the phonon localization is mainly contributed by modes in the moderate- and/or high-frequency range, as to be uncovered later. To ensure the population of these localized phonon modes, the system temperature cannot be too low. On the other hand, the system temperature shall be sufficiently low to suppress the anharmonic phonon–phonon scattering, which will destroy the phonon coherence. To sum up, to observe phonon localization, the system temperature can be neither too low nor too high, which is similar to the condition in the periodic SLs with embedded nanodots.^34,41 To estimate the effect of anharmonic scattering missing in our ballistic NEGF simulation, we calculate the phonon mean free path (MFP) distribution of bulk Si at different temperatures by the first-principles method.⁵⁹ As shown by the normalized thermal conductivity accumulation function vs MFP in Fig. 4(c), the phonons with MFP < 100 nm contribute to only ∼10% and nearly 0% of the total thermal conductivity at 200 and 100 K, respectively. Considering that the thermal conductivity minimum takes place at a system length of <100 nm, the temperature range of 100–200 K with negligible anharmonic scattering shall be good for the observation of phonon localization. Throughout the present work, we mainly consider and discuss the simulation results at 200 K.

Finally, to quantify the effect of phonon localization, we provide a comparison of the thermal conductance of graded SLs to that of the reference periodic SLs with p = 1 uc in Table II. The length of the graded SLs is set to Np = 4 where the thermal conductivity minimum takes place. The length of the periodic SLs is the same as that of graded SLs with each specific number of SL units. The thermal conductance of the periodic SLs is shown to be almost independent of the total length in the present range, which is consistent with the trend in Fig. 2(a). In contrast, the thermal conductance of graded SLs is reduced a lot due to phonon localization. The reduction ratio increases with the number of SL units and reaches as high as ∼96% at seven SL units. Therefore, the graded SL is a good system to achieve low thermal conductivity. Although its performance might be slightly lower than that of the random multilayer structures³⁶ or optimized aperiodic SLs,⁵⁶ the graded SLs seem to be more easily fabricated experimentally due to the deterministic mathematical rule.

To illustrate the phonon localization, we show the length-dependent PR, weighted PR, and transmission of the graded SLs with six SL units in Fig. 5. The PR, weighted PR, and transmission all decrease with increasing length and converge after L ∼ 84 uc (Np ∼ 4). In comparison to the reference periodic SLs with p = 1 uc, the PR of the moderate-frequency (3–7 THz) and high-frequency (>∼10 THz) phonons in graded SLs is reduced considerably, which is the evidence of localization. The reduction in PR of very low-frequency phonons (<∼1 THz) is relatively smaller, indicating ballistic heat transport. The convergence length of Np ∼ 4 is well consistent with the length where the thermal conductance converges and the thermal conductivity minimum takes place in Sec. III A. The convergence length of Np ∼ 4 represents the number of periods that is required in each SL unit to well form the corresponding folded eigenmodes. Once the folded modes in each SL unit are well formed and interfere destructively, the extent of localization will converge at Np > 4–5. In the SL units of the graded SLs, there are confined (localized) high-frequency optical modes due to the spectrum mismatch.⁶⁰ To distinguish this kind of localization and that induced by the destructive interference of folded eigenmodes among different SL units, we compare the PR, weighted PR, and transmission of six periodic SLs with each period length (p = 1–6 uc) in the graded SLs to those of the graded SLs in Fig. 6. The length of the periodic SLs is fixed to 60 uc, which contains a sufficient number of periods to ensure well-formed band-folding in all the cases. With increasing period length in the periodic SLs, the high-frequency phonons are more strongly localized due to the confinement and become almost fully localized at p = 5 uc and p = 6 uc, as shown in Figs. 6(a) and 6(b), which shall be the dominant cause of localization in the high-frequency range in the graded SLs. However, for the moderate-frequency phonons, the PR and weighted PR in the periodic SLs are more or less independent of the period length. Once the different SL units are connected to constitute the graded SLs, the PR and weighted PR are much reduced and smaller than any of those of the periodic SLs. The reduction shall come from the interference of folded modes in different SL units, which introduces extremely strong extinction of transmission in the moderate-frequency range, as shown in Fig. 6(c). This can also be inferred from the results of overall DOS in Fig. 7, where the DOS of graded SLs is quite close to those of periodic SLs in the moderate-frequency range. The large reduction of transmission in this range in the graded SLs shall be caused by the interference of the modes in different SL units. In the high-frequency range, the overall DOS is only significant in periodic SLs with p = 1 uc and p = 2 uc and is negligibly small (except at very few frequency points) in the ones with larger period length (p = 3–6 uc). This is mainly caused by the fact that the high-frequency modes cannot be supported in the Ge segment when it is closer to the bulk state as its thickness increases. As a result, these high-frequency modes in graded SLs can only stay in the first two SL units with short period lengths (p = 1–2 uc). In this perspective, these modes could be termed “gradons,” which represent confined modes in some parts of the graded systems due to the nonexistence of them in other parts.^61,62 The phonon number density profiles of the “gradons” at three typical frequencies within 11–13 THz are shown in Fig. 8. The “gradons” within 13–15 THz are similar and not shown here. The “gradons” are nearly uniformly excited inside the confined region and sharply decay at the boundary of this region, as consistent with the previously reported features.⁶¹ As the “gradons” are essentially confined modes similar to the impurity-induced localized modes,⁶¹ they almost do not contribute to heat transport. Therefore, we are mainly concerned with the localization of moderate-frequency phonons due to coherent backscattering and interference. The partial phonon localization in this work represents a weak-localization, i.e., a weak version of the strong localization in Anderson sense.²² 

The trends of both PR and transmission vs length in the graded SLs with other configurations in Table I are rather similar. The PR, weighted PR, and transmission at the convergence length (Np = 4) for different configurations are shown in Figs. 9(a)–9(c), respectively. With increasing number of SL units and level of disorder, the PR, weighted PR, and transmission all decrease in the entire spectrum, especially in the moderate-frequency range where the phonon localization is the strongest. The normalized thermal conductivity accumulation function vs frequency in Fig. 9(d) provides a corroboration to phonon localization from another point of view. With the increasing level of disorder, the slope in the range of 4–6 THz (the main portion of the moderate-frequency range) gradually declines, which represents decreasing contribution to heat transport from the partially localized states. The localization becomes strong enough to make the thermal conductivity decrease with length at six and seven SL units, where two small plateaus also appear around 1 and 6.5 THz.

Finally, we provide the length-dependent transmission of graded SLs with six SL units for several frequency points in Fig. 10. The frequency points are chosen from those with low PR in Fig. 5(a). A clear exponential decay of the transmission with length is obtained, which has been widely adopted as the evidence of phonon localization.^{32–35,37–39} However, one salient difference is that the transmission decays to a non-zero constant instead of zero in the strong localization regime.^34,35,38 This is a feature of partial localization, which could be caused by either the decoherence of some phonons as highlighted in previous works^37,38 or the insufficient strength of disorder in the present work. Attributed to the moderate disorder, some portion of phonons at specific frequency could still pass through the graded SLs even if the system length increases continuously. The non-zero constant phonon transmission makes the thermal conductivity of graded SLs increase with length at Np ≥ 4–5 in Sec. III A. In a recent contribution,³⁸ a correction term considering the phonon decoherence has been added to the conventional exponential decay formula of transmission, which fits well their direct MD simulation results. We do not attempt to fit the length-dependent phonon transmission in Fig. 10 by the phenomenological formula³⁸ as the partial localization is caused by a different mechanism here. A relevant theoretical model is desired to fit the present results in the full range of length and will be the focus of our future work. Comparing to the electron localization that is easier to achieve,⁴² phonon localization in heat conduction is usually more difficult to observe due to its broadband feature. Thus, further theoretical study is required on the partial phonon localization, especially its dependence on the strength of disorder.

To provide a more intuitive picture of phonon localization, we display the phonon number density distribution computed by Eq. (7) for different frequencies in the graded Si/Ge SLs with six SL units and Np = 4 (L = 84 uc) at 200 K in Fig. 11. The phonon number density distribution has a similar profile to that of the local DOS distribution shown elsewhere;⁵⁶ however, it also includes the effect of phonon occupation number dependent on the system temperature. For consistent discussion, the considered phonon frequency points are the same as those in Fig. 10. For all the considered frequencies except 5.28 THz, the phonon number density is mainly concentrated around a finite spatial region corresponding to a SL unit, as shown by the shaded area in Fig. 11. Two concentrated regions appear at 5.28 THz, as shown in Fig. 11(c). On the other hand, there are also some small peaks outside the concentrated regions in the number density profile. Those are clear features of localized states [7.04 THz in Fig. 11(e)] or partially localized states (other frequencies). In addition, the size of the localization region is roughly inversely proportional to the phonon frequency. The lower-frequency (1.44, 2.88, and 5.28 THz) phonons are mainly localized in the SL units with longer period lengths (p = 4 uc, 3 uc, and 5 uc, respectively), while the higher-frequency (5.76, 7.04, and 8.0 THz) phonons are mainly localized in the SL units with shorter period lengths (p = 1 uc and 2 uc). This could be explained by the fact that the coherence length of phonon wave-packets generally decreases with frequency.^38,48 The coherence lengths of the phonon modes of Si/Ge superlattices within 2–4 and 5–8 THz have been estimated to be around 10 and 5 nm, respectively,³⁸ which are comparable to the length of the SL units with p = 4 uc and p = 2 uc separately. In other words, the SL units with different period lengths in the graded SLs will localize phonons at different frequencies, which introduces significant localization in the broadband spectrum, especially around the moderate-frequency range.

As a comparison, we display the phonon number density distribution for the same set of frequencies in periodic SLs with p = 1 uc and the same total length as that of the graded SLs in Fig. 12. The phonon number density is almost uniformly distributed throughout the SLs with fluctuations, which is a clear feature of extended states. As an intermediate case, the graded SLs with three SL units and Np = 4 (L = 24 uc) have a phonon number density profile combining the features of extended states and localized states, as shown in Fig. 13. With the increasing number of SL units, the localized states will become more dominant in the graded SLs. However, fully localized states are difficult to reach due to the short-range order in such a kind of system, where the partial localization is the main pattern.

Finally, to dig into the localization pattern at the atomic sites, we demonstrate in Fig. 14 the contour of the normalized phonon number density distribution within each concentrated region from Fig. 11. Within each localization region, there are very few bright zones [for instance, at positions of 30 uc and 32 uc for ω = 1.44 THz in Fig. 14(a)] with high phonon number density, which correspond to the sharp peaks of the profiles in Fig. 11. These highly localized vibrational states are almost inside the heavier Ge atomic layers [except two cases for ω = 5.28 THz within the considered frequencies around the positions of 4.5 uc and 57 uc in Fig. 14(c)]. Furthermore, they are located quite close to the Si/Ge interface, which may be caused by the destructive interference of incoming wave-packets and reflected ones from the interface. The present real-space distribution of phonon number density gives a direct picture of phonon localization in the graded Si/Ge SLs.

In summary, a novel scheme is presented to demonstrate the phonon localization in designed graded Si/Ge superlattices. With a sufficient level of long-range disorder, a thermal conductivity minimum with system length appears due to the partial localization of moderate-frequency phonons. We illustrate the clear evidence of localization via the length-dependent participation ratio and transmission, as well as via the thermal conductivity accumulation function vs frequency. In the partial localization regime, the phonon transmission decays exponentially with length to a non-zero constant. We also provide an intuitive picture of localized states by mapping the concentrated distribution of phonon number density. This work will contribute to a new perspective of phonon localization and also a new avenue to engineering heat conduction via the wave nature of phonons.

This work was supported by the Postdoctoral Fellowship of Japan Society for the Promotion of Science (Grant No. P19353), the CREST Japan Science and Technology Agency (Grant Nos. JPMJCR19I1 and JPMJCR19Q3), and the Grant-in-Aid for Scientific Research S (Grant No. 20H05649). This research used the computational resource of the Oakforest-PACS supercomputer system, the University of Tokyo.

The data that support the findings of this study are available from the corresponding author upon reasonable request.