*Article* **Band Tunability of Coupled Elastic Waves along Thickness in Laminated Anisotropic Piezoelectric Phononic Crystals**

#### **Qiangqiang Li 1, Yongqiang Guo 2,\*, Yajun Wang <sup>2</sup> and Haibo Zhang <sup>1</sup>**


Received: 22 May 2019; Accepted: 13 August 2019; Published: 16 August 2019

**Abstract:** Although the passively adjusting and actively tuning of pure longitudinal (primary (P-)) and pure transverse (secondary or shear (S-)) waves band structures in periodically laminated piezoelectric composites have been studied, the actively tuning of coupled elastic waves (such as P-SV, P-SH, SV-SH, and P-SV-SH waves), particularly as the coupling of wave modes is attributed to the material anisotropy, in these phononic crystals remains an untouched topic. This paper presents the analytical matrix method for solving the dispersion characteristics of coupled elastic waves along the thickness direction in periodically multilayered piezoelectric composites consisting of arbitrarily anisotropic materials and applied by four kinds of electrical boundaries. By switching among these four electrical boundaries—the electric-open, the external capacitance, the electric-short, and the external feedback control—and by altering the capacitance/gain coefficient in cases of the external capacitance/feedback-voltage boundaries, the tunability of the band properties of the coupled elastic waves along layering thickness in the concerned phononic multilayered crystals are investigated. First, the state space formalism is introduced to describe the three-dimensional elastodynamics of arbitrarily anisotropic elastic and piezoelectric layers. Second, based on the traveling wave solutions to the state vectors of all constituent layers in the unit cell, the transfer matrix method is used to derive the dispersion equation of characteristic coupled elastic waves in the whole periodically laminated anisotropic piezoelectric composites. Finally, the numerical examples are provided to demonstrate the dispersion properties of the coupled elastic waves, with their dependence on the anisotropy of piezoelectric constituent layers being emphasized. The influences of the electrical boundaries and the electrode thickness on the band structures of various kinds of coupled elastic waves are also studied through numerical examples. One main finding is that the frequencies corresponding to *qH* = *n*π (with *qH* the dimensionless characteristic wavenumber) are not always the demarcation between pass-bands and stop-bands for coupled elastic waves, although they are definitely the demarcation for pure P- and S-waves. The other main finding is that the coupled elastic waves are more sensitive to, if they are affected by, the electrical boundaries than the pure P- and S-wave modes, so that higher tunability efficiency should be achieved if coupled elastic waves instead of pure waves are exploited.

**Keywords:** coupled elastic waves; laminated piezoelectric phononic crystals; arbitrarily anisotropic materials; band tunability; electrical boundaries; dispersion curves

#### **1. Introduction**

Periodically multilayered composite structures [1–3] are constituted by periodically arranged unit cell with multilayered configuration. The different constituent layers in the unit cell have disparate material parameters such as material density and elastic stiffness constants. When elastic waves propagate in the periodically multilayered composites, their abovementioned particular construction form leads to the frequency band property resulting from the periodic scattering and further interference phenomenon of partial waves due to impedance mismatch at the interfaces between alternating constituent layers. This property, caused by the well-known Bragg scattering mechanism [2,3], refers to the elastic wave with frequency in pass-bands that propagates without attenuation, and the elastic wave with frequency in stop-bands that attenuates without propagation. During the recent three decades, the Bragg bands in periodically multilayered composites have been subsequently extended by the innovative concepts such as superlattices [4,5], phononic crystals [2,3,6,7], and metamaterials [2,6,8–10]. Besides the frequency bands, other novel elastic wave phenomena in these periodically multilayers such as the negative refraction, beam steering, and mode switching have also been revealed. On account of this progress, elastic wave (also referred to Floquet/Bloch waves) propagation in periodically multilayered composites (also known as laminated phononic crystals) becomes an even more attractive research topic.

Among the extensive studies on elastic waves in laminated phononic crystals, the very earliest object of study is the periodically multilayered composites made of only elastic materials. The incipient motivation for studying elastic wave propagation in periodically elastic multilayers is to develop nondestructive evaluation strategy for composite materials [2]. To achieve this goal, various elastic waves in the corresponding periodically elastic multilayers have been investigated, i.e., the bulk waves including the pure longitudinal (primary (P-)), pure transverse (secondary or shear (S-)), and their coupling modes in infinitely periodic multilayered elastic media [10–15], the surface/interface waves including the Rayleigh, Love, and surface transverse waves in semi-infinitely periodic multilayered elastic half-spaces [15–17], and the guided waves such as the Lamb and SH-type guided wave in periodically laminated elastic slab of finite thickness [10,15]. Besides the frequency bands of these wave modes, particular attentions have also been paid on the dispersion properties of various waves influenced by material anisotropy [12,13] as well as material elastic constants [13] and on the anisotropy of the characteristic waves like shear horizontal (abbr. SH) wave [14]. With the deepening understanding of the Floquet/Bloch waves in periodically multilayered elastic composites and the pressing demand for acoustic wave devices, researchers realized as early as 1980s that the piezoelectric materials can be introduced periodically into multilayers to form periodically laminated piezoelectric composites [4,5,18], and developed high-performed filters, guiders, and splitters, or delicately control mechanical waves through electricity. Since the piezoelectric material couples the mechanical field with the electrical field that is relatively easier to excite, detect and control as compared to other physical fields, the piezoelectric material is the most important and extensively-used intelligent material, especially in many functional devices like the bulk acoustic wave (BAW) and surface acoustic wave (SAW) devices [19]. Therefore, the combination of the piezoelectric effect and the elastic wave bands in the laminated piezoelectric phononic crystals is a very natural advance to improve the performance of and to add new function for acoustic wave devices [18]. It is also a promising strategy to control mechanical waves in a feasible and relatively easy way.

To push forward the development and design of these wave devices and the control strategy for various elastic waves, multifarious laminated piezoelectric periodic structures (phononic crystals) have been conceived. Moreover, the properties of the corresponding elastic waves thereof have been studied [4,5,18,20–71]. Nevertheless, the studies on elastic waves in laminated piezoelectric phononic crystals may be difficult on four counts. The first aspect, which appears definitely, is that the coupling between mechanical and electrical fields exists in piezoelectric phononic crystals, which represents the most essentially different property when compared to the purely elastic periodic composites. This electromechanical coupling causes the dependence of some or all the mechanical waves on the electrical field, besides results in the addition of a new wave mode (electric potential wave) mainly dominated by the electric field. Accordingly, these waves are sensitive to the electrical boundary conditions, which should be specified during the study of the wave properties and can be sorted into two main classes. One class is called as the passive-type electrical boundary, which remains fixed after applied like either

the electric-open or the electric-short one as the most common boundary condition. The other class is referred to as the active-type electrical boundary, which is able to be switched between diverse passive electrical boundaries or adjusted through connection to external electric circuits. The second aspect, which may possibly happen, is that the constituent materials in multilayered piezoelectric phononic crystals are anisotropic with enough low crystal symmetry to bring about the coupling among the original pure horizontally-transverse (shear horizontal (SH)), pure longitudinal (P), and pure vertically transverse (shear vertical (SV)) modes even for the normally-propagating elastic waves. The third aspect, which depends on the utilized intention of the multilayered piezoelectric phononic crystals, is that the propagation direction of the elastic waves may affect their properties, even as the constituent materials have weak anisotropy (i.e., high symmetry). The fourth aspect is that the multiple reflection and transmission together with possible mode conversion of elastic waves at the surfaces/interfaces of constituent layers give rise to the complex bulk wave modes in infinite composites, or the surface/interface waves including the surface SH, Rayleigh, Love, and Stoneley waves in semi-infinite half-spaces, or the guided waves such as Lamb wave in finite thickness slabs. Note that the first aspect is exclusive of piezoelectric multilayered phononic crystals, while the latter three aspects are all common to multilayered phononic crystals of any materials. In addition, the second aspect, i.e., the material anisotropy actually affects the coupling and the electric-field dependence of wave modes in combination with the third and fourth aspects, i.e., the propagation direction and surface/interface property. For example, as the constituent piezoelectric and elastic layers in the multilayered piezoelectric phononic crystals all have enough high crystal symmetries, the mechanical waves thereof will be pure P-, pure SV-, and pure SH-waves, among which only one is influenced by the electric field and boundary. Nonetheless, the uncoupling requirement to the crystal symmetry may be different for disparate propagation direction. When the constituent materials entail crystal symmetries that do not satisfy the uncoupling condition of mechanical waves in a specific direction, some or all wave modes propagating in that direction will be coupled. The uncoupling or coupling of mechanical waves also influences the complexity of bulk wave and the formation of surface, interface, and guided waves.

Taking the above four aspects into account, we can review the research progresses on elastic waves in laminated piezoelectric phononic crystals by classification in two levels including the electrical boundary and wave type as follows.

Firstly, the investigations considering passive-type electrical boundaries (actually nearly all literature so far, except those specified otherwise concerning the electric open condition although not clearly pointed out), will be surveyed according to the wave type studied. For all kinds of bulk waves, including the pure SH wave, the uncoupled P and SV waves, and the coupled P-SV wave, propagating both normally and obliquely in infinite media with transversely isotropic (hexagonal crystal) constituent materials: Sapriel & Djafari-Rouhani [4] and Nougaoui & Djafari-Rouhani [5] summarized some research achievements before 1990 involving the dispersion curves, the effective constants and the analysis methods. Li and Wang [20,22] and Li et al. [21,23] studied the localization factor and length in randomly disordered piezoceramic–polymer composites by the transfer matrix method (TMM) with particular attention on the disorders of the thicknesses and thickness ratio of constituent layers and the piezoelectric/elastic constants of the piezoelectric layer. The effects of the piezoelectricity, the piezocomposite sort, the propagation direction, the nondimensional wavenumber, and the electrical potential on wave band and localization were discussed. Guo and Wei [24], considering initial stresses and their effects on constitutive equations, governing equations and boundary conditions, analyzed by TMM the dispersion curves in phononic crystal composed of two piezoelectric materials, whose dependences on the normal and shear initial stresses and the corresponding boundary conditions were discussed based on the numerical results. Golub et al. [25] and Fomenko et al. [26] analyzed by TMM the dispersion properties (such as dispersion curves and transmission/reflection coefficients), the localization factor and the classification of pass-bands and band gaps in phononic crystals composed of a specific number of periodically arranged unit cells with homogenous or functionally-graded interlayers

and elastic half-spaces on both sides, whose dependences on the incident angle and the gradation and geometrical properties of interlayers were also discussed. Very recently, this functionally-graded model was extended by Fomenko et al. [27] to two cases, i.e., the infinite layered phononic crystals and finite counterparts between two isotropic half-spaces. The unit cell of both cases was composed of four piezoelectric sublayers with two being homogeneous and two being functionally graded. A semi-analytical method based on the transfer matrix of a unit cell was proposed to analyze the dispersion curves (phase and attenuation coefficients), the energy transmission coefficients, the localization factor and the classification of pass-bands and band gaps, whose dependences on the number of unit cells, the angle and type of incident waves, the thickness and material properties of the functionally graded sublayers, and the geometrical and material properties of the homogeneous layers were also discussed. Chen et al. [28] and Yan et al. [29], based on the nonlocal piezoelectricity continuum theory, analyzed the dispersion curves and the localization factor by TMM and the transmission/reflection spectra by the stiffness matrix method (SMM), respectively, in nanoscale phononic crystals consisting of two piezoelectric materials. The dependences of these wave propagation behaviors together with the found cutoff frequency on the ratio of internal to external characteristic lengths (R), the piezoelectric constant, the impedance ratio, and the incident angle were discussed with referring to numerical results. The influence of R on the mode conversion and the influence of the mode conversion on band gaps were also analyzed. Only concerning the pure SH waves in infinite media with transversely isotropic (hexagonal crystal) constituent materials: Zinchuk et al. [30] proposed a matrizant method to analyze the dispersion characteristics with particular attention on the effects of the piezoelectricity, the unit cell configuration, and the relative thickness ratio in numerical examples. This analysis was later extended by Zinchuk et al. [31] to periodic medium made of alternating metal and piezoelectric layers. Alshits and Shuvalov [32] analyzed the reflection/transmission of inclined SH waves in periodic composites having finite number of unit cells consisting of identical piezoelectric layers with metallized interfaces of fixed electropotential at alternating distances and semi-infinite substrates on its both sides, in which the existence of Bragg resonances was revealed. Qian et al. [33] obtained the phase velocity equations of normally and parallelly propagating waves in piezoelectric-elastic composites, and discussed the basic wave properties such as the filter effect and the effects of thickness and shear modulus ratios of the piezoelectric layer to elastic layer on phase velocity. Lan and Wei [34] studied the influence of imperfect interfaces (modeled by the mass-spring parameters) on dispersion curves and band gaps of both normally and obliquely propagating waves in piezoceramic–polymer composites with incidental attention on the piezoelectricity and the thickness-ratio effects. Only regarding the pure P waves, Faidi and Nayfeh [35] developed an improved continuum mixture model for analyzing the dispersion curves (phase velocity spectra) of the parallelly propagating longitudinal waves in periodic bi-laminated orthotropic piezoelectric media. As far as only the coupled P-SV waves are concerned, Geng and Zhang [36] analyzed the dispersion curves of parallelly propagating coupled P-SV waves in periodic piezoceramic–polymer composites by the method of partial wave expansion with special attention on the effects of the volume fraction and the polymer properties, which are partially validated through the experimentally measured thickness resonance (with polarization parallel to interface) and lateral resonance (with polarization along periodic direction) spectra. Regarding the three-dimensional (3D) coupled waves, besides some occasional discussions in Fomenko et al. [27], Podlipenets [37] presented without validation a Hamiltonian system formalism to analyze the dispersion equations of bulk, surface, and plate waves in respectively the infinite, semi-infinite, and finite phononic crystals with constituent materials of *mm2* or higher symmetry crystal. With respect to the surface/interface and guided waves, the laminated semi-infinite and finite transversely isotropic piezoelectric phononic crystals, respectively, with bounding plane either parallel or perpendicular to the layering plane have all been considered. Initially in the parallel case, Zinchuk et al. [38] analyzed the dispersion curves of the SH-type surface, Stoneley, and guided waves and discussed the state variables distribution of corresponding modes in periodic piezoelectric and metal composites under metallized mechanically free condition. The distinctive/interrelation characteristics of dispersion spectra for guided (normal) wave in even-layered

periodic finite thickness plate with those for the surface wave in semi-infinite half-space were also discussed. Yan et al. [39] investigated the dispersion curves and the mechanical displacements and electrical potential variations of the symmetrical Lamb waves in nanoscale periodic piezoelectric composites based on the nonlocal piezoelectricity continuum theory. The influences of the piezoelectric effect, the ratio of internal to external characteristic lengths (i.e., R representing nanoscale size-effect), and the volume fractions on these wave behaviors particularly like the found mode conversion and cut-off frequency were analyzed in details. Later in the perpendicular case, Wang et al. [40] analyzed using TMM the localization factor of Rayleigh waves in periodic piezoceramic–polymer half-space due to disorders in polymer thickness or piezoceramic material constant, and discussed its influence by the thickness ratios of constituent layers. Alippi et al. [41] proposed an approximate theoretical model by neglecting the mode coupling to analyze the dispersion curves of Lamb wave below the frequency of thickness resonance. By experimentally exciting the band edge resonances, the analysis was validated and the fractional volume dependence of stop-bands was discussed. This work was later improved by Craciun et al. [42] who experimentally and theoretically studied the resonance spectra of the lateral and thickness resonances with polarizations along and perpendicular to the periodic direction, respectively, based on the Kronig-Penney and effective medium models, when considering the anisotropy, piezoelectricity, and volume fraction factors, which together with the coupling of two resonant modes were verified by the corresponding dispersion curves of pure elastic P and electroelastic SV waves. Note from the above literatures considering passive electrical boundaries that although the passive-type electrical boundaries are easy to realize and already provide the laminated piezoelectric phononic crystals with abundant elastic wave properties such as the frequency bands, but these wave properties are fixed or work at fixed frequency ranges once the laminated piezoelectric phononic crystals are fabricated. This passive feature obviously limits the application of laminated piezoelectric phononic crystals with passive electrical boundaries in working occasions where the external dynamic excitations have varying frequencies. Therefore, active-type electrical boundaries are proposed to tune the elastic wave properties such as the frequency bands anytime for adapting them to the external excitations. Two approaches have been proposed so far to actively adjust the electrical boundaries: One is to switch between diverse passive electrical boundaries and the other is to connect to external electrical network with alterable parameters such as capacitance, inductance, resistance, voltage, etc. In what follows, the studies related to the two adjusting manners will be reviewed successively.

Secondly, the investigations simultaneously considering many passive electrical boundaries for possibly switching between will be surveyed according to the wave type studied. For all kinds of bulk waves, including the pure SH wave, the uncoupled P and SV waves, and the coupled P-SV wave, propagating both normally and obliquely, Guo et al. [43] used TMM to analyze the dispersion curves of the phononic crystals composed of two transversely isotropic piezoelectric materials, considering four kinds of dielectric imperfect interfaces including weak conducting, high conducting, low dielectric, and grounded metallized interfaces together with four kinds of mechanical imperfect interfaces such as normal compliant, tangent compliant, tangent fixed, and tangent slippery interfaces. The influences of the piezoelectricity and these mechanically and dielectrically imperfect interfaces on the dispersion curves were discussed based on the numerical results. As far as only the pure SH waves in infinitely laminated piezoelectric phononic crystals with transversely isotropic (hexagonal symmetric) constituent materials are concerned, Ghazaryan and Piliposyan [44] investigated the effects of three kinds of interfaces, including electrically shorted with mechanically continuous condition, magnetically closed with mechanically continuous condition, and electrically open with mechanically smooth contacts, on dispersion properties such as band structures and bandgap width of obliquely-propagating Bloch–Floquet waves, with special attention on the factors like piezoelectricity, incident wavenumber, and filling fraction. In cases of electric-open and electric-short conditions, Zhao et al. [45] computed the dispersion curves, the transmission coefficients, and the reflective spectrum by the global transfer matrix method with attention on the effects of the piezoelectricity, the volume fraction and the propagation direction. With respect to coupled P-SV waves in infinite media, Zhao et al. [46], also in cases of

electric-open and electric-short conditions and by the modified transfer matrix recursive algorithm, theoretically studied the dispersion curves (band gaps) and the transmission/reflection spectra in periodic structures containing fluid layer and orthotropic piezoelectric layer with special attention on the influences of piezoelectricity and propagation direction. As regard to the surface/interface or guided waves, the usually considered laminated piezoelectric semi-infinite or finite phononic crystals have bounding plane parallel to the layering plane. For this kind of phononic crystals, Sapriel & Djafari-Rouhani [4] and Nougaoui & Djafari-Rouhani [5] summarized research results on the dispersion curve of SH-type surface wave in cases of metallized or non-metallized surfaces before 1990s, with special attention on the influence of the nature of the film at the surface. Considering the same electrical boundaries, Zinchuk et al. [47] and Zinchuk & Podlipenets [48] analyzed the dispersion curves of SH-type surface and guided waves, the corresponding mode distributions of state variables by a matrizant method based on periodic Hamiltonian system, as the constituent materials belong to *6mm* crystal class. The effects of the electrical boundaries, the unit cell configuration and relative thickness ratio on the dispersion properties were studied by numerical examples. Zinchuk and Podlipenets [49] introduced their previous method to the analysis of dispersion curves of Rayleigh waves also with the *6mm* materials assumption for three kinds of electrical boundaries including electric-short, non-metalized contact with a vacuum and electric-open conditions, whose influence by the piezoelectricity was examined in the case of free mechanical boundary. Alami et al. [50], on basis of Green's function method, derived and validated the dispersion relation and state density of the SH-type surface waves in a semi-infinite superlattice composed of *6mm* class piezoelectric–metallic layers and capped with a piezoelectric layer with open- or short-circuited surfaces. The interaction between this SH-type surface wave and the possibly-existing interface, guided and pseudo-guided (leaky) waves induced by the cap layer and the dependences of electromechanical coupling coefficient on the cap layer thickness and the metallic layer property were also investigated. In addition, the laminated piezoelectric finite phononic crystals with bounding planes perpendicular to the layering plane have also been considered. The model studied so far is the plane-strain plate consisting of piezoceramic–polymer or two piezoelectric layers with hexagonal (usually *6mm*) symmetry and with the length in the periodic direction and the thickness parallel to the interfaces. In cases of electric open and short boundaries, Otero et al. [51] computed the dispersion curves of Lamb waves along periodic direction using the effective coefficients approximated by the first order asymptotic homogenization method considering four piezoelectric volume fractions. The behaviors of the mechanical displacement and the electric potential for different wave modes were illustrated. Considering the same electric boundaries and Lamb wave, Zou et al. [52] studied, using the plane wave expansion (PWE) method, the band properties such as the dispersion curves, widths and starting frequencies of the first bandgaps and their influences by the filling fraction, the thickness to lattice pitch ratio, and the polarizations directions. In cases of non-piezoelectricity, open-circuit, and short-circuit conditions, Zhu et al. [53] analyzed and compared the dispersion curves, transmission spectra, eigenmode displacements of different modes by the finite element method (FEM) on COMSOL Multiphysics software, based on which the piezoelectric-sensitive mode was defined and its physical mechanism such as its dependences on the piezoelectric constants, the filling ratio, and the ratio of thickness to lattice pitch were revealed. Considering both metallized and non-metallized interfaces, Piliposyan et al. [54] analyzed by TMM the dispersion curves and the reflection/transmission of inclined SH-type guided waves in the infinite periodic composites or in the finite counterparts with a defect layer as the mirror symmetry center and two piezoelectric half-spaces on both sides as substrates. The dependences of these wave behaviors on the ratio of the unit cell's length to the waveguide's height, the piezoelectric material properties, the boundary condition distribution on the lower and upper walls, and the presence of defect layer were discussed, with special attention on the Bragg resonances and the presence of trapped modes and slow waves. Notice from the above surveyed research works that the switching method between diverse passive electrical boundaries can adjust the frequency bands and other wave properties of laminated piezoelectric phononic crystals among several discrete states. However, in practical applications there

