B Z Cunha

M Ichchou

C Droz

A.-M Zine

S Foulard

Interpretable and Physics-Supported Machine Learning Model for Sound Transmission Loss Analysis

Lately, there has been a growing interest in applying Machine Learning and Digital Twins for the speed-up of acoustic simulations. However, the lack of interpretability and physics foundation inhibit the widespread usage of these black-box models by the scientific community. In this article, global sensitivity analysis and feature engineering techniques are leveraged to improve the interpretability and physical consistency of MLbased simulations of the Sound Transmission Loss problem for a variety of plate materials. Computationally efficient sensitivity analysis is obtained via the Mean Decrease in Impurity, which is the byproduct of the training of the Random Forest surrogate models. The resulting sensitivity indices were shown to be similar to the traditional Sobol indices and more accurate than Fourier amplitude sensitivity testing for small datasets. Moreover, introducing basic expert knowledge into the ML inputs helped reduce the surrogate prediction error and interpret the physical meaning of sensitivity indices throughout the frequency spectrum.

Introduction

Recently Machine Learning (ML) and Digital Twins have enhanced vibroacoustic analyses, being able to, e.g., detect, diagnose and predict system failures with high accuracy [START_REF] Stender | Deep learning for brake squeal: Brake noise detection, characterization and prediction[END_REF] and speed up costly simulations [START_REF] Ibrahim | Surrogate-based acoustic noise prediction of electric motors[END_REF]. These methods, however, are classified as "black-box", meaning that the causal relationship between the model's inputs and outputs is not explicitly determined. As understanding the underlying physics of the model is crucial in product design to support decision-making and to advance innovation, the adoption of ML in vibroacoustics has been lagging.

One of the most used techniques to improve the interpretability of a model or phenomenon is Sensitivity Analysis (SA). SA methods evaluate the importance of each input parameter to determine the output variability of a system [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]. In the local sensitivity analysis, commonly used in the vibroacoustic domain [START_REF] Zhang | Sea and contribution analysis for interior noise of a high speed train[END_REF][START_REF] Kim | Shape sensitivity analysis of sequential structural-acoustic problems using fem and bem[END_REF][START_REF] Scarpa | Parametric sensitivity analysis of coupled acoustic-structural systems[END_REF], the inputs are varied individually around a given point in the design space to assess their effects on the model output. On the other hand, Global Sensitivity Analysis (GSA) methods, such as Fourier amplitude sensitivity test (FAST) [START_REF] Cukier | Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. i theory[END_REF] and Sobol's method [START_REF] Sobol | On sensitivity estimation for nonlinear mathematical models[END_REF], provide sensitivities indices valid for a broad input variation range. Moreover, GSA can capture input interaction effects and nonlinearities when higher-order sensitivity terms are considered [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]. The drawback of GSA is that it usually demands a specific sampling method and numerous evaluations [START_REF] Sarrazin | Global sensitivity analysis of environmental models: Convergence and validation[END_REF]. Complex vibroacoustic phenomena have been investigated with GSA, e.g., the importance of micro-perforated panel absorbers parameters on sound absorption levels [START_REF] Laly | Sensitivity analysis of micro-perforated panel absorber models at high sound pressure levels[END_REF], the effects of several vehicle components on its interior noise [START_REF] Lee | Application of global sensitivity analysis to statistical energy analysis: Vehicle model development and transmission path contribution[END_REF], and the sensitivity of Sound Transmission Loss evaluated with analytical methods [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF] and with Statistical Energy Analysis [START_REF] Christen | Global sensitivity analysis and uncertainties in sea models of vibroacoustic systems[END_REF].

SA has also been employed on ML-based methods to render more meaningful and explainable models [START_REF] Humphrey | Improved validation framework and r-package for artificial neural network models[END_REF][START_REF] Razavi | The future of sensitivity analysis: An essential discipline for systems modeling and policy support[END_REF]. Apart from classic SA methods, the so-called permutation feature importance or Mean Decrease in Accuracy (MDA) [START_REF] Breiman | Random forests[END_REF] can also be applied to any ML model to assess its global first-order sensitivities. In MDA, each feature's importance is evaluated as the drop in the model accuracy given by randomly permuting the feature values, i.e., defining it as a noise input. Chai et al. [START_REF] Chai | Comparison of uncertainty quantification process using statistical and data mining algorithms[END_REF] explored the convenience of evaluating the MDA method at the out-of-bag samples of an RF and applied it to the STL problem of a honeycomb sandwich plate. Besides that, some algorithms have sensitivity indices as a post-processing step of the training, as in polynomial chaos expansions (PCE) [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF], gaussian process regressors [START_REF] Marrel | Calculations of sobol indices for the gaussian process metamodel[END_REF] and decision treebased methods [START_REF] Louppe | Understanding variable importances in forests of randomized trees[END_REF]. The Mean Decrease in Impurity (MDI) is an example of a sensitivity method that comes as a by-product of Random Forest (RF) and which sensitivity indices are similar to the Sobol's total sensitivity indices, being more flexible to implement [START_REF] Jaxa-Rozen | Tree-based ensemble methods for sensitivity analysis of environmental models: A performance comparison with sobol and morris techniques[END_REF].

Recent research also integrates physical knowledge into the ML model to increase its accuracy and interpretability for applications in diverse scientific fields. Willard [START_REF] Willard | Integrating scientific knowledge with machine learning for engineering and environmental systems[END_REF] surveyed the methods applied with this purpose including, e.g., penalizing physical inconsistencies during the ML training and defining ML architectures that encode the symmetries and variables dependencies of the physical problem. A straightforward method to embed domain expertise in the ML is to include physics-guided features through feature engineering, as recently applied in thermal [START_REF] Zobeiry | A physics-informed machine learning approach for solving heat transfer equation in advanced manufacturing and engineering applications[END_REF], fluid flows [START_REF] Vaddireddy | Feature engineering and symbolic regression methods for detecting hidden physics from sparse sensor observation data[END_REF], and material science problems [START_REF] Liu | Physics-informed machine learning for composition-process-property design: Shape memory alloy demonstration[END_REF].

This study aims to apply physics-guided features and the MDI global sensitivity analysis to develop an interpretable and explainable RF-based surrogate model for the classical vibroacoustic problem of STL with no additional computational cost. The present paper extends the work in [START_REF] Cunha | On machine learning-driven surrogates for sound transmission loss simulations[END_REF] by performing a comparison of MDI indices with Sobol and FAST indices regarding their formulation, accuracy, and time elapsed, and by encompassing different incident fields and material properties, namely isotropic, orthotropic, and sandwich materials. The performance and the results obtained are compared with the STL sensitivity analyses in [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF][START_REF] Chai | Comparison of uncertainty quantification process using statistical and data mining algorithms[END_REF]. The methodology presented is extendable to other engineering problems and experimental setups, where it may enhance system identification and scientific discoveries.

The paper is organized as follows. Section 2 briefly presents the theory of the classical Sobol's and FAST methods. Section 3 introduces RF and MDI concepts and formulation. Section 4 shows the analytical formulations for the STL for the different plate materials analyzed. The results of the comparison between Sobol, FAST, and MDI are presented in Section 5. The MDI-based sensitivity results in Section 6 bring light to different aspects of the coincidence phenomenon in the STL problem. The influence of including physics-guided features on the surrogate accuracy and the total sensitivity indices are presented in Section 7. Dimension reduction has also been applied following GSA. Section 8 contains the conclusions.

Variance Based Sensitivity Analysis

In the variance-based sensitivity analysis, the importance of each input and input interaction is assessed by their contribution to the output variance. The Sobol's decomposition of variance states that the total variance of the output of a model Y = f (x) = f (x 1 , x 2 , . . . , x D ) can be represented as the sum of the individual effects of the each input x i , denoted as V i , and the higher-order effects from the interactions between a set of inputs as long as some conditions are met [START_REF] Sobol | On sensitivity estimation for nonlinear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF]. That is:

