Multi-resonator metamaterials as multi-band metastructures

Abstract Introducing multi-resonator microstructure into phononic metamaterials provides more flexibility in bandgap manipulation. In this work, 3D-acoustic metamaterials of the body- and face-centered cubic lattice systems encompassing nodal isotropic multivibrators are investigated. Our main results are: (1) the number of bandgaps equals the number, n, of internal masses as each bandgap is a result of the classical analog of the quantum level-repulsion mechanism between internal and external oscillations, and (2) the upper boundary frequencies, ωupper2i, i = 1, 2, ⋯, n, of the gaps formed coincide with eigen-frequencies, ωint;i2 ≠ 0, of the isolated multivibrator, ωupper2;i = ωint; i2, and the lower boundary frequencies, ωlower2,i2, are in good agreement with estimations as ω lower , i 2 ≈ ω ^ int ; i 2 ( ω lower , i 2 ω ^ int ; i 2 ), where ω ^ int ; i 2 represent the eigen-frequencies of the multivibrator when its external shell is motionless. The morphologies of the set of dispersion surfaces, ωm2(k), m = 1, 2, …, 6, in the corresponding passbands are similar to each other and to that of the set of dispersion surfaces, ωext; m2(k), obtained through the exclusion of internal masses. Thus, the problem of analyzing the acoustic properties of the complicated system is reduced to the study of two simple sets {ωint; i2} and ω ^ int ; i 2 , along with {ωext; m2(k)}, the morphology of which depends only on the type of lattice symmetry. This splitting renders controlled phononic bandgaps formation in homogeneous multi-resonator metamaterials feasible.


