Ultra-low and ultra-broad-band nonlinear acoustic metamaterials

Linear acoustic metamaterials (LAMs) are widely used to manipulate sound; however, it is challenging to obtain bandgaps with a generalized width (ratio of the bandgap width to its start frequency) >1 through linear mechanisms. Here we adopt both theoretical and experimental approaches to describe the nonlinear chaotic mechanism in both one-dimensional (1D) and two-dimensional (2D) nonlinear acoustic metamaterials (NAMs). This mechanism enables NAMs to reduce wave transmissions by as much as 20–40 dB in an ultra-low and ultra-broad band that consists of bandgaps and chaotic bands. With subwavelength cells, the generalized width reaches 21 in a 1D NAM and it goes up to 39 in a 2D NAM, which overcomes the bandwidth limit for wave suppression in current LAMs. This work enables further progress in elucidating the dynamics of NAMs and opens new avenues in double-ultra acoustic manipulation.


Supplementary note 1. Method for equivalent parameters
As shown in Fig. 2b, the flexural motion of the entire attachment in the meta-cell can be modeled as a torsional oscillating system. To determine the equivalent parameters J 0 , J r and k T , the finite element method (FEM) is used to calculate the low-frequency responses of the entire attachment, and the displacement of the terminal points are shown in Supple-

Supplementary note 2. Dispersion theory based on FEM FEM of the nonlinear oscillators
The finite elements of the nonlinear oscillators derive from their motion equations. In case N1, that is, neglecting the influence of the vibro-impact oscillator, the node vector is v=[w 0 , θ 0 , r w , θ r ] T , where r w denotes the transverse dis- If we take into consideration of the vibro-impact oscillator, as shown in Fig. 2c, the motion equation of the coupled torsion-vibro-impact system in the NAM beam is

FEM of a periodic cell of the NAM beam
By a standard FEM process, the motion equation of all the nodes in a cell is 3 , ] ,  f f 0 f 0 (16) in which M , K and N denote the mass, the stiffness and the nonlinear stiffness coefficient matrices, respectively; u denotes the displacement vector of all the nodes; f is the vector of the node forces; the subscripts, L, I, R and a, symbolize the nodes at the left boundary, inside the primary beam, at the right boundary and of the attachments, respectively;

FEM of a periodic cell of the NAM plate
For the plate, the displacement vector of every node is . With a proper transform and a standard FEM procedure, the motion equation of a cell reads where the displacement vector of the whole system is T T RT I  r   { } {  , , , , , , , , , }  q q q q q q q q q q q (24)

T T T T T T T T T T LB B RB L R LT
As shown in Supplementary Fig. 2, q j denotes the set of node displacements, for example, q L represents the displacement of all the nodes in the left boundary of a cell. q r =[ r w , θ rx , θ ry ] T or q r =[ r w , θ rx , θ ry , ϕ rx , ϕ ry ] T , which depends on the system N1 or N2. Correspondingly, the node vector of the force is in which I j denotes the identity matrix whose dimension is equal to the corresponding vector q j . With the relationships of interaction forces between neighbor cells, we can prove that 6 At last, the motion equation of a cell is transformed to (22) and (37) have a same form.

Methods to calculate the dispersion relations
The analytical dispersion solutions of the linearized system are obtained with the eigen function However, only approximate dispersion solutions can be found for nonlinear metamaterials.
For the weakly nonlinear metamaterial, perturbation approach (PA) is used to solve the dispersion relations, as elaborated in the Refs. [1,2]. The nonlinear dispersion solution derives from PA is where 0 U is the generalized eigenvector belonging to the eigenfrequency ω 0 of the linearized system, that is, In our manuscript, PA can only be used to in the case NAM-N1. For NAM-N2, the nonlinearity in vibro-impact oscillator is so strong that PA would present wrong results.
For the strongly nonlinear acoustic metamaterial, harmonic balance method (HBM) can present the dispersion relationships more accurately. Here the wave solution is assumed as ( )= sin tt  uU . Balancing the coefficient of sin t  , we obtain the first-order harmonic balance solution In both PA and HBM, by specifying the initial amplitude 01 U of u 1 , we can solve the dispersion frequency ω and other parameters in U.
In the FEM process, a cell of the primary beam is meshed by three conformal elements with same length. And a 2D NAM plate cell is meshed by 2×2 elements. These meshes are accurate enough for low-frequency studies (lower than 1500 Hz). For NAM beam-N2, analytical solutions can be found for the 11-dimensional (11D) system of algebraic equations in Eq. (40). However, for the 2D NAM plate-N2, the system of equations is 16D. It is a high-dimensional system so that it becomes difficult to find analytical solutions with HBM. Therefore only the dispersion solutions of the NAM 7 plate-N1 are considered and solved by PA.