V ar[Y ] = V = D i=1 V i + D i=1 D j>i V ij + • • • + V 1,2,...,D . (1)
Following [START_REF] Sobol | On sensitivity estimation for nonlinear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF], the first-order Sobol sensitivity indices or main effect indices of the input variable x i can be expressed in terms of conditional expected values as

S i = V i V = E x ∼i (Y | x i ) V , (2)
where x ∼i denotes the set of all variables but x i and the numerator of S i is equivalent to the expected average reduction of the output variance when x i is fixed [START_REF] Borgonovo | Sensitivity analysis: a review of recent advances[END_REF]. The total contribution of the variable x i to the variance output, including the effects caused by its interactions with other variables, is measured by the total sensitivity index [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]:

S T i = 1 - V ar x ∼i (E x i (Y |x ∼i )) V . (3)
The total sensitivity indices can be interpreted as the remaining variance when all terms but x i are fixed [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF].

The sum of the total sensitivity indices is equal to or greater than one. When equals to one, the function is additive, i.e., only the first-order terms contribute to the variance, and there is no effect of the interaction between any parameter, what read as S T i = S i .

Evaluating variance-based sensitivity indices involves sampling and uncertainty quantification stages usually based on Monte Carlo or quasi-Monte Carlo integral methods, which are computationally intensive. A popular and less costly technique to evaluate variance-based sensitivity indices is the Fourier amplitude sensitivity test (FAST), developed in 1973 by Cukier et al. [START_REF] Cukier | Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. i theory[END_REF] to evaluate the main effect indices. The FAST method was extended by Saltelli et al. [START_REF] Saltelli | A quantitative model-independent method for global sensitivity analysis of model output[END_REF] to include the computation of the total effects and is now integrated with Sobol's indices. FAST uses a particular period sampling approach to define the multiple Fourier series expansion that approximates the function f (x). As the Fourier series expansion can be decomposed as an Analysis Of Variance (ANOVA)-like decomposition, the components of the model variance can be modeled in terms of the Fourier coefficients [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF]. The FAST and Sobol methods assume uncorrelated and independent variables, otherwise, they lead to biased results if specific methodologies for correlated and dependent variables are not applied [START_REF] Razavi | The future of sensitivity analysis: An essential discipline for systems modeling and policy support[END_REF].

Random Forest Surrogates and Mean Decrease in Impurity SA

Surrogate models or metamodels are approximations of costly functions that are constructed based on inferred information from a set of input-output pairs (q, y) sampled with the true function. An increasingly popular way to construct these surrogates is to use ML regressors to make the data inference. Most ML models are constructed based on a structure/representation that defines the relation between inputs and outputs depending on its parameters. These parameters are optimized during the ML training to minimize the cost function, which defines the error between the predicted ML output and the true output.

Many ML methods use decision trees as its basic structure to define the data relations. The decision tree splits the input feature space using simple rules at each of its branches [START_REF] Breiman | Random forests[END_REF]. Each split at a node n is defined by a threshold t n applied to an input feature q in and divides the space at the current branch C into two new subspaces (two new branches):

C -(q in , t n ) = {(q, y) ∈ C : q in ≤ t n }, C + (q in , t n ) = {(q, y) ∈ C : q in > t n }.
Each node split aims to aggregate samples with similar outputs. Therefore, during the training process the pair (q in , t n ) that minimizes the splitting error at node n, also called impurity criterion G n (q in , t n ), is chosen to split the domain. A common impurity criterion used in regression problems is based on the Mean Squared Error (MSE) of each subspace, defined as

G n (q in , t n ) = L - L   1 L - y∈C - (y -µ -) 2   + L + L   1 L + y∈C + (y -µ + ) 2   , (4)
where µ -and µ + denote the mean output value of the subspaces C -and C + , respectively, and L, L -and L + are the number of samples in subspaces C, C -and C + .

The binary divisions of the decision tree continue successively until an established criterion is met, e.g., the decision tree reaches its maximum depth. At the terminal nodes, a scalar is assigned as output prediction for each final subspace. Individual decision trees can be combined in an ensemble model using bagging to create a Random Forest (RF) predictor [START_REF] Breiman | Random forests[END_REF]. The bagging procedure randomly selects distinct sub-datasets to train each tree, ensuring different trees. The final RF prediction is the average of the decision tree predictions in the forest and, therefore, it has lower variance and generalizes better than individual predictors.

Besides being efficient and able to approximate complex models, RFs are notable for being interpretable ML models. They provide different ways to access feature importances that relate to sensitivity indices. One of them is the Mean Decrease in Impurity (MDI), a by-product of the RF training. The principle of MDI is that the contribution of each input feature to reduce the impurity error is a measure of its importance. Given an RF with Φ κ trees, κ = 1, . . . , T , the MDI importance can be defined as [START_REF] Breiman | Random forests[END_REF][START_REF] Louppe | Understanding variable importances in forests of randomized trees[END_REF]:

MDI(q i ) = 1 T T κ=1 n∈Φκ 1(q in = q i ) L L tot ∆G n (q in , t n ), (5)
where 1(q in = q i ) is one if q in = q i and zero otherwise. The term ∆G n (q in , t n ) denotes the changes in impurity at node n ∈ Φ κ and it is weighted by the fraction of samples at that node L Ltot . Therefore, similarly to how the total Sobol indices can be interpreted as the effect of one input variable in the output variance decrease, the MDI indices are the effect of the input feature in the RF impurity decrease. Indeed, Scornet [START_REF] Scornet | Trees, forests, and impurity-based variable importance[END_REF] proved that, for the case of additive functions with independent variables, the MDI of a decision tree provides a decomposition of the output variance. In the presence of interactions, consistent MDI can be obtained via RF, i.e., by averaging the MDI of several decision trees [START_REF] Scornet | Trees, forests, and impurity-based variable importance[END_REF]. Scornet also showed that the MDI is biased towards positively correlated variables.

MDA is also often applied to evaluate the SA of RF models. Jaxa-Rozen and Kwakkel [START_REF] Jaxa-Rozen | Tree-based ensemble methods for sensitivity analysis of environmental models: A performance comparison with sobol and morris techniques[END_REF] compared MDI and MDA indices with Sobol and Moris metrics for different study cases and highlighted that both MDI and MDA meet the criteria of an "ideal" sensitivity metric, i.e., "suitable for global sampling projects, independent of model structure, relatively easy to implement numerically, stable across all sample sizes and bootstrap resamples and applicable with generic input sampling and heterogeneous input types". Nonetheless, MDI and MDA do not provide the direct effects of inputs on output variance, but their relative importances. The direct effects can be assessed by applying traditional SA approaches to the fitted RF-based surrogate.

Analytical Sound Transmission Loss

The sound insulation capacity is an important characteristic of a structure, which is deeply studied in acoustic projects. The assessment of the sound insulation capacity is usually done by the sound transmission loss (STL) analyses, which evaluate the relation between the sound power reaching the structure I I and the sound power transmitted through the structure I T . Following [START_REF] Fahy | Sound and structural vibration: radiation, transmission and response[END_REF][START_REF] Cremer | Structure-borne sound: structural vibrations and sound radiation at audio frequencies[END_REF], the STL of a plane wave is evaluated by

ST L = -log 10 τ (ϕ, θ, ω) = -log 10 I T I I = -log 10 P T P I 2 , ( 6
)
where τ is the acoustic transparency of a harmonic plane wave of frequency ω, azimuth angle ϕ and incident angle θ, and P I and P T are the incident and transmitted sound fields, respectively.

The sound pressure field can be generically defined as

P * = p * e i(ωt-kxx-kyy±kzz) . (7)
Here, P * can refer to either the incident, transmitted, or reflected sound pressure field, and p * is the associated pressure amplitude. The components of the wave number are defined as k = k x k y k z T = ω c 0 sin θ cos ϕ sin θsinϕ cos θ T , where c 0 is the sound speed through the fluid. Both the incident and the transmitted sound pressure P I and P T have a z component in the positive direction according to the convention in Figure 1, while the sound field reflected by the plate P R has a negative component in z. Hereinafter, the term e i(ωt-kxx-kyy) will be omitted by convenience.

