Experimental. **Femtosecond pulsed laser.** The experimental setup to obtain the pulse time delay resulting from the interaction of light with a sample consisted of a femtosecond pulsed laser (Libra Ti:sapphire from Coherent) with 850 nm, 450mW and 1KHz of frequency, and a Single Shot Autocorrelator (SSA - From Coherent). The SSA has two intern paths: a fixed and a moving path. We placed the wings in both of them, as before the SSA. The positions for placing the samples are schematically presented in Fig. 8(a). The SSA was assembled at the laser exit, with a diaphragm for selecting the laser center beam. These measurements were made in order to obtain a signal of autocorrelation between the split laser beam, to measure the pulse width, and the correlation between one original laser beam and the other one that crossed the sample. When the laser beam crosses the sample, it changes the optical pathway, changing the pulse width (Fig. 8). This characteristic brings the opportunity to determine the measurement of time delay due to the interaction of light with a natural photonic structure. Due to the high precision required on this measurement, we developed a tested analysis, with different samples like glass blades. The processing of the signal is done with the intensity of the pulse autocorrelation. These data are modeled with a Gaussian fit, to obtain the FWHM of the signal. This is used to measure the pulse width, we used it to establish the calculation to determine the time delay.

As identified in Fig. 8(a), the original SSA configuration was maintained. We inserted the samples on the laser path before it entered into the SSA (C1) as well as the two paths inside it (C2 and C3). Figure 8(a) presents the three spots (C1, C2 and C3) where the samples were placed and 8(b) and 8(c) presents how the samples were inserted in the laser path. On pathway C1 (Fig. 8(a)), the sample is placed before the SSA. It means that the two split pulses that will interact with the autocorrelation measurement are identical, and there is no alteration or delay on it. Placing the sample on pathway C2, we are making the laser pulse cross the sample only on the moving pathway inside the SSA. The one used to delay or advance the pulse on the correlation measurement. When we place the sample on the C3 pathway, it is placed on the fixed pathway. The process for measuring the autocorrelation was based on turning the laser on in front of the SSA and moving the moving pathway with the micrometric screw to find the maximum autocorrelation intensity. With this point established, the screw was turned from this point until the minimum value for the autocorrelation, for both positive and negative sides of the maximum.

The autocorrelation intensity was collected for the experimental apparatus with and without the samples, so we were able to compare the effects of the presence or absence of samples on the autocorrelation profile, besides using the glass blades to compare. The curves obtained were modeled with a Gaussian profile. We used the full width and half maximum (FWHM) for determining the pulse width, following the manual instructions from SSA for calibrating [31]. The pulse duration was measured by the equation:

$$\frac{\varDelta x}{2} =\frac{c}{\varDelta t}\to \varDelta t=\frac{\varDelta x/2}{c}, \left(1\right)$$

Where \(\varDelta x\) is the FWHM of the Gaussian-modeled curve obtained by shifting the moving path, \(c\) is the light speed and \(\varDelta t\)is the pulse width in femtoseconds. After obtaining the pulse width for each case, it was needed to understand how the presence or absence of a sample would interfere with the autocorrelation profile. The pulse width without sample was subtracted from the measurements with the samples, so the effect of the presence of the sample on the autocorrelation profile was removed. Then, we compared this result with the obtained from placing the sample on C1 (where the sample was placed before the SSA and the laser pulses would autocorrelated with the same interference), to measure the pulse delay due to only the sample.

With the placement of the glass blades on the laser path we obtained a shift of the laser autocorrelation profile for more (placing on C3) and for less (placing on C2), without intensity reduction, but an intensity reduction on the placement of the blade at C1 (Fig. 2(a)). The pulse width of each case is presented, as is the difference between the pulse width of each laser path and without a sample (*ws*) on the third column of Tables 1–3. The fourth column is the absolute difference on the laser pulse width on paths C2 and C3 from the result of placing the sample at C1 (where the pulse, after interacting with the sample, will pass through the other two paths).

**Interferometry measurements.** The interferometry measurements were performed using a Michelson interferometer setup. It was composed by a diode laser (OZ Optics FOSS-21-3S-4/125-635-S-1) with emission centered at 635 nm and 1 mW, a cube beam splitter (Thorlabs BS029), and a CCD camera (Firefly MUFNVU-03MTM). The setup was made as presented in Fig. 9, where M1 to M4 are mirrors used to level the laser path and BS is the cube beam splitter. Samples were individually placed in the beam path, in order to have its transmitted light information interfering with the reference beam at the image sensor position.

These interference images were taken from the wings of specimens of *Q. gigas* (Brazilian cicada) and *G. oto* (glasswing butterfly). A frequency domain filter was used to remove some undesired frequencies from the original interference images. Then, the filtered images were analyzed on spatial domain by intensity profiles. To compare with the other measurements, a reference line position was defined on the interference pattern without sample (with uniform fringes) as a central line perpendicular to the fringes. Therefore, the intensity profile was observed for the other interference patterns at the same position. For each interference pattern of a sample, its intensity profile was plotted together with the reference measurement (Fig. 4).

Before performing the image analysis, an image process was implemented in order to filter some frequency noises from the original data. The main interference pattern (Fig. 10(a)) presented some secondary fringes, which were subtracted from its Fourier transformation (Fig. 10(b)- black dots) and we obtained the pattern on Fig. 10(c). The same procedure was performed with the original interference pattern obtained for each sample *Q. gigas* wing (Fig. 10(d-f)) and *G. oto* wing (Fig. 10(g-i)).