is a trend toward the flexible usage of the frequency spectrum since the external dynamic excitations in working environments usually have continuously-varying frequency components rather than discrete-varying ones. This trend requires continuously-tunable frequency bands and other wave properties, which can be realized via installing external electrical circuits with alterable electrical parameters on the laminated piezoelectric phononic crystals that serves as the other approach to actively adjust the electrical boundaries. Next, the studies on elastic waves in laminated piezoelectric phononic crystals with connecting to external electrical network will be reviewed.

Thirdly, by the category of studied wave type and the alterable parameters such as capacitance, inductance, resistance, and voltage in the connected electrical circuits, the investigations on laminated piezoelectric phononic crystals possessing continuously-varying wave behaviors will be surveyed. For all kinds of uncoupled waves including pure P and S (SH or SV) modes, Li et al. [55] studied the dispersion curves of waves along the thickness direction in infinite media with orthotropic materials (the lowest symmetry crystal guaranteeing the decoupling of three wave modes) in cases of four electrical boundaries including the electric-open, applied electric capacitance, electric-short, and applied feedback voltage that are capable of both discretely and continuously tuning the frequency bands. Particular attentions were also on the influence of electrode thickness and on the characteristics of wave dispersion. Regarding to pure P waves in infinite phononic crystals with the faces of the piezoelectric layers being electroded and connected to electrical capacitor: Ponge et al. [56,57] briefly reviewed the dispersion properties induced by the periodic electrical boundary conditions (with open-circuit and short-circuit as reference) in homogeneous piezoelectric material, and then accordingly designed a Fabry-Perot cavity whose length might be tuned by a spatial shift of electrical boundaries along with the position of the transducer driven electrically by a voltage. The optimum performance of the cavity device achieved through compromise among the series/parallel resonance frequencies, the band gaps of phononic crystal (influenced by number and length of unit cells), and the coupling coefficient was corroborated by the 1D analytical analysis using TMM, the numerical simulation using finite element method (FEM), and the experiment. Kutsenko et al. [58,59] analyzed using TMM the dispersion spectra of normally-propagating wave in periodic structure of identical transversely isotropic piezoelectric layers or of alternating elastic and piezoelectric layers. Special attentions were on the unusual features of dispersion spectra in the case of negative capacitance value, such as the quasistatic stop-band, the poles of attenuation constant spectra corresponding to jumps of phase constant from 0 to π, and the occurrence of quasiflat dispersion branches and dispersion curves with infinitely growing group velocity [58,59]. Kutsenko et al. [60] extended the model to a general case where a unit cell may consist of several piezoelectric or elastic–piezoelectric multilayers (possibly functionally graded), and studied the effective properties characterizing the homogenized medium such as effective density, elastic constant, as well as Willis coupling constant, and their tunability by the capacitance. Recently, Kutsenko et al. [61] further extended the electrical boundary to 2D semi-infinite periodic network whose unit cell contains two capacitors in parallel and in series, and studied the dispersion spectra and wave fields by deriving explicit equations. The control of the dispersion spectrum, which can come out either as a discrete set of curves, or as a continuous band, or as a superposition of both, by the signs and values of the two capacitances in the unit cell was studied with considering examples and the limiting cases of open and short circuits. As far as pure P waves in infinite media with electrical circuits other than simple capacitor are concerned, Mansoura et al. [62] investigated theoretically the dispersion curves of normally-propagating waves in phononic piezoelectric-elastic crystals, paying special attention to the newly-opening hybridization gap due to a coupling of electrical resonance with the wave propagation as the electrical impendence was inductances. The results were verified experimentally by the measured transmission. Zhu et al. [63] proposed an active acoustic metamaterial consisting of periodic layers of steel, polyuria, and piezoelectric ceramic transducer (PZT) with the PZT layer shunted by an inductor. Its band structure and transmission coefficient were calculated by the TMM, and its effective material parameters were computed by homogenization method. The extremely narrow stop band induced by the resonance circuit was able to be controlled through the inductor and

the corresponding effective parameters behaved negative. Parra et al. [64] calculated using TMM the dispersions, impedances, and displacements in phononic crystals with periodically distributed local and interconnected LC (inductance-capacitance) shunts, which were experimentally validated. The ability of local or interconnected shunts to control the width or depth of bandgap was also discussed. As regard to the coupled P-SV (or identical P-SH) waves, Fomenko et al. [65] analyzed using TMM the dispersion curves, transmission coefficients, and localization factor of the obliquely-propagating waves in infinite phononic crystals with the unit cell consisting of two piezoelectric interlayers and with or without electric potentials at the metallic interfaces of two kinds of configurations, i.e., an infinite periodic structure and a periodic structure surrounded by two half-spaces. The influences of the electrical potential and the incident angle on these wave properties were also investigated. As for the pure longitudinal waves and vibrations along the thickness in finite thickness plates have bounding plane parallel to the layering plane: Mansoura et al. [66,67] theoretically and experimentally investigated the electrical impedance of one piezoelectric layer within phononic crystals made of alternating piezoelectric and elastic layers while the other piezoelectric layers were connected to the external circuit including electric-open, electric-short, and external capacitance conditions, which is related to the band structures and transmissions together with the electromechanical coupling. Mansoura et al. [68] further verified, by comparing the experimentally measured transmissions with theoretically analyzed dispersion curves, the wave control in phononic piezoelectric elastic crystals through the negative capacitance connected on the electroded faces of piezoelectric layers that broadening band gaps. Darinskii et al. [69] theoretically studied using TMM the reflection/transmission along the 6-fold symmetry axis in periodic piezoelectric structure constructed by inserting infinitesimally thin metallic electrodes into a homogeneous piezoelectric material of *6mm* symmetry class and then interconnecting each two successive electrodes by an external capacitor (with electric-open and electric-short conditions corresponds to zero and infinity capacitance). Allam et al. [70] proposed an active acoustic metamaterial (AMM) with the unit cell consisting of transversely isotropic piezoelectric and isotropic brass layers clamped in air, whose effective density in real-time might be controlled by the closed feedback control loop connected between the piezoelectric layers. The stability, characterization, and behavior of the AMM were predicted by vibro-acoustic analytic model and were verified experimentally. The adaptive and programmable control of AMM's effective density through three types of controllers was also studied. Wang et al. [71] theoretically studied using TMM the transmission and the pass-band in piezoelectric laminated phononic crystals inserted with a 0.2 mol% Fe-doped relaxor-based ferroelectric 0.62 Pb(Mg1/3Nb2/3)O3–0.38PbTiO3 (PMN–0.38PT) single crystal defect layer, whose dependences on the thickness/strain adjusted nonlinearly by the external electrical voltage and on the acoustic impedance of constituent materials were analyzed. As for the vibrations associated with coupled P-SV waves propagation in finite thickness plates that have bounding plane perpendicular to the layering plane, Geng and Zhang [36] studied theoretically and experimentally the vibration properties of thickness and lateral resonances (such as the surface vibration profile, the electric impedance distribution, and the electromechanical coupling factor) under an external driving AC electric field in both air and water media and under external incident pressure wave from the water, respectively, with particular attention on the influence of the aspect ratio.

In summary of the above literature review on various types of elastic waves in laminated piezoelectric phononic crystals with both passive and active electrical boundary conditions, we propose the following five insufficiencies that exist in the investigations so far:

1. Although the electric-open, the electric-short, and the applied electrical capacitance boundaries have all received plenty attentions, but the feedback control condition, which plays an important role in vibration control of structures in many engineering fields [72,73], is nearly not concerned except that in Allam et al. [70] the feedback strategy was adopted to control the real-time effective density of an unit cell of the proposed active acoustic metamaterial involving pure P wave propagation in laminated piezoelectric phononic crystals of finite thickness. Hence, the authors Li et al. [55] extended the external feedback control voltage boundary into the infinite media and investigated its contrast with the electric-open, the electric-short, and the applied electrical capacitance boundaries when considering the pure P- and S-waves propagation. However, when the coupled mechanical waves in laminated piezoelectric phononic crystals are investigated, the applied feedback control electrical condition has never been considered so far.


In order to make up for the above deficiencies, this paper analyzes, by introducing the TMM [76] based on the state space formalism [77], the coupled elastic waves, including P-SH, P-SV, SV-SH, and P-SV-SH modes, along thickness direction in infinitely laminated piezoelectric phononic crystals with the unit cell consisting of any number of arbitrarily anisotropic piezoelectric and elastic layers and having four electrical boundaries, such as the electric-open, the applied electrical capacitance, the electric-short, and the applied feedback control conditions. The electrodes on the surfaces of the piezoelectric interlayers are all modeled as elastic layers along with the inserted elastic materials themselves, whose mechanical and electrical functions are both considered in the analysis. The configuration of the analysis model hereof renders the presence of coupled waves without introducing the effect of inclined propagation direction. Consequently, the forming of coupled waves is solely caused by the material anisotropy, the individual effect of which on the multifarious frequency related dispersion curves in cases of four electrical boundaries is the main motivation of our research. Compared with our previous similar work [55], this paper innovatively considers the coupled-mode waves due to material anisotropy as an extension to the previously concerned pure P and pure S waves and as a remedy to deficiencies (2) and (3), in spite that a similar analysis process is adopted.

This paper is organized as follows. Following the research background and motivation in Section 1, the basic model of the general periodically laminated arbitrarily-anisotropic piezoelectric composites is

described in Section 2. Section 3 provides the state space formalism for describing the elastodynamics of the constituent layers. Based on the wave solutions to the state variables of constituent layers, the formulation of transfer matrix method (TMM) is presented in Section 4 to establish the governing relations of the unit cell, which are combined with the Floquet–Bloch theorem [2,78,79] to bring out the dispersion equation of the whole system. The numerical examples for phononic crystals with various layouts of piezoelectric materials and four kinds of electrical boundaries are given in Section 5, with special attention on the general features of diversified frequency related dispersion curves and their dependences on the electrode thickness and electrical boundaries. Some findings from this research are drawn as conclusions in Section 6.

#### **2. Basic Model**

Consider the elastic waves propagating along the thickness direction, i.e., perpendicular to the interface, in infinitely periodically laminated piezoelectric composite structures whose unit cell, as shown in Figure 1, consists of any number of arbitrarily anisotropic (triclinic) piezoelectric and elastic layers with the piezoelectric interlayers having anyone of four electrical boundaries such as the electric-open, the applied electrical capacitance, the electric-short, and the applied feedback control conditions. These electrical boundaries are applied through the electrodes coated on the surfaces of piezoelectric layers. Here in this paper, the electrodes are also modeled as elastic layers along with the inserted elastic materials themselves, whose mechanical and electrical functions are both considered in the analysis, so that their effect on the wave propagation can be revealed. According to the Floquet–Bloch theorem [2,78,79], the unit cell model in Figure 1 together with the periodic boundary condition is adequate for the analysis of characteristic waves with wavenumber *q* and circular frequency ω. Since arbitrarily anisotropic (triclinic) materials are considered in the model, whose crystal symmetry is far lower than orthotropy, the lowest crystal symmetry requirement for decoupling the elastic waves along thickness direction [55]. Therefore, coupled-mode waves are actually analyzed via the model. In the unit cell shown by Figure 1, altogether any number (*m*) of interlayers, involving the shaded elastic layers, gray piezoelectric layers, and black electrode layers, and correspondingly *N* (*N* = *m* + 1) interfaces are assumed for the sake of modeling all sorts of structural configurations.

As shown in Figure 1, the constituent layers, with finite thicknesses and being unbounded on their layering planes, in the unit cell are labeled in sequence from top to bottom as (1),(2), ··· ,(*j*), ··· ,(*m*), while correspondingly those surfaces/interfaces are denoted as 1, 2, ··· , *J*, ··· , *N*, for the convenience of description. The thickness of a typical layer, say (*j*), is denoted as *h*(*j*), shown both in Figure 1 and in the enlarged view of layer (*j*) in Figure 2b, while the thickness of the unit cell is *<sup>H</sup>* = #(*m*) (*j*)=(1) *h*(*j*). For the piezoelectric layer (*j*) as depicted by Figure 2b, a local coordinate system *x*(*j*), *y*(*j*), *z*(*j*) is established with its origin on the top surface of the layer for facilitating the description of physical quantities. *ux*(*j*), *uy*(*j*), and *uz*(*j*) signify the displacements along *x*, *y*, and *z* axes at any position on a plane within layer (*j*) parallel to the layer surfaces, while τ*zx*(*j*), τ*zy*(*j*), and σ*z*(*j*) denote the corresponding stresses on that plane. These mechanical quantities, if the linear theory of piezoelectricity [80–82] are resorted, are known to be related to the six partial mechanical waves whose wavenumbers are the components of vector **Λ**(*j*). The six partial waves can be divided into two groups propagating respectively in downward and upward directions and into three pairs representing respectively the three modes. Because of the piezoelectric effect in the piezoelectric layer, these partial waves can also be tuned by the electrical boundaries applied on the electrode layers, which are coated on the piezoelectric layer surfaces and are connected to the four external circuits including the electric-open, the applied electric capacitance, the electric-short, and the applied feedback control conditions. Both switching between the four electrical boundaries and adjusting the applied capacitance *C*(*j*) in the case of connecting external capacitor or the gain coefficient *Kg*(*j*) in the connecting feedback control condition can actively change the electric charge on the electrode *Q*(*j*) and the voltage (electric potential difference) between electrodes *V*(*j*), and further alter the thickness distribution of the electric potential

φ(*j*), the electrical displacement *Dz*(*j*) along the thickness direction, and the intensity of electrical field *Ez*(*j*) along the thickness direction. Note that in the case of applied feedback control boundary, *V*(*j*) is related to *Kg*(*j*) through *V*(*j*) = −*Kg*(*j*)[*uz*(*j*)(*h*(*j*)) − *uz*(*j*)(0)], where *uz*(*j*)(0) and *uz*(*j*)(*h*(*j*)) represent the displacements along the thickness direction on the top and bottom surfaces of layer (*j*) at the sampling position, respectively.

The descriptions of the local coordinates and the pertaining physical quantities of other piezoelectric layers and of elastic layers denoting either the electrodes or the inserted elastic materials are exactly the same as those for piezoelectric layer (*j*), except that the quantities related to electrical field of elastic layers are not provided as can be noticed from Figure 2a for elastic layer (*i*), because they will not be involved in later formulation. Moreover, no electrical boundaries are applied on the elastic layers, as also noted from Figure 2a. Here and after, all the quantities pertaining to a layer are denoted by a subscript of layer label, while those pertaining to a surface/interface will be signified by a superscript of surface/interface label. All the vector-type physical variables are deemed as positive when their directions are coincident with the positive coordinate axes, while the scalar electric quantities are deemed as positive when they are corresponding to the positive electric charge, and vice versa.

**Figure 1.** The unit cell of the general periodically laminated piezoelectric composite structures with arbitrarily anisotropic constituent interlayers and the description about the propagation characteristics of coupled elastic waves along thickness direction.

(**a**)

**Figure 2.** The descriptions of the local coordinate systems and the physical quantities of typical elastic/piezoelectric layers in the unit cell. (**a**) A typical elastic layer. (**b**) A typical piezoelectric layer covered by the electrode layers with four electrical boundaries: electric-open, applied capacitance, electric-short, and applied feedback control.

#### **3. State Space Formalism**

Since the considered piezoelectric phononic crystals and their unit cells are laminated, the state space formalism [77], which is essentially a displacement–stress hybrid method suitable for structures with unidirectional configuration, is introduced to describe the elastodynamics of the constituent layers. For any elastic or piezoelectric constituent layer, by using the method of separation of variables and by choosing the displacements and stresses on the plane parallel to the surfaces as the components in the state vector, the state equation governing the spectral state vector in frequency–wavenumber domain can be derived from all the fundamental equations of three-dimensional (3D) elasticity and piezoelectricity, respectively. The solution to the state equation can then be obtained by virtue of the theory of first order ordinary differential equations. It should be pointed out that all the equations and their solutions in this section are given in the local coordinate system, but the subscript indicating the pertinent layer and coordinates is omitted for brevity.

#### *3.1. The State Equation of a Layer Derived from the 3D Elastodynamics*

The derivation of the state equation for elastic and piezoelectric layers will be provided successively, although their route is generally similar except for the introduction of electrical boundaries to the piezoelectric layer. For the piezoelectric layer, the differences of the deriving process as compared to that for the elastic layer will be emphasized, and the similarities will be abbreviated, as indicated in Section 3.1.2.

#### 3.1.1. The State Equation of An Elastic Layer

According to the 3D elasticity, the fundamental equations governing the elastodynamics of a typical arbitrarily anisotropic elastic layer (*i*) as seen in Figure 2a can be written based on mechanical and electrical considerations [82]. The mechanical equations include the equation of motion (without body force here), the strain–displacement relation, and the elastic constitutive equation, which are written in fully matrix notation as

$$\rho \left(\nabla\_s\right)^T \mathbf{T} = \rho \frac{\partial^2 \mathbf{u}}{\partial t^2}, \mathbf{S} = \nabla\_s \mathbf{u}, \ \mathbf{T} = \mathbf{c} \mathbf{S},\tag{1}$$

where the superscript "T" denotes the transposition of a matrix or a vector here and after; **u** = [*ux*, *uy*, *uz*] T, **T** = [*Txx*, *Tyy*, *Tzz*, *Tyz*, *Tzx*, *Txy*] T, and **S** = [*Sxx*, *Syy*, *Szz*, 2*Syz*, 2*Szx*, 2*Sxy*] <sup>T</sup> are the displacement, stress, and strain vectors, respectively; ρ denotes the mass density; and ∇*<sup>s</sup>* and **c** are the 6 × 3 operator matrix and the 6 × 6 symmetric elastic stiffness constant matrix, respectively, with the following components

$$
\nabla\_{\mathbf{s}} = \begin{bmatrix}
\partial / \partial \mathbf{x} & \mathbf{0} & \mathbf{0} \\
0 & \partial / \partial \mathbf{y} & \mathbf{0} \\
0 & 0 & \partial / \partial \mathbf{z} \\
0 & \partial / \partial \mathbf{z} & \partial / \partial \mathbf{y} \\
\partial / \partial \mathbf{z} & \mathbf{0} & \partial / \partial \mathbf{x} \\
\partial / \partial \mathbf{y} & \partial / \partial \mathbf{x} & \mathbf{0}
\end{bmatrix}, \mathbf{c} = \begin{bmatrix}
\boldsymbol{c}\_{11} & \mathbf{c}\_{12} & \mathbf{c}\_{13} & \mathbf{c}\_{14} & \mathbf{c}\_{15} & \mathbf{c}\_{16} \\
\boldsymbol{c}\_{21} & \mathbf{c}\_{22} & \mathbf{c}\_{23} & \mathbf{c}\_{24} & \mathbf{c}\_{25} & \mathbf{c}\_{26} \\
\boldsymbol{c}\_{31} & \mathbf{c}\_{32} & \mathbf{c}\_{33} & \mathbf{c}\_{34} & \mathbf{c}\_{35} & \mathbf{c}\_{36} \\
\boldsymbol{c}\_{41} & \mathbf{c}\_{42} & \mathbf{c}\_{43} & \mathbf{c}\_{44} & \mathbf{c}\_{45} & \mathbf{c}\_{46} \\
\boldsymbol{c}\_{51} & \mathbf{c}\_{52} & \mathbf{c}\_{53} & \mathbf{c}\_{54} & \mathbf{c}\_{55} & \mathbf{c}\_{56} \\
\boldsymbol{c}\_{61} & \mathbf{c}\_{62} & \mathbf{c}\_{63} & \mathbf{c}\_{64} & \mathbf{c}\_{65} & \mathbf{c}\_{66}
\end{bmatrix}.
\tag{2}
$$