The continuity between the speed of the fluid and the plate at the interface can be stated as

ik z (p R -p I ) = ρ 0 ω 2 w, for z = 0 at the incident side, ik z p T = -ρ 0 ω 2 w, for z = 0 at the transmission side, (8)
where ρ 0 is the fluid density and w the plate displacement. From Eq. 8, it is deduced that p R = p I -p T , and, therefore, the fluid pressure force acting in the plate is given by F = (p

I + p R ) -p T = 2(p I -p T ).
The plate constitutive equation can be formulated as a function of the plate impedance Z as

iZωw = F = 2(p I -p T ). (9)
Combining Eq. 8 and Eq. 9, the ratio between incident and transmitted pressures is given by

p I p T = 1 + Z cos (θ) 2ρ 0 c 0 . ( 10
)
This ratio can be used to calculate the acoustic transparency of a plane wave τ and the STL for τ (Eq. 6). Moreover, the acoustic transparency for a diffuse field τ d (ω) is the weighted average transmission coefficient over all possible angles θ ∈ [0, π/2] and ϕ ∈ [0, 2π]:

τ d (ω) = 2π 0 π/2 0 τ (ϕ, θ, ω) cos θ sin θ dθdϕ 2π 0 π/2 0 cos θ sin θ dθdϕ . (11)
The sound transmission reaches its maximum value whenever the wavelength of the acoustic wave projected onto the plate, also called the trace wavelength, coincides with the wavelength of the natural propagating waves of the plate, which is known as the coincidence phenomenon. Thus, the coincidence frequency ω coin is the frequency with which the projected acoustic wavenumber k t = ω c 0 sin θ is equal to the plate natural propagating wavenumber k p (ω). In the case of diffuse fields, the minimum frequency at which coincidence takes place (i.e., when sin θ = 1) is called critical frequency ω crit .

Isotropic Model

The structural impedance of a Kirchhoff-Love isotropic plate with thickness h, density ρ, Young's modulus E, Poisson ratio v, structural damping η and plate bending stiffness

D = Eh 3 12(1-v 2 ) (1 + iη)
impinged by an acoustic plane wave is given by

Z(ω, θ) = 1 -ω 2 D ρhc 4 0 sin 4 θ iωρh. ( 12
)
Figure 1: MDI indices for the STL analysis of isotropic plates for a plane wave with θ = 60 • .

By replacing the impedance expression in Eq. 10 with Eq. 12, one can evaluate the acoustic transparency and the STL for the isotropic plate. As observed in [START_REF] Fahy | Sound and structural vibration: radiation, transmission and response[END_REF], while the first term of Z(ω, θ) in Eq. 12 is dominated by the plate inertia, the second term is stiffness-controlled. This characteristic of isotropic plates is later observed in this paper in STL curves.

The natural propagating wavenumber of the thin isotropic plate is the bending wavenumber k p = ω 2 ρh D 1/4 [START_REF] Fahy | Sound and structural vibration: radiation, transmission and response[END_REF]. Therefore, its coincidence frequency occurs at

ω coin (θ) = c 2 0 sin 2 θ ρh D . (13)

Orthotropic Plates

Following [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF], the structural impedance of an infinite orthotropic plate in the STL problem is given by

Z(ω, θ, ϕ) = 1 -ω 2 D(ϕ) ρhc 4 0 sin 4 θ iωρh, (14)
where D(ϕ) is the heading bending stiffness of the plate and is defined as

D(ϕ) = D x cos 4 ϕ + 2D xy sin 2 ϕ cos 2 ϕ + D y sin 4 ϕ = E x cos 4 ϕ + 2(E y v xy + 2G xy ) sin 2 ϕ cos 2 ϕ + E y sin 4 ϕ h 3 (1 + iη) 12 1 - Ey Ex v 2 xy . (15)
In Eq. 15 the sub-indices correspond to the direction relative to each material property, and G xy is the inplane shear modulus.

The behavior of the plate bending wave on an orthotropic plate varies according to the wave's direction. Thus, the coincidence frequency occurs at different azimuth angle ϕ, and is given by

ω coin (ϕ, θ) = c 2 0 sin 2 θ ρh D(ϕ) . ( 16
)
Two special cases of Eq. 16 are important in the STL analysis of orthotropic plates, namely the critical frequencies corresponding to the directions x and y, i.e. when ϕ = 0 • and ϕ = 90 • , respectively, and are calculated as

ω crit,x =

Sandwich Panels

Sandwich structures are composite structures composed of two thin and stiff plates, the skin layers, interspersed with a thick layer of lightweight material, the core layer, and are widely used due to their good stiffness-weight ratio. Sandwich material can be designed to shift the critical frequency to higher frequencies and increase the transmission loss at coincidence due to the dampening effects of the core layer. Narayanan and Shanbhag [START_REF] Narayanan | Sound transmission through a damped sandwich panel[END_REF] proposed an analytical implementation of the STL for infinite damped sandwich panels, considering that the deformations of the skin and core can be modeled respectively as purely bending and shear deformations and that the skin material is non-dissipative.

According to [START_REF] Narayanan | Sound transmission through a damped sandwich panel[END_REF], the structural impedance of the infinite sandwich panels approximates as

Z(ω, θ) = k 6 t DsgβY ω(k 2 t +g) 2 +g 2 β 2 + i ωm - Dsk 4 t ω 1+gY (k 2 t +g+gβ 2 ) (k 2 t +g) 2 +g 2 β 2 , ( 18
)
where β is the core loss factor, Y = 3(1 + h c /h s ) 2 is the geometric parameter, m is the mass density of the entire structure, D s is the skin bending stiffness,

g = 2Gc(1-v 2 s )
Eshshc is the shear parameter of the core and G c is the shear modulus of the core. The subindices c and s relate respectively to core and skin. The coincidence frequency equation for the sandwich panels can be found in [START_REF] Narayanan | Sound transmission through a damped sandwich panel[END_REF].

Comparison of Variance-Based and MDI Sensitivity Analyses

In this section, the total importances obtained with the variance-based SA using Sobol and FAST are compared with the MDI results regarding accuracy and computing time. All analyses were carried out on a desktop with a hexa-core 3.1 GHz processor and 32 GB of RAM using Python as the primary programming language and Matlab functions to solve the analytical STL. The Salib library [START_REF] Herman | SALib: An open-source python library for sensitivity analysis[END_REF] was used to obtain the sensitivity indices with Sobol and FAST approaches, while the RF models were instantiated and trained using sklearn. One RF with 100 trees is trained for each of the N ω outputs and 15% of the database is used to assess the RFs accuracy. The same desktop and RF architecture is used in the remaining sections of this paper. The following comparison was made for the STL of isotropic infinite plates impinged by a diffuse field. The five-dimensional design space considered is presented in Table 1 and encloses a broad range of material properties, including the aluminum grades. The MDI indices are comparable with the variance-based total effect indices S T but on different scales, as the MDI indices sum up to one, while i S T i >= 1. Thus, to enable a prompt comparison between MDI indices and S T , the importances will be scaled as

ST i = S T i i S T i
. MDI and FAST accuracies were evaluated having as a reference model R the Sobol analyses with 24576 evaluations, which was determined after a convergence test. The accuracy error of the sensitivity indices from analysis A in relation to the reference R is evaluated by

ε(A, R) = 1 N ω D Nω ω D i | SA T i -SR T i |. (19)
Figure 2 shows the comparison between the performance of FAST and RF. The MDI indices are more accurate with few evaluations, while the FAST method is slightly more accurate for large databases. The elapsed time for the same number of samples is similar for both methods since sampling is the most timeconsuming stage of both approaches. Recently Chai et al. [START_REF] Chai | Comparison of uncertainty quantification process using statistical and data mining algorithms[END_REF] showed good agreement between FAST and MDA importances obtained with RF-based surrogates for an STL analysis, however, the approach did not prove to be advantageous due to the extra cost of MDA and its lack of mathematical foundation.

