Mimicking bio-mechanical principles in photonic metamaterials for giant broadband nonlinearity

Microscopic structuring can change the effective properties of a material by several orders of magnitude. An example of this is animal bone, which has an effective elastic modulus that is more than 1,000 times larger than that of the constituent proteins. Here, we propose a broadband-enhancement principle of photonic nonlinearity that has a similar mathematical origin as the bone example. The proposed staggered array metamaterials violate the standard Miller’s rule in nonlinear optics and can enhance the third-order nonlinearity by more than a thousand to a billion times, depending on target operation frequencies. This metamaterial principle also enables manipulation of the individual components of the linear and nonlinear susceptibility tensors. Our biomimetic approach overcomes the fundamental speed-efficiency trade-off in current resonant enhancement schemes, making faster and more efficient all-optical devices possible for 1.55 μm wavelength. The principle is also applicable to ionic diffusion, heat conduction, or other transport problems. Biological materials such as bone show enhanced mechanical properties due to their specific structures, which can inspire new biomimetic materials. Here, a broadband metamaterial exhibiting giant optical nonlinearity is proposed, who’s properties share a mathematical origin with mechanically enhanced biomaterials.

N onlinear photonic processes, such as harmonic generation, phase modulation, and multi-photon absorption, have diverse applications including laser sources and optical communications 1-3 , three-dimensional lithography 4 , bioimaging 5 , and many others [6][7][8][9][10] . However, light-matter interaction is predominantly linear for common transparent photonic materials, thus necessitating the injection of high-intensity artificial light sources along the optimal crystallographic directions 11 . Therefore, numerous studies have been conducted with the aim to attain photonic nonlinearity enhancement. Most of these studies have focused on the introduction of tailored electromagnetic resonances near the operation wavelength to boost light-matter interactions [12][13][14][15][16][17][18][19] . Both quantum-electrodynamic approaches involving the engineering of intersubband transition energy levels 12,13 , and classical electrodynamic enhancement methods-mostly utilizing photonic cavities or resonant structures [13][14][15][16][17][18][19] -have been proposed. However, these resonant-enhancement schemes exhibit an inherent trade-off relationship between the degree of enhancement and the frequency bandwidth, thereby hindering their practical applications. Even for applications that do not require large bandwidths, large resonance-induced energy losses and a narrow fabrication error margin represent substantial obstacles. Moreover, only a few studies have investigated the controllability of the tensor form of nonlinear susceptibilities 10,13,18 , despite the fact that it may have important practical consequences.
Here, we achieve both extreme enhancement and tensor form controllability of the linear and nonlinear susceptibilities of a nonlinear staggered array metamaterial (SAMM) over a broad bandwidth. The SAMM represents a periodic composite system of nonlinear materials and conductors with a deep subwavelength period. We analytically and numerically verify that the thirdorder nonlinear susceptibility can be enhanced by more than 10 9 times at terahertz frequencies and by 10 4 times at optical telecommunication frequencies. In addition, we show that 18 different crystal classes with second-order nonlinearity can be realized based on simple variations of the spatial configuration of a given nonlinear material and a conductor. The proposed broadband-enhancement principle is mathematically analogous to that of naturally occurring mechanical SAMMs, such as animal bones and bone-like materials (e.g., nacre and sea shell), which have enhanced elastic modulus by more than three orders of magnitude compared to their constituent soft proteins 20,21 .