Supplementary note 3. Nonlinear vibration of NAM beam
In this section, based on the nonlinear FEM, we take into account of four nonlinear sources in the NAM beams: the inertial nonlinearity, the geometrical nonlinearity, the Duffing oscillator and the vibro-impact oscillator. HBM and Harmonic average approach (HAA) are introduced. By studying the influences of different nonlinear sources, we confirm that the Duffing and the vibro-impact oscillators are the most important nonlinear factors. Moreover, the strong nonlinearity of the vibro-impact oscillator is analyzed with a comparative way. At last, we show the experimental results on a LAM beam.

Nonlinear FEM of a pure beam
The 2D theory for the geometrically exact beam element is based on the assumption of the plane across section. It was derived by Reissner and is valid for finite deflections, rotations and strains 3,4 . The associate strain-deflection relations are where ε, γ, κ are the axial strain, the shear strain and the curvature, respectively; u is the displacement in axial direction, w is the deflection and ψ is the rotation, as shown in Supplementary Fig. 3; (·) x =d/dx. The Eula-Bernoulli beam assumes γ=0. By using this shear constrain, ε is obtained depending upon u and w By some algebra to solve ψ and then using a derivation, one yields For the free-end beam, the axial deformation can be neglect, that is, ε=0. Using this constrain to Eq. (44) yields Therefore u x ≤0, which indicates the beam generates shorten inertial effect. Through an integral, we obtain the displacement induced by the shorten effect at point With the Taylor series Using that Taylor series, we obtain Based on these displacement and stain expressions, the kinetic energy T and strain energy U of an element within length s are where (·) x =d/dx and ρ is the mass of unit length. The second term in T denotes the inertial nonlinearity and the second term in U represents the geometrical nonlinearity. In this expression, the geometrical nonlinearity depends both on the transverse strain and the curvature.
At present, most studies on the nonlinear beams use the Galerkin method by assuming that the motion depends on the superposition of several linear modals. This method reduces much DoFs in analyzing for convenience. However, this 9 method is proper for simple cases. For the nonlinear metamaterial beam coupled with many nonlinear oscillators, modal superposition is improper for broadband analyses. Here, we adopt the nonlinear FEM. In general, we can use the expression of the linear FEM to approximate the displacement field in nonlinear FEM, that is, where () t η is the nodal vector, is the vector of shape functions. If we adopt the finite element coordinate x=lξ, and length of an element is 2l, the shape function of the conformal element is where   in which the subscript 'i' symbolizes the corresponding matrix of the ith element; ()'=d/dξ; the superscript 'T' denotes the matrix transposition. These formulas indicate the nonlinear elemental matrices consis t of quadratic terms and they are depends both upon the shape function and the nodal coordinates.

Nonlinear FEM of the entire NAM beam
With the elemental matrices of the coupled nonlinear oscillators and the nonlinear beam, the motion equation of the entire system can be formulated with a standard FEM procedure 5 . It is where M, K, C are the mass, the stiffness and the damping matrices of the system, respectively. The subscripts , L, I, G, r, denote the variables corresponding to the linear case, nonlinear inertia, geometry and oscillator. r () Kq can be decoupled, but I () Mq and G () Kq are coupled. 10 We adopt the dynamic proportional damping model, where c 0 is a constant damping coefficient. In the main text and in this document, c 0 =0.001 in the bifurcation analyses based on dimensional reduction method. In other method, the damping effect is neglect.