The scaled total sensitivity importances calculated with the three approaches are shown in Figure 3 for visual comparison. Both MDI and FAST importances were computed using 1280 evaluations, being that the root mean squared error (RMSE) of the surrogate evaluated in the test dataset was 0.40 dB. Figure 2a and Figure 3 indicate that there is a good agreement among the sensitivity indices calculated by all approaches and that the MDI indices obtained with RF-based surrogates are suitable indicators of input importances when few samples are available.

Besides that, the MDI presents other advantages. First, the RF can be trained based on a generic sampling strategy, while FAST requires a particular sampling scheme and can not reuse previous evaluations of the model. Furthermore, once the RF-based surrogate is trained, it can be used to perform a multitude of forwarding analyses in a fraction of the time, being suitable for uncertainty propagation, optimizations, and classical sensitivity analysis. Figure 4 illustrates an example in which the RF-based surrogate was used to replace the true STL model in the calculation of Sobol's first order indices with 24576 evaluations. While the elapsed time using the true model was 1867s, the analysis using the surrogate took 23s, that is, more than 80 times faster. At a previous stage, the sampling and training of this surrogate took more 111s. The higher-order sensitivity indices can be readily obtained similarly.

MDI Results for Isotropic, Orthotropic and Sandwich Plates

MDI analysis can be used to approximate expensive GSA analysis but also to bring interpretability to the surrogate model, a constant demand from the scientific community regarding black-boxes. This section analyzes the STL for different materials and incident acoustic fields through the use of surrogate models and the MDI sensitivity analysis. It should be noted that the MDI results depend on the selected design space, as is the case for every GSA, thus, the results cannot be generalized to other input ranges. In this work, design spaces with a broad variation range were selected to train surrogates suitable for design exploration.

Isotropic Plates

The same design space and database from Section 5 is used in this Section. According to the theory of STL for infinite isotropic plates [START_REF] Fahy | Sound and structural vibration: radiation, transmission and response[END_REF], it is known that two regions are expected: the inertia controlled region, in which the STL can be approximated by the mass law, and the stiffness controlled region. These regions correspond to the region where the first and second terms of the structural impedance of the plate in Eq. 12 predominate, respectively.

The MDI sensitivity analysis performed for a plane wave with θ = 60 • (Figure 5) clearly shows both regions, the inertia controlled, where the variance in ρ and h defines all the output variance and goes up to around 2000 Hz, and the stiffness controlled region where E and h play major roles. The transition frequency range between these regions matches the range of the coincidence frequency for this space design, which ranges from approximately 1680 Hz to 2680 Hz. According to STL physics, the amplitude of the STL valley at the coincidence frequency is controlled by the plate damping η, which importance slightly rises at this frequency range in Figure 5. In the case of an incident diffuse field (Figure 3), the coincidence phenomenon occurs at the frequencies above the critical frequency, which varies approximately between 1260 and 2000 Hz. It is observed that in this frequency range, the importance of the inputs E and h, related to the stiffness, predominates, while the damping effect is predominant for higher frequencies. This behavior was also remarked in [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF], which concluded that for plane waves the location of the STL dip has a greater effect on the output variance than its amplitude, explaining why the damping importance near the coincidence is smaller than expected. For the diffuse field analysis, the dip location is important at the critical frequency range, but the STL amplitude is more important in the coincidence region, where damping predominates. These results are summarized in Table 2 and corroborate that the RF-based surrogate is physically consistent and interpretable.

Table 2: Main MDI importances at different regions of STL of isotropic plates.

ω < ω crit At ω ≈ ω crit ω > ω crit Plane Wave ρ, h ρ, h, E + η peak E, h Diffuse Field ρ, h E,h η

Orthotropic Plates

The investigations of the STL for an orthotropic plate were performed considering a unidirectional fiberreinforced composite plate in which carbon fibers are oriented along the x axis with resin as matrix material, resulting in different Young Modulus in the in-plane directions. The material properties ranges are presented in Table 3, from where 1792 points were sampled for each analysis. The complexity of the STL analysis increases for the orthotropic plates as it has two more dimensions than isotropic plates and the coincidence phenomena occur in multiple directions. The incident plane waves cases were all evaluated with θ = 80 • , close to grazing angle, and with azimuth angle ϕ equals 0 • , 45 • , and 90 • , to evaluate the behavior of the coincidence phenomenon and STL curve for different wave directions. The MDI indices evaluated for STL of each plane wave are depicted in Figures 6a to 6c. It is noted that, regardless of the direction of the incident wave, the low-frequency region is controlled by plate inertia. Nonetheless, the behavior of the MDI indices in the high-frequency region differs for each wave direction.

The results for ϕ = 0 • are similar with the ones for the isotropic plate (Figure 5), where h and E x are the most important variables. Besides that, the small peaks of damping and stiffness importances coincide with the range of ω crit,x , which is from approximately 1140 to 1740Hz, according to Eq. 16. Similar peaks are identified for the results with ϕ = 45 • at higher frequencies, indicating where the coincidence is probably occurring. The in-plane shear modulus G xy has considerable effects when ϕ = 45 • , which is in line with the expected physical behavior for a highly orthotropic material. Finally, for ϕ = 90 • , the bending stiffness in the y direction becomes the most important variable at high frequencies and the coincidence peak also coincides with the ω crit,y range, which goes from around 5800 to 10660 Hz. Figure 6d presents the MDI indices for the diffuse field analyses, where there is a remarkable increase in E x and E y importances in the regions of ω crit,x and ω crit,y ranges, respectively. In both regions, h also had increased importance. In the entire coincidence region, the damping importance is primordial. The findings for the MDI indices of the orthotropic plates are summarized in Table 4. In general, for STL with an incident plane wave, the STL at ω > ω coinc is controlled by the corresponding stiffness. For STL with the diffuse field, the importance of variables related to stiffness increases at the critical frequencies, while the coincidence region is damping-controlled.

Table 4: Main MDI importances at different regions of STL for orthotropic plates. Results were obtained for plane wave (PW) cases with incident angle θ = 80 • and variable azimuth angle ϕ, and for the diffuse field case.

ω < ω crit,x ω ≈ ω crit,x ω crit,x < ω < ω crit,y ω ≈ ω crit,y ω > ω crit,y PW: ϕ = 0 • ρ, h ρ, h,E x + η peak E x , h E x ,h E x ,h PW: ϕ = 45 • ρ, h ρ, h ρ, h,G xy + η peak E x ,G xy , h E x ,G xy , h PW: ϕ = 90 • ρ, h ρ, h ρ, h ρ, h, E y + η peak E y , h Diffuse Field ρ, h E x , h η G xy , η η

Sandwich Panels

The RF-based surrogates for sandwich plates were trained with 4096 evaluations from the design space in Table 5 considering eight material parameters. The MDI indices were evaluated for the STL with an incident plane wave with θ = 80 • , Figure 7a, and STL with a diffuse field, Figure 7b. The MDI results for a sandwich plate show that the core shear modulus also has considerable importance at low frequencies, different from the materials from previous sections in which only inertia controlled the low-frequency region. According to the coincidence frequency equation presented by [START_REF] Narayanan | Sound transmission through a damped sandwich panel[END_REF], the expected range would be from around 4900 to 11400 Hz, which coincides with the peak in β at Figure 7a. For higher frequencies, the STL is controlled by thickness and Young Modulus, variables that compose the skin bending stiffness. Similar observations can be made for the results at low frequency with a diffuse field, while above the critical frequency the sandwich plate STL with a diffuse field is controlled by β, as shown in Figure 7b.

Physics-guided features leveraged by surrogates

Previous knowledge of the physics of the problem can be built into surrogate models through the addition of physics-guided features by simply adding pre-processing steps that rearrange the system variables in terms of physical meaning. This section evaluates the impact of adding these physics-guided features on the accuracy of the surrogate and on the sensitivity indices. The analyses were carried out for isotropic, orthotropic, and sandwich plates with an incident diffuse field. It is worth mentioning that the RFs are retrained with the same databases used in previous sections, therefore, the time elapsed to perform the following analyses is negligible. To keep only the most relevant features, the features with an average importance of less than 2.5% were excluded from the analyses.

