Dielectric elastomers are active materials that undergo large deformations and change their instantaneous moduli when they are actuated by electric ﬁelds. By virtue of these features, composites made of soft dielectrics can ﬁlter waves across frequency bands that are electrostatically tunable . To date, to improve the performance of these adaptive phononic crystals, such as the width of these bands at the actuated state, metaheuristics-based topology optimization was used. However, the design freedom offered by this approach is limited because the number of function evaluations increases exponentially with the number of design variables. Here, we go beyond the limitations of this approach, by developing an efﬁcient gradient-based topology optimization method for maximizing the width of the band gaps in an exemplary case study. We employ a ﬁnite element formulation of the governing equations, and use the properties of each element as the design variables. In order to iteratively update the design variables, we employ gradient-based optimization, namely the Method of Moving Asymptotes. We carry out and implement fully analytical sensitivity analysis for computing the gradient of the objective function with respect to each one of the design variables. The numerical results of the method developed here demonstrate prohibited frequency bands that are indeed wider that those that were generated using metaheuristics-based topology optimization, while the computational cost to identify them is reduced by orders of magnitude. Keywords: Dielectric elastomers, nonlinear electroelasticity, wave propagation, band gaps, composites, Bloch-Floquet waves, gradient-based topology optimization 1. Introduction 1 Dielectric elastomers (DEs) are active materials that undergo large deformations when they 2 are actuated by electric ﬁelds (Pelrine et al., 2000; Hajiesmaili & Clarke, 2021). Moreover, the 3 constitutive relations of DEs are nonlinear, such that their instantaneous moduli vary as functions 4 of the electromechanical loads. These features have led to the design of various applications, such 5 as haptic devices (Kim et al., 2016), soft robots (Gu et al., 2018), actuators (Hajiesmaili & Clarke, 6 2021; Sharma et al., 2018, 2019; Su et al., 2018, 2019; Kashyap et al., 2020; Sharma et al., 2017; 7 Email addresses: atulksharma@iitj.ac.in (Atul Kumar Sharma), majd-costa@campus.technion.ac.il (Majd Kosta), meshmuel@technion.ac.il (Gal Shmuel), odedamir@technion.ac.il (Oded Amir) October 1, 2021 Ortigosa et al., 2021; Ortigosa & Martínez-Frutos, 2021; Khurana et al., 2021) and artiﬁcial muscles 8 (Lu et al., 2016), that are based on DEs. 9 More recently, the potential of DEs to function as adaptive waveguides has been explored (Zhu 10 et al., 2020; Mohajer et al., 2021; Jandron & Henann, 2018; Ziser & Shmuel, 2017; Wang et al., 11 2020; Zhu et al., 2020). These studies are based on the understanding that the physical properties 12 that govern elastic waves are electrostatically tunable. Following the pioneering work of Gei et al. 13 (2011) on ﬂexural waves, in a series of papers, Shmuel and collaborators (Shmuel & Debotton, 14 2012; Shmuel, 2013; Shmuel & Pernas-Salomón, 2016; Getz et al., 2017; Getz & Shmuel, 2017; 15 Bortot & Shmuel, 2017) have integrated this concept of material tunability with the theory of 16 acoustic band structure, in order to design adaptive wave suppressors. Band theory predicts that 17 waves in elastic composites (also termed phononic crystals) cannot propagate at certain frequencies 18 (band gaps) owing to Bragg reﬂections, that are functions of the periodicity and the impedance 19 mismatch of the phases (Kushwaha et al., 1993, 1994; Yang & Kim, 2018; Zhang & Gao, 2018; 20 Ren et al., 2020). In addition to this wave isolation characteristic in the band gap frequency region, 21 the phononic crystals with local resonators also exhibit the fracture resistance and crack-arrest 22 characteristics which enhance the strength of the structure (Huang et al., 2020, 2021). The idea to 23 the band gaps using DEs is simple: by applying a bias electromechanical ﬁeld, the periodicity and 24 impedance mismatch between the DE phases is tuned, which in turn change the Bragg reﬂections 25 and the desired band gaps (Shmuel & Band, 2016; Lustig & Shmuel, 2018). 26 To achieve a signiﬁcant tunability, aforementioned works have shown that large electric ﬁelds 27 are needed—a requirement that hinders the practical utility of DEs as tunable phononic crystals. 28 Since this tunability is a function of the microstructure of the composite, Bortot et al. (2018) 29 have developed a scheme to identify unit cells that exhibit better performance based on topology 30 optimization —a numerical method that systematically searches for optimal material distribution 31 that under speciﬁed design and response constraints. As explained in more detail later, here we 32 go beyond this meta-heuristic optimization approach, by developing a gradient-based optimization 33 scheme for the problem considered by Bortot et al. (2018), building upon the model, governing 34 equations and analysis by Shmuel (2013). 35 Research on topology optimization began with the seminal paper of Bendsøe & Kikuchi (1988) 36 and recent advancements are discussed in several review articles (Sigmund & Maute, 2013; Deaton 37 & Grandhi, 2014). For computational implementation, the design domain—a square unit cell in 38 the current context—is discretized into a ﬁnite number of squares and a design variable is attached 39 to each of their centroids. The value of the design variable then indicates which material is placed 40 at the corresponding position in space. This allows for signiﬁcant design freedom as there is no 41 assumption regarding the topology of the unit cell. Consequently, vast improvement in performance 42 can be achieved in comparison to unit cells of prescribed topology, e.g. a circular ﬁber within a 43 surrounding matrix. 44 The use of topology optimization for maximizing phononic band gaps was pioneered by Sigmund 45 & Søndergaard Jensen (2003). Their work was followed by many studies employing optimization 46 techniques for synthesizing phononic crystals such that desired band gap characteristics are obtained 47 (Halkjær & Jensen, 2006; Gazonas et al., 2006; Hussein et al., 2007; Bilal & Hussein, 2011; Liu 48 et al., 2014, 2016a; Li et al., 2016a; Liu et al., 2016b; Xie et al., 2017; Lu et al., 2017; Xie et al., 49 2017; Chen et al., 2018; Yi et al., 2019; Liu et al., 2020; Quinteros et al., 2021; Zhang et al., 50 2021). Detailed reviews on topology optimization of phononic crystals were provided by Yi & 51 Youn (2016) and Li et al. (2019). 52 2 The studies mentioned above aimed at optimizing the linear elastic response, while optimization 53 techniques have been employed only recently for band gaps that depend on nonlinear elastic deformations 54 (Hedayatrasa et al., 2016; Bortot et al., 2018; De Pascalis et al., 2020). 55 The limitation of these three studies, from the optimal design point of view, is that they are based 56 on metaheuristics, speciﬁcally Genetic Algorithms (GA). These methods utilize only function 57 evaluations and no gradient information is used to progress towards an optimum. Their main 58 drawback is that the number of function evaluations that is needed increases exponentially with 59 the number of design variables. Each evaluation consists of solving many generalized eigenvalue 60 problems for generating the band diagram. This computational bottleneck limits the design resolution, 61 i.e., the number of design variables that represent the material distribution, as can be seen in the 62 respective results (Hedayatrasa et al., 2016; Bortot et al., 2018; De Pascalis et al., 2020). 63 The full potential of topology optimization can only be accessed if gradient-based optimization 64 is used, so that the computational effort can be reduced by orders of magnitude and ﬁne design 65 resolution can be accommodated. Such efﬁcient approach is developed here for the case study 66 analyzed by Shmuel (2013), to which Bortot et al. (2018) applied GA optimization. The medium 67 that they analyzed comprises soft dielectric ﬁbers embedded periodically in a matrix made of a 68 different deformable dielectric. The composite is quasistatically actuated by application of axial 69 electric ﬁeld. Under the assumption of incompressible phases, Shmuel (2013) has found a closed- 70 form solution for the resultant deformation in terms of the bias electric ﬁeld and the properties 71 of the composite. Using a plane-wave expansion approach, Shmuel (2013) determined the band 72 structure of incremental anti-plane waves propagating in the deformed conﬁguration, the gaps of 73 which were optimized by Bortot et al. (2018) using a GA method. 74 Here, we introduce a ﬁnite element formulation of the incremental problem, and use the properties 75 of each element as the design variables, the objective of which is to maximize the width of the gap in 76 the audible frequency range. In order to iteratively update the design variables, we employ gradient- 77 based optimization, namely the Method of Moving Asymptotes—MMA (Svanberg, 1987). We 78 carry out and implement fully analytical sensitivity analysis for computing the gradient of the 79 objective function with respect to each one of the design variables. We generate optimized designs 80 with superior performance compared to those obtained by Bortot et al. (2018), while the computational 81 cost—in terms of the number of function evaluations—is smaller by orders of magnitude. 82 Before we proceed with details on the structure of this paper, it is important to note that the 83 constitutive model upon which we develop our equations neglects viscosity (Destrade & Saccomandi, 84 2004; Hong, 2011; Chiang Foo et al., 2012), to establish a starting point for gradient-based optimization 85 in a simpler setting. While our assumption is reasonable for silicons, it is not appropriate for 86 acrylics (Kornbluh & Pelrine, 2008), there viscosity is much more pronounced, and future work 87 will require the integration of viscoelastic theory into (see, e.g., the way in which Mohajer et al. 88 (2021) extended the work of Shmuel et al. (2012)). We also emphasize that the dynamic problem 89 considered here focuses exclusively on the anti-plane motion, i.e., assumes that the more complicated 90 in-plane motion is not excited. This assumption signiﬁcantly simpliﬁes the analysis; a more 91 complete investigation, which allows motions, is left for future work, and would require the framework 92 of Getz et al. (2017). 93 The paper is organized as follows. In Sec. 2 we provide the governing equations of electroelastic 94 composites (Dorfmann & Ogden, 2005; Suo et al., 2008; Gei et al., 2013; DeBotton et al., 2007; 95 Dorfmann & Ogden, 2010), upon which the forthcoming analysis is based. In Sec. 3 we recall the 96 solution of Shmuel (2013) to the quasistatic deformation of the periodic DE composite owing to 97 3 axial electric ﬁelds, and the equations developed therein for the incremental problem of anti-plane 98 shear waves. In Sec. 4, we present the ﬁnite element formulation for the solution of the incremental 99 problem, and speciﬁcally for the calculation of the band diagram of the deformed DE composite. 100 The contribution of this work is detailed in Sec. 5 and Sec. 6, where we present the analytical design 101 sensitivity analysis for computing the gradient of the objective function with respect to each one 102 of the design variables; implement the method using in-house MatLab code; and present numerical 103 examples. We conclude this paper with a summary of our results in Sec. 7. 104 2. Governing equations 105 To provide a self-contained report, here we brieﬂy summarize the governing equations pertaining 106 to nonlinear electro-mechanical deformation of dielectric elastomers following the theory of nonlinear 107 electroelasticity (Dorfmann & Ogden, 2005; Suo et al., 2008) and the associated linearized incremental 108 theory (Dorfmann & Ogden, 2010). 109 2.1. Nonlinear electroelasticity theory 110 Consider an arbitrary deformable dielectric composite body Bm[Bf=Bmade up of two 111 incompressible homogeneous dielectric phases BmandBf, and surrounded by vacuum. In the 112 undeformed conﬁguration, the dielectric composite occupies a region Bm 0[Bf 0=B0with the 113 boundary@B0. An arbitrary material particle in this conﬁguration is denoted by its position vector 114 X. When subjected to electromechanical loadings, the dielectric composite body is deformed and 115 occupies a regionBm t[Bf t=Btwith boundary @Bt. The position of a material particle in the 116 deformed conﬁguration is denoted by x, which is related to Xthrough the nonlinear deformation 117 mapx= (X;t). The components of deformation gradient are given by FiJ=@ i(X;t) @XJ. The 118 determinant of the deformation gradient relates the volume ratio of an inﬁnitesimal element in the 119 deformed state, dv, and its counterpart in the undeformed state, dV, such thatJ= det( F)>0. 120 The components of the right Cauchy-green tensor are CIJ=FiIFiJ. 121 For quasi-electrostatic processes, in the absence of free charges and currents, Gauss’s law and Faraday’s law are written as @di @xi= 0;@ej @xkijk= 0i; (1) respectively, where dandeare the electric displacement and the electric ﬁeld vectors in the 122 deformed conﬁguration, respectively; 0represents the zero vector and represents the Levi-Civita 123 symbol. In addition, the mechanical equilibrium equation in the deformed conﬁguration, in the 124 absence of mechanical body forces, is 125 @ij @xj= 0i; (2) where is the symmetric total Cauchy stress tensor consisting of both mechanical and electrical 126 contributions. 127 In this paper, we consider inﬁnite, periodic, dielectric composites, which posses internal, mechanical 128 traction and free charge interfaces between the two dielectric phases mandfwith jump boundary 129 conditions 130 [ [ij] ]nj= 0i;[ [di] ]ni= 0;[ [ei] ]njijk= 0k; (3) 4 where [ [] ] = ()m()frepresents the jump of quantities between the two phases, and nrepresents 131 the unit normal vector on the deformed interface and directed from phase mto phasef. 132 In Lagrangian description that uses Xas the independent variable, the total ﬁrst Piola–Kirchhoff stress tensor is related to the total Cauchy stress tensor by PiJ=JikFT kJand satisﬁes the Lagrangian form of the equilibrium equation @PiJ @XJ= 0i; (4) in the undeformed conﬁguration B0. The Lagrangian electric displacement and electric ﬁeld are given asDJ=JF1 JidiandEJ=FiJei, respectively. They satisfy the Lagrangian form of Gauss’s law and Faraday’s law @DJ @XJ= 0;@EJ @XKIJK= 0I; (5) respectively. 133 We consider non-dissipative materials for which the ﬁrst Piola-Kirchoff stress tensor Pand 134 Lagrangian electric ﬁeld Eare obtained constitutively in terms of the deformation gradient F 135 and Lagrangian electric displacement Dthrough an augmented energy density function (F,D) 136 as (Dorfmann & Ogden, 2005) 137 PiJ=@ @FiJp0FT iJ; EJ=@ @DJ; (6) wherep0is a Lagrange multiplier that accounts for the incompressibility constraint. 138 2.2. The linearized incremental theory 139 Following Dorfmann & Ogden (2010), we present the governing equations pertaining to the time dependent inﬁnitesimal increment in both elastic and electric displacement ﬁelds _x=_ (X;t), _D(X;t), respectively, superimposed on the static deformed conﬁguration characterized by the ﬁelds (X;t)andD(X;t). Henceforth, the inﬁnitesimal incremental quantities are denoted with the superposed dot. The push-forwards of increments in the total ﬁrst Piola–Kirchhoff stress, the Lagrangian electric displacement and the Lagrangian electric ﬁelds are ij=1 J_PiKFjK; di=1 JFiJ_DJ; ei=F1 Ji_EJ; (7) respectively, satisfying the incremental forms of Gauss’s, Faraday’s laws, and momentum balance 140 equation, 141 @di @xi= 0;@ej @xkijk= 0i;@ij @xj=@2_xi @t2; (8) respectively, where _x(x;t) := _ (X;t). 142 For an incompressible dielectric elastomer material, the linearized incremental constitutive laws 143 in terms of the increments andeare obtained using Eq. (6) and the push-forward operation as 144 ij=Cijklhkl+p0hji_p0ij+Bijkdk; ei=Blki_hkl+Aijdj;(9) 5 x1 x2 x3 a =A a m f f f f  X1 X2 X3 A A m f f f f (a) (b) Figure 1: Dielectric elastomer ﬁber composite in the (a) undeformed conﬁguration, and (b) deformed conﬁguration, when subjected to an axial electric ﬁeld. where Cijkl=1 JFjI@2 @FiI@FkJFlJ;Bijk=FjI@2 @FiI@DJF1 Jk; Aij=JF1 Ii@2 @DI@DJF1 Jj;(10) are the spatial constitutive tangent moduli, and hij=@_xi @xjis subjected to incompressibility constraint 145 hkk= 0. 146 3. Quasistatic ﬁnite deformation and incremental anti-plane shear waves in DE composites 147 As mentioned, using the above theory, Shmuel (2013) determined the quasistatic deformation of 148 the periodic DE composite owing to axial electric ﬁelds [see also the study by Sharma et al. (2021)], 149 and developed the equations that govern incremental anti-plane shear waves propagating through 150 the deformed composite. In order to provide a self-contained report, we concisely recapitulate the 151 quasistatic solution and incremental equations, to which we will employ a ﬁnite element solution 152 and develop the optimization process. 153 We consider a dielectric elastomer composite that is inﬁnitely periodic in the (x1;x2)plane. The 154 composite is made of incompressible DE ﬁbers along the x3direction with arbitrary cross-section 155 (phasef), that are embedded in a different incompressible DE matrix (phase m). The constitutive 156 response of the phases is given by the augmented Gentian free energy density function (Zhao et al., 157 2007; Sharma, 2020): 158 (FiJ;DP) =Jm 2ln 1FiJFiJ3 Jm +1 2"DPFiPFiQDQ; (11) where"is the dielectric permittivity, is the shear modulus of the dielectric elastomer, Jmis a dimensionless material parameter that accounts for the strain stiffening that arises from the limited extensibility of the elastomers. Figure 1a shows the schematic of a square DE composite unit cell 6 of lengthAin the undeformed conﬁguration. The composite is subjected to an average electric ﬁeld of magnitude e3in thex3direction, and is free to expand in the transverse (x1;x2)plane. From an experimental viewpoint, such actuation is achieved by the application of a voltage drop between two compliant coated electrodes on the top and bottom faces of the composite, located far enough one of each other to avoid edge effects. Shmuel (2013) has shown that the resultant deformation and electric ﬁeld are homogeneous, in the form F(m)=F(f)=2 40 0 00 0 023 5;e(m)=e(f)=2 40 0 e33 5; (12) where the in-plane stretch ratio is given by the solution of the nonlinear equation 159  24 = "e2 3; (13) where (p)=(p)J(p) m J(p) m224+ 3,() = (m)() + (f)()and (p)is the volume fraction of 160 phasep. Thus, up to the point of bifurcation, for any given value of applied electric ﬁeld e3and 161 properties of the composite, the solution of Eq. (13) delivers the resultant deformation. 162 On the top of deformed conﬁguration, we consider anti-plane waves propagating in the (x1;x2) 163 plane of the composite. The size of the composite in this plane is assumed large enough to be 164 analyzed as an inﬁnite periodic medium. Figure 2a shows the schematic of the DE ﬁber composite 165 unit cell in the deformed (x1;x2)plane. Let _x3(x;t)and_ (x;t)denote the anti-plane incremental 166 displacement and the incremental electric potential, respectively; the equations that govern these 167 two unknown ﬁelds are 168 @ @xj ~(x)@_x3(x;t) @xj~d(x)@_ (x;t) @xj! =(x)@2_x3 @t2; (14) @ @xj ~d(x)@_x3(x;t) @xj"(x)@_ (x;t) @xj! = 0; (15) respectively, where j= 1;2,~(x) = (x)2"(x)e2 3, and ~d(x) ="(x)e3. 169 Eqs. (14)-(15) are the starting point of our analysis: in the next Sec. we develop a ﬁnite element 170 formulation to solve them, and in turn obtain the band diagram of incremental anti-plane shear 171 waves in the DE composite. Later on, based on the ﬁnite element formulation, we will develop a 172 gradient-based topology optimization method, the objective of which is to identify optimal unit- 173 cells that generate the widest band gaps at designated electric ﬁelds. 174 4. Finite element formulation 175 We introduce next a ﬁnite element formulation with Bloch-Floquet periodic boundary conditions 176 (Kittel et al., 1996) to solve Eqs. (14)-(15). With _x3and_ denoting virtual incremental anti-plane 177 displacement and virtual incremental electric potential, respectively, the weak form statements 178 7 (a) (b) x1 x2 a a m f f f f 2 1 3 4 k1 k2     Figure 2: (a) Schematic description of the twice-symmetric square unit cell in the periodicity plane (x1;x2), and (b) the corresponding ﬁrst Brillouin zone with line segments to be scanned for determining the gaps extrema. corresponding to Eq. (14) and Eq. (15) are written as 179  Bt ~@_x3 @xj~d@_ @xj! @_x3 @xjds= Bt_x3@2_x3 @t2ds; (16)  Bt ~d@_x3 @xj"@_ @xj! @_ @xjds= 0; (17) respectively, where sis the area of the unit cell in the deformed conﬁguration. 180 The DE composite unit cell is discretized using ﬁnite elements, Bt=[Be t, and ﬁnite element 181 approximation for the incremental displacement and electric potential ﬁelds inside the element Be t182 are written using the V oigt notation as 183 _x3=Na_xa3;_ =Na_ a; (18) whereNarepresents the shape function associated with node aof the element, and _xa3and_ a 184 are the incremental anti-plane displacement and the incremental electric potential values at node 185 a. Substituting above ﬁnite approximations and the linearized constitutive relations of Eq. (9) into 186 Eqs. (16)-(17), the elemental level equations corresponding to momentum balance and Gauss law 187 are written as 188 [Kmm]f_xb3g+ [Kme]n _ bo = [M]@2_xb3 @t2 ; (19) [Kem]f_xb3g+ [Kee]n _ bo = 0; (20) whereKmmis the mechanical stiffness matrix, KmeandKemare the electro-mechanical coupled 189 8 stiffness matrices, Keeis the electric stiffness matrix, and Mis the mass matrix. The expressions 190 for these elemental stiffness matrices are given as 191 Kmm= Be t~@Na @xj@Nb @xjds; (21) Kme= Be t~d@Na @xj@Nb @xjds; (22) Kem= Be t~d@Na @xj@Nb @xjds; (23) Kee= Be t"@Na @xj@Nb @xjds; (24) M= Be tNaNbds: (25) Further, by assuming both the nodal incremental displacement and the incremental electric 192 potential to be time harmonic as 193 _xb3= ~xb(x)ei!t; _ b=~ b(x)ei!t; (26) and inserting these harmonic expressions into Eqs. (19)-(20), we get the incremental elemental 194 level equations governing the electro-elastic wave motion with angular frequency !as 195 [Kmm]f~xbg+ [Kme]n ~ bo !2[M]f~xbg= 0; (27) [Kem]f~xbg+ [Kee]n ~ bo = 0; (28) where ~xband~ bdenote the spatially dependent nodal incremental displacement and electric potential 196 ﬁelds, respectively. Eliminating ~ bfrom Eq. (27) by substitutingn ~ bo =[G] [Kem]f~xbgfrom 197 Eq. (28), we ﬁnally obtain the following eigenvalue problem 198 ([KmmKmeGKem]!2[M])f~xbg = ([K]!2[M])f~xbg= 0;(29) whereK=KmmKmeGKemandGdenotes the Moore–Penrose pseudoinverse of the electric 199 stiffness matrix Kee. 200 Using the Bloch theorem (Kittel et al., 1996), we specify the periodic boundary conditions for 201 analyzing the anti-plane wave propagation through the inﬁnite, periodic DE composite. For the DE 202 composite with square lattice shown in Fig. 2a, the Bloch periodic boundary conditions for the 203 9 nodal displacements are expressed as (Vatanabe et al., 2014; Jandron & Henann, 2018) 204 ~xb(x2) = exp(ik1a)~xb(x1); ~xb(x4) = exp(ik2a)~xb(x3):(30) Similarly, the Bloch periodic boundary conditions for the nodal electric potential ﬁelds are expressed 205 as 206 ~ b(x2) = exp(ik1a)~ b(x1); ~ b(x4) = exp(ik2a)~ b(x3);(31) wherek1andk2are the component of real Bloch wave vector k=kieialongx1andx2directions, 207 as shown in Fig. 2b. In order to implement these complex-valued boundary conditions, we used 208 the augmented penalty method (Felippa, 2001; Oliveira et al., 2010). 209 The eigenvalue problem stated in Eq. (29) along with the Bloch periodic boundary conditions 210 of Eqs. (30)-(31) is solved for computing the band diagram, i.e, !as a function of the wave vector 211 k. Due to the quarter symmetry of the problem considered in this paper (see Fig. 2a), it is sufﬁcient 212 to consider only the wave vector along the edges [XMK]of the ﬁrst Brillouin 213 zone, marked by red lines in Fig. 2b (Sigmund & Søndergaard Jensen, 2003; Vatanabe et al., 2014; 214 Meng et al., 2017; Yi & Youn, 2016). The coordinates of various points along these edges are as 215  = (0;0),X = (=a;0),M = (=a;=a ), and K = (0;=a). 216 Having at hand a ﬁnite-element formulation to calculate the band diagram, we proceed next to 217 develop the gradient-based topology optimization method, whose aim is to ﬁnd the microstructure 218 that results with the widest band gaps at prescribed electric ﬁelds. 219 5. Topology optimization 220 In the design of phononic band gap structures, we need to ﬁnd the unit cell that yields the 221 desired band gap characteristics. In this section, we introduce the topology optimization method 222 for obtaining the optimized distribution of the matrix and ﬁber phases in the DE composite unit cell 223 that maximizes the band gap size. 224 5.1. Material property interpolation for topology optimization 225 In this work, the design domain is a unit cell that is discretized by ﬁnite elements. To interpolate 226 the material properties for each ﬁnite element used to discretize the design domain, artiﬁcial design 227 variablese2[0;1](e= 1;2;:::;N ,Nbeing the total number of ﬁnite elements in the design 228 domain) are assigned to all ﬁnite elements as shown in Fig. 3. 229 For example, e= 0 represents an element ethat is composed of matrix phase, while e= 1 230 represents an element consisting of ﬁber phase. As discussed by Sigmund & Søndergaard Jensen 231 (2003) and Meng et al. (2017), the material properties for each ﬁnite element in the design domain 232 can be interpolated linearly between the matrix and ﬁber phases as 233 e= (1e)(m)+e(f); (32) e= (1e)(m)+e(f); (33) 10 Figure 3: Discretization of the unit cell into Nﬁnite elements and assignment of a design variable to each ﬁnite element. "e= (1e)"(m)+e"(f); (34) Je m= (1e)J(m) m+eJ(f) m: (35) 5.2. Optimization problem 234 The topology optimization problem aims to maximize the band gap width between the two adjacent bands in the audible frequency range 0-10 kHz. For two adjacent bands, assuming the eigenfrequencies corresponding to overlying band and underlying band to be !j+1and!j, respectively, the objective function may be written as f=8 >>>>< >>>>:10kHzmax!j(;k)=2 ifmin!j+1(;k)=210kHz; (min!j+1(;k)max!j(;k))=2 ifmin!j+1(;k)=210kHz;(36) where, max!j(;k)andmin!j(;k)represent the maximum and minimum of the jtheigenfrequency 235 !j, respectively, over the kvector for a given design of the unit cell domain, . 236 We note that this objective function is not strictly differentiable because the location of maximum 237 and minimum eigenvalues over the k-vector may change from one design iteration to another. The 238 problem can be regularized using smooth maximum and minimum functions, as described in the 239 appendix. The results we achieved with the smooth formulation are essentially identical to those 240 achieved with the objective of Eq. (36) and therefore we limit the discussion to a brief comparison 241 in the appendix. It should be noted that in most gradient-based topology optimization studies on 242 band gaps, the issue of differentiability is not discussed and we presume that the results reported 243 in the literature were obtained with standard maximum and minimum operators. One exception we 244 found is the work of Qian & Sigmund (2011) where a smooth bound formulation is employed. 245 11 The topology optimization problem may be stated as 246 max f s:t:: K(k)!2M ~xb= 0;k2[XMK] (37) 0e1; e = 1;:::::;N: 5.3. Design sensitivity analysis 247 In order to iteratively update the design variables, we utilize gradient-based optimization, speciﬁcally 248 MMA (Svanberg, 1987). This requires the evaluation of the gradient of the objective function with 249 respect to any individual design variable, e. The derivative of the objective function fin Eq. (36) 250 with respect to the design variable eis obtained as 251 @f @e=8 >>>>< >>>>:1 2@(max!n(;k)) @e;if min!n+1(;k)=210kHz; 1 2@(min!n+1(;k)) @e@(max!n(;k)) @e ;if min!n+1(;k)=210kHz: (38) For evaluating the objective function derivative in Eq. (38), the derivative of the eigenfrequency 252 with respect to the design variable must be evaluated. For a given wave vector k, the derivative of 253 an eigenfrequency !nwith respect to the design variable is computed by differentiating Eq. (29) 254 with respect to the design variable e, as follows 255 @(!n(;k)) @e=1 2!n~xT b@K(k) @e!2 n@M @e ~x b; (39) where, ~x bis the global mass normalized eigenvector, and the derivative of the condensed stiffness 256 matrixKwith respect to the design variable is obtained as 257 @K @e=@Kmm @e@Kme @e(Kee)1Kem+ Kme(Kee)2Kem@Kee @eKme(Kee)1@Kem @e:(40) The derivatives of the stiffness and the mass matrices@Kmm @e,@Kme @e,@Kem @e,@Kee @e, and@M @e258 in Eqs. (39)-(40), are evaluated using Eqs. (21)-(25) as follows 259 @Kmm @e= Be t@~e @e@Na @xj@Nb @xjds; (41) @Kme @e= Be t@~de @e@Na @xj@Nb @xjds; (42) 12 @Kem @e= Be t@~de @e@Na @xj@Nb @xjds; (43) @Kee @e= Be t@"e @e@Na @xj@Nb @xjds; (44) @M @e= Be t@e @eNaNbds; (45) in which the calculation of the derivatives of material parameters e,eand~dewith regard to eis 260 straightforward from Eqs. (32) and (34):@e @e=(f)(m),@"e @e="(f)"(m), and@~de @e=e3@"e @e, 261 while the evaluation of the derivative of ~ewith regard to eis more involved, as follows1262 @~e @e= 2eJe m@ @e+2 e@Je m @e+@e @eJe m  Je m224+ 3 eJe m2@Je m @e4@ @e+ 45@ @e  Je m224+ 32@"e @ee2 3;(46) where@Je m @e=J(f) mJ(m) m, and the derivative of the stretch parameter with respect of design 263 variable@ @eis evaluated by differentiating Eq. (13) with respect to eand expressed as 264 @ @e="  "(f)"(m) e2 3(24) e(J(f) mJ(m) m) +Je m((f)(m)) (Je m224+ 3)+eJe m(J(f) mJ(m) m) (Je m224+ 3)2!#  (22+ 45)NP e=1e + (24) (445)NP e=1eJe m (Je m224+ 3)2: (47) 1Note that the expression (46) corresponds to ~eandeof the same ﬁnite element. However, in the evaluation of ~ethere exists a weak coupling to other elements due to the solution of that is a global state variable. Hence, for two distinct elements denoted e1ande2, the derivative is @~e1 @e2=2e1Je1m@ @e2 (Je1m224+ 3)e1Je1m2h 4@ @e2+ 45@ @e2i (Je1m224+ 3)2: These coupling terms are several orders of magnitude smaller than the terms of (46) and were found to have no inﬂuence on the optimized design. Hence they are omitted from subsequent derivations. 13 Initial material distribution and design of unit cell Finite deformation of infinite periodic DE composite Finite element analysis for the evaluation of band gap Evaluation of the objective function Convergence check Design sensitivity analysis Solution of the optimization problem New design of unit cell Yes End No Figure 4: A schematic representation of the optimization process. 5.4. Implementation 265 The topology optimization framework presented in this paper is implemented by developing 266 an in-house MATLAB code. The topology optimization process starts with generating one quarter 267 of a twice-symmetric structure of the unit cell. In order to generate a symmetric unit cell, we 268 randomly distribute design variables in one quarter of the unit cell, and assign these design variables 269 symmetrically in the other three quarters. We note that this initial unit cell structure may not have 270 gaps in its spectrum. Then, based on the ﬁnite element framework introduced in Section 3, we 271 obtain the band diagram and the corresponding eigenvectors. Accordingly, the objective function 272 is evaluated. Subsequently, we evaluate the gradient of the objective with respect to the design 273 variables using Eq. (38). According to the evaluated objective function gradients, MMA proposes 274 a new vector of design variables by solving a convex subproblem. This yields the formation of a 275 14 Table 1: Material properties. Material Silicone CF19-2186 Polyurethane PT6100S Shear modulus (kPa) 333 5667 Locking parameter Jm 46.3 6.67 Relative permittivity r 2.8 7 Density(kg=m3) 1100 1200 new structure with the updated design variables. Thereafter, we repeat the ﬁnite element analysis 276 for the evaluation of band gaps and MMA optimization scheme for updating the design variables. 277 Such an iterative process continues until the maximum change in any design variable is less than 278 106. Figure 4 depicts schematic representation of the iterative optimization process. 279 6. Numerical results and discussions 280 In this section, we provide the numerical results for the phononic band gaps in dielectric 281 elastomers designed by using the topology optimization framework presented in the previous section. 282 The dielectric elastomer composite is assumed to be composed of polyurethane PT6100S ﬁbers 283 embedded in the silicone CF19-2186 matrix (Getz et al., 2017). The properties of these materials 284 are listed in Table 1. In the undeformed conﬁguration, we set the lattice parameter or the size of 285 the square unit cell to be A= 6:3mm. In order to discretize the unit cell, we used the four nodes, 286 linear quadrilateral ﬁnite elements with two degrees of freedom (one anti-plane displacement and 287 one electric potential) per node. 288 6.1. Validation of the ﬁnite element framework 289 First, for validating the ﬁnite element framework presented in Section 4 and investigating the 290 inﬂuence of the bias electric ﬁeld on the band structure, we compute the performance of a standard 291 unit cell composed of circular ﬁbers embedded in the matrix, as shown in Fig. 5a. The volume 292 fraction of ﬁbers in the unit cell is taken to be 0.5. The unit cell is discretized by 1139 quadrilateral 293 ﬁnite elements as shown in Fig. 5a (left). In all the plots in this paper, the yellow color denotes the 294 ﬁber and the blue color denotes the matrix. For investigating the effect of the bias electric ﬁeld on 295 the band diagram, we perform two cases, one without bias electric ﬁeld ( e3= 0) and the other with 296 bias electric ﬁeld e3= 300 MV/m. For both cases, without and with bias electric ﬁeld, Figs. 5b and 297 5c, respectively, compare the band structure computed using the ﬁnite element framework presented 298 in this work and the plane wave expansion method (PWE) (Shmuel, 2013). The plot shows a very 299 good agreement between FEM and PWE, validating our ﬁnite element framework. From the band 300 diagrams, it is evident that the width of the widest gap at e3= 0, i.e.,( !=2= 0:7626 kHz), is 301 smaller than that at e3= 300 MV/m, i.e., ( !=2= 0:8073 kHz), denoting an increase in the band 302 gap width of 5.86%. This example suggests that the bias electric ﬁeld has a favorable impact as it 303 improves the maximization of band gap widths. 304 6.2. Topology optimization results 305 In this section, ﬁrst, we present the topologies for the maximization of the band gap width 306 between the ﬁrst four adjacent bands (i.e., 1stand2ndbands, 2ndand3rdbands, 3rdand4thbands, 307 15 (a) (b) (c) 0.01.53.04.56.07.5/s119/2π (kHz)k ΓΧ Μ Γ PWE FEM 0.01.53.04.56.07.5/s119/2π (kHz)k ΓΧ Μ Γ PWE FEM Figure 5: Inﬂuence of the bias electric ﬁeld on the band structure of a DE composite with circular ﬁber and comparison of FEM results (solid lines) with the plane wave expansion method (PWE) results (symbols): (a) Unit cell (left) and periodic matrix (right), band diagrams (b) when bias electric ﬁeld e3= 0, and (c) when bias electric ﬁeld e3= 300 MV/m. &4thand5thbands) and identify the two adjacent bands that have the largest optimized band gap 308 in the audible frequency range. For these four optimization cases, we consider two values of the 309 applied electric ﬁeld, i.e., e3= 150 MV/m ande3= 300 MV/m, and discretize the unit cell by 310 2020, bi-linear quadrilateral elements. The initial material distribution of the unit cells is taken 311 to be random. The resulting designs of the unit cells that optimize the band gap widths between 312 the ﬁrst four adjacent bands at applied electric ﬁeld e3= 150 MV/m ande3= 300 MV/m are 313 depicted in the left panels of Figs. 6 and 7, respectively. In both ﬁgures, the middle panels represent 314 the periodic matrices composed of four optimized unit cells, while the right panels show the band 315 diagrams corresponding to the optimized unit cells depicted in the left panels. The estimates of 316 the optimized band gaps between the ﬁrst four adjacent bands and the corresponding ﬁber volume 317 fractions for both the applied electric ﬁeld values are listed in Table 2. From Table 2 and the left 318 panels of Figs. 6 and 7, we observe the maximum band gap width when the band gap is optimized 319 between the 4thand5thbands. Ate3= 150 MV/m, the value of the maximal band gap width in 320 16 Table 2: Optimized band gaps and ﬁber volume fractions, when maximizing gaps between two different adjacent bands. Bias electric ﬁeld ( e3)! 150 MV/m 300 MV/m Band gap Fiber fraction Band gap Fiber fraction 1stand 2ndbands 1.58 0.65 1.71 0.66 2ndand 3rdbands 1.94 0.58 2.25 0.57 3rdand 4thbands 1.73 0.57 2.08 0.58 4thand 5thbands 4.31 0.67 4.61 0.68 the audible frequency range is 4.31kHZ and the corresponding ﬁber volume fraction is 0.67. The 321 maximal band gap width achieved in the audible frequency range at e3= 300 MV/m, is 4.61kHZ 322 with 0.68 ﬁber volume fraction. 323 Next, for the identiﬁed adjacent bands (i.e., 4thand5thbands) representing maximum optimized 324 band gap in the audible frequency range, we present the optimized designs for different initial 325 layouts of the unit cells. We consider three initial layouts of the unit cells (fully random, single 326 central square ﬁber inclusion, and ﬁber inclusions at multiple locations) as shown in Table 3. For 327 these three initial layouts, the iteration history of the objective function and the evolution of the 328 optimal topology of the unit cells at applied electric ﬁelds e3= 150 MV/m ande3= 300 MV/m 329 are represented in Figs. 8 and 9, respectively. In both ﬁgures, the left columns show the iteration 330 history of the objective function and the evolution of the unit cell structure for a random initial 331 material distribution; the middle columns show the same data for an initial design consisting of 332 a central square ﬁber inclusion; and the right columns show the same data for an initial design 333 with ﬁber inclusions at multiple locations. From Figs. 8 and 9, we observe that the optimized 334 topologies depend on the initial layouts. This phenomenon of dependency on the initial layouts for 335 the unconstrained band gap maximization problem has been previously described in the literature 336 (Yi et al., 2019; Li et al., 2016b) and is somewhat expected in gradient-based optimization. From 337 the plots of the iteration histories of the objective function, we observe that the objective function 338 increases rapidly at the starting stage of the topology optimization process and then the convergence 339 rate becomes signiﬁcantly slowed down over the iterations. The estimates of the optimized band 340 gaps and the corresponding ﬁber volume fractions for all three initial layouts and both the values of 341 the applied electric ﬁeld are listed in Table 3. As evident from Table 3, the third initial guess (ﬁber 342 inclusions at multiple locations) for the unit cell layout is best among all the three initial guesses, 343 as it shows the maximum optimized band gap for both values of the applied electric ﬁelds. The 344 maximum band gap achieved at e3= 150 MV/m is 4.37 kHz and that at e3= 300 MV/m is 345 equal to 4.61 kHz. The ﬁber volume fractions corresponding to optimized layouts are the same for 346 both applied electric ﬁeld and are equal to 0.68. As evident from the convergence history plots, 347 for the third initial layout of the unit cell, the topology optimization process converges faster in 348 comparison to the other two initial layouts. This is expected because it utilizes an educated initial 349 guess. 350 Next, we investigate the dependence of the optimized band gaps and the unit cell designs 351 on the FE mesh resolution, that affects the design resolution. For this purpose, we study the 352 maximization of band gap width between the 4thand5thbands (i.e., having the highest band gaps in 353 the audible frequency range). For mesh dependency analyses, we discretize the unit cell by 20 20, 354 17 (a) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (b) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (c) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (d) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ Figure 6: Unit cells (left), periodic matrices composed of 4 4 unit cells (middle), and band diagrams (right) of the optimized results when maximizing the band gap between (a) 1stand2ndbands, (b) 2ndand3rdbands, (c)3rdand4thbands, and (d) 4thand5thbands, when bias electric ﬁeld e3= 150 MV/m. 3030, and 4040 bi-linear square ﬁnite elements, and perform the topology optimizations for 355 maximizing the band gaps at the aforementioned two values of the applied electric ﬁeld. For all 356 18 (a) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (b) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (c) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (d) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ Figure 7: Unit cells (left), periodic matrices composed of 4 4 unit cells (middle), and band diagrams (right) of the optimized results when maximizing the band gap between (a) 1stand2ndbands, (b) 2ndand3rdbands, (c)3rdand4thbands, and (d) 4thand5thbands, when bias electric ﬁeld e3= 300 MV/m. three mesh densities, we consider the same initial layout of the unit cell, i.e., the third layout of the 357 aforementioned discussion for which we obtained the maximum at both the applied electric ﬁelds. 358 19 Table 3: Effect of initial material layout on the optimized band gaps and ﬁber volume fractions for maximization of band gap between 4thand5thbands. Initial material layout ! Bias electric ﬁeld ( e3)Band gap Fiber fraction Band gap Fiber fraction Band gap Fiber fraction 150 MV/m 4.31 kHz 0.67 4.30 kHz 0.68 4.37 kHz 0.68 300 MV/m 4.61 kHz 0.68 4.51 kHz 0.68 4.61 kHz 0.68 Table 4: Effect of mesh size on the optimized band gaps and ﬁber volume fractions for maximization of band gap between 4thand5thbands. Mesh density! 2020 30 30 40 40 Bias electric ﬁeld ( e3)Band gap Fiber fraction Band gap Fiber fraction Band gap Fiber fraction 150 MV/m 4.37 kHz 0.680 4.33 kHz 0.715 4.30 kHz 0.680 300 MV/m 4.61 kHz 0.680 4.56 kHz 0.706 4.53 kHz 0.680 The optimized designs of the unit cell and the corresponding band diagrams for 20 20, 3030, 359 and 4040 mesh densities and at e3= 150 MV/m are compared in the left, middle and right panels 360 of Fig. 10a, respectively. For both values of bias electric ﬁeld, the estimates of the optimized band 361 gaps and the corresponding ﬁber volume fractions at the aforementioned three mesh densities are 362 listed in Table 4. From the table, it is evident that the value of the optimized band gap decreases 363 with an increase in the mesh density for both the values of the bias electric ﬁeld. This is attributed 364 to the effect of FE mesh reﬁnement on the evaluation of the band gap. Interestingly, at 20 20 and 365 4040 mesh density, we obtain the same designs of the unit cell and ﬁber volume fractions for both 366 values of the applied electric ﬁeld. 367 Further, for a one-to-one comparison of the optimized unit cell structures obtained at different 368 mesh sizes, we evaluate the band gaps for all three optimized structures on a 120 120 mesh 369 resolution. This resolution is chosen such that all optimized designs can be projected precisely 370 to a common mesh. Table 5 enlists the band gap values evaluated on a 120 120 mesh resolution 371 for the optimized unit cells obtained at the aforementioned three mesh densities and two values 372 of bias electric ﬁeld. The band diagrams evaluated at 120 120 mesh resolution of the unit cells 373 optimized at three mesh sizes are shown in Fig. 11. From Table 5, we can see that as the design 374 resolution increases, so does the width of the optimized band gap. This demonstrates the efﬁcacy 375 of the topology optimization method: better performance can be achieved by increasing the design 376 resolution, hence expanding the design space. The unit cell layouts obtained from the optimization 377 at 4040 mesh density show the maximum band gap for both bias electric ﬁelds. For 120 120 378 mesh resolution, the value of maximum band gap achieved at e3= 150 MV/m is 4.24 kHz with 379 0.68 ﬁber volume fraction and that at e3= 300 MV/m is equal to 4.44 kHz with the same ﬁber 380 volume fraction of 0.68. 381 Finally, we discuss the advantages of the gradient-based topology optimization framework 382 developed in this paper over the method based on PWE and a genetic algorithm (GA) for maximization 383 20 (a) (b) (c) 081624324048-1012345482 05 Objective functionI teration number0 0102030405058-101234558205 Objective functionI teration number0 06 12182427-10123455 272 00 Objective functionI teration number Initial state (Iteration=0) Initial state (Iteration=0) Initial state (Iteration=0) Intermediate state (Iteration=5) Intermediate state (Iteration=5) Intermediate state (Iteration=5) Intermediate state (Iteration=20) Intermediate state (Iteration=20) Intermediate state (Iteration=20) Final state (Iteration=48) Final state (Iteration=58) Final state (Iteration=27) Figure 8: Optimization history and the layout evolution history of the unit cells for maximizing the band gap between 4thand5thbands at bias electric ﬁeld e3= 150 MV/m, for various initial layouts of the unit cells: (a) fully random; (b) central square ﬁber inclusion; and (c) ﬁber inclusions at multiple locations. of band gaps in dielectric elastomers (Bortot et al., 2018). For this purpose, we compare the 384 optimized designs of the unit cells and the corresponding optimized band gaps obtained using 385 21 (a) (b) (c) 091827364554-10123455 4205 Objective functionI teration number0 081624324044-10123454 42 05 Objective functionI teration number0 08 1 62 43 2-10123455 322 00 Objective functionI teration number Initial state (Iteration=0) Initial state (Iteration=0) Initial state (Iteration=0) Intermediate state (Iteration=5) Intermediate state (Iteration=5) Intermediate state (Iteration=5) Intermediate state (Iteration=20) Intermediate state (Iteration=20) Intermediate state (Iteration=20) Final state (Iteration=54) Final state (Iteration=44) Final state (Iteration=32) Figure 9: Optimization history and the layout evolution history of the unit cells for maximizing the band gap between 4thand5thbands at bias electric ﬁeld e3= 300 MV/m, for various initial layouts of the unit cells: (a) fully random; (b) central square ﬁber inclusion; and (c) ﬁber inclusions at multiple locations. both methods of optimization. For one-to-one comparison, we consider the design structures 386 optimized at the same mesh density i.e., 20 20. Figures 12a and 12b show the designs of the 387 22 2020 30 30 40 40 (a) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ (b) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ 0246810/s119/2π (kHz)k ΓΜ Κ ΓΧ Figure 10: Effect of mesh size (left, 20 20; middle, 3030; right, 4040) on the optimized layouts of the unit cells for maximizing the band gap between 4thand5thbands, when bias electric ﬁeld (a) e3= 150 MV/m, and (b)e3= 300 MV/m. optimized unit cell obtained using gradient-based and GA-based topology optimization methods, 388 respectively, at bias electric ﬁeld e3= 150. The ﬁber volume fraction for the unit cell optimized 389 using a gradient-based optimization scheme is 0.68 and while for that optimized using GA is equal 390 to 0.64. The band diagrams for optimized unit cells in Figs. 12a and 12b are shown in Fig. 12c. 391 From the band diagrams, we observe the larger band gap for the unit cell optimized using the 392 23 (a) ( b) 0246810/s119/2π (kHz)k Optimization mesh 20X20 Optimization mesh 30X30 Optimization mesh 40X40Γ Χ Μ Κ Γ 0246810/s119/2π (kHz)k Optimization mesh 20X20 Optimization mesh 30X30 Optimization mesh 40X40Γ Χ Μ Κ Γ Figure 11: Comparison of the band diagrams for 120 120 mesh resolution of the unit cells optimized at 2020, 3030, and 4040 mesh densities, (a) when bias electric ﬁeld e3= 150, and (b) when bias electric ﬁelde3= 300 MV/m. Table 5: Band gap between 4thand5thbands at 120120 mesh resolution for the unit cells optimized on 2020, 3030 and 4040 mesh densities. Optimization mesh density ! 2020 30 30 40 40 Bias electric ﬁeld ( e3)! 150 MV/m 300 MV/m 150 MV/m 300 MV/m 150 MV/m 300 MV/m Optimized unit cells Band gap at 120120 mesh resolution 4.17 kHz 4.38 kHz 4.22 kHz 4.43 kHz 4.24 kHz 4.44 kHz gradient-based formulation (4.37 kHz) in comparison to the one optimized using GA (4.13 kHz). 393 The present approach also holds a signiﬁcant computational advantage: the computational cost of 394 the GA is higher by several orders of magnitude, and this limits the achievable design resolution. 395 For example, Bortot et al. (2018) reported that 12,500 function evaluations were required for GA 396 to converge in a problem with 100 design variables. Herein, these problems are solved in less 397 than 60 design cycles, where each cycle requires a single function evaluation and the computation 398 of gradients. The cost of computing gradients can be estimated as roughly the same as a single 399 function evaluation. As shown above, larger design freedom can yield better band gap performance, 400 and this can only be achieved with the gradient-based approach that easily scales to ﬁner mesh 401 resolutions. 402 7. Summary 403 When subjected to electric ﬁelds, dielectric elastomers undergo ﬁnite deformations and their 404 instantaneous moduli change. These features can be exploited to electrostatically tune the width of 405 band gaps, across which waves cannot propagate in dielectric elastomer composites. 406 24 (a) (b) (c) 0246810/s119/2π (kHz)k Gradient based optimization [Our results] GA based optimization [Bortot et al. 2018] Γ Χ Μ Κ Γ 4 .37 kHz4 .13 kHz Figure 12: Designs of the unit cells optimized using (a) the gradient based optimization formulation presented in this paper and (b) a GA-based optimization method (Bortot et al., 2018). (c) Comparison of the corresponding band diagrams. The width of these gaps in the actuated state is a function of the microstructure of the composite. 407 To date, to identify optimal microstructures with wider gaps in the actuated state, topology optimization 408 methods based on metaheuristics, speciﬁcally genetic algorithms, were used. These methods, 409 which utilize only function evaluations without any gradient information, are limited by the computational 410 cost that increases exponentially with the number of design variables. 411 In this work, we have developed a gradient-based topology optimization method, which utilizes 412 information from the gradient of the objective function in order to iteratively update the design 413 variables (Svanberg, 1987). The case study to which we have developed the method is of incremental 414 anti-plane shear waves, propagating in a planar dielectric elastomer composite that is subjected to 415 25 axial electric ﬁelds (Shmuel, 2013; Getz & Shmuel, 2017; Bortot & Shmuel, 2017). To this end, 416 we have employed a ﬁnite element formulation to obtain the eigenmodes of the linearized problem, 417 and extract the band diagram (Benchabane et al., 2006; Li et al., 2012; Veres et al., 2013). Based on 418 ﬁnite element formulation, we have derived and implemented a fully analytical sensitivity analysis 419 for computing the gradient of the objective function: the width of the gap in the audible range, at 420 two prescribed electric ﬁelds. 421 In comparison with the genetic algorithm-based method (Bortot et al., 2018), our gradient-based 422 method not only identiﬁed microstructures with wider gaps: it also reached these microstructures 423 in an orders-of-magnitude cheaper computational cost. Thus, the gradient-based approach easily 424 scales to ﬁner mesh resolutions. 425 As a potential future work, the present topology optimization framework can extended to 426 optimize band gaps associated with in-plane motions (Getz et al., 2017); and account for the 427 viscosity of the elastomers by the integration of viscoelastic theory (Mohajer et al., 2021; Destrade 428 & Saccomandi, 2004; Chiang Foo et al., 2012; Hong, 2011). 429 Acknowledgement 430 The authors are thankful to the anonymous reviewers for their valuable input. This research 431 was supported by the Israel Science Foundation, funded by the Israel Academy of Sciences and 432 Humanities (Grant no. 2061/20), the United States-Israel Binational Science Foundation (Grant 433 no. 2014358), and Ministry of Science and Technology (grant no. 880011). 434 Appendix: Smooth formulation with approximate maximum and minimum eigenfrequencies 435 along the wave vector k 436 The objective function of the optimization problem in Eq. (36) is not strictly differentiable, 437 because the k-vector(s) corresponding to the maximum and minimum eigenfrequencies may change 438 during the optimization iterations. Here, we present a topology optimization formulation in which 439 the maximum and minimum values of eigenfrequencies along the wave vector kare approximated 440 using ap-norm function, and compare the results for one test case to those in Section 5. We examine 441 the maximization of the band gap between the fourth and ﬁfth bands which have a maximum 442 optimized band gap in the audible frequency range. 443 Utilizing the standard p-norm measure (Duysinx & Sigmund, 1998; Le et al., 2010), the maximum 444 and minimum values of the eigenfrequencies corresponding to lower and upper bands, respectively, 445 are approximated as 446 !PN max= nX i=1(!i)p!1 p ;!PN min=1 nP i=1 1 !ip1 p; (A.1) where,pis thep-norm parameter, !irepresents the value of the eigenfrequency corresponding to 447 theithpoint along the wave vector k. In this work, we took p= 80 andn= 40 points along the 448 wave vector k. 449 26 (a) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ 06121824303641-10123455 302 00 Objective functionI teration number41 (b) 0246810/s119/2π (kHz)k ΓΜ Κ Γ Χ 011223344556677-10123455 77200 Objective functionI teration number40 Figure 13: Unit cells (left), band diagrams (middle), and objective function history (right) when maximizing the band gap between 4thand5thbands using the p-norm formulation, (a) when bias electric ﬁeld e3= 150 MV/m, and (b) when bias electric ﬁeld e3= 300 MV/m, The modiﬁed version of the objective function is written as 450 fPN=8 >>>>< >>>>:10kHz!PN max=2 ifmin!j+1(;k)=210kHz; (!PN min!PN max)=2 ifmin!j+1(;k)=210kHz:(A.2) and the sensitivities of the maximum and minimum frequencies approximated using p-norm with 451 respect to design variable eare expressed as 452 @!PN max @e= nX i=1(!i)p!1 p1 nX i=1(!i)p1@!i @e! ; @!PN min @e= nX i=11 !ip!11 p nX i=11 !ip+1@!i @e! ;(A.3) in which the derivative@!i @eis evaluated using Eq. (39). Further, utilizing Eq. (A.3), we evaluated 453 the gradient of the modiﬁed version of the objective function@fPN @eused for updating the design 454 variables. 455 27 (a) (b) 03691215182124270510152025303540Location of /s119max and /s119min a long the wave vector k I teration number/s119max/s119 min 0481216202428320510152025303540Location of /s119max and /s119min a long the wave vector k I teration number/s119max/s119 min Figure 14: Variation of the location of the !maxand!minalong the wave vector k, for the optimzation of the band gap between 4thand5thbands using discrete formulation (Section 5), at (a) bias electric ﬁeld e3= 150 MV/m and (c) bias electric ﬁeld e3= 300 MV/m. The optimized layouts of the unit cell, the corresponding band diagrams, and objective function 456 history obtained using the aforementioned p-norm formulation, for the maximization of the band 457 gap between the 4thand5thbands are shown in Figs. 13a and 13b, when bias electric ﬁeld is equal 458 toe3= 150 MV/m and e3= 300 MV/m, respectively. From the optimized unit cell designs, we 459 observe that the designs are the same as those obtained using the original formulation (Figs. 8c and 460 9c) for both values of bias electric ﬁeld. In Fig. 14, we show that the location of !maxand!minalong 461 the wave vector kis almost constant same for all iterations of the non-smooth formulation. This 462 explains why the same designs are obtained using both formulations. However, from the right panel 463 of Figs. 13a and 13b, we observe that the optimization with p-norms takes 41 (for e3= 150 MV/m) 464 and 77 (fore3= 300 MV/m) iterations to converge, which consumes 14 (for e3= 150 MV/m) and 465 39 (fore3= 300 MV/m) more iterations in comparison to the non-smooth formulation (Figs. 8c and 466 9c). Further, the computational cost of optimization with the p-norm formulation is signiﬁcantly 467 higher in comparison to the non-smooth formulation, due to the evaluation of design sensitivities 468 or the derivatives of eigenvalues corresponding to different values of wave-vector kalong the ﬁrst 469 Brillouin zone. 470