Results
Bespoke nonlinear susceptibility tensor of the proposed SAMM. The proposed SAMM structure is composed of alternating metal and nonlinear dielectric layers (Fig. 1a). The metal layers are composed of insulated metal patches in a planar lattice, and the nearest metal layers are staggered with half-unit-cell shifts in both the x-and y-directions (denoted as "A" and "B" layers in Fig. 1a). The metal and nonlinear dielectric layers ('I') are stacked in the sequence A-I-B-I-. The lateral unit-cell dimensions in the x-and y-directions, the lateral gap sizes in the same directions, and the thicknesses of the metal and dielectric layers are denoted as a x and a y , g x and g y , and h m and h d , respectively. Figure 1b depicts one unit cell of the SAMM.
The proposed microscopic (subwavelength) structure closely resembles the mechanical microscopic structure of mineralized collagen fibrils in animal bone and other hard biological materials (Fig. 1c) 20,21 . In a simplified model, the fibrils in the bones are an array of isolated hard mineral plates embedded in a soft protein matrix, and the mineral plates are in a staggered formation as in Fig. 1a. Figure 1d schematically shows the microscopic deformation of the mechanical SAMM under macroscopic tensile strain. The mathematical parallelism between the enhancement of mechanical and electromagnetic properties is presented in the discussion section.
From the point of view of materials, SAMMs can be considered homogeneous, crystalline photonic materials because their periods are much smaller than the wavelength-as those of natural crystals-and their effective material properties, such as the linear and nonlinear susceptibilities and refractive indices, are well defined and remain constant regardless of the thickness of the sample [22][23][24] . For these metamaterials, their linear (n = 1) and nonlinear (n > 1) homogenized electric susceptibilitiesχ ðnÞ eff can be calculated based on the reciprocity theorem 13,19,25  Using this approach, one can analytically derive the homogenized nonlinear susceptibility tensor of the proposed structure for second-harmonic generation (SHG). In this section, we assume for the sake of simplicity, the use of a metal with perfect electric conductivity and an insulator with frequencyindependent nonlinear coefficients, thus resulting in the Kleinman symmetry condition 11,26 These restrictions are eliminated in the rigorous analysis of the actual enhancements and numerical demonstrations presented in the following sections. For the proposed structure, the local electric field is mainly oriented in the z-direction inside the dielectric layers owing to the parallelplate capacitor-like geometry of the metal plates above and below the dielectric. Subsequently, the relevant field-enhancement factors are only F zx , F zy , and F zz , and other factors are being negligibly small. The z-directional local electric fieldenhancement factors under an x-directional macroscopic electric field, F zx , of a specific structure (a x = 4800 nm, a y = 800 nm, g x,y = 0.2a x,y , and h d = h m = 20 nm) is plotted in Fig. 1e-f (see Methods). The actual profiles of F zx and F zy are very simple with a large and almost spatially uniform value (around 120 in the aforementioned case) throughout the space between the top and bottom metal plates, and with a negligible value elsewhere. These maximum field-enhancement values are approximately |F zx | = 0.5a x /h d and |F zy | = 0.5a y /h d 27 . For the nonlinear dielectric, we assume that the microscopic pole directions are either the +zor −z-direction or unpoled, thus resulting in the intrinsic secondorder nonlinear susceptibility of the form χ 33 pðrÞ with p(r) = +1, −1, or 0, respectively. As a representative example, we consider a specific spatial pole configuration of p(r) = sign [F zy (r)]. Figure 1g-h show z-directional local second-order nonlinear enhancement factor under an x-directional macroscopic electric field, χ 33 j, of this example. The nonlinear enhancement factor is very large (around 14,000) and spatially nearly uniform throughout the aforementioned high-field region between plates. The contracted notation 11,26,28 of the effective nonlinear susceptibility simplifies tõ if g/a x and g/a y are small (see Supplementary Note 2 for detailed derivation). The tensor form in Eq. (1) corresponds to the monoclinic system with a single two-fold rotational symmetry. In a similar manner, one can realize various other nonlinear crystal systems by appropriately designing the spatial symmetry of p(r). If we also include the crystal designs in which the sides of the metal plates are intentionally arranged in a non-parallel orientation to any of the lattice vectors to break mirror symmetries, all of the noncentrosymmetric triclinic, monoclinic, orthorhombic, and tetragonal crystal classes can be realized (Fig. 2). If we consider the use of hexagonal plates in a hexagonal lattice, all of the noncentrosymmetric trigonal and hexagonal crystal classes are also covered (see Supplementary Fig. 1). Note that by tuning both the microscopic structure and the pole configuration in this simple manner, all noncentrosymmetric crystal classes-with the exception of the cubic classes-can be realized (a total of 18 point groups) and that only χ ð2Þ zzz (χ ð2Þ 33 in the contracted notation) of the constituent dielectric is utilized in all cases.
Extreme enhancement of nonlinear susceptibilities at various frequency ranges. The local electric field-enhancement factors of the proposed SAMM show three noticeable characteristics. First, the magnitude of the enhancement factors can become very large for high-aspect-ratio structures ( Fig. 1e-f, for example). The resulting nonlinear enhancements are even larger because they are products of multiple field-enhancement factors, depending on the order of the nonlinearity ( Fig. 1g-h, compared to Fig. 1e-f). Additionally, these enhancement factors are nearly frequencyindependent in the long-wavelength regime. A simple model, assuming perfect electric conductors and dispersion-less dielectrics, exhibits no frequency dispersion at all (F zx = 0.5a x /h d and F zy = 0.5a y /h d ), and a rigorous model with realistic materials also exhibits a small frequency dependence, except for the cases with extremely high aspect ratios. Third, the enhanced local field appears over a significant portion (50% or even more) of the total volume of the effective medium, as opposed to small "hot spots", as in previous tip-enhanced nonlinear structures. In structures with a hot spot near a sharp tip, or a narrow gap between metallic inclusions, the volume of the hot spots becomes smaller if one increases the degree of field enhancement by reducing the tip's radius of curvature or the gap size. Subsequently, after homogenization with volume integration, the effective nonlinear susceptibility of these structures is not as drastically enhanced as their fields. These are the key factors that differentiate this nonresonant, polarization accumulation scheme of extreme nonlinearity enhancement via space-filling of local electric dipoles from many previous approaches. The conceptual explanation regarding to the polarization accumulation and the large susceptibility enhancement are in Supplementary Note 3 and Supplementary Fig. 2.
The effective third-order nonlinear electric susceptibility tensor,χ ð3Þ eff , has an overall enhancement factor which is proportional to the fourth power of F zx or F zy . The longwavelength value ofχ ð3Þ eff for a SAMM with a x = a y = a becomes based on the simple model, and it satisfies the tetragonal symmetry (4mm) (see Supplementary Note 2 for detailed derivation). In the following section, we focus on χ eff;11 , which leads to an extremely large third-order nonlinear polarization along the x-direction under an incident x-directional electric field.
For more accurate calculations in the rest of the paper, we now adopt a rigorous analytical model, considering the actual complex permittivity of both the metal and the dielectric. This model is used to calculate the effective refractive index and the nonlinear susceptibilities (see Supplementary Note 4 for detailed   [29][30][31][32][33] . Because the magnitude of the metal permittivity is reduced as the frequency increases, the resonance frequency can be increased by lowering the aspect ratio. Therefore, one should design the structural parameters by considering the target frequency range, the amount of To reduce the computational load, we consider structures that are uniform in the y-direction (g y = 0) in the following calculations. χ ð3Þ eff ;11 for non-zero g y is also of the same order of magnitude as long as g y /a y ≪ 1. Figure 3b shows the enhancement factor (jχ ð3Þ eff ;11 =χ ð3Þ 33 j) for two different aspect ratios calculated by a nonlinear transfer-matrixbased parameter extraction method using nonlinear finitedifference time-domain (FDTD) simulation results. The values are compared to the analytical model predictions. For a moderate aspect ratio (a = 900 nm, g = a/5, h d = h m = 10 nm, i.e. a/h d = 90), we observed a million-fold enhancement of the effective nonlinear susceptibility compared to that of the constituent nonlinear dielectric (Fig. 3b, blue line), which is already two orders of magnitude higher than that of a previously proposed theoretical resonant structure 17 . Moreover, the enhancement is almost uniform from a near-zero frequency up to 2 THz, which is much lower than the electrostatic resonance frequency (5 THz). A SAMM with a larger aspect ratio (a = 6000 nm, g = a/5, h d = h m = 10 nm, i.e. a/h d = 600) shows more than a billion-fold enhancement (Fig. 3b, green line). Although there is an electrostatic resonance at a sub-terahertz frequency for the structure with this high aspect ratio, the enhancement remains extremely large from the zero to the terahertz frequency range. The enhancement factors from the FDTD simulations and the rigorous analytical model are in excellent agreement. The retrieved χ ð3Þ eff ;11 value remained constant regardless of the thickness and the input intensity (see Supplementary Fig. 3). Therefore, the extracted χ ð3Þ eff ;11 is well defined as an effective bulk material property of the SAMMs.
At higher frequencies, including the infrared and visible, the highest aspect ratio that one can use for broadband enhancement becomes lower owing to the smaller magnitude of the metal permittivities. Figure 3c shows the enhanced third-order nonlinearity of the proposed structure with a = 75 nm, g = 15 nm, and h d = h m = 5 nm with realistic materials at 200 THz, which is close to an important optical-fiber communication band. The electrostatic resonances occur at 67 THz and 200 THz. The χ ð3Þ eff ;11 enhancement is larger than 4×10 4 at 200 THz. In addition, the enhancement in the long-wavelength limit is larger than 10 3 , which is still a very high enhancement.
Thin and efficient nonlinear reflectors utilizing SAMMs incorporating nonlinear quantum wells. The proposed microscopic structure of the staggered conductors functions as a universal magnifier and symmetry-modifier of the intrinsic nonlinear susceptibility tensor of the constituent nonlinear material of choice, regardless of the latter. Accordingly, the quantitative enhancement factors shown in Eqs. (1) and (2) as well as the symmetry-changing capabilities demonstrated in Fig. 2 remain the same when we change the nonlinear materials. Thus, one can utilize quantum wells (whose sub-band energy levels are carefully designed to possess resonant-enhanced nonlinear susceptibilities that are several orders of magnitude larger than those of typical nonlinear materials in their off-resonance regime) as the nonlinear material between the conductors and the microscopic structure will further enhance these large nonlinearities, by another several orders of magnitude. This allows the use of metasurface types of nonlinear reflectors for massively parallel, pixel-by-pixel nonlinear optical operations (Fig. 4), rather than using waveguides or crystal-type devices with long propagation paths.
The nonlinear reflectors in Fig. 4a comprise a thin SAMM slab with GaN/AlN quantum well 34 , dielectric impedance-matching layers above and below 35,36 , and a metallic mirror (see Methods). The electromagnetic properties of the GaN/AlN quantum well is summarized in Supplementary Fig. 4. Figure 3b shows the switching performance of the optimized nonlinear reflector (a = 57 nm, g = 20 nm, h d = 5.1 nm, h m /2 = 10.2 nm, h u = 357 nm, and h l = 57 nm) for pulses which are 30 fs with a center wavelength of 1.55 μm. The thickness of the SAMM and the total thickness including dielectric layers are only of the order of 1/60 and 1/3 of the vacuum wavelength, respectively. While resonant-cavity-based all-optical switching approaches sacrifice speed and bandwidth in exchange for lower switching power [2], our non-resonant-enhancement scheme does not suffer from this limitation, which is apparent from the 26 fs output pulses (red curve in Fig. 4b inset). The recovery speed is ultimately limited by the relaxation time of the intersubband transitions in quantum wells, which can be of the order of 100 fs or faster 34,37 . The system exhibits a modulation in excess of 10 dB with a pulse energy utilization of only 380 fJ, while the reflectance is more than 50% at the on-state mode of operation. This 10 dB modulation pulse energy is a record-low value for all-optical switches with similar modulation speed [2], and is 250 times COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-020-0352-0 ARTICLE COMMUNICATIONS PHYSICS | (2020) 3:79 | https://doi.org/10.1038/s42005-020-0352-0 | www.nature.com/commsphys smaller than previous quantum-well-based all-optical switches 37 . We visually compare switching performance of the result in Fig. 4 with several state-of-the-art all-optical switches [37][38][39][40][41][42][43] in Supplementary Fig. 5 using a Pareto optimality plot based on the switching time and the required switching energy. The maximum local electric field magnitudes are of the order of 5 × 10 8 V m −1 for the 30 fs, 380 fJ pulses, which are below the damage threshold for many solid dielectric materials 44 . Although the effective susceptibilities are now dispersive owing to the resonant nature of the quantum wells, the bandwidth is still large enough to cover a significant part of the optical communication bands or to allow femtosecond-pulse operations. In addition, we verify that the designed device maintains more than 5 dB modulation depth with same pulse energy (380 fJ) under up to 30% (±10 nm) lateral miss-alignment (zero miss-alignment being the perfectly staggered case) between two metal layers of the SAMM slab, relative to the lateral unit cell size (Fig. 4c).

Discussion
Nonlinear susceptibilities of many natural nonlinear crystals can be estimated using Miller's rule, which states that the nth-order nonlinear susceptibilities (e.g., χ ð2Þ ð2ω; ω; ωÞ) can be approximated to an (n + 1)th-order linear susceptibilities (e.g., χ ð1Þ ð2ωÞχ ð1Þ ðωÞχ ð1Þ ðωÞ) with a constant factor, δ, known as the Miller's delta 11,19,45 . However, according to O'Brien et al. 19 , the nonlinear response predicted by Miller's rule can be several times larger or smaller than the measured values for composite media such as metamaterials. In case of the proposed SAMMs, the actual nonlinear susceptibilities (Eqs. (1) and (2) . It is confirmed that the modified Miller's rule with δ 0 ¼ δ½M d À3 ½M s À1 estimates third-order nonlinear susceptibilities of SAMMs well. While large effective oscillator strength enhancement factor M s can be found in many previous studies 16,17,31,33,46,47 , large effective oscillator density enhancement factor M d is a unique characteristic to SAMMs and induces a dramatic deviation from the conventional Miller's rule. In fact, our SAMMs can be considered as an example of a new category of electromagnetic media with designable oscillator densities (see Supplementary Note 5 for details).
The enhancement of the electric susceptibilities demonstrated herein can be viewed as a special case of a general enhancement principle of effective material properties that occur in two-phase composite materials. For example, the enhancement of the elastic   4 Record-efficient all-optical switch with a subwavelength thickness. a Concept of the functionality of metasurface type all-optical switch with staggered-array metamaterial (SAMM). Two-dimensional optical computation (e.g., AND operation) can be performed. A cross sectional view is also shown (not drawn to scale). b Intensity-dependent reflection is shown for a designed nonlinear reflector. Temporal profile of the input and output pulses are shown in the inset. c Dependence of the modulation depth on the lateral miss-alignment between two metallic layers. The input pulse energy is fixed at 380 fJ. modulus of the mineralized collagen fibrils in bone (Fig. 1c, d) can be conceptually explained as follows: First, the macroscopic tensile strain induces a shear strain enhancement by a factor of 0.5a/h d inside the protein layer, which results in an enhanced shear stress by a factor of 0.5a/h d . This enhanced shear stress accumulates over the distance of 0.5a and appears as tensional stress in the mineral plates with a thickness of h m , thereby inducing another enhancement factor equal to 0.5a/h m . Rigorous derivations of the effective mechanical properties can be found in previous reports 20,21 , and the expression of the effective elastic modulus is similar to Eq. (3) in Methods with n = 1 (linear case) and becomes M eff ≈ G p (a/h d ) 2 /8, where G p is the shear modulus of the protein matrix. This shows that the system can have an effective elastic modulus that is several orders of magnitude larger than that of the host material, for cases involving high aspect ratios (a/h d ≫ 1). This mechanical enhancement principle has been previously adopted for biomimetic hard and tough materials made of ceramics and polymers 48 . We note that this mechanism for enhanced elastic modulus is analogous to the explanation in Supplementary Note 3.
In general, we argue that the same principle can be encountered in various homogenization problems (including the examples already discussed) with a constituent equation of the form F = σE, where F is the relevant flux density (e.g., heat or mass), E is the conservative driving force (e.g., the temperature gradient or mass density gradient), and σ is the conductivity. E is an irrotational vector field owing to the conservative nature, and F becomes divergence-free in a steady state. Subject to these conditions, the effective conductivity for a two-phase composite material with isolated high-conductivity phases embedded in a continuous lowconductivity phase host material in a staggered formation, as in our SAMM, is σ eff = σ(a/h d ) 2 /8 (see Supplementary Note 6 and Supplementary Fig. 6 for detailed derivation). Thus, the effective properties of the composite materials can be many orders of magnitude higher than those of the continuous host, which may occupy more than half of the total volume, while the enhancement factor is purely geometrically controlled. Moreover, the nonlinearity of the conductivity is enhanced by even higher orders of magnitude in a manner identical to that of the enhancement of the nonlinear susceptibility, which might find use in nonlinear control of heat and mass flow, for example.
It has been shown that the symmetry and strength of individual components of linear and nonlinear susceptibility tensors can be tuned by the appropriate design of microscopic composites of nonlinear dielectrics and metals and that 18 nonlinear crystal classes can be realized. A near dispersion-less enhancement of χ (3) of the order of 10 6 was demonstrated for broad terahertz frequency ranges, and resonant enhancements of the orders of 10 9 and 10 4 were demonstrated for terahertz and near-infrared frequencies, respectively. Ultrafast, low-power, pixel-by-pixel alloptical modulation schemes at optical communication wavelengths was presented using a structure with a subwavelength thickness (λ/3). These SAMMs may serve as a useful platform for scientific studies of photonic nonlinearity with weak light or for high-density integration of nonlinear photonic devices. Moreover, the enhancement principle may find applications in other branches of science, while synergistic enhancements of multi-physical interactions are also imaginable.

Methods
Linear and nonlinear effective electric susceptibility tensors of metamaterials. The general form of the effective electric susceptibility tensors for a threedimensional metamaterial can be formulated in a similar manner to that of the two-dimensional case (metasurface), which was previously derived from the Lorentz reciprocity theorem 13,19,25 . In general, the nth-order nonlinear susceptibilities are (n + 1)th-order tensors. However, for nth-harmonic generation problems, the corresponding nonlinear susceptibilities can be represented with a more concise second-order tensor form using a contracted notation 11,26,28 , Local electric-field-enhancement factor plot in Fig. 1. A lossy Drude model is assumed for gold, which is the metallic material (plasma frequency: 2.18 PHz, damping frequency: 6.45 THz) 49 . Crystalline zinc oxide is assumed to be the dielectric 50 . The enhancement factor values are obtained from FDTD simulations (We used commercial FDTD software of Lumerical Inc. for all FDTD simulations in this paper). A three-unit-cell-thick SAMM slab is used in the simulations. The enhancement factor is measured at 2 THz, which is a sufficiently low frequency and can be considered to be in the long-wavelength regime.
Effective nonlinear electric susceptibilities calculation in Fig. 3. χ eff ; 11 is retrieved from the simulations using a nonlinear transfer matrix method. The complex amplitude of the third-harmonic wave is obtained based on nonlinear FDTD simulations, whereas the linear electromagnetic property of the SAMM is obtained from separate linear FDTD simulations. A dispersion-less nonlinear dielectric (ε 11 = ε 33 = 6, χ ð3Þ 33 = 10 −21 m 2 V −2 for 10 GHz to 100 THz) and a lossy Drude model for gold 49 are assumed as the constituent materials for Fig. 3a. Zinc oxide 50 (ε 11 = ε 22 = 7.8, ε 33 = 8.9, χ ð3Þ 33 = 10 −18 m 2 V −2 for 0.2 to 6 THz frequency) and the lossy Drude model for gold are assumed as the constituent materials, and a SAMM slab with two-unit cell thickness is used for the FDTD simulation for Fig. 3b. Zinc oxide 51 (ε 11 = ε 22 = ε 33 = 4, χ ð3Þ 33 = 10 −18 m 2 V −2 for 20 to 750 THz frequency) and aluminum 52 are assumed as the constituent materials, and a SAMM slab with a half-unit cell thickness is used for the FDTD simulations for Fig. 3c.
All-optical switch demonstration in Fig. 4. Owing to the limitation of the computational load, we obtained the results based on two-dimensional FDTD simulations (i.e. uniform in the y-direction in Fig. 4a). The pulse energy was calculated for a three-dimensional beam whose electric field profile in the radial direction is the same as that of the two-dimensional beam used in the simulations. The optical property of the GaN/AlN quantum well is calculated from two-level electron model 11,34 whose carrier density, excitation energy, transition dipole length, relaxation time, dephasing time, and background relative permittivity are 5 × 10 19 cm −3 , 0.8 eV, 0.5 nm, 100 fs, 10 fs, 5.37, respectively 34,37 . Calculated intensity-dependent permittivity of the GaN/AlN quantum well is presented in Supplementary Fig. 4. The dynamics of the material property of the quantum well is treated as two-level electron model which incorporate equation of motion of electron population in sublevels 53 in nonlinear FDTD simulations. Aluminum 52 and silver 52 are assumed to be the materials used for the metallic part of the SAMM and for the back-reflecting mirror, respectively. Silica is used for the upper and lower impedance-matching dielectric layers (relative permittivity: 2.074). A rigorous analytical model for the electric susceptibilities of SAMMs, a transfer matrix method, and a particle swam optimization are used for the structural parameter optimization.