Table 6 shows the physics-guided features added for each analysis and underlines the selected relevant features. While the parametric study in [START_REF] Narayanan | Sound transmission through a damped sandwich panel[END_REF] concluded that the core shear parameter g was the most influential parameter for the STL of a sandwich plate, the current analysis indicates that the importance of g is negligible. This can be explained as here the output variance is better explained by other terms correlated to g than by g itself.

Table 6: Physics-guided features included and important variables selected underlined.

STL Variables Physics-guided features

Isotropic Plates h, E, ν, ρ, η m, D Orthotropic Plates h, E x , E y , G xy , ν, ρ, η m, D x , D y , D xy Sandwich Plates h s , E s , ν s , ρ s , h c , G c , ρ c , β m, D s , g, Y
Figure 8 shows the MDI indices for the plate materials considered. Three general aspects are noted:

• In the mass-controlled region the variance is entirely defined by m, except for the sandwich plate results where G c also has a significant contribution.

• The regions of critical frequencies are preceded by an increase in the importance of densities with a subsequent increase in the importance of the stiffness that controls the vibration mode respective to the coincidence taking place.

• In the coincidence regions and STL variability is composed mainly of the effects of damping. The results show that the variable effects can be better represented within the added physical-guided terms. For example, the original analysis of the isotropic case (Figure 3c) showed a contribution of h both at low and high frequencies, while in Figure 8a these contributions are encompassed by m and D at low and high frequencies, respectively. However, it should be stressed that the physics-guided features depend on or [START_REF] Scornet | Trees, forests, and impurity-based variable importance[END_REF]. Therefore, the MDI indices of Figure 8a can be helpful in identifying the contributions of complex physical terms in the system output, but they should not quantitatively determine their importance. In addition to assisting the interpretation of MDI analysis results, the addition of the engineered features guides the RF model to follow the known physics. In this way, the physics-guided features increase the physical consistency of the black box and improve its accuracy, as seen in Table 7.

Conclusion

Interpretable and physics-guided RF-based surrogate models were implemented for the STL analysis through MDI sensitivity analyses and feature engineering. The STL problems were analytically modeled for infinite plates with isotropic, orthotropic, and sandwich materials, considering plane wave and diffuse field as incident acoustic fields.

The MDI-based sensitivity indices showed good accuracy related to the reference Sobol's total sensitivity indices even for a small number of function evaluations. The MDI method proves to be potentially advantageous over FAST since it does not depend on a minimum number of evaluation points and specific sampling methods, furthermore, its performance was superior to or close to FAST performance for the same number of evaluations. However, it should be highlighted that the MDI accuracy depends on the RF's accuracy. The inclusion of physics-guided features was beneficial for surrogate accuracy. Furthermore, the combined use of MDI with feature engineering allows exploring the importance of user-defined terms, clarifying the physical relationships of the analyses.

These capabilities were illustrated for the STL analyses performed, which elucidated different aspects of the STL, enriching the discussions in [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF][START_REF] Chai | Comparison of uncertainty quantification process using statistical and data mining algorithms[END_REF][START_REF] Cunha | On machine learning-driven surrogates for sound transmission loss simulations[END_REF]. Therefore, physics-based interpretations enabled by the methods implemented can increase the level of confidence in the surrogate and allow an in-depth comprehension of the physical phenomena with the potential to boost scientific discovery supported by data.

  Sensitivity indice error ε.

Figure 2 :

 2 Figure 2: Comparison of Sensitivity Analyses approaches.

  Sobol total sensitivity indices.

  FAST total sensitivity indices.

Figure 3 :

 3 Figure 3: Total sensitivity indexes of the STL analyses for isotropic plates.

  True STL model is used to evaluate the outputs.

  RF-based surrogate is used to predict outputs.

Figure 4 :

 4 Figure 4: First-order sensitivity indices with Sobol algorithm.

Figure 5 :

 5 Figure 5: MDI indices of isotropic plates STL for a plane wave with θ = 60 • . Surrogate RMSE: 0.97 dB.

Figure 6 :

 6 Figure 6: MDI indices of the STL analyses of orthotropic plates evaluated with RF-based surrogate.

Figure 7 :

 7 Figure 7: MDI indices of the STL analyses of sandwich plates evaluated with RF-based surrogate.

Figure 8 :

 8 Figure 8: MDI indices of STL analyses with addition of physics-guided features and carrying out features selection.

Table 1 :

 1 Design space considered for the isotropic plate analyses.

	Inputs	h	E	ν	ρ	η
	Units	mm GP a	-	kg /m 3	%
	Lower Bound	7	55	0.20 2300 0.1
	Upper Bound	8	85	0.35 3300	1

Table 3 :

 3 Design space considered for the orthotropic plate analyses.

	Inputs	h	E x	ν	ρ	η	E y	G xy
	Units	mm GP a	-	kg /m 3	% GP a GP a
	Lower Bound 3.25 190 0.25 1300 0.1	5	50
	Upper Bound	4	260 0.35 1800	1	10	65

  Plane wave with θ = 80 • and ϕ = 45 • . Surrogate RMSE: 0.40 dB. Plane wave with θ = 80 • and ϕ = 90 • . Surrogate RMSE: 0.81 dB.

		1		1
	Sensitivity Index	0.4 0.6 0.8	Sensitivity Index	0.4 0.6 0.8
		0.2		0.2
		0		0
		3 2 6 4 5 7 6 4 2 9 0 0 1 2 6 3 1 7 7 1 2 4 8 4 3 4 8 4 4 8 8 6 6 8 5 2 9 6 1 0 1 3 4 7 8 1 8 9 0 3		3 2 6 4 5 7 6 4 2 9 0 0 1 2 6 3 1 7 7 1 2 4 8 4 3 4 8 4 4 8 8 6 6 8 5 2 9 6 1 0 1 3 4 7 8 1 8 9 0 3
		Frequency [Hz]		Frequency [Hz]
	Sensitivity Index	(a) Plane wave with θ = 80 • and ϕ = 0 • . Surrogate RMSE: 0.27 dB. 0 0.2 0.4 0.6 0.8 1 (b) 3 2 6 0 0.2 0.4 0.6 0.8 1 Sensitivity Index 4 5 7 6 4 2 9 0 0 1 2 6 3 1 7 7 1 2 4 8 4 3 4 8 4 4 8 8 6 6 8 5 2 9 6 1 0 1 3 4 7 8 1 8 9 0 3 Frequency [Hz] 4 5 7 6 4 2 9 0 0 1 2 6 3 1 7 7 1 2 4 8 4 3 4 8 4 4 8 8 6 6 8 5 2 9 6 1 0 1 3 4 7 8 1 8 9 0 3 Frequency [Hz] (c) 3 2 6 (d) Diffuse Field. Surrogate RMSE:
				0.19 dB.

Table 5 :

 5 Design space considered for the sandwich plate analyses.

	Inputs	h s	E s	ν s	ρ s	β	h c	G c	ρ c
	Units	mm GP a	-	kg /m 3	% mm GP a kg /m 3
	Lower Bound 1.5	55	0.20 2300 0.1 25	15	200
	Upper Bound	2	85	0.35 3300	1	35	25	300

Table 7 :

 7 Accuracy of the RF-based surrogate with and without physics-guided features (PGF) in terms of the RMSE. input variables and the MDI indices may be biased towards them

		Accuracy without PGF Accuracy with PGF
	Isotropic Plates	0.40 dB	0.22 dB
	Orthotropic Plates	0.19 dB	0.16 dB
	Sandwich Plates	1.02 dB	0.71 dB
	correlate with other

Acknowledgements

This work received financial support of the European Union's Horizon 2020 research and innovation program under Marie-Curie grant agreement No 860243 to the LIVE-I project.


Aurelien Manchon
email: manchon@cinam.univ-mrs.fr

Fengjun Zhuo
email: fengjun.zhuo@kaust.edu.sa