Harmonic balance method (HBM)
The wave propagation depends on the periodic solutions. A frequency response is also a periodic solution. If the stability of a solution is unnecessary, HBM can be used. A first-order approximation of the solution is ( ) sin tt   qQ . Subscribing it into (61) and balancing the harmonic, we obtain the solution for the cubic nonlinear system, When the nonlinear matrices contain higher-order terms, the approach is still valid except for the different coefficients.
However, there is larger round-off error for higher order nonlinear terms. In this case, we need to take into consideration of high-order harmonic waves, for example, supposing the solution is The Newton iteration numerical algorithm is used to find the solutions of the vector Q with a specified force F at point E (see Fig. 2f,g). In this process, each cell of the NAM beam is meshed into two conformal elements, which is accurate in low-frequency range. Therefore, there are 78 DoFs in total for the NAM beam-N1 and the NAM beam-N2 has 90 DoFs.
However, HBM doesn't provide a regime to analyze the stability of a solution, so it cannot be directly used for bifurcation studies, even for simple cases. To reach this problem, we adopt the harmonic average approach (HAA) 6 .

Harmonic average approach (HAA)
To analyze the stability of periodic solutions, we need to solve   1 to derive the Jacobian matrix of the system. However, it is difficult at present to solve a high-dimensional   1

Mq
because it contains many unknown variables. Therefore we neglect the inertial nonlinearity in this method and our analyses will address that this neglecting does not introduce large errors.
We consider the equation Comparing the expressions of q in Eq s. (65) and (67), we obtain ( cos sin ) The further substitution of (65) and (67) In practice, for the cubic nonlinear system, we can reduce the order of triangle functions in N k and rewrite it as ( , ,cos ,sin ) ( , )cos ( , )cos3 ( , )sin ( , )sin 3 Then, assuming that u and v are constant, by calculating (70) The steady solutions correspond to the condition  u v 0 and their expressions are The amplitude of the response is 22  Y u v . Generally, we cannot obtain an analytical solution of these algebraic equations. Newton method and continuation are used to find the solution.
In fact, HBM can obtain the same result as (74), but with the differential formula, we can analyze the stability of a solution. We rewrite (72) Then the Jacobian matrix of the system described by (76) is The stability of a solution is determined by the real parts of the eigenvalues of J: if there is an eigenvalue having a pos itive real part, the periodic solution is unstable; otherwise, it is stable.
However, though this stability theory of the periodic solution is rigorous in mathematics, for a high-dimensional sys-

Influences of nonlinear sources
The theoretical influences of different nonlinear sources on the NAM beam-N1 are illustrated in Supplementary Fig. 4.
The influences of the vibro-impact oscillator are explained in Fig. 3c and not repeated here. As illustrated in Supplemen- By comparing displacements at point A s A under the three cases in Supplementary Fig. 4b, we find that: in 0-20 Hz (it is in 0-LR1), the geometrical nonlinearity plays a main effect; however, in the frequency range 20 Hz-LR2, especially near LR1, the Duffing oscillator plays a main effect to significantly suppress the elastic waves; the nonlinear effect deriving from the Duffing oscillator is broadband but it becomes weaker for the frequencies over LR2, meanwhile, the nonlinear effect from geometrical nonlinearity becomes a little stronger than the Duffing oscillator. Therefore in the 13 broadband analyses, the geometrical and inertial nonlinearities can be neglected when we had considered the nonlinear oscillators.