The electrical equations under the quasi-static assumption to electric field include the charge equation (Gauss's law) of electrostatics (with free charge here), the electric field-electric potential relation, and the electric constitutive equation, which are written in fully matrix notation as

$$(\nabla)^{\mathsf{T}}\mathbf{D} = 0,\ \mathbf{E} = -\nabla\phi\_{\prime}\ \mathbf{D} = \mathbf{e}\mathbf{E},\tag{3}$$

where φ is the scalar electric potential; **D** = [*Dx*, *Dy*, *Dz*] <sup>T</sup> and **E** = [*Ex*, *Ey*, *Ez*] <sup>T</sup> denote the electric displacement and electric field intensity vectors, respectively; and ∇ and ε are the 3 × 1 Hamilton operator vector and the 3 × 3 dielectric constant matrix, respectively, having the following components

$$
\nabla = \begin{bmatrix}
\partial / \partial \mathbf{x} \\
\partial / \partial \mathbf{y} \\
\partial / \partial \mathbf{z}
\end{bmatrix}, \mathbf{c} = \begin{bmatrix}
\varepsilon\_{11} & \varepsilon\_{12} & \varepsilon\_{13} \\
\varepsilon\_{21} & \varepsilon\_{22} & \varepsilon\_{23} \\
\varepsilon\_{31} & \varepsilon\_{32} & \varepsilon\_{33}
\end{bmatrix}.
\tag{4}
$$

According to the method of separation of variables, the harmonic solutions to the physical quantities can be given as

$$
\Gamma(\mathbf{x}, y, z, t) = \hat{\Gamma}(k\_x; k\_y; z; \omega) \mathbf{e}^{i(\omega t - k\_x x - k\_y y)} = \hat{\Gamma}(z; \omega) \mathbf{e}^{i\omega t} (\Gamma = \mathbf{u}, \mathbf{T}, \mathbf{S}, \phi, \mathbf{D}, \mathbf{E}),
\tag{5}
$$

where <sup>i</sup> <sup>=</sup> <sup>√</sup> −1 is the imaginary unit; ω is the circular frequency; *kx* and *ky* are the wavenumbers along the local coordinates *x* and *y*, respectively; and the superscript "ˆ" means corresponding physical quantities in the *kx* −*ky* − ω domain. In Equation (5), the vanishing of wavenumbers in the layering plane *kx* = *ky* = 0 has been taken into account, because only the elastic waves propagating along the thickness direction are considered in this paper. Besides, since the considered layer is infinite in the layering plane and has finite thickness, then only the displacements and stresses on the surfaces will appear in the mechanical boundaries. Consequently, we choose the displacements and stresses on the plane parallel to the surfaces [77] as the components of the state vector **v**ˆ = [(**v**ˆ *<sup>u</sup>*) T,(**v**ˆ*T*) T] <sup>T</sup> = [*u*ˆ*x*, *<sup>u</sup>*ˆ*y*, *<sup>u</sup>*ˆ*z*, *<sup>T</sup>*<sup>ˆ</sup> *zx*, *<sup>T</sup>*<sup>ˆ</sup> *zy*, *<sup>T</sup>*<sup>ˆ</sup> *zz*] T , with **v**ˆ *<sup>u</sup>* = **u**ˆ = [*u*ˆ*x*, *u*ˆ*y*, *u*ˆ*z*] <sup>T</sup> and **v**ˆ*<sup>T</sup>* = [*T*ˆ *zx*, *T*ˆ *zy*, *T*ˆ *zz*] <sup>T</sup> being referred to as the displacements and stresses state vectors, respectively.

Note from Equations (1) and (3) that the mechanical and electrical fields in the elastic layer are mutually independent each other. Since here in this paper the mechanical field rather than electrical field is continuous across the adjacent layers, only Equation (1) is used to derive the state equation of the elastic layer. Substitution of Equation (5) into the first formula in Equation (1) leads to

$$\frac{d\hat{\mathbf{v}}\_T}{d\mathbf{z}} = -\rho a^2 \hat{\mathbf{v}}\_{\mathbf{u}}.\tag{6}$$

Substituting the second formula into the third formula in Equation (1), then introducing Equation (5) into the resulting equation **T** = **c**∇*s***u**, and finally selecting from the consequent equation out those pertaining to the stress state vector, we can obtain

$$
\hat{\mathbf{v}}\_T = \mathbf{G} \frac{d\hat{\mathbf{v}}\_u}{dz},
\tag{7}
$$

which further leads, through matrix inversion, to

$$\frac{d\hat{\mathbf{v}}\_{\mathcal{U}}}{dz} = \mathbf{G}^{-1}\hat{\mathbf{v}}\_{T\_{\prime}}\tag{8}$$

The 3 × 3 coefficient matrix **G** is composed of

$$\mathbf{G} = \begin{bmatrix} c\_{55} & c\_{54} & c\_{53} \\ c\_{45} & c\_{44} & c\_{43} \\ c\_{35} & c\_{34} & c\_{33} \end{bmatrix} . \tag{9}$$

The combination of Equations (6) and (8) results in the state equation of the elastic layer

$$\frac{d\hat{\mathbf{v}}}{dz} = \mathbf{M}\hat{\mathbf{v}},\tag{10}$$

where the 6 × 6 coefficient matrix **M** of the state equation are formed as

$$\mathbf{M} = \begin{pmatrix} 0 & \mathbf{G}^{-1} \\ -\rho \omega^2 \mathbf{I}\_3 & 0 \end{pmatrix} \tag{11}$$

with **I**<sup>3</sup> the third order identity matrix.

#### 3.1.2. The State Equation of a Piezoelectric Layer

According to the 3D piezoelectricity, the fundamental elastodynamic equations of a typical arbitrarily anisotropic piezoelectric layer (*j*), as seen in Figure 2b, can be given. From both mechanical and electrical considerations [80–83], all the equations are identical to those of elastic layer except that the constitutive equation currently should be written in fully matrix notation as

$$\mathbf{T} = \mathbf{c}\mathbf{S} - \mathbf{e}^{\mathrm{T}}\mathbf{E}, \ \mathbf{D} = \mathbf{e}\mathbf{S} + \mathbf{e}\mathbf{E},\tag{12}$$

where **e** is the 3 × 6 piezoelectric constant matrix composed of

$$\mathbf{e} = \begin{bmatrix} \mathbf{e}\_{11} & \mathbf{e}\_{12} & \mathbf{e}\_{13} & \mathbf{e}\_{14} & \mathbf{e}\_{15} & \mathbf{e}\_{16} \\ \mathbf{e}\_{21} & \mathbf{e}\_{22} & \mathbf{e}\_{23} & \mathbf{e}\_{24} & \mathbf{e}\_{25} & \mathbf{e}\_{26} \\ \mathbf{e}\_{31} & \mathbf{e}\_{32} & \mathbf{e}\_{33} & \mathbf{e}\_{34} & \mathbf{e}\_{35} & \mathbf{e}\_{36} \end{bmatrix}. \tag{13}$$

Equation (12) indicates that the mechanical and electrical fields in piezoelectric layer are coupled through the constitutive relation. Consequently, although here in this paper only the mechanical field is continuous across the adjacent layers, the effect of the electrical field on the mechanical field needs to be considered during the derivation of the state equation for the piezoelectric layer.

*Crystals* **2019**, *9*, 426

Bear in mind that the solutions to the mechanical and electrical quantities of the piezoelectric layer can be still expressed as Equation (5), and the displacements, stresses, and whole state vectors of the piezoelectric layer are formed identically as those of the elastic layer shown before. Therefore, Equation (6) is still right for piezoelectric layer when Equation (5) is substituted into the equation of motion without body force. In addition, by substituting the strain–displacement relation as the second formula in Equation (1), the electric field-electric potential relation as the second formula in Equation (3), and Equation (5) into the constitutive Equation (12), through the similar process as before in Section 3.1.1, we gain

$$\boldsymbol{\hat{\Psi}}\_{T} = \mathbf{G} \frac{d\boldsymbol{\Psi}\_{u}}{d\boldsymbol{z}} + \mathbf{F}^{T} \frac{d\boldsymbol{\hat{\phi}}}{d\boldsymbol{z}},\\\boldsymbol{\mathcal{D}}\_{z} = \mathbf{F} \frac{d\boldsymbol{\Psi}\_{u}}{d\boldsymbol{z}} - \boldsymbol{\varepsilon}\_{33} \frac{d\boldsymbol{\hat{\phi}}}{d\boldsymbol{z}}.\tag{14}$$

where **G** is the same as Equation (9)

$$\mathbf{F} = \begin{bmatrix} \ \ c\_{35} & \ c\_{34} & \ c\_{33} \ \end{bmatrix} \tag{15}$$

and only the expression of *D*ˆ *<sup>z</sup>* among the three components of electric displacement vector is sorted out since only it can be used, after the electrical boundaries of the piezoelectric layer have been introduced, to represent dφˆ/d*z* by d ˆ**v***u*/d*z* as

$$\frac{d\hat{\phi}}{dz} = \frac{\mathbf{F}}{\varepsilon\_{33}} \frac{d\hat{\mathbf{v}}\_{u}}{dz} - \mathbf{R}[\hat{\mathbf{v}}\_{u}(h) - \hat{\mathbf{v}}\_{u}(0)].\tag{16}$$

The 1 × 3 coefficient matrix **R** is related to **F** diversely in cases of four different electrical boundaries as shown from Table A1 in Appendix A, where the derivation of Equation (16) is provided in detail. Substituting Equation (16) into the first formula of Equation (14), it is obtained by further simplification that

$$\frac{d\boldsymbol{\Psi}\_{\boldsymbol{u}}}{d\boldsymbol{z}} = \left(\boldsymbol{\mathbf{G}} + \frac{\mathbf{F}^{\mathrm{T}}\mathbf{F}}{\varepsilon\_{33}}\right)^{-1} \boldsymbol{\hat{\Psi}}\_{\boldsymbol{T}} + \left(\boldsymbol{\mathbf{G}} + \frac{\mathbf{F}^{\mathrm{T}}\mathbf{F}}{\varepsilon\_{33}}\right)^{-1} \mathbf{F}^{\mathrm{T}} \mathbf{R} [\boldsymbol{\hat{\Psi}}\_{\boldsymbol{u}}(\boldsymbol{h}) - \boldsymbol{\hat{\Psi}}\_{\boldsymbol{u}}(\boldsymbol{0})].\tag{17}$$

The combination of Equations (6) and (17) gives the state equation of the piezoelectric layer as

$$\frac{d\hat{\mathbf{v}}}{dz} = \mathbf{M}\hat{\mathbf{v}} + \mathbf{N}[\hat{\mathbf{v}}(h) - \hat{\mathbf{v}}(0)],\tag{18}$$

where the 6 × 6 coefficient matrix **M** and the 6 × 6 inhomogeneous term matrix **N** are given by

$$\mathbf{M} = \begin{pmatrix} 0 & \left(\mathbf{G} + \mathbf{F}^{\mathrm{T}} \mathbf{F} / \varepsilon\_{33}\right)^{-1} \\ -\rho\omega^{2}\mathbf{I}\_{3} & 0 \end{pmatrix}, \mathbf{N} = \begin{pmatrix} \left(\mathbf{G} + \mathbf{F}^{\mathrm{T}} \mathbf{F} / \varepsilon\_{33}\right)^{-1} \mathbf{F}^{\mathrm{T}} \mathbf{R} & 0 \\ 0 & 0 \end{pmatrix}. \tag{19}$$

Notice that the state Equation (18) of the piezoelectric layer is naturally degenerated to that of the elastic layer in Equation (10), when the piezoelectric constants are set to zero, i.e., **N** is a matrix of zeros due to **F** = 01×3, with 1 × 3 signifying the numbers of row and column of the zero matrix 0 here and after.

#### *3.2. The Traveling Wave Solution to the State Equation of a Layer*

The state equation of an elastic layer is a set of one order homogeneous linear ordinary differential equations as shown in Equation (10), whose general solution can be directly written according to the mathematical theory of this kind of equations as

$$
\dot{\mathbf{v}} = \Phi \mathbf{e}^{\Lambda z} \mathbf{w}\_{\prime} \tag{20}
$$

where **w** is the 6 order column vector composed of undetermined amplitudes; and **Λ** and **Φ** are the 6 × 6 matrices composed of eigenvalues and eigenvectors of coefficient matrix **M**, respectively. The <sup>6</sup> <sup>×</sup> 6 diagonal matrix <sup>e</sup>**Λ***<sup>z</sup>* is formed by placing <sup>e</sup>γ*iz* on the main diagonal with <sup>γ</sup>*<sup>i</sup>* the *<sup>i</sup>* th (*i* = 1, 2, 3, 4, 5, 6) eigenvalue.

Nevertheless, the state equation of a piezoelectric layer is a set of one order *inhomogeneous* linear ordinary differential equations as shown in Equation (18), whose complete solution is comprised of the general solution to the corresponding homogeneous equations and the particular solution to the inhomogeneous equations and is expressed as

$$
\hat{\mathbf{v}} = \Phi(\mathbf{e}^{\mathbf{A}z} - \mathbf{P})\mathbf{w},
\tag{21}
$$

where **<sup>P</sup>** is the 6 <sup>×</sup> 6 matrix related to matrix **<sup>N</sup>** via **<sup>P</sup>** = **<sup>Λ</sup>**−1**Φ**−1**NΦ**(e**Λ***<sup>h</sup>* <sup>−</sup> **<sup>I</sup>**6) with **<sup>I</sup>**<sup>6</sup> denoting the sixth order identity matrix here and after.

Notice that the solution to the state equation of the elastic layer in Equation (20) can also be achieved through degeneration to that of the piezoelectric layer when the piezoelectric coefficients equal to zero, i.e., **P** is a matrix of zeros due to **N** = 06×6.

#### **4. Transfer Matrix Method**

After the state equations and their solutions determining the state vectors of the elastic and piezoelectric constituent layers have been obtained, the classical transfer matrix method (TMM) [76] is further introduced to derive the equation governing the elastodynamics of the whole unit cell. This equation is represented by a relation between the state vector on the top surface and that on the bottom surface in the TMM formulation. The derived transfer relation of the unit cell will then be combined with the Floquet–Bloch theorem [2,78,79] to bring out the dispersion equation governing the dispersion characteristics of coupled elastic waves in the laminated arbitrarily anisotropic piezoelectric phononic crystals.

#### *4.1. Transfer Relation of an Elastic Layer*

Take a typical elastic layer (*i*) (as depicted in Figure 2a for illustration), the state vectors of its top and bottom surfaces, when referring to Equation (20), should be expressed respectively as

$$
\hat{\mathbf{v}}\_{(i)}(0) = \boldsymbol{\Phi}\_{(i)} \mathbf{w}\_{(i)} \; \hat{\mathbf{v}}\_{(i)}(h\_{(i)}) = \boldsymbol{\Phi}\_{(i)} \mathbf{e}^{\mathbf{A}\_{(i)}h\_{(i)}} \mathbf{w}\_{(i)}.\tag{22}
$$

It is obtained from the former formula in Equation (22) that **w**(*i*) = **Φ**−<sup>1</sup> (*i*) **v**ˆ(*i*)(0), which is substituted into the latter formula in Equation (22) to provide the transfer relation of elastic layer (*i*)

$$
\hat{\mathbf{v}}\_{(i)}(h\_{(i)}) = \boldsymbol{\Phi}\_{(i)} \mathbf{e}^{\mathbf{A}\_{(i)}h\_{(i)}} \boldsymbol{\Phi}\_{(i)}^{-1} \boldsymbol{\Psi}\_{(i)}(0) = \mathbf{B}\_{(i)} \boldsymbol{\Psi}\_{(i)}(0),\tag{23}
$$

where **<sup>B</sup>**(*i*) = **<sup>Φ</sup>**(*i*)e**Λ**(*i*)*h*(*i*)**Φ**−<sup>1</sup> (*i*) is referred to as the transfer matrix of elastic layer (*i*).

#### *4.2. Transfer Relation of a Piezoelectric Layer*

A typical piezoelectric layer (*j*) is accounted for demonstration. By referring to Equation (21) one can express the state vectors at the top and bottom surfaces of layer (*j*) as

$$
\hat{\mathbf{v}}\_{(j)}(0) = \boldsymbol{\Phi}\_{(j)}(\mathbf{I}\_6 - \mathbf{P}\_{(j)}) \mathbf{w}\_{(j)}, \\
\hat{\mathbf{v}}\_{(j)}(h\_{(j)}) = \boldsymbol{\Phi}\_{(j)}(\mathbf{e}^{\mathbf{A}\_{(j)}h\_{(j)}} - \mathbf{P}\_{(j)}) \mathbf{w}\_{(j)}.\tag{24}
$$

When the expression **w**(*j*) = (**I**<sup>6</sup> − **P**(*j*)) −1 **Φ**−<sup>1</sup> (*j*) **v**ˆ(*j*)(0) gotten from the former formula in Equation (24) is introduced into its latter formula, the transfer relation is easily obtained by this process of eliminating **w**(*j*) as

$$
\hat{\mathbf{v}}\_{(j)}(h\_{(j)}) = \boldsymbol{\Phi}\_{(j)}(\mathbf{e}^{\mathbf{A}\_{(j)}h\_{(j)}} - \mathbf{P}\_{(j)}) (\mathbf{I}\_6 - \mathbf{P}\_{(j)})^{-1} \boldsymbol{\Phi}\_{(j)}^{-1} \hat{\mathbf{v}}\_{(j)}(0) = \mathbf{B}\_{(j)} \hat{\mathbf{v}}\_{(j)}(0),\tag{25}
$$

where **<sup>B</sup>**(*j*) <sup>=</sup> **<sup>Φ</sup>**(*j*)(e**Λ**(*j*)*h*(*j*) <sup>−</sup> **<sup>P</sup>**(*j*))(**I**<sup>6</sup> <sup>−</sup> **<sup>P</sup>**(*j*)) −1 **Φ**−<sup>1</sup> (*j*) is referred to as the transfer matrix of piezoelectric layer (*j*).

#### *4.3. Transfer Relation of an Interface*

Consider a typical interface *J* between adjacent layers (*i*) and (*j*). According to the compatibility of displacements and equilibrium of stresses, the state vector at the top surface of layer (*j*) is related to that at the bottom surface of layer (*i*) by

$$
\hat{\mathbf{v}}\_{({\}}(0) = \hat{\mathbf{v}}\_{({\}}(h\_{({\}})) = \mathbf{B}^{\mathsf{J}}\hat{\mathbf{v}}\_{({\}}(h\_{({\}}))}.\tag{26}
$$

where **B***<sup>J</sup>* = **I**<sup>6</sup> is referred to as the transfer matrix of interface *J*.

#### *4.4. Transfer Relation of the Unit Cell*

For the whole unit cell, by recurring the transfer relations of layers and interfaces alternately from bottom to top, the state vector at the bottom surface of the bottommost layer (*m*) can be expressed by the state vector at the top surface of the topmost layer (1) as

$$\begin{array}{lcl}\mathfrak{v}\_{(m)}(h\_{(m)}) &= \mathbf{B}\_{(m)}\hat{\mathfrak{v}}\_{(m)}(0) = \mathbf{B}\_{(m)}\mathbf{B}^{M}\hat{\mathfrak{v}}\_{(m-1)}(h\_{(m-1)}) = \cdots\\ &= \mathbf{B}\_{(m)}\mathbf{B}^{M}\mathbf{B}\_{(m-1)}\mathbf{B}^{(M-1)}\cdots \mathbf{B}\_{(j)}\mathbf{B}^{J}\cdots \mathbf{B}\_{(2)}\mathbf{B}^{2}\mathbf{B}\_{(1)}\hat{\mathfrak{v}}\_{(1)}(0) = \mathbf{B}\hat{\mathfrak{v}}\_{1}(0)\end{array} \tag{27}$$

which is referred to as the transfer relation of the unit cell. The 6 × 6 matrix **B** = **<sup>B</sup>**(*m*)**B***M***B**(*m*−1)**B**(*M*−1) ··· **<sup>B</sup>**(*j*)**B***<sup>J</sup>* ··· **<sup>B</sup>**(2)**B**2**B**(1) is the transfer matrix of the unit cell.

#### *4.5. Dispersion Relation of the Laminated Arbitrarily anisotropic Piezoelectric Phononic Crystals*

Since the unit cells with arbitrarily anisotropic elastic and piezoelectric constituent layers are periodically arranged in the considered laminated phononic crystals, in view of the Floquet–Bloch theorem [2,78,79] for periodic structures, the state vector **v**ˆ(*m*)(*h*(*m*)) should also be related to **v**ˆ(1)(0) through

$$
\hat{\mathbf{v}}\_{(m)}(h\_{(m)}) = \mathbf{e}^{\mathrm{iqM}} \hat{\mathbf{v}}\_{(1)}(0),\tag{28}
$$

where *q* and *H* are the wavenumber of the characteristic coupled elastic waves in and the height of the unit cell, respectively, as can be seen in Figure 1.

Combining Equation (28) with Equation (27) and eliminating the state vector **v**ˆ(*m*)(*h*(*m*)), one gets

$$\mathbf{B}\dot{\mathbf{v}}\_{(1)}(0) = \mathbf{e}^{\mathrm{i}qH}\dot{\mathbf{v}}\_{(1)}(0) \text{ or } (\mathbf{B} - \mathbf{e}^{\mathrm{i}qH}\mathbf{I}\_{6})\dot{\mathbf{v}}\_{(1)}(0) = 0. \tag{29}$$

The former formula in Equation (29) comes down to the matrix eigenvalue problem, while the latter indicates that the determinant of the coefficient matrix (**<sup>B</sup>** <sup>−</sup> <sup>e</sup>i*qH***I**6) should be zero to give nontrivial solution to ˆ**v**1(0), i.e.,

$$\mathbf{e}^{\text{iq}\mathbb{H}} = \text{Eigenvalues}(\mathbf{B}) = \mu \text{ or } \det(\mathbf{B} - \mathbf{e}^{\text{iq}\mathbb{H}}\mathbf{I}\_6) = 0. \tag{30}$$

These formulas are the dispersion equations governing the characteristics of the coupled elastic waves along the thickness direction in the laminated arbitrarily anisotropic piezoelectric phononic crystals, which actually specify the relation between the wavenumber *q* and the frequency ω. Notice that the phase velocity and the wavelength are expressed as *v* = ω/*q* and λ = 2π/*q*, respectively. When the frequency ω is given within a range at specified increment, the frequency-related dispersion curves including the eigenvalue, wavenumber, wavelength and phase velocity spectra, can all be obtained after the former formula in Equation (30) has been solved numerically by direct eigenvalue operation. Otherwise, as anyone among ω, *q* (or λ) and *v* is specified, the other two can be obtained by numerically solving the latter formula in Equation (30), so as to provide the comprehensive dispersion curves

including the frequency-related, wavenumber-related, wavelength-related, and phase velocity-related dispersion curves. Regardless, in the following section, only the frequency-related dispersion curves are illustrated since they can already describe the dispersion characteristics of coupled elastic waves in a relatively complete way.

#### **5. Numerical Examples**

The above derived analysis method will be utilized in this section to calculate the comprehensive frequency-related dispersion curves of elastic waves along thickness direction, including the eigenvalue amplitude spectra, the wavenumber spectra, the wavelength spectra, and the phase velocity spectra in laminated piezoelectric phononic crystals with four exemplified unit cell configurations and with four kinds of electrical boundaries as stated in Section 2. The four unit cell configurations guarantee the presence of four representative patterns about the elastic wave coupling and decoupling, which are composed of the Glass as the inserted elastic layer, the Brass as the electrode layer, and anyone or two layers among the three types piezoelectric layers labeled as "LiNbO3 A", "LiNbO3 B", and "LiNbO3 C", respectively. These three piezoelectric layers are actually all acquired from the piezoelectric material LiNbO3 crystal with *3m* symmetry but are arranged in three different directions. The layers "LiNbO3 A", "LiNbO3 B", and "LiNbO3 C" are formed, respectively, when the threefold symmetry axis *3* are parallel to the local coordinate axes *z*, *y*, and *x*, and correspondingly the three layers are interpreted as *3*//*z*, *3*//*y*, and *3*//*x*, respectively. Because "LiNbO3 B" and "LiNbO3 C" have lower crystal symmetry in the layering plane than the orthotropic crystal discussed by the authors of [55], the coupling of normally-propagating elastic waves exists when they are the constituent layers of the unit cell. By contrast, "LiNbO3 A" layer still has high crystal symmetry in the layering plane. Thus, when the unit cell is Glass-Brass-LiNbO3 A-Brass, which is the first analysis model, all the elastic waves along thickness direction are decoupled. The computed dispersion curves of the pure mode waves in this configuration will be later taken as the reference to show the differences of dispersion characteristics of the coupled mode waves in the other three configurations, i.e., Glass-Brass-LiNbO3 B-Brass, Glass-Brass-LiNbO3 C-Brass, and Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass. In all the four considered configurations, the thicknesses of the inserted elastic Glass layer and the piezoelectric layers are all 10 mm, and the thickness of the Brass electrodes is 0.025 mm, except that it is varied as stated in Section 5.3 where the effect of the electrode thickness on the dispersion curves is discussed. The material parameters of these constituent layers are listed in Table 1. Note that the parameters for "LiNbO3 A" are exactly excerpted from Auld [82], while those for "LiNbO3 B" and "LiNbO3 C" are obtained from further coordinate transformation. The (*c*11)*Glass* = (*c*22)*Glass* = (*c*33)*Glass* = 83.34 GPa and (ρ)*Glass* <sup>=</sup> <sup>2540</sup> kg/m3 for Glass are quoted from Kutsenko et al. [60] and the other parameters are computed from the Poisson's ratio 0.2163 and the Young's modulus 73.39 GPa determined by rough average of the many values in material handbook [84]. The stiffness constants of Brass are all computed from the Poisson's ratio 0.337 and the Young's modulus 106.80 GPa that are roughly determined together with the mass density (ρ)*Brass* <sup>=</sup> <sup>8320</sup> kg/m<sup>3</sup> by referring material handbook [84]. The material parameters of Glass and Brass have also been provided in our previous studies [55,85].

For the convenience of drawing the dispersion curves, the dimensionless frequency ω*H*/(π*v*), the dimensionless wavenumber *qH*/π, the dimensionless wavelength λ/*H*, and the dimensionless phase velocity *v*/*v* are employed with *v* = \$ (*c*55)*LiNbO*<sup>3</sup> **<sup>A</sup>** /(ρ)*LiNbO*<sup>3</sup> the velocity of shear wave in LiNbO3 A.

Since in previous literatures the coupled waves and the electrode thickness are not considered, in order to validate our proposed analysis method, these factors have to be left out. Thus, we further calculated the band structures (frequency–wavenumber spectra) of pure P and S waves along thickness in periodically multilayered Glass-LiNbO3 <sup>A</sup> composites with electric-open boundary, which can be compared with the results computed by analytical formula in Galich et al. [86] or identically in Shen & Cao [87]. All the results are provided in Appendix B.


**Table 1.** Material parameters of the constituent layers in the unit cell of exemplified laminated phononic crystals.

#### *5.1. Dispersion Properties of Coupled Elastic Waves*

To search the general features of the comprehensive frequency-related dispersion curves, the electric-short boundary is adopted in the calculation, because the spectra associated with this electrical boundary serve as the limiting curves for both the applied capacitance and the applied feedback control boundaries as will be shown in Section 5.2. Figures 3–6 give the results of these comprehensive frequency-related dispersion curves for the periodic structures consisting of Glass-Brass-LiNbO3 A-Brass, Glass-Brass-LiNbO3 B-Brass, Glass-Brass-LiNbO3 C-Brass, and Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass configurations, respectively; subfigures (a–d) in these figures represent the relations between the frequency and the eigenvalue amplitude (! ! !μ ! ! !), the wavenumber (*qH*/π), the wavelength (λ/*H*), and the phase velocity (*v*/*v*), respectively. The corresponding wave modes of all the spectra in these figures are marked out.

electric-short boundary: (**a**) eigenvalue amplitude !!!μ!!! spectra; (**b**) wavenumber spectra; (**c**) wavelength spectra; (**d**) phase velocity spectra.

The subfigures in Figures 4–6 indicate, as compared with their counterparts in Figure 3, that the coupled mechanical waves in laminated phononic crystals have some dispersion features identical to those of the uncoupling pure waves, which can be found in [55] but some of which will mention briefly as follows. For example, (1) any kind of frequency-related dispersion curve can show the frequency bands clearly. In frequency ranges where the spectra associated with the real wavenumber (or wavelength/phase velocity) appear, the corresponding wave modes propagate without attenuation, while in frequency ranges where the spectra associated with the imaginary wavenumber (or non-unit eigenvalue amplitude) emerge, the corresponding wave modes attenuate without propagation. These frequency ranges are called as the pass-bands and stop-bands, respectively. (2) When the wavelength tends to infinity, the first order spectrum of each mode asymptotically tend to zero frequency, and the corresponding phase velocity gradually tends to a cut-off value; the higher order spectra of each mode asymptotically tend to the bounding frequencies with phase 2*n*π (*n* is an integer), and the corresponding phase velocity also asymptotically tends to infinity. These features are actually common to all kinds of waves in diverse periodic structures such as the longitudinal waves in periodic elastic rods [88] and periodic piezoelectric rods [85]. In addition, as a complement to the summarization in previous studies, it is also found that the eigenvalue amplitude spectra appear in reciprocal pairs.

Beside these same features, the subfigures in Figures 4–6 also exhibit three dispersion features of coupled elastic waves that are different from those of the pure mode waves. The first remarkable one is that a narrow stop-band can be observed near the intersections of the dispersion curves of different coupled and/or uncoupled waves, which are shown in the green ellipse of Figure 4b, Figure 5b, and Figure 6b. The reason for forming these bands may involve strong coupling and interaction between different wave modes near specific frequency as a result of the material anisotropy. The second feature is that the pass-bands and stop-bands of the coupled wave of one mode no longer appear alternately, which is exactly opposite to those of the pure wave of one mode shown in Figure 3b. The third feature is that for some pass-bands/stop-bands, the frequency lines corresponding to phase *qRH* = 2*n*π and *qRH* = (2*n* + 1)π, i.e., the Brillouin zone boundaries, are no longer the demarcation line between a pass-band and a stop-band. Thus the dispersion curves of these pass-bands/stop-bands are not entirely constrained in-between the Brillouin zone boundaries.

#### *5.2. Tuning the Dispersion Characteristics of Coupled Elastic Waves by the Electrical Boundaries*

In order to show the influence of the electrical boundaries on the dispersion characteristics of coupled elastic waves, the band structures (wavenumber spectra) of the periodic anisotropic piezoelectric composites with the four configurations of the unit cell are computed in cases of four electrical boundaries, including the electric-open, the applied electrical capacitance, the electric-short, and the applied feedback control conditions. The wavenumber spectra are selected from the four kinds of frequency-related dispersion curves because they contain the most information and thus they are the representative dispersion curves. For the sake of expressing the influences more clearly, we always take the electric-open and electric-short boundaries as reference and individually consider the applied positive capacitance, the applied negative capacitance, and the applied feedback control conditions. Figure 7 shows the band structures in cases of electric-open, electric-short, and applied capacitances *<sup>C</sup>*/*<sup>A</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>8</sup> F/m2 and *<sup>C</sup>*/*<sup>A</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>7</sup> F/m2 , with the subfigures (a), (b), (c), and (d) denoting the results for the unit cell configurations of Glass-Brass-LiNbO3 A-Brass, Glass-Brass-LiNbO3 B-Brass, Glass-Brass-LiNbO3 C-Brass, and Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass, respectively. Figure 8a–d provides the corresponding wavenumber spectra in cases of the applied capacitance with the values *<sup>C</sup>*/*<sup>A</sup>* <sup>=</sup> <sup>−</sup>2.5 <sup>×</sup> <sup>10</sup>−<sup>8</sup> F/m2 and *<sup>C</sup>*/*<sup>A</sup>* <sup>=</sup> <sup>−</sup>2.6 <sup>×</sup> <sup>10</sup>−<sup>8</sup> F/m2 for periodic Glass-Brass-LiNbO3 A-Brass composite as well as*C*/*<sup>A</sup>* <sup>=</sup> <sup>−</sup>3.8×10−<sup>8</sup> F/m2 and*C*/*<sup>A</sup>* <sup>=</sup> <sup>−</sup>4.0×10−<sup>8</sup> F/m2 for the other configurations, together with the electric-open (*C*/*A* = 0) and electric-short (*C*/*A* = ∞) condition, while Figure 9a–d provides the corresponding ones in the case of applied feedback control with the gain coefficient being *Kg* = 5.0 <sup>×</sup> <sup>10</sup><sup>10</sup> V/m and *Kg* = 1.0 <sup>×</sup> <sup>10</sup><sup>12</sup> V/m together with the electric-open and electric-short (*Kg* = 0.0 V/m) conditions. Note that the values of the positive capacitance and the gain coefficient are

chosen so as to most clearly show their influences on the band structures and their correlation with the electric-open and electric-short conditions. The negative capacitance values are determined so as to closely near but not at the singular point specified by ε<sup>2</sup> <sup>33</sup>*S*/*C* + ε33*h* = 0 that can be obtained through the vanishing of the denominator of the expression for vector **R** as given in Table A1. In addition, notice also from the expression for vector **R** in Table A1 that *C*/*A* = +∞ and *C*/*A* = −∞ actually lead to the same **R** and accordingly the same computed band structures. In Figure 10, the band structures associated with the negative capacitances are compared with those associated with a very big gain coefficient *Kg* = 1.0 <sup>×</sup> 1012 V/m, so that the correlation between applied negative capacitance and applied feedback control can be revealed. In Figures 7–10, the band structures corresponding to the pure wave that is independent on the electric field are marked by solid arrows pointing to the direction of increasing frequency, these are trivial curves for our investigation. In the following interpretations, we only focus on the elastic waves that are dependent on electric field and thus the electrical boundary. Each unit cell configuration exactly corresponds to one elastic wave. As also noticed in Section 5.1, these waves are pure P, coupled P-SV, coupled SV-SH, and coupled P-SV-SH modes in subfigures (a–d), respectively. Common in this numerical example Section, the pure P wave in the periodic Glass-Brass-LiNbO3 A-Brass structure is used as the object of comparison, so that the differences between the coupled waves and the pure modes waves can be clearly exhibited.

From all the subfigures in Figure 7, it is seen that as the positive capacitance varies from *C*/*A* = 0 to *C*/*A* = ∞, the band structures gradually changes from those associated with the electric-open condition to those associated with the electric-short condition. Generally, the central frequency of bandgaps moves downward, and the bandgap widths becomes narrower. Nevertheless, the entire alteration is not very significant, since the band structures associated with the electric-short condition are actually relatively near the corresponding curves associated with the electric-open condition. Comparing Figure 7b–d with Figure 7a, it is seen that the dispersion curves of coupled waves are more obviously influenced than those of the pure wave, so that the coupled waves are more sensitive to the applied positive capacitance.

From all the subfigures in Figure 8, it is seen that for any waves related to the electric field, either uncoupled or coupled modes, sharp variation of the attenuation constants (imaginary wavenumbers) to infinity exists at some frequency ranges, where the phase constants also sharply and corresponding changed. These phenomena are similarly found in Refs. [55,58] for pure P and S waves. As the negative capacitance varies from *C*/*A* = 0 to *C*/*A* = −∞, the band structures does not change in a monotonous way, which is attributed to the singular property as the applied capacitance corresponds to the zero denominator of the expression for vector **R**. The comparison between Figure 8a–d further indicates that the coupled P-SV-SH wave has more complex variation pattern of the attenuation and phase constants spectra near the singular range.

It is seen from Figure 9 that the gain coefficient of the applied feedback control has very obvious influence on the band structures of both uncoupled and coupled waves except those of the coupled SV-SH wave. The insensitivity of coupled SV-SH wave to the applied feedback control condition is resulting from the *e*<sup>33</sup> = 0 for LiNbO3 <sup>C</sup> as shown in Table 1, which further leads to the vanishing of the last component of vector **F** and accordingly the vanishing of the components at the last row of **F**T**R**. As the gain coefficient *Kg* increases from zero, the band structures of coupled P-SV and P-SV-SH waves associated with the applied feedback control condition, as compared to those associated with the electric-short condition, changes substantially including forming new phase/attenuation constants spectra as shown in subfigures (b,d). However, those of pure P wave simply move downwards to the smaller frequency side.

Figure 10 indicates, except the coupled SV-SH wave, all other waves can realize the singular pattern of dispersion by either applying proper negative capacitance boundary or applying feedback control condition with relatively big gain coefficient. For coupled SV-SH wave in periodic Glass-Brass-LiNbO3 C-Brass composites discussed here, only the proper negative capacitance can be applied.

Moreover, from all the above discussion it is known that the electric-short boundary is both the upper value limiting of the applied positive/negative capacitance boundary and the lower value limiting of the applied feedback control boundary. This is the reason that the electric-short boundaries are utilized in the previous Section 5.1 and the following Section 5.3.

**Figure 7.** The influence of the applied capacitance on the band structures of the uncoupled and coupled elastic waves in periodic anisotropic piezoelectric composites with the unit cell consisting of (**a**) Glass-Brass-LiNbO3 A-Brass, (**b**) Glass-Brass-LiNbO3 B-Brass, (**c**) Glass-Brass-LiNbO3 C-Brass, and (**d**) Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass.

**Figure 8.** The variation of the band structures of uncoupled and coupled elastic waves when the applied negative capacitance is at the vicinity of the singular point in periodic anisotropic piezoelectric composites with the unit cell consisting of (**a**) Glass-Brass-LiNbO3 A-Brass, (**b**) Glass-Brass-LiNbO3 B-Brass, (**c**) Glass-Brass-LiNbO3 C-Brass, and (**d**) Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass.

**Figure 9.** The influence of the gain coefficient in the applied feedback control condition on the band structures of the uncoupled and coupled elastic waves in periodic anisotropic piezoelectric composites with the unit cell consisting of (**a**) Glass-Brass-LiNbO3 A-Brass, (**b**) Glass-Brass-LiNbO3 B-Brass, (**c**) Glass-Brass-LiNbO3 C-Brass, and (**d**) Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass.

**Figure 10.** The comparison of the band structures associated with the applied feedback control condition of big gain coefficient (*Kg* = 1.0 <sup>×</sup> <sup>10</sup><sup>12</sup> V/m) and those associated with the applied negative capacitance near the singular point in periodic anisotropic piezoelectric composites as the unit cell consisting of (**a**) Glass-Brass-LiNbO3 A-Brass, (**b**) Glass-Brass-LiNbO3 B-Brass, (**c**) Glass-Brass-LiNbO3 C-Brass, and (**d**) Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass.

#### *5.3. Electrode Thickness Influence on the Dispersion Characteristics of Coupled Elastic Waves*

For the sake of knowing the mechanical effect of electrode layers on coupled elastic waves besides their electrical function, the band structures (wavenumber spectra) of the periodic anisotropic piezoelectric composites with the above-mentioned four configurations of unit cell are computed by specifying the electric-short boundary and four values of the electrode thickness, i.e., 0 mm, 0.025 mm, 0.25 mm, and 2.5 mm. The results are provided in Figure 11, with subfigures (a–d) corresponding to the configuration of Glass-Brass-LiNbO3 A-Brass, Glass-Brass-LiNbO3 B-Brass, Glass-Brass-LiNbO3 C-Brass, and Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass, respectively. Adopting the electric-short boundary is because the associated band structures serve exactly as the limiting curves for both the applied capacitance and the applied feedback control boundaries as obtained from Section 5.2, and thus they are typical. In Figure 11a–d, for any kind of wave, when the ratio of the electrode thickness to the thickness of the piezoelectric layer (called as the relative electrode thickness) is less than 0.25%(= 0.025/10), the mechanical effect of the electrode is too weak to affect the band structures. When it reaches to 2.5%(= 0.25/10), the band structures of the relatively high order modes are distinguishably affected by the mechanical effect of the electrode, though the band structures of low order modes are nearly not influenced. Also inferred from these figures that, when the relative electrode thickness is higher than 2.5%, with the increasing of electrode thickness, the band structures will be more and more obviously affected by the mechanical effect of the electrode. The variation is characterized by the decreasing of the central frequencies of the pass-bands/stop-bands and the changes of the bandgap widths without uniform rule. From the comparisons of Figure 11b–d with Figure 11a, it can be noticed that the band structures of coupled wave modes associated with the narrow bandgaps discussed in Section 5.1 are more sensitive to the electrode thickness than those of the other coupled wave modes and uncoupled waves. With the increasing of electrode thickness, these sensitively-influenced bandgaps gradually move down, but their widths remain almost unchanged.

**Figure 11.** The influence of electrode thickness on the band structures of the uncoupled and coupled elastic waves in periodic anisotropic piezoelectric composites with electric-short boundary and with the unit cell consisting of (**a**) Glass-Brass-LiNbO3 A-Brass, (**b**) Glass-Brass-LiNbO3 B-Brass, (**c**) Glass-Brass-LiNbO3 C-Brass, and (**d**) Glass-Brass-LiNbO3 B-Brass-LiNbO3 C-Brass.

#### **6. Conclusions**

This paper studies the active tuning of coupled elastic waves propagating perpendicular to the interface in general laminated piezoelectric phononic crystals composed of arbitrarily anisotropic piezoelectric and elastic materials, where the wave coupling is attributed only to the material anisotropy. The active tuning is realized through the continuous adjustment of the electrical parameters in the external electric circuit connected to the piezoelectric layers or the switching between different electric circuits. The considered electrical boundaries include the electric-open, the electric-short, the applied capacitance, and the applied feedback control conditions. The analytical matrix method is proposed,

on the bases of the state space formalism for modeling the elastodynamics of constituent layers, the transfer matrix method for modeling the relations between constituent layers, and the Floquet–Bloch theorem for modeling the periodic conditions of the unit cell, to derive the dispersion equation of general periodic laminated piezoelectric composites. By validating and adopting the analysis method in numerical examples, the dispersion characteristics of coupled elastic waves as well as the influences of electrical boundaries and electrode thickness on the band structures of coupled waves were investigated. The following conclusions may be drawn from these studies for the coupled waves induced by only material anisotropy other than inclined propagation.


All these findings will push forward the design and nondestructive evaluation of laminated anisotropic piezoelectric phononic crystals that may be applied in the smart structures with active wave control function.

**Author Contributions:** Conceptualization, Y.Q.G.; methodology, Y.Q.G.; software, Y.Q.G. and Q.Q.L.; validation, Q.Q.L., Y.J.W., and H.B.Z.; formal analysis, Q.Q.L., Y.J.W., and H.B.Z.; investigation, Y.Q.G.; resources, Y.Q.G.; data curation, Q.Q.L.; writing—original draft preparation, Q.Q.L., Y.J.W., and H.B.Z.; writing—review and editing, Y.Q.G. and Q.Q.L.; visualization, Q.Q.L.; supervision, Y.Q.G.; project administration, Y.Q.G. and Q.Q.L.; funding acquisition, Y.Q.G. and Q.Q.L.

**Funding:** This work was funded by the National Natural Science Foundation of China, grant number 11372119, and the Youth Innovation Funds in Shanxi Agricultural University, grant number 2018003.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**

In this section, the derivation of the relation between dφˆ/d*z* and d ˆ**v***u*/d*z* will be given in detail, based on all the electrical considerations including the electrical boundary conditions and their influences on mechanical field.

First, substituting Equation (5) into the charge equation (Gauss's law) of electrostatics with free charge as the first formula in Equation (3), one gets d*D*ˆ *<sup>z</sup>*/d*z* = 0, which implies that *D*ˆ *<sup>z</sup>* is a constant. The relation between *D*ˆ *<sup>z</sup>* and the electric charge on electrodes *Q*ˆ is *D*ˆ *<sup>z</sup>* = *Q*ˆ /*A* with *A* the area of electrodes, which further indicates that *Q*ˆ is a constant. It is then substituted into the latter formula of Equation (14) to express dφˆ/d*z* as

$$\frac{d\hat{\phi}}{dz} = \frac{\mathbf{F}}{\varepsilon\_{33}} \frac{d\psi\_u}{dz} - \frac{\hat{Q}}{\varepsilon\_{33}A} \tag{A1}$$

Second, the voltage (potential difference) *V*ˆ between the electrodes can be computed from the electric potential by

$$\begin{split} \hat{V} \quad = \hat{\phi}(h) - \hat{\phi}(0) = \int\_{0}^{h} \frac{d\hat{\phi}}{dz} dz = \int\_{0}^{h} \left( \frac{\mathbf{F}}{\varepsilon\_{33}} \frac{d\psi\_{\mathbf{z}}}{dz} - \frac{\hat{Q}}{\varepsilon\_{33}\mathbf{A}} \right) \mathrm{d}z \\ = \frac{\mathbf{F}}{\varepsilon\_{33}} [\hat{\Psi}\_{\mathrm{il}}(h) - \hat{\Psi}\_{\mathrm{il}}(0)] - \frac{\hat{Q}h}{\varepsilon\_{33}\mathbf{A}} \end{split} \tag{A2}$$

where Equation (A1) has been involved. Equation (A2) provides the relation with respect to three quantities, i.e., the voltage *<sup>V</sup>*<sup>ˆ</sup> , the electric charge *<sup>Q</sup>*ˆ, and the displacement difference [**v**<sup>ˆ</sup> *<sup>u</sup>*(*h*) <sup>−</sup> **<sup>v</sup>**<sup>ˆ</sup> *<sup>u</sup>*(0)]. On the basis of Equation (A2) and the mathematical relation specifying either *V*ˆ or *Q*ˆ or a relation between *V*ˆ and *Q*ˆ corresponding to the four considered electrical boundaries, the expressions of *V*ˆ and *Q*ˆ can be further obtained as listed in Table A1.


**Table A1.** Expressions of relative electrical quantities in cases of four electrical boundaries.

Finally, substituting the expressions of *Q*ˆ into Equation (A1) gives,

$$\frac{d\boldsymbol{\hat{\phi}}}{d\boldsymbol{z}} = \frac{\mathbf{F}}{\boldsymbol{\varepsilon}\_{33}} \frac{d\boldsymbol{\Psi}\_{\boldsymbol{\mu}}}{d\boldsymbol{z}} - \mathbf{R}[\boldsymbol{\Psi}\_{\boldsymbol{\mu}}(\boldsymbol{h}) - \boldsymbol{\Psi}\_{\boldsymbol{\mu}}(\boldsymbol{0})] \tag{A3}$$

i.e., Equation (16), where **R** is naturally obtained from *Q*ˆ and represented by **F** in cases of four electrical boundaries as also listed in Table A1.

#### **Appendix B**

In this section, the frequency–wavenumber spectra (band structures) of elastic waves along thickness in laminated piezoelectric phononic crystals with the unit cell consisting of Glass-LiNbO3 A layers and the piezoelectric layer being in a state of electric open are analyzed by the formulation described above, so that the proposed method and the corresponding computer program can be validated in this degeneration situation. The elastic waves in this case are decoupled into pure SH, SV and P waves, so the computed results can be compared with their counterparts obtained by the explicit formula in Galich et al. [86] or identically in Shen & Cao [87], as provided in Figure A1. It is worth noting that the equivalent elastic constant *c*<sup>33</sup> + *e*<sup>2</sup> <sup>33</sup>/ε<sup>33</sup> is used for the LiNbO3 <sup>A</sup> layer in the explicit formula so as to take its piezoelectric effect into account. From Figure A1, it is seen that the band structures gained by the two methods are exactly coincident, which proves our proposed method in this paper. Besides, note that the wavenumber spectra of the SH wave are superposed exactly on those of the SV wave due to the same wave speed.

**Figure A1.** Comparisons of the wavenumber spectra for the periodic Glass-LiNbO3 <sup>A</sup> multilayers with electric-open boundary computed by the proposed method and the explicit dispersion relation from Galich et al. [86].

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An Optimization of Two-Dimensional Photonic Crystals at Low Refractive Index Material**

#### **Thanh-Phuong Nguyen 1, Tran Quoc Tien <sup>2</sup> and Quang Cong Tong 2,\* and Ngoc Diep Lai 3,\***


Received: 1 August 2019; Accepted: 21 August 2019; Published: 24 August 2019

**Abstract:** Photonic crystal (PC) is usually realized in materials with high refractive indices contrast to achieve a photonic bandgap (PBG). In this work, we demonstrated an optimization of two-dimensional PCs using a low refractive index polymer material. An original idea of assembly of polymeric multiple rings in a hexagonal configuration allowed us to obtain a circular-like structure with higher symmetry, resulting in a larger PBG at a low refractive index of 1.6. The optical properties of such newly proposed structure are numerically calculated by using finite-difference time-domain (FDTD) method. The proposed structures were realized experimentally by using a direct laser writing technique based on low one-photon absorption method.

**Keywords:** photonic crystals; photonic bandgaps; polymer materials; direct laser writing

#### **1. Introduction**

Photonic crystal (PC), an artificial material in which the refractive index is modulated at wavelength scale, offers presently many interesting applications in different domains [1,2]. The most important property of the PCs is the existence of a so-called photonic bandgap (PBG). In case of one-dimensional (1D) PCs, the PBG can be easily obtained by assembling multiple layers of two different materials with very low refractive indices contrast [3–6]. However, in case of two- and three-dimensional (2D and 3D) PCs, it requires that the refractive indices contrast should be as high as possible. Some 2D and 3D PCs are made by polystyrenes or SiO2 nanoparticles [7–10], which are fabricated by a simple fabrication method, but these PCs do not have a true PBG, or in other word, they have a so-called pseudo or partial PBG, due to their weak refractive index. Therefore, most 2D and 3D PCs are made by semiconductor materials with high refractive indices, which however makes the fabrication procedure complicated and expensive. Different optimizations have been proposed to obtain the PBG as large as possible and with a material with a smaller refractive index. By increasing the symmetry of the PCs, for example in the case of 2D PCs, from square (four-fold) to hexagon (six-fold), and to quasi-periodic PCs, such as Penrose structure [11,12], the PBG can be opened at lower refractive index materials, such as polymers [13,14].

The most difficult task of PCs investigations is the fabrication, which should satisfy several conditions, such as simple, rapid, inexpensive, and particularly flexible. Today, many fabrication techniques have been proposed to realize PCs, for example, e-beam lithography, self-assembly of opals, laser interference, and direct laser writing. The choice of fabrication technique depends on desired applications. Among available fabrication technologies, the direct laser writing (DLW) is the most powerful one allowing fabrication of any desired PCs, in multi-dimension ((1D, 2D, 3D) PCs) [15,16]. The working principle of the DLW is to focus tightly a laser beam into a small region of a photoresist to induce a local polymerization or depolymerization effect. By moving the focusing spot following a well-designed structure, a corresponding polymeric pattern can be obtained.

In this work, we demonstrated an optimization of the 2D PCs to obtain a large PBG even with low refractive index materials, such as polymer. We proposed an original idea of high symmetry 2D PCs, and a new experimental way to realize these PCs by the DLW technique. The PBGs of a standard honeycomb structure and of a newly proposed structure are calculated and compared by using a Finite-difference time-domain (FDTD) method. These structures are then fabricated by a new scanning method using the one-photon absorption-based DLW technique [17].

#### **2. Theoretical Calculation**

For calculation of the PBG of different PCs, we have used the FDTD method (Lumerical Software) with Bloch boundary conditions. We assumed that the refractive index of material is 1.6, a typical value of polymer without considering the dispersion effect of the material. This choice is justified by the values experimentally measured in our case, i.e.; the refractive index of the used polymer is in between 1.55 and 1.62 for the wavelengths ranging from 400 nm to 800 nm. For simulations of 2D PCs, we distinguished two different modes, TM and TE, which are the transverse magnetic or transverse electric modes where the magnetic or electric field is parallel to the PC plane, respectively.

It is well-known that the honeycomb lattice is the best 2D PC, which provides a largest PBG, as compared to the square or hexagonal lattice. We first calculated the PBG of this structure at low refractive index material, by considering two types: one is the honeycomb lattice comprising of polymeric cylinders (*n* = 1.6) in an air background (*n* = 1) and the another is an inverse geometry consisting of air-holes in the polymeric background. Figure 1a,b show a honeycomb lattice and its reciprocal space, where a is lattice constant, and r is the radius of polymeric cylinders or air-holes. For the case of polymeric cylinders in the air background, the PBG only reveals a small gap for TE mode. Figure 2a shows a band structure for TE mode, obtained with an optimum radius of r/a = 0.3. Figure 2b shows gap map obtained with various filling factors (different ratios of r/a). We see that the largest PBG for TE mode is about 0.02 at r = 0.3a. When considering a honeycomb lattice made of air cylinders in the polymeric background, we found that the PBG only reveals for TM mode. A band structure for TM mode (r/a = 0.35) is shown in Figure 3a. The PBG (for TM mode) as a function of the r/a ratio was also calculated, and a gap map is presented in the Figure 3b. We can see that the PBG is increased as compared to the one obtained in previous case, and reaches 0.05 at r = 0.35a. These results are somehow similar to results obtained with higher refractive index [1,2]. But at low refractive index material, the PBG exists but not complete, i.e.; no common PBG for both TE and TM modes. Table 1 summarizes the PBG of the honeycomb lattice at low refractive index material for both configurations.

**Figure 1.** (**a**) A two-dimensional honeycomb lattice with two main lattice vectors: **<sup>a</sup>**<sup>1</sup> = [√3a/2, a/2], **<sup>a</sup>**<sup>2</sup> = [√3a/2, <sup>−</sup> a/2]. The artificial atome has a cylindrical form with a radius r. (**b**) The reciprocal space of corresponding two-dimensional honeycomb lattice with the primitive vectors: **<sup>b</sup>**<sup>1</sup> = 2*π*/a (1/3√3, 1/3) and **<sup>b</sup>**<sup>2</sup> = 2*π*/a (1/3√3, <sup>−</sup>1)/3). The reduced first Brillouin zone is identified by <sup>Γ</sup>, K and M points.

**Figure 2.** Photonic band structureof a 2D honeycomb lattice of polymeric cylinders in an air background (see insert of b). (**a**) Photonic bandgap for TE mode with r/a = 0.3. (**b**) Map of photonic bandgap as a function of r/a (shaded blue region).

**Figure 3.** Photonic band structure of a 2D honeycomb lattice of air cylinders in a polymeric background (see insert of b). (**a**) Photonic bandgap for TM mode with r/a =0.35. (**b**) Map of photonic bandgap as a function of r/a (shaded blue region).

**Table 1.** Summary of the best photonic bandgaps of honeycomb lattices obtained with a polymer material with a refractive index of *n* = 1.6.


It has been well-known that the larger PBG can be obtained by increasing the structure symmetry [18]. Different structures have been proposed [11–14], and larger PBG are obtained, even at low refractive index materials. In particular, it is important to emphasis on the most symmetric one, called circular photonic crystal (CPC) [19,20], which possesses a perfect symmetry in a 2D plane. This CPC structure enables particularly an isotropic PBG, and consequently a PBG at low refractive index. However, the CPC is a non-periodic structure and usually developed around a symmetric point. Therefore, its photonic property is only characterized by a transmission calculation or measurement but not its PBG.

In this work, we proposed a novel kind of equivalent 2D PC by arranging multi-rings in a periodic configuration, such as a hexagonal lattice. This novel structure allowed us to combine the advantages of both CPC structure, thanks to the perfect symmetry of the rings (dielectric or air), and the periodic organization of these rings in a well-known hexagonal lattice, which can be numerically simulated to obtain a true PBG. As it will be shown in the experimental part, the proposed structure is also matched to the capacity of the fabrication technology.

Figure 4a presents an assembled multi-rings lattice and Figure 4b illustrates its corresponding reciprocal space, where a is the lattice constant, r is the radius of the rings, and d is the width of rings. Similar to honeycomb PC, the PBG of the assembled multi-rings lattice can be calculated for two possible configurations: the lattice containing multi-rings of air (*n* = 1) in the polymeric background (*n* = 1.6) and the another is an inverse geometry consisting of multi-rings of polymer in the air background. In practices, these two configurations can be optically realized by using a negative photoresist or a positive photoresist, or simply by using a positive photoresist with controllable laser power [17]. We demonstrated that the multi-rings assembly allowed obtaining different 2D PCs with various filling factors. That can be done by changing the r/a ratio. It is clear that if a > 2r, all rings will be separated, and the structure is a simple hexagonal PC with a ring as an "artificial atom". However, if a < 2r, multiple rings will overlap to each other, resulting in a complicated 2D PC, with a perfect symmetry determined by the ring.

**Figure 4.** (**a**) Proposal of a 2D assembled-multirings structure arranged in a hexagonal configuration. The main lattice vectors are **<sup>a</sup>**<sup>1</sup> = a 1/2, <sup>√</sup>3/2 and **<sup>a</sup>**<sup>2</sup> = a 1/2, − <sup>√</sup>3/2 ; r is the radius of rings; d is the width of ring-walls; a is the lattice constant. (**b**) The corresponding 2D reciprocal space off the 2D assembled-multirings structure. **<sup>b</sup>**<sup>1</sup> <sup>=</sup> <sup>2</sup>*π*/a 1, 1/√<sup>3</sup> and **<sup>b</sup>**<sup>2</sup> <sup>=</sup> <sup>2</sup>*π*/a 1, <sup>−</sup>1/√<sup>3</sup> : the primitive vectors of the reciprocal lattice. The Brillouin zone is identified by a yellow hexagon.

Figure 5 shows the simulation results for the case of assembled polymeric multi-rings (*n* = 1.6) in the air background. Similar to the case of the honeycomb, it has seen that the PBG reveals only for TE mode. For the simulation shown in Figure 5a, we assumed that the ring width is 300 nm, the ring separation a = 1.5 μm, and the ring radius r/a = 0.51, which is suitable to the fabrication capacity. A larger PBG is obtained, as compared to similar case of the honeycomb lattice shown in Figure 2a. By changing the radius of rings, a gap map (shaded blue region) as a function of r/a for TE mode has been carried out, and shown in Figure 5b. It can be seen that the largest PBG is obtained with a value of aΔ*ω*/2*π*c = 0.035 with r = 0.51a.

For the case of a 2D lattice consisting of air rings (*n* = 1) in the polymeric background, in contrast, the PBG reveals only for TM mode. Figure 6a shows its PBG obtained with a radius ratio of r/a = 0.78 with a = 1.5 μm. Again, the PBG is larger, as compared to that obtained by the honeycomb lattice shown in Figure 3a. A map of PBGs for TM mode versus the ratio of r/a is shown in Figure 6b. It can be seen that multiple gaps can be obtained with this structure, even with low refractive index contrast. The PBG can reach a value of aΔ*ω*/2*π*c = 0.089 at r = 0.78a, which is larger than 0.05, obtained by the honeycomb PC. A summary of PBGs of multi-rings lattice is shown in Table 2. This novel PC still does not produce a complete PBG, but the individual PBG (for TE or TM mode) is better than those obtained by a best standard 2D PC (e.g., honeycomb).

**Figure 5.** (**a**) Photonic band structure (TE mode) of an assembled-multirings lattice consisted of polymer-rings in an air background. (**b**) Map of photonic bandgap (TE mode) as a function of r/a (shaded blue region). The simulations are realized with r/a = 0.51, a = 1.5 μm, d = 300 nm, and *n* = 1.6.

**Figure 6.** (**a**) Photonic band structure (TM mode) of an assembled-multirings lattice consisted of air rings in the polymer background. (**b**) Map of photonic bandgap (TM mode) as a function of r/a (shaded blue region). The simulations are realized with r/a = 0.78, a = 1.5 μm, d = 300 nm, and *n* = 1.6.

**Table 2.** Summary of photonic bandgaps of assembled multi-rings structures. The calculations were realized with a ring width of d = 300 nm and the polymer material has a refractive index of *n* = 1.6.


We note that in a standard case of square or hexagonal 2D PCs, the rod-type has a TM-bandgap and the hole-type has a TE-bandgap, respectively [1,2]. However, when working with particular structures, such as honeycomb or multi-ring PCs, the situation changes. In these cases, the existence of TM or TE bandgaps varies as a function of the ratio between the rod size and the lattice periodicity. With low refractive index, we can see that the rod-type favors a TE-bandgap while the hole-type favors a TM-bandgap, respectively.

#### **3. Fabrication of Proposed Structure by Direct Laser Writing Method**

To fabricate these proposed PCs, two possible methods can be used, such as e-beam lithography and DLW method. In this work, we have demonstrated the fabrication by using one-photon absorption-based DLW method, which is simple, low cost, and it can even produce 3D PCs [16,17,21]. The experimental setup is shown in Figure 7 and can be explained as follows. A simple and low-cost continuous-wave (cw) 532 nm laser was used to write the desired structures. The laser power was adjusted by using a combination of a half-wave plate and a polarizing beam splitter. The laser beam is

sent through an electronic shutter (S), which allows controlling the exposure time. The quarter-wave plate was used to obtain a circular polarization of the excitation beam, just before focusing. The laser beam was sent on a 50/50 beam-splitting cube (BS), after which one part goes to a powermeter and the other is directed to an oil-immersion objective lens (OL) with a numerical aperture (NA) of 1.3, which focuses the laser beam on the photoresist sample. The last one was placed on a piezoelectric translation stage (PZT), which allows moving the sample in 3D with a resolution of 0.1 nm. Thanks to the high NA OL, the laser is tightly focused onto the sample with a focusing spot of about 300 nm. To determine the focusing position, the same OL is used to collect the fluorescence photons, emitted from the photoresist sample, which then propagate in opposite direction with respect to the excitation beam. The fluorescence beam is passed through a long-pass filter (F) with a cut-off wavelength of 580 nm, to remove the excitation laser light, and directed towards an avalanche photodiode (APD). It has been demonstrated that this technique can fabricate structures with a thickness as large as 100 micrometers, depending on the absorption coefficient of the photoresist [16,21], with a periodicity down to 400 nm and the smallest size of pattern can be as small as 57 nm [17].

**Figure 7.** Experimental setup of the low one-photon absorption-based direct laser writing method used to fabricate desired polymeric photonic crystals. HWP: half-wave plate, PBS: polarizing beam splitter, S: electronic shutter, M: mirror, QWP: quarter-wave plate, BS: beam splitter, OL: objective lens, PZT: piezoelectric translator, F: long-pass filter, L1,2: lenses, PH: pinhole, APD: avalanche photodiode.

In this work, we proposed a new way to realize desired structures. Indeed, in e-beam lithography or DLW, the structures are usually created by moving the focusing spot following a design. We proposed to assembly multiple lines or rings to achieve desired structures. Various structures can be created with the same lines or rings assembly, which simplifies a lot for fabrication procedure. For demonstration, we used a positive S1818 photoresist with a thickness of 1 μm. With low excitation laser power, the focusing spot produced an air-hole in polymer thin film (after development), and by moving this focusing, we obtained an assembly of air multi-rings lattice. Please note that with the same photoresist and with higher excitation laser power, we can obtain a polymeric multi-ring thanks to the optically induced thermal effect [17].

Figure 8a shows a new proposed trajectory of the focusing laser beam to draw a honeycomb lattice. By this method, we can obtain different structures, starting from a standard air-hole hexagonal structure to a honeycomb structure with controllable filling factors and dots shape. Figure 8b–e show scanning electron microscope (SEM) images of different structures (period of a = 1.5 μm) fabricated by a constant laser power of 10 μW and by different moving velocities from 9 μm/s to 6 μm/s. We can see that with velocities of 9 μm/s and 8 μm/s, typical hexagons with air-holes are obtained, while with velocities of 7 μm/s and 6 μm/s (higher dose), honeycombs with different filling factors are demonstrated. Thus, by using the simple scanning scheme shown in Figure 8a, and by controlling the scanning speed, various types of honeycomb lattices can be obtained.

**Figure 8.** Fabrication of a honeycomb lattice by a simple direct laser writing method. (**a**) Illustration of a writing model in which the focused laser beam is scanned line by line with a period a, following a triangular configuration. The line width, d, is controlled by the writing parameters, such as laser power, writing speed, and focusing lens. The final structures vary from air-holes hexagonal structure to polymeric-dots honeycomb, depending on the writing doses. (**b**–**e**) SEM images of fabricated 2D structures (period, a = 1.5 μm), obtained with different doses. For these structures, the laser power was fixed at 10 μW and the scanning speeds were 9 μm/s (**b**), 8 μm/s (**c**), 7 μm/s (**d**), and 6 μm/s (**e**), respectively.

By the same method, assembled multi-rings lattices were fabricated by scanning the focusing laser beam following the trajectory shown in Figure 9a. Different structures have been demonstrated by changing the structures parameters, such as the ring width, d, the lattice constant, a, as well as the ring radius, r. Figure 9b–e show SEM images of four structures realized by keeping the same laser power, 10 μW while changing the moving velocities of laser beam from 9 μm/s to 6 μm/s. The lattice constant and the ring radius were chosen as 2 μm and 1.51 μm, respectively. When varying the scanning velocities, the ring width changes from 150 nm (Figure 9b), 200 nm (Figure 9c), 300 nm (Figure 9d), to 400 nm (Figure 9e), resulting in different 2D structures with different filling factors and structure compositions. This perfectly agrees with the theoretical prediction, and their PBGs are also demonstrated to be better than those obtained by the standard 2D PCs.

**Figure 9.** Fabrication of various high symmetry 2D photonic crystals by assembling multiple rings. (**a**) Illustration of the assembly of multi-rings, created by moving the focused laser beam, following a hexagonal configuration. The ring width, d, is controlled by the writing dose, and the rings are separated, from center to center, by a lattice constant a. By changing the exposure doses (laser power and writing speed) as well as the lattice constant, a, the final structures can be varied from a Penrose-like structure (**b**) to a hexagonal structure (**e**) consisting of different polymeric dots. (**b**–**e**) SEM images of 2D assembled-multi-rings lattices with a lattice constant a = 2 μm and fabricated by the same laser power of 10 μW. Different structures were obtained, as predicted by the simulations, with different ring widths d = 150 nm (**b**), 200 nm (**c**), 300 nm (**d**), and 400 nm (**e**), realized by different scanning speeds, 9 μm/s (**b**), 8 μm/s (**c**), 7 μm/s (**d**), and 6 μm/s (**e**), respectively.

Please note that it is not possible at this moment to experimentally characterize the PBGs of these fabricated structures, because the structures sizes are limited to 100 μm × 100 μm due to the scanning system of the DLW technique. Also, with a lattice constant of 2 μm, an assembled multi-rings PC produces a PBG in the range of 2.85–3.15 μm, which is in the mid-infrared domain. The limited PC thickness (1 μm) also makes the direct PBG characterization difficult. Further reduction of the lattice constant (a) and of the ring radius (r), and increase of the PCs thickness are therefore necessary to obtain a PBG in visible or near IR range, which can be very promising for many applications. A direct application of such fabricated structures is to make PC-based organic laser, with controllable emission wavelength and desired beam shape.

#### **4. Conclusions**

In conclusion, we have demonstrated an optimized two-dimensional photonic crystal with a better photonic bandgap at low refractive index contrast. The numerical calculations of the photonic bandgaps of honeycombs and assembled multi-rings lattices at low refractive index (*n* = 1.6) have been done by using FDTD method. The results showed that the TM mode has larger PBGs, as compared to the TE mode, and the largest PBG is approximately 0.089 (= aΔ*ω*/2*π*c) at r = 0.78a, obtained by the newly proposed assembled multi-rings photonic crystal. This is due to the perfect symmetry at local area of the photonic crystal, and the PBG calculation is possible thanks to the organization of multiple rings in a periodic structure. Both types of lattices, honeycomb and assembled multi-rings, have then been fabricated by using a low one-photon absorption-based direct laser writing technique, in a simple way. By controlling the exposure doses (laser power and exposure time), the lattice parameters, such as structure filling factor and lattice periodicity, were well controlled, allowing an optimization of the photonic bandgap.

**Author Contributions:** Q.C.T., T.Q.T. and N.D.L. proposed the idea and designed the experiments. T.-P.N. performed simulations and Q.C.T. realized experiments and analyzed the data. All authors wrote, reviewed and approved the final version of manuscript.

**Funding:** This research is supported by a public grant of Ministry of Science and Technology of Vietnam under the project [ÐTÐLCN.01/2017] and by a public grant overseen by the French National Research Agency in the frame of GRATEOM project [ANR-17-CE09-0047-01].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Bubbly Water as a Natural Metamaterial of Negative Bulk-Modulus**

#### **Pi-Gang Luan**

Department of Optics and Photonics, National Central University, Jhongli District, Taoyuan City 32001, Taiwan; pgluan@dop.ncu.edu.tw

Received: 18 August 2019; Accepted: 30 August 2019; Published: 1 September 2019

**Abstract:** In this study, an oscillator model of bubble-in-water is proposed to analyze the effective modulus of low-concentration bubbly water. We show that in a wide range of wave frequency the bubbly water acquires a negative effective modulus, while the effective density of the medium is still positive. These two properties imply the existence of a wide acoustic gap in which the propagation of acoustic waves in this medium is prohibited. The dispersion relation for the acoustic modes in this medium follows Lorentz type dispersion, which is of the same form as that of the phonon-polariton in an ionic crystal. Numerical results of the gap edge frequencies and the dispersion relation in the long-wavelength regime based on this effective theory are consistent with the sonic band results calculated with the plane-wave expansion method (PWEM). Our theory provides a simple mechanism for explaining the long-wavelength behavior of the bubbly water medium. Therefore, phenomena such as the high attenuation rate of sound or acoustic Anderson localization in bubbly water can be understood more intuitively. The effects of damping are also briefly discussed. This effective modulus theory may be generalized and applied to other bubble-in-soft-medium type sonic systems.

**Keywords:** acoustic metamaterial; effective medium; bubble resonance; negative modulus

#### **1. Introduction**

Metamaterials are usually manmade structures designed to have some exotic properties that are rarely or never seen in natural materials. The split-ring resonators (SRRs) array is a well-studied structure that responds to incident electromagnetic waves in a range of frequencies like a negative-permeability material [1]. Similarly, the structure of periodically arranged rubber-shelled metal spheres embedded in an epoxy background is an acoustic metamaterial with negative mass density for low frequency sound [2]. In both cases, the exotic wave properties are caused by the coupling of the propagating wave in the background medium to the local resonance mechanism in each cell. The negativity of the material parameters in these two cases also implies the existence of a frequency gap [3]. No wave with frequency in the gap can propagate in this metamaterial if it occupies the whole space. If the metamaterial is distributed in a finite region, waves get attenuated before they penetrate through the medium.

Frequency gap and wave attenuation effects can also be achieved using the Bragg scattering [3] mechanism instead of local resonance. The bandgap of a usual photonic [4] or sonic crystal [5,6] is caused by the destructive interference between different 'component waves' in the Fourier series expansion of the Bloch mode. To make such a bandgap large, a high contrast of material parameters between the periodically distributed inclusions and the background material is required [7]. In addition, the center of the gap usually corresponds to a wavelength of about two to three times the lattice constant (spatial period) of the structure. This fact means that in order to isolate waves using non-resonant structures, the entire structure has to be bulky. Therefore, resonant structure is the better choice in this respect.

Recently, some researchers have noticed that bubble phononic crystals can attenuate sound very effectively, so they designed bubble phononic screens to isolate acoustic waves [8–11]. Similar to the negative-density metamaterial mentioned above, the high attenuation ability of the bubble phononic crystal is caused by the resonance coupling of elastic waves and bubble pulsation. In fact, a similar effect of bubble pulsation resonance in water is a well-known phenomenon and has been studied for decades [12]. The sound attenuation ability of bubbly water has also been discovered theoretically for about two decades [13,14]. To the best of our knowledge, the sound attenuation effect has been explained so far by the multiple scattering theory of waves when the bubble positions and bubble sizes are randomly distributed, whereas the concept of phononic bandgap was used if the bubbles were periodically arranged. We will show in this paper that bubbly water can be treated as a natural metamaterial, and the resonance coupling of the pressure wave with the bubble pulsation results in a negative modulus of the medium, giving a low frequency gap.

Acoustic gaps and dispersion effects in a short wavelength regime for sound waves propagating in liquid have also been studied recently and attracted much attention [15,16]. The phononic gaps at elevated frequencies emerge due to symmetry breaking in phonon interactions. These theoretical predictions bear a general character shared among soft and biological (non-crystalline) materials/metamaterials, as have been verified through inelastic X-ray scattering (IXS) experiments.

The metamaterial theory of bubbly water we developed in this paper is based on an oscillator model for a pulsating bubble in water. According to this model, the pulsating bubble can be regarded as a radially vibrating oscillator. The spherical bubble (the air region) acts as a spring, and a shell of water around the bubble is the mass connected to the spring. To derive the effective bulk modulus of the bubbly water, it is assumed that all the bubbles have the same size and are periodically distributed in the water background. For the sake of simplicity, we only consider simple cubic (SC) structures. Therefore, a cubic unit cell contains three regions: the air sphere, the mass shell (water), and the water outside the shell. We treat the material in each region as a lump material. We then derive the dynamical modulus by considering the radial oscillator motion and the pressure–volume (P–V) relation of the air region.

The purpose of this paper is to provide a simple theory to study the acoustic properties of bubbly water at low filling fraction and long wavelengths. The scattering function of a single air bubble in water and the corresponding oscillator model are derived in Section 2. We then use this model to derive the bulk modulus of the bubbly water in Section 3 and briefly discuss the effects of absorption and the possibility of generalizing this theory. Numerical results are presented in Section 4. Finally, we conclude this paper in Section 5.

#### **2. Scattering Function and the Oscillator Model for a Single Bubble**

In this section, we derive the acoustic scattering function *fs* of a small air bubble for the incident plane wave. The resonance frequency ω<sup>0</sup> of the radial vibrating motion of the bubble can be read out directly from the scattering function. On the other hand, the radial stiffness constant *C* of the bubble can be derived by analyze the relation between the inner pressure and the radial displacement of the bubble. We then use the relation *C* = *m*ω<sup>2</sup> to derive the effective mass *m* of the oscillator. The result obtained in this section will be used in the next section for constructing the dynamic model we mentioned in the last section.

#### *2.1. Scattering Function of a Bubble*

Suppose a spherical air bubble of radius *ra* is surrounded by water, and the pulsation of the bubble is driven by a plane incident wave *pinc* of wavelength much larger than the bubble size (*kra* 1); here *k* = 2π/λ is the wave number. In this limit the bubble scatters the incident wave isotropically, leading to the scattering wave of the form

$$p\_{\text{scat}} = f\_{\text{\\$}} p\_{\text{inc}} G(r). \tag{1}$$

Here, *G*(*r*) = *eikr*/*r* represents the spherical wave solution of the Helmholtz equation (with source) <sup>∇</sup><sup>2</sup> + *<sup>k</sup>*<sup>2</sup> *<sup>G</sup>*(*r*) <sup>=</sup> <sup>−</sup>4πδ(*r*), and *fs* is the scattering function to be determined later.

Although the bubble is assumed to be small, we also assume that the acoustic process is fast enough to avoid heat conduction between the bubble and the environment, thus the pulsation of the bubble can be approximated as an adiabatic process in the thermodynamics sense, satisfying the adiabatic *P*–*V* relation *PV*<sup>γ</sup> = constant. Here *P* is the interior pressure of the bubble, *V* is the bubble volume, and γ = 1.4 is the specific heat ratio for air. Under this assumption, the radial motion of the bubble satisfies

$$\frac{\Delta P}{P} = -\frac{\gamma \Delta V}{V} = -\frac{3\gamma \Delta R}{R},\tag{2}$$

where Δ*P*, Δ*V*, and Δ*R* are the variation of the pressure, the volume, and the radius of the bubble, respectively, and *R* is the instantaneous radius of the bubble, respectively. Denote Δ*P* as *p*, Δ*R* as ξ, and assume both of them are small comparing to the static value of *P* and *R*. We can therefore rewrite Equation (2) as

$$p = -\frac{\Im \gamma P\_0}{r\_a} \xi = -\frac{\Im B\_d}{r\_a} \xi. \tag{3}$$

Here *P*<sup>0</sup> is the hydrostatic ambient pressure when the incident wave is absent, and *Ba* = γ*P*<sup>0</sup> is the bulk modulus of the air.

The pressure across the bubble surface must be continuous, therefore *p* = *pinc* + *pscat*. In addition, both *p* and ξ get a time factor *e*−*i*ω*<sup>t</sup>* under the influence of the incident wave. Differentiate Equation (3) twice with respect to time, we get

$$
\omega^2 (p\_{\rm inc} + p\_{\rm scat}) = \frac{3\gamma P\_0}{r\_a} \ddot{\xi}.\tag{4}
$$

However, the radial acceleration .. ξ is related to the pressure gradient via

$$
\rho\_w \ddot{\xi} = -\frac{\partial p}{\partial R} = -\frac{\partial p\_{scat}}{\partial R}.\tag{5}
$$

Here ρ*<sup>w</sup>* is density of the water, and the second equality is approximately correct under the long wavelength assumption. Combine Equations (1), (4), and (5) and use the fact *kra* 1, we find

$$f\_s = \frac{r\_d}{\alpha\_0^2/\alpha^2 - 1 - ikr\_d}.\tag{6}$$

Here ω<sup>0</sup> is the Minnaert frequency, which is the resonance frequency of the 'bubble oscillator', given by

$$
\omega\_0 = \frac{1}{r\_a} \sqrt{\frac{3\gamma P\_0}{\rho\_w}}.\tag{7}
$$

As can be easily observed, the strongest scattering happens at a frequency very close to the resonance frequency ω0. Therefore, the effective medium properties might be modified by the resonance. We will derive the effective bulk modulus in the next section.

#### *2.2. The Sti*ff*ness Constant and the Mass of the Oscillator Model*

The mechanism of this low frequency resonance can be understood intuitively. When the wavelength is much longer than the bubble radius, the pressure difference across the bubble itself is negligible, thus we can think of the bubble as a perfect spring. However, even under the long wavelength condition, the scattered wave rapidly drops in the radial direction, so there is a pressure difference between the bubble surface and an imaginary "outer surface" that surrounds the bubble and a water shell around the bubble. These two effects together cause the bubble to vibrate radially like an oscillator that is driven by a driving force. The bubble itself acts as a spring, the water shell is the mass, and the radial pressure difference is the driving force. The detailed derivation of this oscillator model can be found in [12]. Here we just derive the stiffness constant and the mass of the water shell for later use.

We denote the surface area 4π*r*<sup>2</sup> *<sup>a</sup>* of the bubble as *Si*, and rewrite the bubble pressure *p* in Equation (3) as *pa,* then according to Equation (3), the radial restoring force *Fr* = *paSi* is

$$F\_r = -\frac{3B\_a S\_i}{r\_a} \xi = -C\xi. \tag{8}$$

The radial stiffness constant *C* is related to the oscillator mass *m* through the relation *C* = *m*ω<sup>2</sup> <sup>0</sup>, so we have

$$m = \frac{3B\_a S\_i}{a r\_0^2 r\_a} = 4\pi r\_a^3 \rho\_w. \tag{9}$$

#### **3. E**ff**ective Bulk Modulus and Dispersion Relation**

In this section, we derive the effective bulk modulus of the bubbly water. For the sake of simplicity, we assume all bubbles are of the same size, and their centers are located periodically on the cites of a simple cubic lattice. We denote the bubble volume and the cell volume as *Va* and *Vc*, and the water volume *Vc* − *Va* in one unit cell as *Vw*. The volume filling fraction of the bubbles is denoted as *f* (do not confuse it with the scattering function *fs*). Thus, we have *f* = *Va*/*Vc* and 1 − *f* = *Vw*/*Vc*. We also denote the bulk modulus for air and water as *Ba* and *Bw*, respectively. Similar notation rules also apply to the densities ρ*a,* ρ*w*, the vibration displacements *ua*, *uw*, and the vibration velocities *<sup>v</sup><sup>a</sup>* <sup>=</sup> . *<sup>u</sup>a*, *<sup>v</sup><sup>w</sup>* <sup>=</sup> . *uw*. In addition, the bubble surface is denoted either as *Sa* or *Si*, and the outer surface of the 'mass shell' is denoted as *So*. Since the volume enclosed by *So* is four times of the bubble volume, and the 'outer radius' *ro* corresponding to *So* is related to the bubble radius *ra* as *r*<sup>0</sup> = 22/3*ra*. Our goal is to derive the effective density of the medium ρ*eff* and the effective bulk modulus *Beff* for use in the following two equations

$$\frac{\partial}{\partial t}(\rho\_{\varepsilon f}\mathbf{v}) = -\nabla p, \quad \nabla \cdot \mathbf{v} = -\frac{\partial}{\partial t} \bigg(\frac{p}{B\_{\varepsilon f}}\Big). \tag{10}$$

The left equation is Newton's second law expressed in density form. Here ρ*eff v* is the effective momentum density of the vibrating medium, and−∇*p* plays the role of the force density. Any component of ∇*p* is calculated from the pressure difference between the corresponding two boundaries of the cell, divided by the lattice constant. The right equation is in fact the time derivative of the equation *<sup>p</sup>* <sup>=</sup> <sup>−</sup>*Beff*∇ · *<sup>u</sup>*, which describes the relation between the volume variation rate <sup>Δ</sup>*Vc Vc* <sup>=</sup> ∇ · *<sup>u</sup>* and the corresponding pressure variation *p* (see Figure 1). With these two equations we can derive the wave equation for monochromatic (single frequency) waves in this sound-bubble coupling system and the dispersion relations of the propagating modes. The results provide simple explanation of some dramatic wave phenomena such as the very high attenuation rate of sound and Anderson localization, both of which have not yet been explained using the simple concept of frequency-dependent effective bulk modulus.

**Figure 1.** Oscillator model for deriving the effective bulk modulus and effective density. The meaning of the notations are explained in Section 3.

#### *3.1. The Pressure–Volume Relations and the E*ff*ective Bulk Modulus*

We begin with the following relation about the volumes:

$$\frac{\Delta V\_{\varepsilon}}{V\_{\varepsilon}} = \frac{\Delta V\_{a} + \Delta V\_{w}}{V\_{\varepsilon}} = f \frac{\Delta V\_{a}}{V\_{a}} + (1 - f) \frac{\Delta V\_{w}}{V\_{w}}.\tag{11}$$

Now replace the volume variation rates by the corresponding divergence of the displacements, and apply the relations ∇ · *<sup>u</sup>* <sup>=</sup> <sup>−</sup>*p*/*B*, for them, we find

$$\nabla \cdot \mathbf{v} = -\frac{\partial}{\partial t} \Big[ f \frac{p\_d}{B\_a} + (1 - f) \frac{p\_w}{B\_w} \Big]. \tag{12}$$

Comparing Equation (12) with Equation (10), and making use of the equality *pw* = *p*, we get

$$\frac{1}{B\_{eff}} = \frac{f}{B\_a} \left(\frac{p\_a}{p}\right) + \frac{1-f}{B\_w}.\tag{13}$$

The ratio *pa*/*p* will be derived in the next subsection using the oscillator model.

#### *3.2. The Dynamics Equation of the Radial Vibrating Oscillator*

We now study the radial motion of the bubble. The mass of the oscillator can be thought of as a water shell of constant density ρ*<sup>w</sup>* between the inner surface *Si* and the outer surface *So*. The total radial force *paSi* − *pSo* drives the radial motion of the mass shell, satisfying the equation of motion:

$$p\_{\mathfrak{a}}\mathbb{S}\_{\mathfrak{i}} - p\mathbb{S}\_{\mathfrak{o}} = m\ddot{\mathbb{S}}.\tag{14}$$

Applying the sinusoidal condition to ξ (ξ has a time factor *e*−iω*<sup>t</sup>* ), we get

$$m\ddot{\xi} = -m\omega^2 r\_a \left(\frac{\xi}{r\_a}\right) = -\frac{ma^2 r\_a}{3} \left(\frac{\Delta V\_a}{V\_a}\right) = \frac{m\omega^2 r\_a}{3} \left(\frac{p\_a}{B\_d}\right). \tag{15}$$

Substitute Equation (15) into Equation (14), we find

$$\frac{p\_a}{p} = \frac{3B\_a S\_o}{mr\_a \left(\omega\_0^2 - \omega^2\right)};\tag{16}$$

here we have used the alternative expression of the Minnaert frequency

$$
\omega\_0 = \sqrt{\frac{3B\_a S\_i}{mr\_a}},
\tag{17}
$$

which can be derived from Equation (7), Equation (9), and the relation *Ba* = γ*P*0.

#### *3.3. The E*ff*ective Bulk Modulus as a Function of Frequency*

Now it is straightforward to derive the frequency-dependent bulk modulus. Substitute Equation (16) into Equation (13), we get

$$\frac{1}{B\_{eff}} = \frac{1-f}{B\_w} + \frac{3fS\_o}{mr\_d(\omega\_0^2 - \omega^2)},\tag{18}$$

which can be written as

$$\frac{1}{B\_{eff}} = \frac{1}{B\_{\infty}} \left( \frac{\omega^2 - \omega\_1^2}{\omega^2 - \omega\_0^2} \right) = \frac{1}{B\_{\infty}} \left( 1 - \frac{F\alpha\_0^2}{\omega^2 - \omega\_0^2} \right) \tag{19}$$

where the parameters *<sup>B</sup>*∞, *<sup>F</sup>*, and <sup>ω</sup><sup>2</sup> <sup>1</sup> are given by

$$B\_{00} = \frac{B\_W}{1 - f}, \quad F = \frac{\omega\_1^2 - \omega\_0^2}{\omega\_0^2}, \quad \omega\_1^2 = \omega\_0^2 + \frac{3f}{1 - f} \left(\frac{S\_0 B\_W}{mr\_a}\right). \tag{20}$$

It is clear from Equation (19) that the inverse bulk modulus *B*−<sup>1</sup> *ef f* as a function of frequency has the same form as the dielectric function in the Lorentz model, or that for the polariton wave in the polar medium. According to this observation we learn immediately that the effective bulk modulus for the bubbly water medium becomes negative in the frequency range ω<sup>0</sup> <ω<ω1. This frequency interval is also the frequency gap or forbidden band of the propagating modes.

The bulk modulus of air is given by *Ba* = ρ*ac*<sup>2</sup> *<sup>a</sup>* , where *ca* is the sound speed in air. Similarly, *Bw* = ρ*wc*<sup>2</sup> *<sup>w</sup>* gives the bulk modulus of the water. Substitute these two relations and *So*/*m* = 24/3/ρ*<sup>w</sup>* = 2.52/ρ*<sup>w</sup>* into Equation (20), we get

$$B\_{\rm co} = \frac{\rho\_{\rm w} c\_w^2}{1 - f}, \quad F = \frac{2.52f}{1 - f} \Big| \frac{\rho\_{\rm w} c\_w^2}{\rho\_{\rm w} c\_a^2} \Big| \quad \omega\_1 = \sqrt{1 + F} \omega\_0. \tag{21}$$

The bulk modulus ratio between water and air is about 1.5 <sup>×</sup> 104, thus according to Equation (21) even a low filling fraction *f* will give a large *F*, corresponding to a wide gap. We will derive in the next section the dispersion relation of the propagating modes and discuss the effect of the frequency gap.

#### *3.4. Dispersion Relation of Acoustic Modes*

To determine the dispersion relation, besides the bulk modulus, we should know the effective mass density. Since we assume the wavelength is long enough, it is reasonable to assume that the bubble moves together (in phase) with the surrounding water in the same cell. This consideration leads us to the simple average result

$$
\rho\_{eff} = f \rho\_a + (1 - f) \rho\_w. \tag{22}
$$

The correctness of this expression will be tested later.

The wave equation for a monochromatic (single frequency) wave can be derived by substituting Equation (19) and Equation (22) into Equation (10). The result is

$$
\nabla^2 p - \frac{p\_{eff}}{B\_{eff}} \frac{\partial^2 p}{\partial t^2} = 0.\tag{23}
$$

Here the ratio ρ*eff* /*Bef f* before the double time derivative is the inverse square speed of sound in the bubbly water medium. Adopting water as the standard medium or reference 'wave vacuum', we can define the refractive index of this medium as

$$n(\omega) = c\_{\overline{w}} \sqrt{\frac{\rho\_{eff}}{B\_{eff}}} = n\_{\circledast} \sqrt{\frac{\alpha^2 - \alpha\_1^2}{\alpha^2 - \alpha\_0^2}}.\tag{24}$$

Here *n*<sup>∞</sup> = *cw* \$ ρ*eff* /*B*<sup>∞</sup> is a reference index at high frequency. This refractive index becomes imaginary in the frequency gap, which indicates that waves cannot propagate in the medium if its frequency is within the gap.

We can deepen our understanding of the refractive index function by studying the dispersion relation of the propagating modes. Substitute the plane wave form *<sup>p</sup>* <sup>=</sup> *<sup>p</sup>*0*ei*(*K*·*r*−ω*t*) into Equation (23), we get

$$K^2 = \frac{n^2(\omega)\omega^2}{c\_\omega^2} = \frac{n\_\infty^2 \omega^2}{C\_\omega^2} \left(\frac{\omega^2 - \omega\_1^2}{\omega^2 - \omega\_0^2}\right). \tag{25}$$

This dispersion relation is of the same form as that for the phonon-polariton waves. According to Equation (25), propagating modes can exist with frequency below ω<sup>0</sup> or above ω1. Within the interval ω<sup>0</sup> <ω<ω1, the wave vector becomes imaginary, thus the wave amplitude decays along the wave vector direction. If the bubbly water medium is filled in a slab region of finite thickness, an incident wave having a frequency within the gap would get attenuated before it penetrates through the slab. Furthermore, if a pulsating source of sound surrounded by the 'bubble cloud' radiates sound waves of frequency in the gap, the sound wave would be trapped by the bubble cloud, leading to the phenomena such like the Anderson localization of acoustic waves.

The negative modulus of the bubbly water is caused by the dynamical coupling of the sound field in water and the local resonance of every bubble oscillator. Similar mechanism for polar crystal also leads to the negative dielectric function in the polariton gap. The main difference between these two systems is that in the bubbly water case the oscillators vibrate radially so they couple with the scalar pressure field, whereas in the polar crystal case the local dipoles (for example, the *Na*<sup>+</sup> <sup>−</sup> *Cl*<sup>−</sup> pairs) are coupled with the electric vector field.

#### *3.5. Absorption E*ff*ect and the Possibility of Generalization*

In the derivation of the effective bulk modulus, we never considered any absorption or loss effects. For real bubbles, the damping effect caused by thermal conductivity (non-adiabatic correction), shear viscosity absorption at the bubble surface, and the radiation loss due to scattering must be considered. These corrections will cause the resonance frequency of the bubble to change slightly, reducing the peak of the scattering function and widening its width. These corrections will also introduce a damping force on the radial motion of the bubble oscillator and change the bulk modulus such that it does not acquire an unphysical infinite value. Since the absorption effect is a very cumbersome effect, and our aim in this study is not to study this effect, we will not further address this issue. We follow the spirit of [13,14] to choose the proper bubble size so that we can avoid the effects of the absorption.

Another question is how to generalize our theory to elastic solid media containing bubbles. As the authors have demonstrated in [8–11], for elastic media with very low shear stress, the primary mechanism of the first band gap is still caused by the resonance of the bubble. Therefore, we believe that our theory can be generalized and applied to such a system without much difficulty.

#### **4. Numerical Results**

In this section we discuss some numerical results based on our theory. Hereafter we assume ρ*<sup>a</sup>* = 1.2 kg/m<sup>3</sup> , ρ*<sup>w</sup>* = 1000 kg/m<sup>3</sup> , *ca* = 343 m/s , and *cw* = 1490 m/s . These parameters give us the bulk modulus *Ba* = ρ*ac*<sup>2</sup> *<sup>a</sup>* = 1.412 <sup>×</sup> 105 Pa and *Bw* = <sup>ρ</sup>*wc*<sup>2</sup> *<sup>w</sup>* = 2.22 <sup>×</sup> 109 Pa. If the radius of a bubble is measured in mm, the Minnaert frequency measured in kHz and we get a simple formula ω<sup>0</sup> = 20.58/*ra* or ν<sup>0</sup> = ω0/2π = 3.28/*ra* . To avoid high absorption, we assume the bubbles have radius larger than 100 μm [12,13].

The scattering function *fs* as a function of dimensionless frequency *kra* = ω*ra*/*cw* is plotted in Figure 2. Here we replaced the *fs* function by its dimensionless correspondence *fs*/*ra* because here the information of the size is irrelevant. According to Figure 1, the absolute value as well as the imaginary part of *fs* goes to the maximum value 72.5*ra* at *kra* = 0.0138, and behaves like an even function with respect to this peak location. The real part of *fs*, on the other hand, first raises to 36.6*ra*, and then drops to zero at this place, and continually goes downwards to −36.5*ra*, then approaches the horizontal axis without crossing the zero value. Resonance of this kind belongs to the Lorentz dispersion, which is responsible for the negative bulk modulus and the low frequency gap of this bubbly water medium.

**Figure 2.** The scattering function of a single bubble. The notations in this figure have been defined in Section 2.

We now study the characteristics of the acoustic band structures of the bubbly-water sonic crystals. Plane wave expansion method (PWEM) were implemented for the calculations. Results for filling fraction *f* = 0.001 and *f* = 0.01 are shown in Figure 3. These results indicate that bubbly-water sonic crystal is an ideal structure to open large gap at low frequency. For the case with filling fraction *f* = 0.001, the gap is opened between ω*Lra*/*cw* = 0.0185 (at X) and ω*Ura*/*cw* = 0.085 (at Γ). Comparing them with the gap boundaries ω0*ra*/*cw* = 0.0138 and ω1*ra*/*cw* = 0.0881 predicted by the effective metamaterial theory, we find excellent agreement between them. The low frequency limit of the gap in Figure 3 is a little higher than the prediction. We believe this discrepancy is due to the bad convergence of the PWEM at low filling fraction. Similarly for the *f* = 0.01 case, we have the gap opened from ω*Lra*/*cw* = 0.0173 (at X) to ω*Ura*/*cw* = 0.271 (at ¡). The gap boundaries provided by the effective theory are ω0*ra*/*cw* = 0.0138 and ω1*ra*/*cw* = 0.277, also in very good agreement with the band structure result. We noticed that the agreement between ω*<sup>L</sup>* and ω<sup>0</sup> in this case is better than in the former case. This fact is consistent with our judgement about the PWEM.

**Figure 3.** The acoustic band structures for filling fraction (**a**) f = 0.001, and (**b**) f = 0.01.

The metamaterial theory developed in this paper not only can provide the range of the gap, but also can give us the attenuation rate of the wave, which is directly related to the imaginary part of the effective wave number that determined by Equation (25). In addition, absorption or loss effect can be easily included by a simple replacement of the frequency: ω → ω(1 + *i*δ) [12]. Here δ is a positive number smaller than 0.1 if we assume the that the bubble radius is larger then 100<sup>1</sup> m. For the purpose of demonstration, in the simulation of lossy media, we take δ = 0.05. Figure 4a,b is the results for the lossless media, while Figure 4c,d is the results for the lossy media. Here the dimensionless variable *Ka* is the product of the wavenumber *K* in Equation (25) times the lattice constant *a* of the SC structure. the lattice constant in the SC structure. The real part of the wave number (wave vector) of the mode waves are represented by the blue solid curves, and the red broken lines represent the imaginary part of the wave vector. These dispersion relations are similar to the dispersion relations of the phonon-polariton that mentioned in the previous section. We noticed that the presence of loss or absorption does not destroy the attenuation effect because it is so large.

**Figure 4.** Dispersion relations for f = 0.001 and f = 0.01, lossless media are shown in (**a**) and (**b**). The corresponding results for lossy media are shown in (**c**) f = 0.001 and (**d**) f = 0.01.

Now we discuss the effective bulk modulus according to Equation (19) and its modified form when loss cannot be ignored. As before, four cases are studied, and their results (*Beff* /*Bw*) are shown in Figure 5. The parameters used in the calculations are the same as those used in Figure 4. The small circles in the four subplots indicate the starting frequencies (changes from positive to negative) of the negative bulk modulus. For the lossless cases, the starting points and the resonance points (infinite effective modulus) are the same as the gap boundaries of the metamaterial. The loss effect seemed not to change the effective bulk modulus much if the frequency is not very close to the resonance point. However, the effective modulus at the resonance frequency no longer goes to the unphysical value of infinity, as expected.

**Figure 5.** Bulk modulus for lossless cases with f = 0.001 and f = 0.01 are shown in (**a**) and (**b**). The corresponding lossy cases are shown in (**c**) f = 0.001 and (**d**) f = 0.01.

Finally, we discuss the valid range of the theory developed in this paper. Based on our numerical simulations, we found that the gap range given by the effective theory is in good agreement with that obtained from the band structure calculation (the discrepancy is less than 3%) if the filling fraction *f* of the bubbles is within the range 0.0005 to 0.02. Beyond this range, we are not sure about the accuracy, and need further investigation. In addition, the bubble radius must not be less than 100 μm, otherwise high absorption effects must be considered [12].

#### **5. Conclusions**

In this paper, we apply the oscillator model of air-bubble-in-water to the bubbly water medium and derive the effective bulk modulus of the medium. We show that in a wide range of frequency the effective modulus becomes negative, while the effective density of the medium is still positive. These properties imply the existence of a wide acoustic gap, consistent with band structure calculations. The dispersion relation for the acoustic modes in this medium has Lorentz type dispersion, like that of the phonon-polariton. The theory proposed in this paper provides a simple mechanism for explaining the long-wavelength behavior of the bubbly water medium. Using this theory, phenomena such as the high attenuation rate through a bubble screen or the acoustic Anderson localization in bubbly water can be understood more intuitively. The effects of damping are also discussed. This effective modulus theory may be generalized and applied to other similar structural sonic medium containing periodically arranged bubbles.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **All-Optical Ultra-Fast Graphene-Photonic Crystal Switch**

#### **Mohammad Reza Jalali Azizpour 1, Mohammad Soroosh 1,\*, Narges Dalvand <sup>2</sup> and Yousef Seifi-Kavian <sup>1</sup>**


Received: 30 June 2019; Accepted: 28 August 2019; Published: 3 September 2019

**Abstract:** In this paper, an all-optical photonic crystal-based switch containing a graphene resonant ring has been presented. The structure has been composed of 15 × 15 silicon rods for a fundamental lattice. Then, a resonant ring including 9 thick silicon rods and 24 graphene-SiO2 rods was placed between two waveguides. The thick rods with a radius of 0.41a in the form of a 3 × 3 lattice were placed at the center of the ring. Graphene-SiO2 rods with a radius of 0.2a were assumed around the thick rods. These rods were made of the graphene monolayers which were separated by SiO2 disks. The size of the structure was about 70 μm2 that was more compact than other works. Furthermore, the rise and fall times were obtained by 0.3 ps and 0.4 ps, respectively, which were less than other reports. Besides, the amount of the contrast ratio (the difference between the margin values for logics 1 and 0) for the proposed structure was calculated by about 82%. The correct switching operation, compactness, and ultra-fast response, as well as the high contrast ratio, make the presented switch for optical integrated circuits.

**Keywords:** graphene; kerr effect; optical switch; photonic band gap; photonic crystal

#### **1. Introduction**

In recent years, increased demand for ultra-fast processing systems has necessitated the development of optical devices. Many attempts have been accomplished to introduce all-optical structures that allow circuits operating at a frequency of terahertz. Achieving a high data transferring rate with a possibility for integration is among the important issues in designing all-optical circuits [1,2].

Photonic crystals (PCs) were initially and simultaneously proposed by Yablonocvitch [3] and John [4] in 1987. PCs are structures composed of periodic layers with different dielectric constants and demonstrate useful characteristics such as photonic band gap (PBG) [3], slow light regime [5–7], and self-collimation [8–11]. Because of their compact size and the aforementioned properties, they are promising candidates in realizing all-optical devices such as optical waveguides [10–12], filters [13–18], demultiplexers [19–25], switches [26–30], logic gates [31–36], encoders [37–42], decoders [43–48], and analog-to-digital converters [49–54].

In photonic integrated circuits (PICs), optical switches play a crucial role in guiding light to the desired directions. Generally, optical switches include ring resonators, input and output ports that drop light with a particular wavelength. PCs-based switches have been accomplished by using the Kerr nonlinear effect [28,55]. Hache et al. developed the first concept using the Kerr nonlinearity effect in PIC. The proposed structure included two materials, Si and SiO2, with a 1.5 μm stop band and 18 Gw/cm2 functional power. Although their structure shed light in the all-optical circuits area, the level of power was not applicable for PICs [56]. Another optical switch suggested by Alipour-Banaei et al. consisted of the chalcogenide material with the nonlinear coefficient 9 <sup>×</sup> 10−<sup>17</sup> m2/W. The threshold power level for switching operation and the footprint were 2 KW/μm2 and 259 μm2, respectively [57].

Serajmohammadi et al. proposed a PC-based NAND gate using an optical switch. The structure was made of one bias signal, two input ports, and one output port. Using the bias signal, the Kerr effect was approached when both input signals were introduced into the switch. They also used the chalcogenide glass to gain the nonlinear properties for a resonant ring. Even though the threshold power and the footprint reduced to 1.5 KW/μm2 and 230 μm2, respectively, the power was so high to employ in integrated circuits [58]. In this way, Alipour-Banaei et al. proposed another optical switch by a focus on the power issue [59]. Although they approached the power level 1.5 KW/μm2, the footprint was notably increased to as large as 360 μm2.

Recently, Daghooghi et al. developed a PC-based switch for a decoding operation which used a 31 × 31 lattice of silicon rods in the air [60]. They employed nanocrystal material in which the nonlinear coefficient was 10−<sup>16</sup> m2/W and succeeded in reducing the threshold power level to 13 W/μm2, while the footprint was considerably increased.

As mentioned, many attempts have been done to reduce the threshold power level for switching, as well as the size of the structure. The possibility of fabrication for the proposed structures is an important challenge that restricts researchers in selecting the different materials and sizes. So, the compatibility with conventional technologies helps researchers achieve new materials with distinct characteristics. Graphene has a substantially more Kerr coefficient than other materials such as chalcogenide and nanocrystal [35,61,62]. Several experiments and calculations have been done for determination of Kerr coefficient of graphene [63–70]. Soh et al. have shown the Kerr coefficient of graphene at a wavelength of 1550 nm is equal to 10−<sup>15</sup> m2/W [62]. They presented a comprehensive analysis for nonlinear optical Kerr effect in graphene including two-photon absorption, Raman transition, self-coupling, and quadratic AC Stark effect. They calculated the absorption rate using the S-matrix element and converted it to nonlinear refractive index coefficients. Then, they obtained the rates of distinct nonlinear processes that contribute to the Kerr nonlinear refractive index coefficient. Due to this fact, using graphene in PICs results in decreasing the functional optical intensity in comparison to other materials such as chalcogenide and nanocrystal. This is the main reason why we are utilizing graphene. Moreover, graphene has an impact on two other factors in PICs fabrication: contrast ratio and footprint.

In this study, a new graphene-based switch is proposed to enhance the power and footprint difficulties. The switch is composed of 9 thick silicon rods at the center of the resonant ring and 24 graphene-SiO2 stacks as the nonlinear rods around them. Each nonlinear rod is made of the graphene monolayers which are separated by SiO2 disks. Using graphene monolayers enhances the light-material interaction and results in decreasing the threshold intensity for switching operation. Furthermore, the more Kerr coefficient of graphene than chalcogenide and nanocrystal used in the previous works assists with enhancing the power difficulty in this work. The obtained results of the simulation present 70 μm2 and 0.235 W/μm<sup>2</sup> for the footprint and the threshold intensity level, respectively. Furthermore, the normalized output power margins, 4% and 86% for logic 0 and 1 illustrate that the structure is potentially a good candidate in PIC applications.

In this paper, three other sections have been organized as follows; in Section 2, the structure will be presented in detail. Then, the simulation results will be provided in Section 3. In Section 4, the evaluation of the device and discussions will be described, along with a comparison of the obtained results with other works. Finally, a conclusion of the work will be presented in Section 5.

#### **2. Materials and Methods**

The primary structure includes a lattice of 15 × 15 rods in air background where the lattice constant is 558 nm. So, the overall size of the structure is about 70 μm2. The refractive index of rods is 3.46, and the radius of rods for the fundamental structure is 0.2a, where a is the lattice constant. To calculate the band structure, the plane wave expansion (PWE) method is used [71]. According to the band structure (as shown in Figure 1), two PBGs at TE mode are obtained for 0.29 ≤ *a*/λ ≤ 0.42 and 0.725 ≤ *a*/λ ≤ 0.74 where λ is the wavelength. So, for 1290 nm ≤ λ ≤ 1990 nm optical waves could

not be allowed for propagation throughout the structure. Including the third optical communication window, the last interval is used for the proposed device.

**Figure 1.** The band diagram of the proposed structure.

The formation of the switch was gained by removing and changing the radius of some rods in *X* and *Z* directions. One of the best features of the proposed structure is that there is not any bias or enable port. It can be seen in Figure 2 that there are two output ports (named as O1 and O2) and one input port (named as IN) that are connected via a resonant ring. A 3 × 3 lattice rod formed by changing in the radius 0.41a was placed in the center of the ring. Around the core, 24 graphene-SiO2 rods (with blue color) were used in which the linear refractive index (*n*1) and nonlinear coefficient (*n*2) of graphene were 2.6 and 10−<sup>15</sup> m2/W, respectively [62]. These rods are made of the graphene monolayers which are sandwiched between SiO2 disks. Berman et al. presented a photonic crystal structure including a periodic array of the graphene-SiO2 stacks [72]. They calculated the transmission coefficient for different frequencies and showed that the structure can be used as frequency filters and waveguides for optical waves. According to their research, we used the graphene monolayers which were separated by SiO2 disks for blue color rods. Using graphene in the ring resonator causes the possibility of achieving the nonlinear Kerr effect. The optical Kerr effect is generally known as the changing in the refractive index (*n*) in response to applied light intensity (*I*) and is defined by *n* = *n*<sup>1</sup> + *n*2*I* [65,66]. The nonlinear coefficient of graphene is defined as [73]:

$$m\_2 = \frac{3\eta}{n\_1^2 \epsilon} \chi^{(3)}$$

where *n*<sup>1</sup> is the linear refractive index, η is the impedance, and χ(3) is the nonlinear susceptibility. The nonlinear susceptibility has been described as follows [65]:

$$\chi^{(3)} = \frac{\sigma\_{\mathfrak{F}}(\alpha)}{\alpha d\_{\mathfrak{F}}}$$

where ω is the frequency and *dg* is the thickness of the graphene monolayer. The dielectric constant is defined as [72]:

$$
\varepsilon(\boldsymbol{\omega}) = \varepsilon\_0 + \frac{4\pi i \sigma\_\circ(\boldsymbol{\omega})}{\boldsymbol{\alpha}d},
$$

where ε<sup>0</sup> is the dielectric constant of SiO2, *d* is the width of SiO2 layer and equal to 1 nm. The dynamical conductivity of the graphene (σ*g*) is defined by the following equation [72]:

$$\sigma\_{\mathcal{S}}(\omega) = \frac{\varepsilon^{2}}{4\hbar} \Big[ \eta (\hbar \omega - 2\mu) + \frac{i}{2\pi} \Big( \frac{16K\_{\text{B}}T}{\hbar \omega} \log \left[ 2\cosh \left( \frac{\mu}{2K\_{\text{B}}T} \right) \right] - \log \frac{\left( \hbar \omega + 2\mu \right)^{2}}{\left( \hbar \omega - 2\mu \right)^{2} + \left( 2K\_{\text{B}}T \right)^{2}} \right) \Big]$$

where *e* is the electron charge, *KB* is the Boltzmann constant, and μ is the chemical potential.

**Figure 2.** The proposed all-optical switch.

Changing the refractive index of the nonlinear rods results in changing the effective index of the ring resonator, and hence, the switching operation could be approached. The more nonlinear coefficient for graphene in comparison with chalcogenide and nanocrystal [60,62] assists in reducing the needed optical power for switching. In this study, with respect to the aforementioned issue, the graphene was used in the ring resonator to propose a low-threshold all-optical switch. Furthermore, the high contrast ratio for digital applications as the difference between the normalized power margins for logic 0 and 1 is another advantage of the presented structure.

#### **3. Results**

To simulate the optical wave propagation throughout the proposed structure, the finite difference time domain (FDTD) method was used [74]. According to the Courant condition, it was assumed as follows [74]:

$$c\Delta t < \frac{1}{\sqrt{(\frac{1}{\Delta x})^2 + (\frac{1}{\Delta z})^2}}\gamma$$

where *c* is the speed of light in vacuum, Δ*t* is the time step, and Δ*x* and Δ*z* are mesh sizes in both *X* and *Z* directions. The second condition is about the grid spacing that should be less than λ/10. The time step was assumed Δ*t* = 0.2 fs and the structure was discretized in which the length of unit cells was Δ*x* = Δ*z* = 100 nm. In the first stage, an optical pulse was applied to the input port and the pulse response was calculated at output ports. Figure 3 shows the frequency harmonic amplitude for O1 and O2 ports. It can be seen that the optical intensity between two output ports at λ = 1547 nm can be approached for switching applications.

**Figure 3.** Response of the switch for different wavelengths.

To evaluate the resonance phenomenon, the transmission ratio for different input intensities was simulated (see in Figure 4). In this study, the intensity value at output port divides to input intensity is defined as the transmission ratio factor. When a signal is introduced to the input port, the output ports will be ON depending on the entrance optical intensity. If the input intensity is less than 0.235 W/μm2, the resonance wavelength of the ring will not change. As a result, O1 and O2 ports will be ON and OFF, respectively. The resonance wavelength can be changed because of the Kerr effect, so the optical waves will considerably be coupled to the O2 port for amounts of more than 0.235 W/μm2. In this case, O1 and O2 ports will be OFF and ON, respectively. It can be concluded that the threshold intensity for switching is around 0.235 W/μm2. This value is less than one in other works [45,57,59,60], so it can be considered an advantage of the proposed switch.

**Figure 4.** Transmission ratio of the structure for different input intensities.

Figure 5 shows the normalized electric field distribution for two optical intensities, 0.1 W/μm2 and 0.45 W/μm2. A color bar has been inserted on the right side in which the red and black colors are the response in +1 and −1 for the domain amplitude of the field. It can be seen that optical waves are dropped from the upper waveguide to the lower waveguide for *I* = 0.45 W/μm2 while the dropping operation is not done for *I* = 0.1 W/μm2.

**Figure 5.** Electric field distribution into the structure for: (**a**) *I* = 0.45 W/μm2; (**b**) *I* = 0.1 W/μm2.

For more evaluation, the pulse response of the presented structure corresponding to Figure 5 has been calculated. As shown in Figure 6, the normalized powers at port O1 and O2 are obtained by 86% and 2% for *I* = 0.1 W/μm2, respectively, while ones are reached to 4% and 95% for *I* = 0.45 W/μm2. Also, as Figure 7 presents, the rise and fall times of the structure are obtained by 0.3 ps and 0.4 ps, respectively. The rise time is defined as the needed time for coupling the optical waves from one waveguide to another. This time is recorded when the power at the output port reaches 90% steady-state value. The fall time is assumed as the during time for decreasing power from a steady-state value to 10% of it. Comparing the obtained times with the ones in other works [75] demonstrates that the proposed graphene-based switch is faster than them. This characteristic is an essential factor in optical switching applications.

**Figure 6.** The normalized output power for: (**a**) *I* = 0.45 W/μm2; (**b**) *I* = 0.1 W/μm2.

**Figure 7.** The pulse response for proposed optical switch.

#### **4. Discussion**

Approaching the fast response of a switch is one of the main challenges for designing different switched-based applications. Furthermore, the small size of the structure results in it being capable of being considered for integrated circuits. Also, the low power and high contrast ratio are other important parameters for designing different switches. Some attempts have been done to improve the characteristics of optical switches [61,75].

In this study, a graphene-based all-optical switch has been presented in Section 2. Figures 5 and 6 demonstrate the correct operation of the switch. To assess this work, the obtained results have been compared with references [57–60] in Table 1.


**Table 1.** Comparing the simulated results with ones in other works.

\* NA: Not Assigned.

As presented in Figure 6, the time analysis shows that the proposed switch is faster than the one in other works [60]. Approaching the amount of 0.3 ps and 0.4 ps for rise and fall times is the main advantage of this work. Because of the large nonlinear coefficient of graphene in comparison with materials such as chalcogenide, silicon, nanocrystals, and gallium arsenide, the needed threshold intensity for switching operation in this study has been reduced to 0.235 W/μm2 in comparison to references [57–60]. This is another advantage of the graphene-based switch. Furthermore, the size of the proposed structure is less than the one in references [57–60]. It is obvious that the compactness of the switch is necessary for use in optical integrated circuits.

Contrast ratio of the output ports is calculated at 82%. The presented switch covers the high contrast ratio while this parameter has not been reported in references [57–60]. Besides all of these advantages, tuneability of some electrical and optical parameters for graphene is a principal characteristic for designing photonic crystal-based devices. Changing the refractive index of graphene in response to electric field affects optical power transmission toward the output ports. Also, the resonant frequency (or wavelength) of the switch could be changed to reach the desired value after the fabrication process.

#### **5. Conclusions**

In this study, a photonic crystal-based structure was proposed as an all-optical switch. Using the graphene-SiO2 stacks as rods in the resonant ring, switching operation was correctly approached. The rise and fall times of the device were obtained by 0.3 ps and 0.4 ps, respectively, that were less than the ones in other reports. The size of the presented structure was just about 70 μm<sup>2</sup> and hence was smaller than the one in the aforementioned works of Section 1. The high nonlinear coefficient of the graphene resulted in a decreasing of the needed optical intensity to 0.235 W/μm2 for switching operation in comparison with other works. Furthermore, simulation results showed the contrast ratio of the switch was about 82%. In respect of the obtained results, it can be concluded that the presented all-optical switch is capable of consideration for optical integrated circuits.

**Author Contributions:** Formal analysis, M.R.J.A. and N.D.; Methodology, M.S. and M.R.J.A.; Writing-Original Draft, M.R.J.A.; Writing-Review and editing, M.S. and Y.S.-K.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no competing financial interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Photonic Crystal Cavity-Based Intensity Modulation for Integrated Optical Frequency Comb Generation**

#### **Henry Francis 1,\*, Si Chen 1, Kai-Jun Che 2, Mark Hopkinson <sup>1</sup> and Chaoyuan Jin 1,3,\***


Received: 30 August 2019; Accepted: 23 September 2019; Published: 25 September 2019

**Abstract:** A simple scheme to generate an integrated, nanoscale optical frequency comb (OFC) is numerically studied. In this study, all optical intensity modulators based on photonic crystal (PhC) cavities are cascaded both in series and parallel. By adjusting the modulation parameters, such as the repetition rate, phase, and coupling efficiency of the modulating wave, it is possible to produce combs with a variety of different characteristics. Unique to PhC intensity modulators, in comparison with standard lithium niobate modulators, is the ability to control the amplitude of the light via a cavity rather than controlling the phase through one arm of a Mach–Zehnder interferometer. This opens up modulation-based OFC generation to new possibilities in both nanoscale operation and cavity-based schemes.

**Keywords:** photonic crystals; microwave photonics; optical frequency combs

#### **1. Introduction**

In recent years, optical frequency combs (OFC) have become an increasing popular research area [1–3]. The rapid adoption of frequency combs into many different research topics has meant the need for diverse functionality and integration [4]. There are many schemes for generating an OFC, which include mode locked lasers [5], electro-optic (EO) modulators [6], and micro toroidal cavities [7]. Mode locked lasers have the ability to generate combs with a fixed phase relationship and stable operation. However, the frequency space between comb lines is governed by the device length due to fundamental operation principles and hence limit their abilities and flexibility to be integrated onto photonic circuits. Another well-established area of OFC generation is in micro-toroidal cavities [8]. This technique relies on nondegenerate four wave mixing within the cavity to generate a broad band of intensity spikes in the frequency domain. Using this technique, on chip integration of an OFC generator has been achieved [9]. This provided a major step forward within the research community. Currently, OFC generation by this method is receiving a lot of attention, a comprehensive review is given in ref. [1]. While the work outlined in ref. [1] shows OFCs with large spectral range and on-chip integration, the work presented here is in parallel and looks towards new technologies for generating OFCs. By generating an OFC from photonic crystal (PhC) cavities and waveguides, the work looks towards improving spectral homogeneity, reducing the device volume and exploring the potential for nanoscale photonic integration of OFCs.

A common method of OFC generation is the use of cascaded EO LiNbO3 intensity modulators (IM) [10]. In this method, light is split into two arms and the phase of the light through one of arms is modulated via a radio frequency (RF) electrical signal. The RF signal induces a refractive index change of the waveguide material, hence changing the phase of the light. The recombination of the two carrier signals with a controllable phase difference causes modulation of the carrier lights intensity. A DC bias is also present through each arm so that a constant phase shift can be applied. The voltage of the RF signal used to modulate the light will determine the modulation depth and the voltage of the DC bias will determine a constant phase shift. To generate an OFC, the light is modulated using a sinusoidal electrical signal. This generates sidebands in the frequency response of the modulated light, and the position of these sidebands, relative to the initial frequency of the carrier light, has a direct relation to the frequency of the electrical signal. The amplitude of the sidebands is determined by the modulation depth and phase shift. By fine control of these parameters, an intensity modulator can generate multiple comb lines in the frequency response of the modulated light [11]. A second IM leads from the output of the first IM, where the voltage of the RF and DC bias are kept the same but the modulation frequency is a fifth of the first IM. This generates a large number of homogeneous comb lines centered around the initial frequency of the input light. This technique provides a simple solution to generating a multiple wavelength source, however this method cannot be implemented with chip-scale components due to the need for relatively large waveguides to achieve deep enough modulation. In the work presented here, the operating principle outlined above is adapted and built upon for the implementation in PhC waveguides and cavities.

The standard modulation process of LiNbO3 modulator is fundamentally very different to the process undergone in photonic crystal cavity modulators. Here, a mismatch between the cavity resonance and the propagating wave will either cause transmission or reflection, depending on the implemented scheme. In this paper, two intensity modulation schemes are analyzed; direct and side coupling, as shown in Figure 1. In the direct coupling scheme, the transmitted power is proportional to the energy inside the cavity, this will decrease the achievable contrast ratio. This scheme will also lead to the transmitted power having a similar pulse shape to the control light. However, in the side coupled scheme, interference between the decaying wave from the cavity and the propagating wave in the waveguide mean a higher contrast ratio and variable pulse shape. In both schemes, an increase in the optical excitation power leads to a larger cavity resonance shift, this has a large effect on the pulse shape of the carrier wave and the sidebands produced from the modulation process. By arranging two or more PhC IMs in either a series or parallel formation, it is possible to increase the number of generated sidebands and hence observe an OFC. In this paper, detailed analysis of the different modulation schemes will be undertaken. This analysis will confirm that an OFC can be produced from a PhC device, thus broadening the application of PhC devices into the mature and prosperous research field of OFC generation.

**Figure 1.** The two all-optical intensity modulator schemes studied in this paper. *a*(*t*) is the cavity field amplitude, *u* is the coupling factor, and *sin* (*t*) and *sout* (*t*) represent the amplitude of the propagating wave in the input and output waveguide, respectively.

#### **2. Materials and Methods**

The scheme proposed here is based on photonic crystal waveguides and cavities to construct all-optical intensity modulators, as shown in Figure 1. This method of all-optical modulation has been studied in great detail by different research groups [12–15], which has led to a rich and common understanding of the fundamental principles behind PhC switches. PhC cavities have a very high quality factor (Q) while maintaining a low mode volume (V). This means that nonlinearity within the cavity can be greatly enhanced, given that the field intensity scales with Q/V. A shift in the cavity resonant wavelength, due to a refractive index change, is therefore possible with a relatively low-power pump light. Intensity modulation of a carrier light is therefore possible when a pump light is used to change the transmission wavelength. The device itself is made up of holes in triangular lattice formation fabricated in an InP membrane with embedded InAs quantum dots (QDs), with broad band emission around 1550 nm. The waveguides are made up of a line of missing holes and the cavity consists of a three hole defect, known as an L3 cavity, as shown in Figure 2. To ensure the cavity has a resonant frequency around 1550 nm the filling factor and lattice constant have been carefully calculated. The lattice constant is calculated to be 480 nm and the filling factor at 0.29. The thickness of the membrane will be 220 nm. The parameters given were calculated in COMSOL by finding the eigenfrequency via the finite element method (FEM). A scanning electron microscope (SEM) image of an L3 cavity along with the FEM simulation results for an L3 cavity with a fundamental mode at 1550 nm is shown in Figure 2.

**Figure 2.** (**a**) An scanning electron microscope (SEM) image of a photonic crystal membrane with an L3 cavity in the center. (**b**) Finite element method (FEM) calculation of the fundamental mode at 1550 nm.

Coupled mode theory (CMT), a standard PhC analysis tool, is used in this paper to simulate the device described above. The equations employed in this paper are based on an extension of CMT models that have been presented previously [13,16,17]. In the simulated device, it is assumed that an active region is present in the form of InAs QDs. This will lead to a dominating third order nonlinearity process based on the saturable absorption (SA) of QDs. Therefore, the refractive index change will be dominated by the Kramers–Kronig relation to the absorption coefficient. This has a direct relation to the field amplitude within the photonic crystal cavity and hence will be the main dynamic variable used throughout. A pump and carrier light are injected into the input waveguide, denoted by *s p,c in* in Figure 1. The dynamical equations for the cavity modes is written as:

$$\frac{da\_{p,\varepsilon}}{dt} = (-i(\omega\_0 + \Lambda \omega\_p(t) - \omega\_{p,\varepsilon}) - \gamma\_{\text{total}})a\_{p,\varepsilon}(t) + \sqrt{2\gamma\_{\varepsilon}}s\_{\text{in}}^{p,\varepsilon}(t) \tag{1}$$

Here, *ap,c* is the amplitude of the cavity mode which is excited by the pump and carrier light, respectively and the energy inside the cavity is represented by |*ap,c(t)|2*. The field of the cavity is represented by *A p,c in (t) = a p,c in (t) exp(-i*ω*p,ct)* and the input light through the waveguide is given by *S p,c in (t) = s p,c in (t) exp(-i*ω*p,ct*), the power of the light is represented by *|s p,c in (t)|2*. It is assumed that *ap* is far higher than *ac*, such that only *ap* will have an effect on the cavity resonance shift. ω*p,c* is the angular frequency of the input pulses, for either the pump or carrier light. In order for the cavity resonance to

change, the amplitude of the pump light must increase. The change is then dependent on the energy in the cavity due to the pump light and a third order non-linearity, denoted as the Kerr constant. This is given below:

$$
\Delta\omega\_p(t) = \mathcal{K}\_{krr} |a(t)|^2 \tag{2}
$$

*Kkerr* represents a simplified value for the nonlinear dynamics in the system [12,17], given by:

$$K\_{krrr} = \frac{\omega\_0 c n\_2}{n\_{eff} V\_{krrr}} \tag{3}$$

*n2* is the enhanced Kerr coefficient, due to the SA of QDs in the membrane structure and the average rate of SA within the cavity is calculated and given by *Vkerr*. These values are greatly enhanced by the PhC cavity due to the small mode volume and high Q attainable in PhC structures. *neff* is the effective refractive index, this due to the membrane thickness, material properties and resonant wavelength of the cavity. The total loss rate in the cavity is split into two parts; the intrinsic loss, γ*int*, and the coupling loss, γ*c*. The intrinsic loss due to coupling into non-relevant modes, predominantly by vertical emission, is dependent on the Q factor of the cavity, such that γ*int =* ω*0/2Q*. The coupling loss is due to the coupling coefficient between the cavity and the waveguide. The light in the L3 cavity will reflect off each side of the cavity in the longitudinal direction which leads to light coupling to the waveguide at each end, hence *2*γ*c*.

In order to induce a change in the transmission through the device, a change in the refractive index is needed. The amount that the refractive index must shift by is dependent on the linewidth of the cavity mode. In the device simulated here, a resonance shift of 1 nm will cause the carrier light to either reflect or transmit through the device. To obtain a shift of 1 nm, it is calculated that a refractive index shift of 0.005 is needed. The relationship between refractive index change and cavity resonance change is defined as <sup>Δ</sup>*<sup>ω</sup> <sup>ω</sup>* <sup>=</sup> <sup>−</sup>Δ*<sup>n</sup> <sup>n</sup>* . Therefore, the profile of the refractive index change is in relation to the profile of the pump light, as given in Equation (1) and shown by the dashed line in Figure 3. Achieving a refractive index change of 0.005 using QD saturable absorption as the refractive index non-linearity is common among nanophotonics research community [18,19].

**Figure 3.** The power of the carrier wave at the output of the device is given by the solid lines. The refractive index change inside the cavity is given by the dotted line.

Using the parameters stated in Table 1, Equations (1)–(7) are solved numerically. This is done using standard routines available for solving differential equations via a self-written computing program.


**Table 1.** Physical parameters used in coupled mode theory (CMT) calculations.

#### **3. Results and Discussion**

As discussed earlier, two common cavity-waveguide all-optical modulator configurations have been considered in this paper: side and direct cavity coupling; these are schematically represented in Figure 1. Given that the waveguide is unbroken in the side coupling scheme, the output will be dependent on both the cavity energy and the energy in the waveguide, giving rise to Equation (4) and the green line in Figure 3.

$$s\_{out}^c(t) = u \frac{a\_\mathcal{\epsilon}(t)}{2} + s\_{in}^c(t) \tag{4}$$

The coupling factor, *<sup>u</sup>*, is defined as *u =* <sup>√</sup>2*γ<sup>c</sup>* [21]. In the direct coupling scheme the output will only rely on the energy in the cavity due to the carrier pump, such that:

$$s\_{out}^c(t) = u \frac{a\_c(t)}{2} \tag{5}$$

The output of this is shown by the solid purple line in Figure 3. The control light has a sinusoidal wave form with a repetition frequency in the GHz range, as shown in Figure 3. This leads to different output dynamics for the carrier light, depending on which scheme is used. In this paper, the PhC-based IMs are cascaded in both series and parallel in order to generate an OFC, as shown in Figure 4. When two IMs are in parallel with each other, it is assumed that the carrier light is split evenly into each arm. Each arm is then coupled to a cavity via the side coupling scheme, as shown in Figure 4 and when they recombine, the output is given by:

$$s\_{out}^{\varepsilon}(t) = u\frac{a\_{\varepsilon1}(t)}{2} + s\_{in}^{\varepsilon1}(t) - \left(u\frac{a\_{\varepsilon2}(t)}{2} + s\_{in}^{\varepsilon2}(t)\right) \tag{6}$$

where *sc1,c2 in (t)* represents the input to the bottom or top cavity and *ac1,c2 in (t)* represents the amplitude of the carrier signal in the bottom or top cavity. When two IMs are put in parallel with each other, as shown in Figure 4a, and the phase of the modulating wave in each arm can be controlled, it displays properties similar to that of a conventional IM. The PhC design, outlined in Figure 4a, consists of a waveguide leading into a Y junction symmetric splitter [22,23], where the pump and carrier light split evenly into each arm. An L3 cavity is coupled to each arm of the device to induce intensity modulation of the split carrier light. The two intensity modulated carrier lights then recombine via a Y junction to produce a dual modulated carrier signal. This technique of waveguide splitting and recombining in PhC waveguides is common when implementing PhC-based Mach–Zehnder interferometers [24,25]. Using this process, the shape of the carrier light will be determined by a number of parameters unique to this type of modulation, these include the phase shift of the pump light in each arm, the coupling efficiency of the light into the cavity and the Q factor of the cavity.

**Figure 4.** Schematics for the proposed photonic crystal (PhC) all-optical optical frequency combs (OFC) generators. (**a**) Two side coupled cavity modulators in parallel. (**b**) Two direct coupled cavity modulators in series.

This method of modulation can generate 3 comb lines with a maximum intensity difference of 1 dB on each side of the central carrier frequency, hence 7 comb lines all together as shown in Figure 5b. Additional comb lines exist at much weaker intensities further away from the initial carrier frequency. The number of generated intensity spikes in the frequency domain is limited by the attainable modulation profile of the carrier light in the time domain. As can be seen, the carrier wavelength has a much stronger intensity than the sidebands produced. Nevertheless, the homogeneity in intensity between the sidebands remains very uniform. The homogeneity of these sidebands is dependent on the parameters stated above and can therefore be optimized to generate a comb, as shown in Figure 5b. During experimental implementation of the proposed scheme, there will be disparity between the optimised parameters given here and the fabricated device. Although this will have an effect on the OFC quality, it is expected that a 7 line OFC can be observed for a broad range of parameter values. Figure 5a shows that for Q factors between 10,000 and 25,000 the side band intensity fluctuates to a maximum of 1.5 dB. This shows that a wide range of Q factor values will still achieve a 7 line OFC. The extinction ratio of the carrier signal is dominated by the coupling factor between the cavity and the waveguide. This will have an effect on the intensity of the generated comb lines relative to the intensity of the central carrier frequency. However, 3 sidebands on either side of the initial carrier frequency can still be generated for a broad range of coupling factors.

**Figure 5.** (**a**) Homogeneity between the sidebands against the cavity Q factor. (**b**) The generated OFC.

Figure 5a shows the effect the Q factor has on the uniformity with an optimum Q factor of 13,000. This is due to the shape of the carrier light through each arm when it recombines to generate the intensity modulated signal. Cavity Q factors in PhC L3 cavities can be greatly enhanced by altering the position of the holes surrounding the L3 cavity [26]. It is therefore possible to design the L3 cavity with the desired Q factor via fine tuning of the hole positions. The 7 comb lines observed have with a maximum fluctuation in intensity of less than 1 dB if the carrier frequency is not taken into consideration. In previous work, the intensity of the initial carrier frequency has also posed a problem in modulation-based combs [27] and micro-cavity-based combs, due to the high energy needed to induce nonlinearity within the micro-cavity [28]. A method to overcome this is to include a cavity-based notch filter. To obtain the required suppression of the carrier signal the coupling coefficient, size and loss coefficient of the cavity can be determined via calculation. The filter will induce an overall loss to the carrier wave signal of around 5 dB, but the flatness of the 7-line OFC will be greatly enhanced.

When two IMs are in series with each other, as shown in Figure 4b, the output from the first IM will be the input for the second. The equation for the energy inside the second cavity is given by:

$$\frac{da\_{\ell 2}}{dt} = (-i(\omega\_0 + \Delta \omega\_{p2}(t) - \omega\_{\varepsilon}) - \gamma\_{\text{total}})a\_{\varepsilon 2}(t) + \sqrt{2\gamma\_{\varepsilon}}a\_{\varepsilon 1}(t) \tag{7}$$

The output of the second cavity is equivalent to the output of the single cavity device, as given in Equation (5). In the parallel scheme outlined above, side coupled cavities were used to generate an OFC because a broad pulse was able to be produced from a sinusoidal control light. The high extinction ratio and pulse shape gained from this scheme meant a short, high intensity pulse could be produced and generate an OFC. However, in order to produce an OFC from two IMs in parallel, the sinusoidal pulse shape at the output of the IM is more desirable. Sidebands are produced in the frequency domain from the initial carrier light by a single IM with a sinusoidal pump light profile [29]. This is shown in Figure 6a where a sideband either side of the initial carrier frequency has been generated. The sideband intensity is within 5 dB of the carrier intensity and their FSR is in direct relation to the modulation frequency. The sidebands produced from the first IM then act as carrier waves leading into the second IM. The second IM has a modulation frequency one third that of the first, meaning that the generated sidebands will have an FSR one third that of the original FSR. Two extra comb lines are generated on each side of the carrier frequency; they are the sum of the generated sidebands from both the initial carrier frequency and the initial sidebands generated from a single IM.

This scheme then generates multiple sidebands that are of equal distance from each other in the frequency domain and of similar intensity. In this scheme, 7 comb lines have been generated, Figure 6b shows the generated OFC, while Figure 6a shows sideband generation from a single IM, in which the FSR is directly related to the repetition rate of the modulating signal. In this case, this is within the microwave range, at around 10 GHz. A microwave signal is then generated from the beat note between the carrier and the sideband, this is a common tool in microwave photonics [30].

**Figure 6.** (**a**) The frequency response of the carrier wave after a single direct coupled intensity modulators (IM). (**b**) The generated OFC from two all-optical PhC direct coupled modulators in series.

#### **4. Conclusions**

To conclude, a new and novel approach to OFC generation using PhC structures is proposed. By integrating multiple PhC cavities, it is possible to observe OFC generation on the nanoscale. A model based on CMT is built upon to simulate the modulation of a carrier light through a PhC modulator using a pump light with a repetition rate in the RF range. By cascading multiple cavities in either series or parallel, sidebands in the frequency response of the carrier wave can be observed. By careful manipulation of the device parameters, these sidebands can be optimised to give a flat spectrum across 7 comb lines. Using the parallel-based scheme, the Q factor of both cavities can be specified to give a flat comb over the generated sidebands. However, the initial carried frequency will have a far greater intensity, which could limit the OFC operation. To overcome this, it is possible to use a specifically designed notch filter. Using a scheme with two cavities in series, it is possible to generate a flatter comb across all lines without the need for a notch filter. In this scheme, two separate RF lines are needed to drive each cavity separately, which will cause complexity within the system.

PhC devices are a mature enough research field that fabrication techniques are commonplace to many institutes; therefore, implementation of this technique experimentally is very promising. PhC devices based on cavities and waveguides are a well-established and rich research field which have the potential to enable on chip photonic networks. To this end, the device proposed here will contribute towards both the field of PhCs as well as expanding OFC generation techniques into a new domain.

**Author Contributions:** Conceptualization, H.F. and C.J.; methodology, K.C. and S.C.; software, H.F., K.C. and S.C.; validation, C.J., M.H. and K.C.; investigation, H.F.; resources, C.J.; writing–original draft preparation, H.F.; writing–review and editing, C.J. and M.H.; visualization, H.F and C.J.; supervision, C.J. and M.H.

**Funding:** This research was funded by EPSRC First Grant, the Royal Society Research Grant (UK) and NSFC Grant No. 61574138 and No. 61974131 (China).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