Hang Li

Aurélien Manchon

Topological thermal Hall effect and magnonic edge states in kagome ferromagnets with bond anisotropy

Topological thermal Hall effect and magnonic edge states in kagome ferromagnets with bond anisotropy

Fengjun Zhuo, Hang Li,

I. INTRODUCTION

Magnonics, which aims to investigate the properties and control the transport of magnons in magnetic insulators, has attracted great interests in modern condensed matter physics [1][2][3]. Due to their nanometer wavelength and low energy consumption, magnons are considered as promising information carriers with great potential for data processing devices in the near future [4][5][6].

Concurrently, understanding the topological phases of matter in condensed matter physics has gained substantial momentum and several paradigms involving various (quasi)particle excitations have been investigated, including electrons [7][8][9][10], photons [11][12][13], phonons [14,15], magnetic solitons [16][17][18] and magnons [19,20]. Specifically in magnonic systems, the concept of topology has also been extended to topological magnon insulators (TMIs) [21][22][23][24], magnonic Dirac and Weyl semimetals [25][26][27][28][29][30][31][32][33][34][35][36], as well as magnonic nodal-line and triple-point semimetals [37][38][39][40]. Meanwhile, a growing number of theoretical and experimental studies has been devoted to investigating various topological features of magnon bands on specific lattice geometries, such as, pyrochlore [26,41], honeycomb [25], Kitaev [42][43][44], triangular [45] and kagome lattice [46][47][48][49][50][51], as well as skyrmion crystals [52,53]. However, in most of them, the nontrivial band topology is brought about by the Dzyaloshinskii-Moriya interaction (DMI) [54,55], which plays a role analogous to the spin-orbit coupling in electronic systems lacking inversion symmetry [56][57][58].

Although a magnon is a charge-neutral quasiparticle immune to the Lorentz force usually experienced by charged particles, it experiences an effective magnetic field in momentum space induced by DMI [59,60]. The magnon thermal Hall effect (THE) is characterized by a transverse ther-mal conductivity in response to a longitudinal temperature gradient, which has been realized experimentally in the pyrochlore ferromagnets [61,62] and the kagome ferromagnets [46,47]. The magnon THE was firstly explained as a result of noncompensated magnon edge currents by Matsumoto et al. [59]. Subsequently, studies have found that the net edge currents are closely related to the nontrivial magnonic edge states, and the sign of the thermal Hall conductivity is associated with the propagation direction of edge states [22,63]. The unique nature of magnon THE has inspired lots of discoveries in phase transitions among magnonic topological phases. Heretofore the topological phase transitions for magnonic kagome lattices have been reported under a varying magnetic field [49,50], magnon-phonon coupling [64][65][66] and next-nearest-neighbor (NNN) exchange interactions [67]. Recently, topological magnons in XXZ kagome lattice antiferromagnets have been widely investigated [32,68,69]. When generalized into a kagome lattice, three different nearestneighbor (NN) hopping parameters could form an alternating pattern, giving rise to a rich variety of magnonic topological phases in two dimensions. From this point of view, it is tempting to ask whether a magnonic analog of the topological insulator exists and, more importantly, what non-trivial transverse transport properties of magnons are shown in such a system.

In this paper, we study a kagome ferromagnet with anisotropic NN ferromagnetic interaction and DMI as illustrated in Fig. 1. By varying the modulation factors η 1 , η 2 and η 3 , we find five different magnonic topological phases and plot the topological phase diagram in η 2 -η 3 parameter plane. We moreover show the topological phase transition characterized by the Chern numbers of magnon bands, which can be detected by the change of magnon thermal Hall conductivity. We also demonstrate the correspondence between the thermal Hall conductivity and the propagation direction of the nontrivial edge states. Furthermore, we find sign changes of the thermal Hall conductivity coupled with the topological phase transition. The studies here will be useful for the future design of magnon-based devices, since the thermal Hall conductivity can be easily tuned via modulation factors.

( ) ( ) ( ) . 2 3 , 2 1 , 2 3 , 2 1 , 0 , 1 3 2 1 a a a - - = - = = a a a   
N e x t -n e a r e s t -n e i g h b o u r v e c t o r s i n t r a c e l l :

( ) ( ) ( ) ,. 2 3 , 2 3 , 2 3 , 2 3 , 3 , 0 3 2 1 a a a - - = - = = b b b   
T h e H e i s e n b e r g H a m i l t o n i a n The paper is organized as follows. In Sec. II we introduce our theoretical model and method. Detailed numerical results are presented in Sec. III, including the topological phase diagram, topological magnon band structures, magnonic edge modes, the thermal Hall conductivity. Finally, we end the paper with a brief discussion on the candidate materials and a summary in Sec. IV.

( ) ( ) ( ) å å å × + × - × - = j i j i ij ij j i j i j i j i ij S S z D S S J S S J H , , 2 , ˆ      n w h e r e D D J J k ij k ij h h = = , 1 w i t h k = 1 f o r ( i , j ) = ( A , B ) , k = 2 f o r ( i , j ) = ( B , C ) , k = 3 f o r ( i , j ) = ( A , C ) .

II. SPIN MODEL AND TOPOLOGICAL MAGNONS

Let us consider a collinear ferromagnet on a kagome lattice in the xy plane as schematically shown in Fig. 1. Its Hamiltonian reads

H = - ij J ij S i • S j -J 2 ij S i • S j + ij D ij ij ẑ • (S i × S j ) , (1)
where S i is the vector of spin operators at site i. The first two terms describe the NN and NNN ferromagnetic interactions (J ij , J 2 > 0), respectively. The third term is the out-of-plane NN DMI, where D ij is the DMI strength and ij = ±1 depends on the chirality of the triangles in the kagome lattice. Specifically, the NN interactions J ij and D ij are modulated around its equilibrium value J 1 and D with different coefficients η i depending on the bond direction as shown in Fig. 1:

J ij (D ij ) = η 1 J 1 (D) along α 1 , J ij (D ij ) = η 2 J 1 (D) along α 2 and J ij (D ij ) = η 3 J 1 (D) along α 3 directions, respec- tively.
The anisotropic NN interactions can be realized by lattice distortion induced by a controllable voltage or mechani- cal strain. Notice that such lattice deformations should change not only the NN exchange interactions but also the NNN ferromagnetic interactions. As discussed in the Appendix, the anisotropy of NNN ferromagnetic interactions does not produce new nontrivial band topology, introducing only restricted modifications to the band dispersion and distribution of Berry curvature. So for simplicity, we will treat the NNN ferromagnetic interactions as isotropic throughout this paper. In addition, the simultaneous changes of NN ferromagnetic interaction and DMI will preserve the collinear ferromagnetic ground state, because the DMI is well below the threshold value for a new spin configuration formation.

In the following, we apply a linear spin-wave theory to obtain the tight-binding magnon Hamiltonian in the momentum space. We express the spin operators in Eq. ( 1) in terms of magnon creation operator b † i and annihilation operator bi by using the Holstein-Primakoff transformations [70]:

S + i = √ 2Sb i , S + i = √ 2Sb † i and S z i = S -b † i b i
, where we introduce the magnon ladder operators S ± i = S x i ± iS y i . We neglect the higher-order terms in the magnon operators from the magnon-magnon interactions, as we consider the low temperature regime here. After a Fourier transformation, Eq. ( 1) becomes

H = S k Ψ † k H (k) Ψ k , where the tight-binding magnon Hamiltonian, H = H 0 + H J1 + H D + H J2 , is a 3 × 3 matrix under the basis Ψ † k = b † A , b † B , b † C . The submatrices of H read H 0 =   m 1 0 0 0 m 2 0 0 0 m 3   , ( 2a
) (e) 3 1  C (c) 1 1   C (d) 1 1   C (a) 1 1  C (b) 1 1  C (f)
H J1 = -2J 1   0 γ 1 γ 3 γ 1 0 γ 2 γ 3 γ 2 0   , (2b)
H D = -2iD   0 γ 1 -γ 3 -γ 1 0 γ 2 γ 3 -γ 2 0   , (2c)
H J2 = -2J 2   0 cos (k • β 1 ) cos (k • β 2 ) cos (k • β 1 ) 0 cos (k • β 3 ) cos (k • β 2 ) cos (k • β 3 ) 0   , (2d) where m 1 = 2J 1 (η 1 + η 3 ) + 4J 2 , m 2 = 2J 1 (η 1 + η 2 ) + 4J 2 , m 3 = 2J 1 (η 2 + η 3 ) + 4J 2 and γ i = η i cos (k • α i ),
respectively. The linking vectors α i and β i are illustrated in Fig. 1.

