Inverse design of a topological phononic beam with interface modes

Inspired by the idea of topological mechanics and geometric phase, the topological phononic beam governed by topological invariants has seen growing research interest due to generation of a topologically protected interface state that can be characterized by geometric Zak phase. The interface mode has maximum amount of wave energy concentration at the interface of topologically variant beams with minimal losses and decaying wave energy fields away from it. The present study has developed a deep learning based autoencoder (AE) to inversely design topological phononic beam with invariants. By applying the transfer matrix method, a rigorous analytical model is developed to solve the wave dispersion relation for longitudinal and bending elastic waves. By determining the phase of the reflected wave, the geometric Zak phase is determined. The developed analytical models are used for input data generation to train the AE. Upon successful training, the network prediction is validated by finite element numerical simulations and experimental test on the manufactured prototype. The developed AE successfully predicts the interface modes for the combination of topologically variant phononic beams. The study findings may provide a new perspective for the inverse design of metamaterial beam and plate structures in solid and computational mechanics. The work is a step towards deep learning networks suitable for the inverse design of phononic crystals and metamaterials enabling design optimization and performance enhancements.


Introduction
Topological physics and mechanics inspired phononic crystals (PnCs) and acoustic metamaterials (AMs) have observed a surge of research studies in unusual control of acoustic and elastic waves. Interface states that are related to topological phase inversion [1] of trivial and non-trivial phases * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. provide robust, lossless, immune to backscattering and sharp edges wave transportation at the interface of topologically distinct PnC and AM [2][3][4][5]. Such interface mode provides potential avenues for wave control including wave energy focusing/harvesting [6], vibration and noise control [7,8], topological insulators [9,10], higher order corner states/modes [11,12], and design of other acoustic devices. For a one-dimensional (1D) periodic system, such an interface mode is initially reported in condensed matter physics using Su-Schrieffer-Heeger model [13,14]. The interface state is observed inside the bandgap (BG) of topologically protected PnC. Later, multiple other studies have reported the interface state generation mechanism, role of geometric phases in the bounding BG edges and procedure to achieve interface modes in 1D [3,5,15], two-dimensional (2D) [16][17][18] and threedimensional [19] periodic structures for both air-born acoustic and solid-born elastic waves.
While discussing topological mechanics and its application in PnCs and AM, one cannot negate the importance of topological/geometric phases [20]. The idea of geometric phase, also refereed as Pancharatnam-Berry phase is first proposed by Pancharatnam in 1956 for investigating light propagation through a series of polarizer. Then this idea is employed by Berry in 1984 for quantum mechanics [21]. Further details can be found in a recently reported review articles [20,[22][23][24]. Chern number is used to characterize the topological invariants for 2D system that is identical to Berry phase enclosing the first Brillouin zone [24][25][26]. For a 1D system, Zak phase (a special type of Berry phase defined along a 1D bulk band) is employed to determine the topological invariants. In such a system, the presence of nontrivial Zak phase indicates the birth of interface or edge states. Therefore, in modern physics topological invariants has a vital role in the discovery of materials. As described in the prior works [2,26], the Zak phase explains the global characteristics of wave functions throughout a whole band of a 1D system and is unaffected by local disturbances and disorders, preserving the system's symmetry. Ma et al [24] applied this concept from condensed matter physics to acoustic and mechanical systems. Hence in the context of topology, Zak phase has a prime importance in identifying the interface states. The interface state is generated when the mode transition frequency is common between the trivial and nontrivial phases bounding the BG edges with symmetric and antisymmetric eigenmode polarization [2,15].
In general, the complex structural geometry and wave responses play a substantial role in the elastic topological states. The works about topological PnC to date reported are mostly based on a forward-design approach where first we precisely define the geometric parameters and subsequently obtain the wave response such as wave dispersion/frequency response. In this process, a set of geometric entities are supplied and wave response is calculated based on trial and error method. Depending upon the complexity of the model, this approach is computationally expensive and time demanding. In contrast, an inverse design mechanism can achieve design parameters from target topological response/wave dispersion property. An example is recently reported by Li et al [27] where machine learning (ML) and deep learning (DL) datadriven methods are employed to inversely design PnC unit cell structures from a large set of generated database based on random geometries. Recently, these data-driven approaches have found a surge of research interest in the field of PnC and AM due to the fact that metamaterial structural features determine the wave performance rather than material composition and/or chemistry [28]. If one maps the data pattern between metamaterial geometric features with wave properties like acoustic/mechanical frequency responses, then by forward and inverse design processes, optimized PnC and AM can be obtained. Prior to applying it for the complex PnC structures that would be promising for real-life applications, this work begins with a simple structure namely a beam model that has simple and accurate analytical solutions.
Therefore, in this study we develop ML and DL data-driven models to inversely design a topological phononic beam to accurately predict topological interface states via both forward and inverse design approaches. The analytical modelling and PnC design strategy is based on results reported in Muhammad et al [2]. Here we extend this study by introducing ML and DL based forward and inverse design approaches to effectively design phononic beam with interface states. The transfer matrix method (TMM) is employed to derive the constitutive wave dispersion relation for longitudinal elastic wave. We employed Euler-Bernoulli beam theorem with modified TMM to obtain the dispersion relation for bending elastic waves. These constitutive analytical equations are used to generate a large database for network training. The choice of network is very much dependent on the input data and desirable output. In this study, an auto-encoder is developed to optimize the phononic elastic beam. Here we transform the band structure data from frequency spectrum to binary digits 0, 1 and −1 where 0 shows passbands and 1 and −1 correspond to BG with different wave energy polarization determined by analytical calculation of Zak phase [26,29]. First by forward design process, for a set of geometric parameters, the band structure is obtained. Here we designate the geometric entities with shape function δ = (2L B − L A ) /L where L is length of beam and subscript shows beam A, B. Then by inverse design process, δ is predicted for the required topological beam band structure governing interface states. These network generated results are validated by FEA based COMSOL Multiphysics and experiment testing on predicted δ values. Overall, an excellent agreement is observed.
The paper is organized as follows. Section 2 briefly explain the modelling strategy and procedure for deriving the constitutive dispersion relation. The ML and DL network architecture are explained in section 3. The results are discussed in section 4. Finally, the conclusion and future prospect of the developed network for metamaterial research are discussed in section 5.

Modelling strategy and governing equations
In this section the TMM and modified TMM for deriving the constitutive relation for longitudinal and bending elastic waves, respectively are briefly discussed. For more details, one can refer to Muhammad et al [2]. The schematic diagram of 1D phononic beam is shown in figure 1. The lattice constant of the periodic beam is a = L = 100 mm. The unit cell of the periodic lattice consists of two topologically distinct beams with varying length and diameter. The beam B with length L B and diameter D B = 8 mm is sandwiching beam A of length L A and diameter d A = 4 mm. The topology of beam is varied by introducing a shape parameter δ = (2L B + L A )/L. We consider aluminium rod with young modulus E = 70 GPa, mass density ρ = 2700 kg m −3 and Poisson ratio v = 0.33 for analytical and numerical modelling, see table 1. Later for   [2,15], bending eigenmodes are present at a lower frequency region compared to the longitudinal counterpart. In numerical codes, such eigenmodes are present in sequential order and need to be separated by postprocessing steps [2,15]. For longitudinal waves, the dispersion relation obeys ω 2 = k 2 E ρ where ω, k are angular frequency and wavenumber respectively. Since the slenderness ratio of the beam is L/D B = 12.5 it is within permissible limits [2], we adopted Euler-Bernoulli beam theory to achieve the analytical dispersion relation for bending elastic wave.

Constitutive equation for longitudinal wave
The constitutive dispersion relation for longitudinal wave is derived as follows. More details are outside scope of this work and can be found in Muhammad et al [2]. In an elastic homogenous rod/bar, the longitudinal wave propagation can be written as where u, ρ and E are displacement, mass density and young's modulus of the beam. The general harmonic solution of equation (1) is Likewise, the constitutive equation for the axial force is where F (x, t) and S are the axial force and cross-section area of the beam. For beam B, the displacement and axial force correlation can be written as According to TMM, the displacement and axial force at the left and right sides of beam B becomes Here A 1 and A 2 are unknown displacement fields representing the eigenmodes. From equations (5) and (6), transfer matrix can be expressed as Likewise, applying the above procedure for beam A, the transfer matrix becomes See appendix for expression of T A and T B . The boundary condition considered between two sub-beams is Combining equations (7)-(9) leads to In order to make the proposed unit cell structure infinitely periodic, Floquet-Bloch condition is applied on equation (10) Further simplification of equation (11) gives where k is Bloch wavenumber in 1D system and I is a 2 × 2 identity matrix. The solution of equation (12) gives the constitutive dispersion relation as [2] cos Equation (13) is so-called dispersion relation for longitudinal wave where for each ω, wavenumber k can be easily calculated. We used equation (13) to developed database for longitudinal wave.

Analytical solution for bending elastic wave
According to Euler-Bernoulli beam theory, the governing equation of motion for bending elastic wave can be expressed as [2] where w (x, t) is lateral displacement, EI Z , ρ and S are bending rigidity, mass density and the cross-sectional area of beam, respectively. For harmonic solution of the form w (x, t) = w (x) e iωt , the expression for out-of-plane displacement filed can be written as where k = ρSω 2 /EI is wavenumber for the flexural waves and A, B, C, and D are unknown parameters. The expression for angular rotation, shear force and bending moment are are the first order, second order and third order derivatives of w (x) with respect to x. The expression for P, Q, R, S are the combination of hyperbolic and triangular functions The equation (16) also satisfy following conditions .
(17) The first, second and third derivatives of w (x) with respect to x are . (18) That leads to Now the unknown parameters A, B, C, D are in the form of displacement w (x), angular rotation φ(x), bending moment M (x) and shear force V (x). The expressions at x = 0 are denoted by w 0 , φ 0 , M 0 and V 0 that can be expressed as Hereafter all the parameters are presented in dimensionless formats. This is required for development of input database. The normalized parameters introduced arẽ (22) Further modification of equation (21) gives In the equation above A is state vector, T j is transfer matrix for the jth section of primitive cell in an infinite beam structure, x L shows the local coordinate relative to the coordinate x in jth section. Here, x L = 0 and x L = L B indicate the left and right ends of the beam B. Similarly, x L = L B and x L = L A + L B indicate the left and right ends of the beam A, respectively. The details for the transfer matrix T is presented in the appendix.
In the nth unit cell of the periodic structure, we have at nL ⩽ x ⩽ nL + L B and 0 ⩽ x L ⩽ L B , specifically we obtain At x = nL + L j and x L = L j , the analogy above gives The continuity of the displacement, angular deflection, bending moment and shear force at the interface of different section yields Equation (28) is the relationship of the initial parameters between the nth and the (n + 1) th unit cell of an infinitely periodic beam. The expression A [(n + 1) L] describes the continuity of displacement, angular rotation, bending moment and shear force between topologically distinct unit cells and their physical quantities at the left end of (n + 1) th and the right end of nth primitive cell, T designates the transfer matrix.
Because of the periodicity condition, the state vector A should also follow the Bloch periodicity condition [30] A [(n + 1) L] = e ikL A (nL) (29) where k is 1D Bloch wavenumber. Substituting equation (28) into equation (29) yields where I is a 4 × 4 identity matrix. Equation (30) indicates that the eigenvalues for the dispersion relation can be obtained by determining the roots of the following matrix determinant The solution of equation (31) shows, for each ω there exists four complex roots. Assuming each root is in the form of x + iy and an associated wavenumber k be in the form of p + iq, these parameters yield the following relationship [31] e i(p+iq)a = x + iy.
Further simplification leads to q = − ln x 2 +y 2 2L and over the Brillouin zone π a , − π a where a = L is the lattice constant. By solving the equation (31) based on the conditions defined in equations (32)-(33), the wavenumber is determined for specified frequency range.

Zak phase
The detailed description and significance of Zak phase is beyond the scope of this work. More details can be found in [26,29]. Xiao et al [26] extended the concept of Zak phase to acoustics. For nth Bloch band, the Zak phase can be expressed as where ξ n,k (x) is the Bloch eigenfunction of a periodic structure for a specified wavenumber k. 1/ 2ρc 2 is a factor for the weight function of elastic medium. In the nth band of the longitudinal and bending waves, the Bloch eigenfunction is [26] U n,k (x) or W n,k (x) = ξ n,k (x) e ikx where U n,k (x) and W n,k (x) denote the longitudinal and bending wave fields, respectively. In this study, the reflection phase of the associated two BGs is determined to calculate the Zak phase sgn (ψ n−1 ) / sgn (ψ n ) = −e iθ Zak n withθ Zak n = 0 or π.
Where ψ ∈ (−π, π) is reflection phase of the nth BG and it satisfies Sgn [ψ] = −1 correspond to ψ n ∈ (−π, 0), Sgn [ψ] = 1 is associated to ψ n ∈ (0, π). Hence the determined value of Zak phase can be either 0 or π. In this study, the unknown displacement fields that represent eigenmodes are determined from equations (4) and (20) for longitudinal and bending Bloch bands respectively. Such information is passed to ML network as binary digits 0, 1 and −1. Further details are given in the preceding sections.

Numerical modelling
The role of numerical modelling is limited to validation of the network findings. We employed COMSOL Multiphysics finite element codes to confirm the presence of interface states from network predicted δ values. The comparison of numerical results with TMM is available in Muhammad et al [2]. The supercell array of topologically predicted δ is built inside COMSOL Multiphysics, solid mechanics physics module and by performing frequency response study, the interface state is highlighted. COMSOL built-in fine tetrahedral mesh is applied to obtain the solution.

Autoencoder (AE) network architecture
A brief explanation on ML and DL network architecture developed in this study is provided. We developed an AE to train and predict the topological properties of phononic beam. An AE is a data compression neural network comprised of two parts i.e. encoder and decoder. The encoder compresses input data into a low-dimensional representation, termed as the bottleneck. Later, the decoder network reconstructs the original input data from the bottleneck, see figure 2. Generally, encoder and decoder networks are symmetrical. However, this is not necessary and they can be asymmetrical. In either case, the dimension of input layer and output layer are identical. Mathematically, the input x ∈ R n is multiplied by the weight matrix W i and summed with a bias b i to obtain the output of the fully connected layer as Where the superscript i refers to ith layer of AE and f i is the activation function. Here in this study we have used ReLU activation function, see table 2. The output a i from the ith hidden layer is used as input for (i + 1)th hidden layer and this process is continued until the last hidden layer of network to reconstruct the input data at the decoder output. The encoder maps the data from input space to latent space i.e. bottleneck layer h as While in an inverse process Γ −1 , the decoder maps the bottleneck back to the input space, see figure 2 The AE is trained in a self-supervised mode, where the network learns to reconstruct its input. The difference between actual input and reconstructed input at the decoder output layer is used to quantify the loss function. This loss function is minimized in order to improve the network performance. Network training is achieved by back-propagating sensitivities (gradient of a loss function with respect to the network parameters) from the output layer of decoder to the input layer of encoder, and subsequently updating the network parameters. Mathematically, the loss function can be defined as where ∥.∥ 2 denotes L 2 norm and m is number of training examples. Therefore, the job of encoder Γ (x) is to extract relevant features from a high dimensional input space to a low dimensional latent space. While the decoder Γ −1 (h) reconstruct the corresponding input from the feature vector h. The input data is generated through analytical formulations based on TMM to train and test the AE, see schematic workflow in figure 1. The topology of a beam is linked with shape parameter −1 ⩽ δ ⩽ 1. The δ is swept with a step size of 0.005 to generate thousands of data samples. Initially, the generated band structure data is in the form of normalized frequencyω and wavenumber k. To distinguish the Zak phase for each Bloch band, the band structure is converted into binary digits 0, 1 and −1 where 0 shows passband and 1 and −1 are associated to symmetric and anti-symmetric Bloch bands. The binary digit band structure data is arranged as label vector α with dimension 200 × 1 (longitudinal wave) and 500 × 1 (bending wave). The network training is performed in two steps: first the encoder is trained to map α to the bottleneck layer which approximates shape parameter δ. Second, the decoder which is connected to the encoder is trained, while the encoder is frozen and runs in inference mode.
The number of neurons in the input and output layers are same as the dimension of α. The input database is spilt into

Results and discussions
This section will discuss the results obtained from AE in two parts, separately for longitudinal and bending elastic waves. A brief summary of network performance in the form of R 2 and mean square error (MSE) are given in table 3. Overall, excellent agreement is observed that indicate AE is learning well the co-relation between band structure and phononic beam topology defined by δ.
First, we will discuss the training and testing of AE for correctly mapping the band structure and δ values for longitudinal wave propagating in the phononic beam. Later, same approach is adopted for bending elastic wave. To train the network, TMM and modified TMM based analytical dispersion relation is formulated, see sections 2.1 and 2.2. Using in-house MAT-LAB code, the δ is swept from −1 to 1 with step size of 0.005 and associated band structures are obtained. The band structure is in the form ofω and k. By using the analytical formulation based on reflection phase of Bloch band, the geometric Zak phase is identified. Since proposed phononic beam is symmetric with respect to central cross-section plane, therefore the Zak phase can have either 0 or π values depending upon symmetry of eigenmodes, see section 2.3. The band structure and Zak phase information is passed to AE input layer in the form of binary digits 0, 1 and −1 via label vector α. The number of neurons in the input layer is same as dimension of α. Figures 3  and 4 show some examples where normalized band structure is transformed into binary digit band structure with distinct geometric phase as 1 (Zak phase is 0 i.e. symmetric) and −1 (Zak phase is π i.e. anti-symmetric) for longitudinal and bending elastic waves propagating in the phononic beams.
The geometric Zak phase provides useful information about the presence of topologically protected interface states. The interface state is generated when mode transition frequency observed from Dirac cone degeneracy, see Muhammad et al  [2] for details, is common in the BG of topologically distinct phononic beams. In order to teach the algorithm about this concept, the Zak phase is presented in the form of 1 and −1. Otherwise from normalized band structure obtained from TMM, it is very difficult to distinguish the geometric phase of Bloch bands. The distinction is only possible through theoretical, numerical and experimental [26] means where symmetry type/reflection phase of Bloch band is investigated.
We generated approximately 4000 data samples for longitudinal and bending elastic wave separately to train and test the AE. The AE is developed and networking training and testing performance is quantified by calculating the R 2 and MSE. Figure 5 shows the learning curve for encoder and decoder for both types of waves. We have used 2500 epochs for longitudinal wave and 5000 epochs for bending wave. Excellent agreement is observed between training and test data. These curves validate that the AE is learning the correlation between binary band structures and δ values.
After AE training and testing, the network is assigned the task for making predictions. The prediction is made for δ values and this newly predicted results were compared with ground truth i.e. δ values swept in the analytical model to obtain normalized band structures. The predicted and target δ values comparison is performed for both longitudinal and bending waves as shown in figure 6. To aid understanding, the error histogram is also shown and one can observe with increase in instances, the error tends to move towards zero.
Next, we asked the network for combination of δ values that has mode transition frequency common in their BGs to predict the interface states. For both longitudinal and bending elastic waves, interestingly the network does output some possible combinations. The obtained δ values are passed to COM-SOL Multiphysics and a supercell phononic beam is constructed. In total 10 unit cell structures with distinct topology are considered to perform the frequency response study. For longitudinal wave, the harmonic excitation is applied at the left   end along the in-plane direction of the supercell array i.e. xdirection, see figure 1. The boundary probe is used to record the output displacement field at the right end. As reported in previous studies [2,5,26], at the topologically protected interface mode frequency, a sharp transmission peak is observed inside the BG. We observed sharp transmission peak designating topologically protected interface mode frequency inside the BG for predicted combination of δ values, see figure 7. The displacement field plots corresponding to these frequencies are also shown. In these plots, one can observe the localization of wave energy at the interface and decaying energy fields away from it.
Likewise, the supercell array consisting of 10 unit cell structures are constructed to validate the presence of topologically protected interface mode frequency in the phononic beams upon bending wave excitation. The harmonic out-ofplane excitation along z-direction (see figure 1) is applied at the left end and response in the form of out-of-plane displacement is recorded at the right end of the supercell array. For predicted δ values, interface state with sharp transmission peak inside the common BG of topologically distinct beam is observed as shown in figure 8. The displacement field plots depicted also show concentration of wave energy at the interface with decaying energy fields away from it. This validates the AE performance for correctly predicting the δ values that governs topologically protected interface states.
To further corroborate AE findings, we converted a 1 mlong steel rod into topologically protected phononic beam with δ 1 = 0.4 and δ 2 = 0.85. We choose steel for the experiment as the mechanical behaviour is closer to linear elasticity and the material strength also allows the manufacture of the prototype in a single piece. Further, steel produces a better signal to noise ratio and a smooth wave transmission spectrum can be observed. In addition, since the results are presented in the normalized frequency, the choice of aluminium or steel should not make much difference for predicted behaviour of interface modes. The length of each beam are as follows: L A1 = 30 mm, L B1 = 35 mm, L A2 =7.5 mm, L B2 = 46.25 mm. The diameter of beam is: d A = 4 mm, D B = 8 mm. The on-site and schematic diagram of experiment setup is shown in figure 9. The out of plane excitation is induced by white noise with frequency range (50-20 000 Hz) using PCB High Frequency Shaker Model 2025E-HF. The shaker is controlled by amplifier (PCB 2100E21-100 Smart Amplifier). The input signal is recorded     by accelerometers (PCB 352A60) in out-of-plane direction. The output signal is captured at the interface of topological beams by a Polytech PDV-100 laser vibrometer. The data is recorded and passed to the computer for post-processing using a NI-USB-4431 24 Bit Data Acquisition USB Module. We used compiled MATLAB programme to postprocess the recorded data.
The obtained experimental result is shown in figure 10. Since the selected beam configuration has two topologically protected interface modes, we observed sharp transmission peaks in response spectrum at 4764 Hz and 11 304 Hz that is indicative of topological protected interface modes. We also scan the beam with laser at multiple points to double check if the selected transmission peaks are interface modes. The obtained signals verify our prediction and selected transmission peaks around noisy signals are observed. This is reasonable and as you move away from the interface, the wave energy decay due to presence of BG. Further, we are not the first who performed an experimental test on topologically protected phononic beams. Prior studies [26,32,33] have also reported experimental results and for solid structures, the transmission signal is always noisy due to material anisotropy and the wide range of frequency spectrum considered. The contribution of the present study includes experimental test on phononic beam with topological properties predicted by a deep learning algorithm. Therefore, in terms of wave transmission peak that is an established approach to observe interface mode, the current numerical and experimental results show excellent agreement.

Conclusion
In summary, this work proposed the application of deep learning methods for the inverse design of topological phononic beams to predict topologically protected interface modes. To achieve this goal, a rigorous analytical model is established using the TMM for both longitudinal and bending elastic waves. The obtained constitutive equations are used to generate input data for deep learning network training and testing. In order to map a correlation between the beam geometry and band structure, the frequency spectrum data is converted into binary digits. The geometric Zak phase provide useful information about the band inversion and variances in symmetry of Bloch eigenmodes. Bloch bands with distinct Zak phase having mode transition frequency common in the BGs of topologically distinct phononic beams induce topologically protected interface modes, evident from the sharp transmission peak in frequency response spectrum. For the set of generated data, the Zak phase is analytically calculated by determining the reflection phase of Bloch eigenmode. The distinct Zak phase is differentiated as 1 (symmetric eigenmode) and −1 (anti-symmetric eigenmode) with reference to central crosssection plane of beam and this information is passed to the autoencoder. The network is trained based on generated data for both longitudinal and bending elastic waves. Upon successful training and testing, the interface mode prediction is made for combination of phononic beams. The network findings are validated by finite element based numerical simulations and experimental test on the manufactured topological beam. Overall, excellent agreement is observed. The study findings may provide a new insight into the application of deep learning methods in solid and computation mechanics for design and optimization of metamaterials based on beam and plate structures. Future studies will focus on applying a more holistic deep learning network for design optimization and performance enhancement of PnCs and metamaterials to control subwavelength acoustic and elastic waves.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.