Nonlinearity of the vibro-impact oscillator
In the theoretical description of the metamaterial beam and plate, the nonlinear restoring force of the vibro-impact os- Fig. 3, we analyses the nonlinear effect by taking n=3 under F=5 N as example in bifurcation analyses. In fact, we can obtain stronger nonlinear effect by increasing n. In Supplementary Fig. 5, we compare the broadband responses under F=2 N for two cases n=3 and n=5. Results indicate that the case n=5 obtains a broader effect that suppresses the elastic wave with a smaller excitation than the case n=3, especially in the third passband. Moreover, the responses near LR1 are approximately equal for the two cases. Therefore we can predict that, further increasing n will make the theoretical features in the frequency responses approach to the experimental results.

Experiments on the LAM and NAM beams
In this section, an experiment on a LAM beam is introduced firstly. We remove all the magnets in the NAM cell in Fig.   2a to fabricate a cell of the LAM. There is a strut and a bolt left, which can be equivalent to a torsional local resonator that produces a LR bandgap. Influences of moderate deformation on the broadband frequency responses are measured.
For this LAM beam, under moderate deformations, there would be inertial and geometrical nonlinearities in the primary beam, but there is not nonlinearity came from the attached oscillators. As illustrated in Supplementary Fig. 6

Supplementary note 4. Dimension-reduction method
Bifurcation properties can help to directly demonstrate the chaotic mechanism in theory. For a low-dimensional model, 16 HAA can present bifurcations of periodic solutions accurately 6 . However, for the high-dimensional model of the NAM beam, HAA would give wrong stability information because of the numerical errors coming from the periodic solutions and the eigenvalues of the Jacobian matrix (eigenvalues are sensitive to periodic solutions). To analyze the bifurcations of the periodic solutions of the NAM beam, we need practical dimension-reduction algorithm 7 in frequency domain to reduce the dimension of the finite element model.
Here we adopt the thought of the post-processed Galerkin method 8 in frequency domain to reduce the dimension of the NAM beam and then we analyze the periodic solutions and their bifurcations based upon the reduced system.

Dimension-reduction algorithm
The dimension d of a system is the number of the second-order differential equations (64). We adopt a rigorous and practical dimension-reduction approach based on the orthogonality of linear modals.
The eigenvalues of the matrix B  , k k k k  C P CP K P K P .
The first formula in (85) is the reduced system. Then, the problem comes to be how to solve the system. The traditional 17 Galerkin approach neglect the influences of high-order modals by specifying y=0, so the reduced system is The stabilities of the periodic solutions of the original system are same with that of the reduced system.
Step 3: With the solved x, we post-process that 11 ( ) ( cos sin ) tt       y x a b and adopt the first-order harmonic approximation.
Step 4: The periodic solution of the original system is solved with Eq. (83).

Bifurcations in the NAM beam-N1
To demonstrate the validity of the dimension-reduction algorithm, we reduce the dimension of the finite element model of NAM beam-N1 from 78 to 20 to study the bifurcations of periodic solutions , as illustrated in Supplementary Fig. 8.
Geometrical nonlinearity is considered here. In the dimension-reduction procedure, we consider the first 20 linear eigenmodes to approximate the responses in 0-120 Hz. In the next section, we will introduce how to expand this frequency range. Moreover, perturbation continuation methom 6 is used to find different branches of solutions. Considering the numerical errors, we used a tiny eigenvalue μ c =1×10 -8 of the Jacobian matrix as the critical value to determine the stability of a solution. In fact, for unstable solutions, the positive real parts of the eigenvalues of Jacobian matrix can be very large; in contrast, they are negative or just in the magnitude 10 -12 for stable solutions. In Fig. 3, we also use this threshold value for NAM beam-N2.
As suggested in Supplementary Fig. 8, near 16 Hz, because the existence of cubic geometrical nonlinearity, the nonlinear resonance consists of 3 branches bending right, among which two branches are stable and another branch is unstable, which is like the nonlinear resonance of a 1DoF nonlinear system. This nonlinear resonance demonstrates the accuracy of the dimension-reduction algorithm. In contrast, near 60 Hz and 85 Hz, nonlinear resonances consist of unstable branches or unstable peaks whose amplitudes are much lower than the linearized responses. Those unstable branches provide high possibilities of chaotic responses, as well established in Ref. [6].