In contrast to electronic systems, the absence of Fermi surface in bosonic systems leads to the ill-defined Chern number of magnons. However, the Chern number of magnons for the nth magnonic bulk band can still be given in a similar way,

C n = 1 2π BZ dk 2 Ω z nk . (3)
and the Berry curvature of magnons Ω z nk is defined as

Ω z nk = -2 m =n Im ψ nk |∂ kx H| ψ mk ψ mk ∂ ky H ψ nk (ε nk -ε mk ) 2 , (4)
where ψ nk and ε nk are the eigenvectors and eigenvalues of H (k) for the nth band, respectively. Another important quantity in topological systems is the winding number ν m in band gap m, which is given by the sum of Chern numbers up to the mth band [71],

ν m = n≤m C n . (5)
The winding number reflects the surface properties of the system, that is to say, we can find |ν m | pairs of topologically nontrivial edge states in the mth band gap and their propagation direction is determined by the sign of ν m . The topological nature of magnonic edge states can be probed by measuring of the magnon thermal Hall effect in magnonic systems. The magnon thermal Hall conductivity is calculated using the following formula [59,60] where

κ xy = - k 2 B T 3 n=1 BZ d 2 k (2π) 2 c 2 ρ B (ε nk ) Ω z nk , (6)
ρ B (ε nk ) = e ε nk /k B T -1 -1 is the Bose-Einstein distribution.
The weighting function is given by

c 2 = (1 + x) ln 2 1+x
x -ln 2 x -2Li 2 (-x), with Li 2 (x) being the polylogarithm function.

III. RESULTS

For the sake of simplicity, the Hamiltonian Eq. ( 1) is normalized to J 1 S, so that energy, temperature and thermal conductivity are expressed in the units of J 1 S, J 1 S/k B and J 1 S (k B / ), respectively. In the following numerical calculations, unless otherwise specified, D = 1 and J 2 = 1. Without loss of generality, we set η 1 = 1 as the scale unit, whereas η 2 and η 3 are tunable parameters. The Chern numbers of magnons are calculated by using the algorithm of Fukui et al. with 1000×1000 points in momentum space [72].

A. Topological phase diagram and band topology

As mentioned above, varying the anisotropic NN exchange interactions provides a handle to tune topological phase transitions in kagome ferromagnets. The topological phase diagram as a function of η 2 and η 3 is depicted in Fig. 2. These phases are characterized by sets of Chern numbers (C 1 , C 2 , C 3 ) of the lower, middle and upper magnon bulk bands, as shown in Fig. 4. Specifically, five different topological phases are found with Chern numbers: phase I (3, -2, -1), phase II (1, -1, 0), phase III (1, 0, -1), phase IV (-1, 2, -1) and phase V (-1, 1, 0). In the following disscusion, we chose {η 2 , η 3 } = {1, 1}, {2, 0.3}, {2, 0.8}, {2, 1.5} and {2, 3} in each topological phase, all of which are marked by the black pentagrams in Fig. 2.

In Figs. . We also note that the Berry curvature exhibits higher densities when the band gap is shifted towards closing, for example, the peaks of Berry curvature along Γ -K line and near K point in Figs.

3(c) and 3(e).