Introduction
Phononic metamaterials, as synthesised latticed materials possessing microstructure, have been the subject of contemporary research due to their beneficial controllable features such as concurrent negative effective stiffness and density [1] and controllable frequency pass-and stop-bands [2]. These features, not observed in natural materials, render them useable in a variety of applications, including blast wave mitigation [3], seismic wave protection [4], enhanced flexural wave sensing [5], and several other applications where wave manipulation or guiding are relevant [6]. The microstructure in phononic metamaterials may comprise a single mass and spring (single resonator) or multiple masses and springs (multiple resonator) or such idealised realisations of other elements when feasible.
Introducing multiple resonators into the microstructure of phononic metamaterials adds new dimensions to bandgap formation, frequency spectral expansion, and passband expansion/contraction as a result of the additional degrees-of-freedom introduced [7,8]. Compared to their single resonator (SR) counterparts ( Fig. 1(a)), multi-resonator (MR) metamaterials ( Fig. 1(b)) possess more degrees-of-freedom per unit cell and, as a consequence, a higher number of dispersion surfaces and multiple bandgaps. SR metamaterials, such as those introduced by Lazarov and Janson [9] or Huang et al. [1], possess only a limited number of (in their case, two) complete bandgaps when acting as phononic filters rendering them unsuitable for some devices in need of multiple gaps. The solution is to increase the number of degrees-of-freedom in the microstructure through the introduction of additional elements. While several related ideas utilising more sophisticated microstructures such as fractals/pseudo-fractals [10], Labyrinthine acoustic metamaterials with space-coiling channels [11], and Hilbert curves [12] are proposed recently, extending the classic SR metamaterials using merely additional masses and springs proves to introduce distinctive features while the system's level of complexity remains undaunting.
The inclusion of microstructure also assists with narrowing down or extending the gap size through correct tuning of model parameters [2,13,14]. This is true of SR as well as MR metamaterials, however, MR metamaterials are more flexible due to an increased number of model parameters involved so the size of bandgaps may be appropriately tuned. It is therefore essential to investigate the primary features of MR metamaterials. While SR metamaterials of different topology and at several length scales are well studied (see e.g. [2,13,[15][16][17][18]) there are a limited number of studies conducted on the characteristics of MR metamaterials (see e.g. [7,8,[19][20][21]).
In the context of MR metamaterials, several forms are conceivable. While Huang and Sun [7] focused on the simple mass-spring-based design, Bonanomi et al. [22] discussed wave propagation in granular chains of directional filtering properties, Li [23] and Li et al. [24] came up with the design for acoustic diodes based on the hybrid tunable sonic crystal-multiphase gel units, Martinez-Sala et al. [25] investigated diffraction structures at macro length scales, and Gorishnyy et al. [26] proposed hypersonic phononic crystals. These studies are very interesting and form the basis of novel designs in the field [6], however, as mentioned before, we concentrate here on simple mass-spring models of MR with distinctive features.
Besides theoretical modelling, experimental research on the topic of design and testing phononic metamaterials has been ongoing ever since the conception of the idea, and some recent studies have focused on the topic. To mention but a few, Liu et al. [27] fabricated and tested sonic crystals encompassing localized resonant units and depicted spectral gaps with a lattice constant two orders of magnitude smaller than the relevant wavelength. They constructed and tested a 2 cm slab of the metamaterial which was shown to break the conventional massdensity law of sound transmission by one or more orders of magnitude at almost half a kilohertz. In another study conducted by Sheng et al. [28], a rigorous theoretical derivation of dynamic mass density yielded frequency-dependent values significantly deviating from the longwavelength (static) value. The dynamic mass density was shown to be negative for certain frequency ranges (bandgaps) and the results were corroborated with experimental data. In a more recent study, Wang et al. [29] reported on the mechanism for simultaneous realization of perfect absorption and broadband insulation by layered acoustic metamaterials. They showed that at a particular frequency (312 Hz) nearly perfect absorption (98.4%) was achieved by testing a sample of 15 mm thickness. Barnhart et al. [19] conducted an experimental program to study dissipative MR metamaterials used for broadband elastic wave attenuation. This broadband wave attenuation was demonstrated through an impact test performed on finite samples where the frequency spectrum of the transmitted amplitude was found to be in close agreement with numerical results. Finally, in a rather recent study, Ma and Sheng [30] provided an overview of the topic as investigated over the one and half decades prior to their work. Most results in this realm hinge upon the simple fact that composites with locally resonant microstructure within their unit cells exhibit simultaneous negative effective density and stiffness at certain frequency ranges. While the experimental study of this class of metamaterials is interesting in itself and essential to adduce evidence in support of the theoretical results obtained also to depict the physical significance of parameters defining the models investigated, such a study falls beyond the scope of the present paper as it merits a great deal of research on its own. This paper is organised in 5 sections as follows: following this introduction, the discrete-parameter multi-resonator metamaterial is formulated in Section 2. Each MR unit comprises mass points of zero rotatory inertia and is connected to adjacent masses by massless linear springs of equal stiffness in each direction (local isotropy). The two kinds of connectivity of MR units considered are those of the BCC and FCC structures. In Section 3, the BCC MR metastructure is analyzed and its results expounded. In Section 4, the FCC lattice is considered. In each case, dispersion surfaces are tuned to form bandgaps and passbands. Finally, the study is concluded in Section 5.

The discrete-parameter multi-resonator (MR) metamaterial
The discrete-parameter MR metamaterial is an idealisation of scenarios where mass and springs stand for structural elements, as in Fig. 2(a). When MR units are embedded in a continuum, as if connected by external springs in a given topology, the MR metamaterial is ensued (See Fig. 2(b)). When feasible assumptions are made regarding the distribution of mass and stiffness in the system, the discrete-parameter model is resulted in and the equations of motion may be cast in a simple form. Considering a three-dimensional metamaterial comprising an infinite arrangement of locally resonant mass-spring MR's connected in a given topology, as shown schematically in Fig. 2(b), the equations of motion for different layers of the mass-spring model ( Fig. 2(a)) are derived as follows: where the layers are numbered from inside outwards; m p is the mass of layer p; and k p q is the stiffness of the spring connecting layer p to layer p + 1 in the direction of q, where q designates a generalized coordinate which may be the x, y, or z components of the displacement vector in the cartesian coordinate system. Equations of motion for the outermost layer depend upon the external topology of MR units in the metamaterial and may be complicated, in general. In this work, two different well-known topologies of the BCC and FCC lattices are considered which render local interactions of each MR unit with immediately adjacent ones possible. In the sequel, the equations of motion are derived for both topologies and subsequent analyses performed which make it possible to derive the dispersion surfaces and bandgap structures for each case, corresponding to acoustic and optic modes of vibration ( Fig. 2(c)).

Analysis of the MR metamaterial
As mentioned before, two topologies, namely the BCC and FCC, are considered here to represent the external connectivity of the MR units. The BCC and FCC topologies are selected as they replicate two of the topologies encountered in nature, for instance in crystalline metals. Fig. 3 shows both of these topologies in 3D along with the nodal microstructure for a four-layered resonator, as an example.

Formulation of the problem
The dispersion equation for the BCC structure under consideration is based on two basic relations. If two locally resonant elements of the system are connected by a spring of stiffness κ, then the second element acts on the first one with a force, f 12 of Eq. (3).
Here x 1 and x 2 are the displacements of the spring ends from the equilibrium position, e 12the unit vector directed along the spring, and the displacements are assumed to be small i.e. |x 2 − x 1 | ≪ l, where l is the unstretched length of the spring.
When a wave propagates in an infinite medium, the displacements, x 1 , x 2 , in the equivalent nodes of the structure with coordinates, r 1 , r 2 , are related by Bloch's equation as follows: where k is the wave vector, and the vector r 2 − r 1 is a translational symmetry vector. The unit cell of the BCC structure includes two lattice nodes (Fig. 4) and the corresponding lattice constant is equal to a. Blue and red nodes form simple cubic lattices which are then inserted into each other to create the model. Nodes A and B (marked with circles) are used as 'reference nodes'. As such, the displacement of any other structural node is associated with the displacement of the corresponding reference node through Eq. (4). For the selected shape of the unit cell (instead of a truncated octahedron, which contains one node only), the Brillouin zone in k-space is a cube of side length 2π/a.
As an infinite lattice with this unit cell, the dynamics of the entire system can be described by equations for the displacements of the reference nodes represented by a column vector, X, composed of 18 entries as in Eq. (5).
In relations (5), X ext represents the displacement vectors of the outer shells for the reference nodes A and B of mass M (see Fig. 4), X int are the displacements of the inner shells of mass m i , and X nuc are the displacements of the inner cores in nodes A and B (the mass of the core is m n ). As mentioned before, nodes A and B have an identical internal structure (local homogeneity) (the stiffness of internal springs is indicated in Fig. 4). The stiffness of external springs connecting the nodes of the same sublattice is χ (in Fig. 4(a), for example, these are the springs between the neighboring blue nodes). The stiffness of springs that connect  adjacent nodes belonging to different sublattices is derived as ϰ ¼ 2χ= ffiffiffi 3 p , provided that the stiffness of each spring is inversely proportional to its length, i.e.: Now, we derive the equations of motion for only the outer shells of reference nodes, i.e. without taking into account their internal microstructure. To this end, we take into consideration all the forces, caused by 6 neighboring blue nodes (see Fig. 4(a)) and 8 red nodes, acting on the reference node A. Similar calculations are performed for reference node B. The basis of these calculations is presented in Eqs. (3) and (4). For example, when calculating the forces exerted by 8 neighboring red nodes and acting on the reference node A, it should be taken into account that their displacements are related by the equation: Here e BA is the vector directed from the reference node B to the reference node A, e j is the unit vector directed from the reference node A to one of the neighboring red nodes (j = 1, 2, 3, …, 8), and X B (red) is the displacement of the reference node B. As a result, we obtain the following system of equations: The third-order matrix, A, describes the interaction of nodes in each of the blue and red sublattices while the matrices B and B * (complex conjugate of B) represent the mutual influence of the sublattices on one another.
Hereinafter, by the wave vector, k, we mean the dimensionless wave vector, i.e. k j ja ¼ 2πa λ ¼ k, where λ is the wavelength, I n is the identity matrix of the n-th order, and the following relations hold: where α i are the functions of wave vector components and related to them as follows: When we take into account the fact that the lattice nodes have an internal microstructure, we obtain a system of equations as: In the sequel, for a qualitative analysis of acoustic properties of the system under consideration, it is useful to know the eigenfrequencies of the isolated lattice nodes i.e. without taking into consideration the influence of external springs. These can be found from the truncated system of Eq. (15) if we single out the equations of motion for the displacements of the outer shell of the node, x e , its inner shell, x i , and its core, x n , at χ = 0 as in Eq. (16).
In the model chosen, it suffices to consider the one-dimensional problem as the excited internal oscillations are independent (due to the geometric linearity approximation) along arbitrarily defined perpendicular directions. For the periodic solutions, i.e. x = x 0 exp (iωt), there are two non-zero eigenfrequencies (the third solution of Eq. (15), ω 2 = 0 corresponds to the free movement of the node as a whole): where.
The frequency ω 0 2 ≡ ω i 2 represents the eigenfrequency of an isolated node in the absence of a core, while the frequency ω n 2 corresponds to the natural oscillations of the isolated inner shell-core pair. The equation for the eigenfrequencies and eigenvectors characterizing the acoustic properties of the entire system is obtained based on Eq. (15), represented in the form X = X 0 exp (iωt) thus: In Eq. (21), the frequency ω is given in units of ω 0 where: Since the matrix of the system of equations in Eq. (21) is a function of the wave vector k, Eq. (20) defines, in the general case, 18 dispersionsurfaces ω(k) in the four-dimensional (4D) space. Four parameters, α, β, γ, and Ω 2 -which determine the acoustic properties of the systemhave a transparent physical sense (including Ω 2 -see the next section) and are easily controlled.

Acoustic properties of the system
The presence of multivibrators inside the lattice nodes leads to the formation of both separate passband regions and bandgaps in the frequency spectrum. Below, we show that the corresponding boundary frequencies can be established as the simple mathematical relations of a clear qualitative interpretation.
We demonstrate the indicated acoustic properties of the system using the results obtained for a specific set of parameters (Figs. 5 and 8: α = 2, β = 1, γ = 2, and Ω 2 = 1.5). Note that the parameter Ω 2 (see Eq. (24) and Fig. 5) corresponds to the frequency of lattice vibrations at k = 0 in the case when internal multivibrators are removed from the system. Such (thrice degenerate) vibrations (according to the shape of the Brillouin zone for the selected unit cell) are equivalent to vibrations with the wave vector k x, y, z = 2π and wavelength λ = 1 (in the units of the lattice constant a). In other words, each of the sublattices of the system is shifted, as a whole, relative to the other sublattice along one of the coordinate axes x, y, or z. In this case, the neighboring planes of lattice sites with an orientation of the type [100] are displaced in the antiphase without deformation of the springs in the sublattices.
Oscillations at the boundary point k(π, π,0) in Fig. 5 (point M 1, 2 in Fig. 12) also correspond to antiphase displacements of neighboring planes of nodes ([110] -type) with frequencies: The first two frequencies correspond to transverse vibrations, and the third to longitudinal. For future estimates, we note that ω M, 3 2 is the maximum frequency in the first Brillouin zone. This follows from the truncated system of Eq. (21), if we exclude from consideration internal multivibrators: For individual values of the wave vector, k, the dispersion equation (Eq. (26)) can be easily solved analytically. At the point k(π, π, π) (point R in Fig. 12), the longitudinal and transverse vibrations have the same frequency ω 2 Fig. 12 The planes outlined in Fig. 5(a) with red and blue contours are corresponding to surfaces: where ω int; 1, 2 2 are the eigenfrequencies of isolated nodes (see Eq. (15)). When defined in terms of dimensionless parameters these are recast as: For the parameters given above ω int; 1 2 = 2/3, ω int; 2 2 = 4/3. The result of the interaction between external and internal elastic forces is shown in Fig. 6(a). The brief qualitative interpretation of this is quite transparent. The internal structure of lattice nodes defines the isotropy of their elastic properties. In this case, a coupling strength between external and internal vibrational modes manifests itself (see Eq. 28) and surfaces ω 2 = ω 2 (k) (see Fig. 5(a)) these modes have the same frequency and the same polarization, so that effects analogous to the quantum effect of energy-level repulsion arise (see [27]): the surfaces ω 2 = ω 2 (k) are essentially 'torn up' by planes ω 2 = ω int; 1, 2 2 . The reconnection of these 'primary' dispersion surfaces is followed by the formation of bandgaps (see Fig. 6(a)). On the adjacent newly formed upper frequency surfaces (ω 2 (k) is higher than the corresponding value ω int; i 2 ), the external elastic forces act on the outer shells antiphase with internal forces. On the lower adjacent surfaces (ω 2 (k) < ω int; i 2 ) these forces act in phase. Thus, the maximal frequency on these surfaces can be achieved when the result of the 'in-phase displacements' vanishes, in other words, the outer shells are stationary. The analogs of the bandgap/passband formation mechanism discussed manifest themselves in other physical systems. The passband for electromagnetic (EM) waves in metals corresponds to the frequencies, ω EM , that is higher than the eigenfrequency of plasma oscillations, ω plasma , in the system of electrons and positive charged nuclei. More complex examples of the formation of passbands/bandgaps were realized by the interaction of EM waves with saturated atomic cesium vapor as shown in previous studies [31].
It is important to establish the quantitative characteristics of the passband regions (alternatively bandgaps) as a function of parameters defining the acoustic system. To this end, we consider a simpler problem when the internal structure of the node is represented by a simple vibrator (see the inset in Fig. 8). Examples of its solution, by analogy with the results of Figs. 5 and 6 are shown in Fig. 7. We draw attention to the fact that the topology of dependencies ω 2 (k xy ) in each passband region, shown in Fig. 6 and Fig. 7(c) and (d), is similar to the structure of the frequency dependencies of the system in the "initial" (without internal structure) state -see Fig. 5 and Fig. 7(a) and (b). Let us consider, in detail, factors determining the result of splitting of this initial state into separate zones.
The interaction between external and internal elastic media does not depend on the nature of the external elastic forces that cause vibrations of the node shell (see the inset in Fig. 8) with a certain angular frequency, ffiffiffiffiffiffiffiffiffiffi ffi χ=M p , at m i = 0. Taking into account the internal structure of the node leads to the formation of hybrid frequencies, which are easy to determine from the system of equations of motion for the displacements of the outer shell, x e , and internal mass, x i . The corresponding dimensionless hybrid frequencies are equal to: where  Fig. 8 is equal to 1). The main result of the solution of the problem is that the lower hybrid frequency has an upper limit (see Fig. 8): lim The limit value 1/(1 + γ) can be obtained from qualitative considerations that lead to a simple method for calculating hybrid frequencies for a multi-level internal structure of nodes.
In the case of free vibrations of an isolated node, the displacements of the shell and the core are out of phase and the displacement of the system's centre of mass is identically zero -Mx e (t) + mx i (t) = 0. The influence of an external field changes the frequency of internal oscillations and shifts the center of mass of the node with time. However, these forced vibrations can be represented as free vibrations of an isolated node with the effective mass of the outer shell M * = − mx i (t)/ x e (t), which ensures the immobility of the center of mass of the node. At the lower hybrid frequency, the displacements x i (t) and x e (t) are in phase and M * < 0. On the other hand, the eigenfrequency of oscillations of the node is ω 2 As the stiffness of the external spring increases, the displacement value x e (t) → 0, and the effective mass M * → − ∞, thereby bringing the value of ω h 2 to the upper limit k i m . Its dimensionless quantity, k i m =ω 2 0 , is equal to the value of 1/(1 + γ) obtained from Eq. (31), and represents the eigenfrequency of internal vibrations under the assumption that the outer shell of the node is motionless.
The result obtained can be generalized to the case of a multilevel metastructure of n−1 internal multivibrators possessing n−1 internal elements/masses (Fig. 2(a)). The total number, n, of elements of an isolated node determine its n eigenfrequencies as: ω 2 int;0 ¼ 0; ω 2 int;1 ; ω 2 int;2 ; …::; ω 2 int;n−1 , which subsequently represent the lower boundaries of the passband regions. The upper boundaries of these regions (or the lower boundaries of the bandgaps) are guaranteed to represent eigenfrequencies: ω 2 h,1 , ω 2 h,2 , ω 2 h,3 , . . . ::, ω 2 h,n−1 of an isolated node under the assumption that its outer shell is motionless. Bandgaps thus correspond to frequency intervals: Data presented in Figs. 5 and 6 confirm the assumption above about the general principles anent the formation of a multiband spectrum of frequencies. In the case under consideration in this section n ¼ 3. The eigenfrequencies of a two-level internal vibrator are defined in Eq. (26). To determine the operator that transforms the dependencies ω 2 = ω 2 (k), shown in Fig. 5, into a series of three passbands, it suffices to consider the one-dimensional oscillations of a single node (taking into account the internal microstructure) based on Eq. (15) taking into consideration external forces, as was done in deriving Eqs. (27) and (28). For the dimensionless parameters we have adopted, the corresponding equations assume the following form: Hence For parameters in Fig. 6, the eigenvalues are ω h, 1 2 = 0.195 and ω h,2 2 = 1.138. Such estimates somewhat narrow the width of the bandgaps and expand the middle passband. The widths of these zones can be determined more accurately if we substitute in Eq. (33) the highest frequency reached in Fig. 5 (see Eq. (25)).
The obtained values of ω h, 1 2 = 0.156 and ω h, 2 2 = 1.087 coincide with the result shown in Fig. 6(a). Data presented in Fig. 7 show the dependencies of ω i 2 (k) only along one straight line in the k-space (i = 1, 2, 3, …, 18 in the increasing order of frequencies). Graphical representation of dispersion hypersurfaces in four-dimensional space (k,ω) is impossible, but it is realistic to depict the alternation of bandgaps and passbands as follows.
Based on Eqs. (28) and (35), it is easy to control the passbands and bandgaps by approximating the surfaces ω (a) 2 , with spheres of radii ω 1, h 2 (k = 0), ω int; 1 2 , ω 2, h 2 , and ω int; 2 2 . Data represented by the surface in Fig. 10 depict the characteristic dependence of frequencies ω int; 1 2 and ω int; 2 2 on the parameters α and β for a fixed value of γ. It is discussed below how the upper bandgap can be adjusted based on these data (see in Fig. 6 the orange-coloured zone referred to in the sequel as the 'orange-zone').
A narrow 'orange-zone' can be formed if the frequency difference 2 ) is small relative to the frequency average value 〈ω 21 2 〉 = (ω int; 2 2 + ω int; 1 2 )/2. According to Fig. 10, one of the suitable parameter sets is: α = 7 and β = 0.5. The corresponding result is shown in Fig. 11(a). The dependencies of frequency on the wave vector are calculated for the trajectory in k-space, presented in Fig. 12. The 'orange-zone' can be significantly expanded (see Fig. 11(b)) for a large gap of Δω 21 2 (α = 0.5 and β = 0.5 -see Fig. 10).
In Fig. 11(a) and (b), the frequency Ω 2 exceeds the value ω int; 2 2 . By an inverse ratio, Ω 2 ≲ ω int; 2 2 , the difference ω 2, h 2 − ω int; 1 2 , which is an estimation of the width of the second passband, can be significantly reduced (recall that all three hybrid frequencies depend on Ω 2 and are easily calculated based on Eq. (33)). Results in Fig. 11(c) show the possibility of forming such a rather narrow passband regions (see the olive passband) with low group velocities of acoustic waves. The results obtained for the BCC lattice are directly related to more complex 3D structures as shown in the next section.