Theory. **Localization of light in modeled natural photonic disordered structure.** Anderson's localization was introduced in the pioneering theory studied in 1958 by Philip Warren Anderson [18]. He studies an electron involved in a potential with random energies caused by the disorder in the lattice, showing that the electron wave function is located at only a few sites at all times, provided the amount of randomness is large enough [10–12, 14–16, 19]. His work was awarded the Nobel Prize in physics.

The study of the photonic system has its bases on solid-state concepts for electronic systems. Concepts like band structures and band gaps were introduced. In photonic crystals, there are concepts of prohibited frequency regions or photonic band gaps (PBGs) [12]. The fundamental concept of photonic crystals was established due to the periodic spatial modulation of the material’s dielectric function ε(**r**). That function conforms to the crystal along with one, two or three spatial directions that define the dimensional characteristic. These arrays affect the propagation of electromagnetic waves in the photonic system [10–12, 14–17, 19, 21, 26, 27, 32]. The theoretical foundation of the study of photonic crystals is involved in the following master equation (equation of eigenvalues) [12]:

$$\overrightarrow{\nabla }\times \left[\frac{1}{\epsilon \left(\overrightarrow{r}\right)}\overrightarrow{\nabla }\times \overrightarrow{H}\left(\overrightarrow{r}\right)\right]={\left(\frac{\omega }{c}\right)}^{2}\overrightarrow{H}\left(\overrightarrow{r}\right), \left(2\right)$$

For the photonic case, if ε(**r**) is periodic, the eigenvalues of Eq. (1) are the allowed real frequencies of the system [12]. This theory is implemented in MIT Photonic Bands code (MPB) [28] and was useful in our calculations of PGBs.

**Kronig-Penney model and transfer matrix method.** In our work we are interested in joining the study and analyses of our analytical results and obtaining a relation of them with the experimental data obtained in the optical characterization of the natural photonics systems, and also can establish a comparison between the both results if possible. In order to develop this work, we introduce the 1D Kronig-Penney model and transfer matrix formalism [30–33]. This is because periodic systems are slightly affected by a disorder model with various kinds of disorder mainly related to numerical methods [33].

The influence of random imperfections is commonly present in real experiments [33], in this way it is important to make an analysis of these properties inside the system studied.. For our case is the real natural photonic structures inside the sample wings studied. In general way, these imperfections are originated from the variations in the medium parameters such as the dielectric constant, magnetic permeability, and barrier widths or heights (33), to studied the photonic system to establish the relation between PGBs and reflectance spectrum was implemented a Kronig-Penney model [33–35] in 1D ordered photonic system model by this mathematical expression:

$$R\left(\omega \right)=\left|\frac{ck\left(\omega \right)/ (\omega -1)}{ck\left(\omega \right)/ (\omega +1)}\right|, \left(3\right)$$

The transmittance calculation to the disordered photonic system model to the *M. cypris* wings starts with the initial structural values for the 1D ordered photonic system model obtained from the electronic microscopic image. Transmittance calculation was performed to an ordered 1D photonic system made of chitin-air with 105 layers. In the model, it is considered that the light comes from the left and is partially reflected in the first layer, while the rest is transmitted through the arrangement, and it is assumed that there are no losses (Fig. 5). All calculations were made with an initial lattice size at ∼0.63µm with the chitin section size dchitin = 0.51 µm and the air section size of dair = 0.12 µm, values obtained from the electronic microscopic image of the wing [8].

The transfer matrix can be generalized for the case of quasi-one-dimensional systems, which are finite in the direction perpendicular to the direction of propagation of the particle [33]. The transmittance calculation has been performed using the formalism of the transfer matrix **M** (Scheme in Fig. 11).

Transmittance calculation included structural random disorder in the width of the layers of chitin-air in the 1D photonic system of the *M. cypris* wings. It was developed as (i) 10% random variation in the air layers, (ii) random variation of the 10% in the chitin layers, (iii) 5% random variation in the air and chitin layers simultaneously, (iv) 10% random variation in the air and chitin layers simultaneously. In the modeled system, the alternate arrangements of both dielectric materials are identified by the refractive index’s values nair=1 and nchitin=1.55, which remain constants in all calculations developed. The transmittance of light through the random arrangement of intercalated layers of chitin-air depends on the natural photonic medium formed by layers in the system studied:

$$M=\frac{1}{2{n}_{2}}\left[({n}_{2}+{n}_{1}){e}^{i\phi } ({n}_{2}-{n}_{1}){e}^{i\phi } ({n}_{2}-{n}_{1}){e}^{-i\phi } ({n}_{2}+{n}_{1}){e}^{i\phi } \right], \left(4\right)$$

The introduction of the disorder in the modeled 1D nature photonic system is done by defining d1 like the width of the air layers and d2 like the width of the chitin layers, the mathematical expression that represents this array is:

$$Layer´s width=\left\{{(d}_{1}\pm {\delta }_{1})+({d}_{2}\pm {\delta }_{2} ) \right\}; -\varDelta b\le {\delta }_{i}\le \varDelta b. \left(5\right)$$