To better visualize the magnonic edge states, we solve the eigenvalue problem of a nanoribbon geometry with open boundary conditions, and then plot the band structures in Fig. 5. The winding numbers in phase II are ν 1 = 1 and ν 2 = 0, so two crossed chiral edge states are clearly found in the lower band gap and two gapped trivial edge states are found in the upper band gap, as shown in Fig. 5(a). In phase III, as the two winding numbers are the same, ν 1 = ν 2 = 1, two equivalent pairs of nontrivial edge states emerge in each band gap in Fig. 5(b). As η 3 is further increasing, the nonzero winding numbers in phase IV are ν 1 = -1 and ν 2 = 1, and we can find two pairs of nontrivial edge states in the two band gaps in Fig. 5(c). Meanwhile, the edge states have opposite slopes v g = ∂ε/∂k (i.e., group velocity): the red one moves from the bottom left to the top right (v g > 0) in the lower gap; instead, the red one moves from the top left to the bottom right (v g < 0) in the upper gap. In other words, they propa- gate in opposite directions. Finally, the system has a pair of nontrivial edge states in the lower gap and a pair of trivial edge states in the upper gap with ν 1 = -1 and ν 2 = 0, as shown in Fig. 5(d). Note that the propagation directions of the nontrivial edge states in the lower gap are opposite in Figs. 5(a

B. Topological thermal Hall effect

In this section, we discuss the topological nature of magnon THE. In magnonic systems, the topological phases and phase transitions can be detected by the measurement of magnon dispersion using inelastic neutron scattering [47]. In Fig. 6, the magnon thermal Hall conductivity κ xy in low temperature regime is plotted as a function of temperature with different modulation factors {η 2 , η 3 } in each topological phases, which all show a monotonically increasing behavior. We also find xy < 0 in phase I, II and III, conversely, κ xy > 0 in phase IV and phase V.

The thermal Hall conductivity κ xy as a function of η 3 when η 2 = 2 for four different temperatures are depicted in Fig. 7. The system undergoes the phase transition from phase II to phase V, as η 3 increases from 0 to 3. In addition, we find a sign reversal of κ xy when the topological phase transits between phases III and IV, which can be explained by connecting the propagation directions of the nontrivial edge states. From Eq. (4) and Eq. ( 6), we see that C n and κ xy only have different weighting functions in the integrand. Hence, the dominant contribution to the thermal Hall conductivity at low temperatures comes from the lowest band, since the lowest band is more occupied due to the c 2 function within the bosonic statis- tics. As a consequence, the sign of κ xy is governed by the sign of ν 1 in low temperature regime as ν 1 = C 1 . As we mentioned above, the propagation directions of edge states in the lower gap are opposite between the topological phases III and IV. Finally, the sign reversal of κ xy is accompanied by the sign reversal of ν 1 , even though the two pairs of nontrivial edge states also have opposite propagation directions in phase IV.

IV. DISCUSSION AND SUMMARY

In this paper, a tunable topological magnon excitation in kagome lattice ferromagnets has been predicted. We shall further briefly comment on the candidate materials for experimental realizations of the above theoretical model. Recent experiments show that the materials Cu(1-3,bdc) [46,47], as well as the mineral haydeeite-α-MgCu 3 (OD) 6 Cl 2 [48], exhibit a stable ferromagnetic phase down to the monolayer limit. In these materials, magnetic ions are arranged on a kagome lattice with their magnetization pointing perpendicular to the plane. Meanwhile, the variations of the NN interactions can be easily achieved from the lattice distortions induced by applying external perturbations, such as mechanical strain or pressure, because the interactions depend sensitively on the inter-atomic distances [73][74][75][76]. The different types of strains studied in the present work can be achieved and tuned by applying a uniaxial mechanical strain field, either along different crystallographic directions. For instance, a uniaxial strain along α 1 direction shall result in η 2 = η 3 = η 1 . Using the parameters in Ref. [48], we can find a thermal Hall conductivity of -3.2×10 -3 W/Km with isotropic NN interaction (η 2 = η 3 = η 1 ) and 3.2 × 10 -3 W/Km with anisotropic NN interaction (η 2 = η 3 = 2η 1 ) at the temperature T = 50 K. Tuning the anisotropy parameters over a wide range to cover the diagram reported on Fig. 2 remains a technical challenge that can be achieved by depositing the kagome ferromagnet on a piezoelectric substrate for instance [77,78]. Detailed first-principles calculations are needed to identify suitable interfaces [79].

Thus far, we have studied the magnon Hall effect in a picture of magnons within the bosonic description. As mentioned in Sec. II, we truncate the bosonic Hamiltonian to quadratic order by neglecting the magnon-magnon interaction in the low temperature regime. However, by using the Holstein-Primakoff transformations, it is unclear whether the quartic terms for the interaction between magnons could be ignored. In this perspective, another theoretical formalism of the thermal Hall effect was developed in the language of spin operators instead of boson operators [80]. In this formalism, the thermal Hall effect in a correlated paramagnetic phase can also be captured with reductions to the Schwingerboson representation of spin operators by using self-consistent mean-field theory. Remarkably, the thermal Hall coefficients obtained from the linear-response theory for noninteracting models [59] and the mean-field Hamiltonian for the interacting spin model are identical. On the other hand, as lowenergy charge-neutral excitations, phonons can also carry angular momenta [81][82][83] and induce a transverse phonon angular momentum flows at finite temperatures, the so-called phonon thermal Hall effect [14]. However, the phonon angular momenta in a phonon Hall current can only create negligible accumulation at the crystal edges in comparison with the magnon spin moment accumulation in a magnon Hall current [84,85]. In addition, since the magnon-phonon interaction may significantly modify the magnon-mediated THE due to symmetry breaking by an in-plane DMI and the magnetoelastic coupling [64,86,87], its role on the nontrivial topology of magnon excitation would be an interesting topic left to future studies [87]. Since magnon-magnon interactions and magnonphonon interactions need to be taken into consideration at high temperatures, the linear spin-wave theory becomes inaccurate in this limit. However, self-consistent Holstein-Primakoff and Schwinger-boson methods within Hartree-Fock meanfield approximation will be applicable to computation of the thermal Hall response functions [80,88], which goes beyond the scope of this paper.

In summary, we have shown that tuning the relative strength of anisotropic NN ferromagnetic interaction and DMI in kagome ferromagnets gives rise to the topological phase transition between different topological phases characterized by the Chern numbers of magnon bands. We describe the band topology of magnons with the magnonic edge states in the band gaps. Furthermore, we have calculated the thermal Hall conductivity κ xy and expound its relevance to the topological invariant ν m . We find that κ xy can easily be controlled by tuning the modulation factors. In addition, a sign reversal of κ xy was found along with topological phase transitions. In this regard, the present studies open thrilling perspectives for experimentalists and can be applied in the field of thermal micro-sensors in the future.

  r e s t -n e i g h b o u r v e c t o r s :

Figure 1 .

 1 Figure 1. Schematics for the kagome ferromagnet with lattice constant a. Spin sublattices A, B, and C are placed at the corners of the triangles. η1, η2 and η3 are the modulation factors along A-B bond (green line), B-C bond (blue line) and A-C bond (purple line) directions, respectively. The NN and second-NN connecting vectors are labeled by αi and βi.

Figure 2 .

 2 Figure 2. Topological phase diagram of the kagome lattice on the η2 -η3 plane with regions characterized by sets of Chern numbers. Five variable sets {η2, η3} we use in later calculations are shown as black pentagrams in each topological phase.

Figure 3 .

 3 Figure 3. Berry curvatures of the lowest band with different modulation factors. (a) η2 = 2 and η3 = 0.3 in phase II. (b) η2 = 2 and η3 = 0.8 in phase III. (c) η2 = 2 and η3 = 1.5 in phase IV. (d) η2 = 2 and η3 = 3 in phase V. (e) η2 = 1 and η3 = 1 in phase I. Corresponding Chern numbers of the lowest band C1 are given in each panel. The dotted red lines denote the edges of the first Brillouin zone. (f) The first Brillouin zone of the reciprocal lattice. The dotted black lines denote the selected k path connecting the high-symmetry points in the magnon bands shown in Fig. 4.

Figure 4 .

 4 Figure 4. Magnon band structures of a kagome ferromagnet for the selected modulation factors in different topological phases. (a) η2 = 1 and η3 = 1 in phase I. (b) η2 = 2 and η3 = 0.3 in phase II. (c) η2 = 2 and η3 = 0.8 in phase III. (d) η2 = 2 and η3 = 1.5 in phase IV. (e) η2 = 2 and η3 = 3 in phase V.

  3(a)-3(e), we plot the distribution of Berry curvature Ω z nk associated with the lower magnon band for five selected sets of modulation parameters {η 2 , η 3 }, that drive the system into five different topological phases. It is remarkable that the Berry curvature exhibits a C 2 rotational symmetry in Figs. 3(a)-3(d) due to the anisotropic NN interactions, instead of a C 6 rotational symmetry when η 1 = η 2 = η 3 shown in Fig. 3(e). The summation of Ω z nk over the first Brillouin zone gives the corresponding Chern number of the lower band C 1 . Then, we turn our attention to the exploration of the band topology in kagome ferromagnets. Fig. 4 shows the magnon bulk bands along the high-symmetry directions (Γ -K -M -Γ) of the Brillouin zone [Fig. 3(f)] in each topological phase. Three well separated magnon bands, i.e., two well defined band gaps, can be found in all phases. The upper band is close to the middle band at M point with increasing η 3 when η 2 = 2 as shown in Figs. 4(b)-4(e)

Figure 5 .

 5 Figure 5. Magnon band structure with coupled magnonic edge states in a kagome ferromagnet with a nanoribbon geometry at different modulation factors. (a) η2 = 2 and η3 = 0.3 in phase II. (b) η2 = 2 and η3 = 0.8 in phase III. (c) η2 = 2 and η3 = 1.5 in phase IV. (d) η2 = 2 and η3 = 3 in phase V. The black lines are the magnon bulk bands. The blue/red lines represent the edge states for the upper/lower edges.

  )-5(b) and Figs. 5(c)-5(d).

Figure 6 .

 6 Figure 6. The magnon thermal Hall conductivity in low temperature regime as a function of temperature for different modulation factors {η2, η3} in each topological phase.

Figure 7 .

 7 Figure 7. Topological thermal Hall conductivity in low temperature regime versus modulation factor η3 at different temperatures with η2 = 2. Four different topological phases labeled with Chern numbers are separated by the critical points shown as vertical dashed magenta lines.

Figure 8 .

 8 Figure 8. Topological phase diagram on the η2 -η3 plane with regions characterized by sets of Chern numbers when taking anisotropic next-nearest-neighbor exchange interactions into consideration. Five different topological phases are labeled with the same Chern numbers as shown in Fig. 2.

ACKNOWLEDGEMENTS F.Z. was supported by King Abdullah University of Science and Technology (KAUST). A.M. acknowledges support from the Excellence Initiative of Aix-Marseille Université-A*Midex, a French "Investissements d'Avenir" program. H.L. acknowledge the support from Henan University (No. CJ3050A0240050) and National Natural Science Foundation of China (No. 11804078).

APPENDIX: ANISOTROPIC NEXT-NEAREST-NEIGHBOR EXCHANGE INTERACTIONS

As described in the main text, we propose that the topological phase transition can be tuned in kagome ferromagnets with anisotropic NN interactions. Generally speaking, such NN bond anisotropy may also influence the NNN ferromagnetic interactions. In this case, the new Hamiltonian reads

where the anisotropic NNN ferromagnetic interactions J ij are given as

Here, all other parameters have the same meanings and values as mentioned in the main text. Fig. 8 shows the topological phase diagram as a function of η 2 and η 3 when taking anisotropic NNN ferromagnetic interactions into consideration. In comparison with the topological phase diagram in Fig. 2, both phase diagrams are divided into five identical phases listed in the main text. In other words, there's no other nontrivial topological phase creating. In addition to that, the anisotropic NNN ferromagnetic interactions redistribute the region of each phase, as the specific magnon band dispersions and distributions of Berry curvature are modified by the difference between the NNN ferromagnetic interactions.