The FCC-multi-resonator mass-in-mass structure
Here we demonstrate the implementation of the physical principles of bandgap formation as discussed above, using the example of an FCC substratum lattice for the MR metamaterial. We construct the mathematical description of the acoustic properties of such a system in a slightly different way than in the previous section.
When considering the BCC structure, we operated with a unit cell in the form of a cube containing two nodes. Accordingly, the entire system was presented in the form of two simple cubic sublattices. This node partitioning made it possible to map the oscillation frequency of the system, Ω 2 , with wave vectors k x , k y , k z = 2π/a to the center of the Brillouin zone (k = 0) and subsequently display hybrid frequencies at the same point (k = 0) for interpretations and estimates of the width of bandgaps. A similar approach applied to the FCC problem (the unit cell contains 4 nodes) will force one to operate with four sublattices, Fig. 9. Representation of bandgaps in the three-dimensional (θ, ϕ, ω 2 )-space. Coordinates θ and ϕ indicate the direction along which the wave propagates (α = 2, β = 1, γ = 2, and Ω 2 = 1.5). which greatly complicates the construction of the dispersion equation. In the presence of internal multivibrators, the formation of bandgaps and passbands in between them obeys the same laws as in the case of the BCC lattice. The equation of motion for the coordinates of the elements of the system below.  Fig. 11. have the form: from which the equation for determining the eigenfrequencies and eigenvectors, X 0 follows: •X 0 ð41Þ Here η ¼ 1 8þ2 ffiffi 2 p , and parameters, α, β,and γ are defined, as before, by Eq. (24). The fourth parameter is the dimensionless eigenvalue defined by the ratio Ω 2 = ω L 2 /ω 0 2 (the definition of ω 0 2 is given in Eq. (18)).
In Fig. 14(a) the frequency characteristics of the system are depicted along the path ΓXWKΓLKWLUX in the k-space (see Fig. 13(b)). The results obtained from Eq. (38) when internal multivibrators are not taken into account and for the value of the parameter Ω 2 = 1.5 are shown in Fig. 14(a). The effect of turning the multivibrators on is shown in Fig. 14(b) with conspicuous gaps in between dispersion curves (projections of dispersion surfaces). It can be seen that the topology of the dependences ω 2 (k) in the three passbands is an analog of the topology of the set of dependencies obtained in the "initial state" -see Fig. 14(a).
The revealed patterns of formation of passbands are applicable for more complex multivibrators (see Fig. 15). Here we present the result  Author statement Vyacheslav N. Gorshkov: Conceptualization, Formal analysis, Supervision, Writing-original draft, Writing-review and editing. Pooya Sareh: Investigation, Visualization, writing-original draft. Navid Navadeh: Methodology, Visualization. Vladimir V. Tereshchuk: Investigation, Software. Arash S. Fallah: Conceptualization, Investigation, Project administration, Supervision, Writing-original draft, Writing-review and editing.