Modeling the Reflectance Changes Induced by Vapor Condensation in Lycaenid Butterfly Wing Scales Colored by Photonic Nanoarchitectures

Gas/vapor sensors based on photonic band gap-type materials are attractive as they allow a quick optical readout. The photonic nanoarchitectures responsible for the coloration of the wing scales of many butterfly species possessing structural color exhibit chemical selectivity, i.e., give vapor-specific optical response signals. Modeling this complex physical-chemical process is very important to be able to exploit the possibilities of these photonic nanoarchitectures. We performed measurements of the ethanol vapor concentration-dependent reflectance spectra of the Albulina metallica butterfly, which exhibits structural color on both the dorsal (blue) and ventral (gold-green) wing sides. Using a numerical analysis of transmission electron microscopy (TEM) images, we revealed the details of the photonic nanoarchitecture inside the wing scales. On both sides, it is a 1D + 2D structure, a stack of layers, where the layers contain a quasi-ordered arrangement of air voids embedded in chitin. Next, we built a parametric simulation model that matched the measured spectra. The reflectance spectra were calculated by ab-initio methods by assuming variable amounts of vapor condensed to liquid in the air voids, as well as vapor concentration-dependent swelling of the chitin. From fitting the simulated results to the measured spectra, we found a similar swelling on both wing surfaces, but more liquid was found to concentrate in the smaller air voids for each vapor concentration value measured.


Introduction
The colors of butterfly wings [1][2][3][4][5][6] are generated by pigmentary (chemical) and structural (physical) factors. Both the pigments and nanostructures responsible for the color are located mostly in the wing scales. Structural colors are produced by constructive and destructive interference of the electromagnetic waves scattered on the photonic nanoarchitecture [7]. Photonic nanoarchitectures are a class of nanocomposites where light propagation obeys special rules. In the case of butterflies, their most important constituents are chitin and air. Optical properties of photonic crystal-type nanoarchitectures were first modeled by Yablonovitch [8] and John [9]. In its most general definition, a photonic crystal is a graded-refractive-index material, where the refractive index is a periodic function of the position in one-, two-, or three-dimensions, giving rise to a stop band [7]. Biological materials, however, always have a certain kind of disorder, and the disorder varies on a large scale among different butterfly species: from rigorously-ordered butterfly scales such as Parides sesostris [10,11] and Cyanophrys remus [12] to the lower end of the disorder scale as represented by seemingly amorphous structures such as those found in Albulina metallica [13].
A reversible change of the color and intensity occurs when the butterfly wing is subjected to vapors [14]. The first experiments were performed on Morpho-type butterflies [14][15][16], but vapor-dependent spectra have also been demonstrated on many other butterflies with different photonic nanoarchitectures [17,18]. The spectral changes are vapor selective, i.e., not only do they depend on the vapor concentration, but also on the kind of vapor. The selectivity was later explained, for the case of the Morpho butterflies, by a polarity gradient along the open Christmas-tree-like photonic structure of this butterfly species [19]. The measured change in the spectrum due to the vapors was small, but not as small as would have occurred by merely changing the contents of the air voids with the vapor, because the refractive index difference between the vapors and air was quite small [20]. One can reproduce, however, the magnitude of the color change by supposing that the vapor is partly condensed into liquid, i.e., the air voids become partly filled with liquid [20]. This phenomenon can be explained by the condensation effects, BET (Brunauer-Emmett-Teller) condensation [21], and the capillary condensation [22][23][24] process, which result in the vapor being condensed into liquid in the small voids below the saturation vapor pressure in open space. The amount of condensed vapor depends not only on the vapor concentration, but also on the size and shape of the air voids and the temperature. As additional proof of the condensation process, we demonstrated that it was possible to create a reversible color change by local cooling and warming of the wing [25]. That study also revealed that the rate of change of the color upon cooling depends on the openness of the structure, which shows that the vapor infiltrates those structures faster, where there are more channels between the air spaces, as seen in transmission electron microscopy (TEM) and scanning electron microscopy (SEM) images. Morpho-type butterflies [14][15][16], are known to have an open photonic nanoarchitecture, where the empty (air-filled) parts of the structure are easily accessible to the vapor. Vapor-sensitive color change was later shown to occur for the case of the closed photonic nanoarchitecture of the Hoplia coerulea beetle [26]. The scales of this beetle have a porous multilayer structure enclosed by a 100 nm-thick encasing envelope, termed as a "photonic cell". In order to infiltrate the air voids, the vapor first has to penetrate the chitin walls. An important difference between the vapor-sensitive spectral response of open and closed structures is that the intensity of the main spectral peak always decreases for the open structures with increasing vapor pressure, but can increase for the case of closed structures. This was explained as an effect of the increased index of refraction of the walls between the air voids, as a result of the liquid filling of the air channels inside the walls [26].
In the present work, we report on simulations of the ethanol vapor concentration-dependent optical spectra on both the dorsal and ventral wing surfaces of the Albulina metallica butterfly and compare these results to measurements [27]. The Albulina metallica butterfly is exceptional because both of its wing surfaces possess structural color and the wing scales of both its dorsal (blue) and ventral (gold-green) sides have a similar photonic nanoarchitecture, but the characteristic sizes-distances and volumes of air voids-are different [28]. This makes it possible to study the effect of these structural parameters on the vapor-dependent spectra. Recently, we tested ten different volatiles in the 0-50% vapor concentration range on whole wing pieces of this butterfly, as well as some of the volatiles on single scales in both reflected and transmitted light [27]. Chemically-selective responses were obtained on both wing surfaces, but the spectral changes were different on the two sides. In our present work, we attempted to create a model to explain these findings. Based on the statistical analysis of the transmission electron microscopy (TEM) images, we built a simple theoretical model of the photonic nanostructure. Subsequently, the vapor-dependent optical spectrum was calculated by exactly solving the 3D Maxwell equations on the model structure utilizing the finite element method (FEM) [29]. We considered three parallel processes: (i) vapor condensation in the air voids, (ii) vapor condensation in the channels inside the chitin walls, and (iii) swelling of the structure. Chitin is known to exhibit a substance-specific swelling in various solvents [30]. First, we calculated the spectrum of a single scale, then by a simple model of the melanin absorption and the statistical distribution of scale directions, we calculated the reflectance spectrum of the whole wing.
The organization of the paper is as follows. Sections 3.1-3.3 contain a detailed analysis of the SEM and TEM images of the Albulina metallica butterfly scales. In Sections 3.4-3.6, we build a model based on the TEM information and calculate the vapor-dependent reflectance spectra. In Section 3.7, we compare the measured and calculated results. Section 4 is devoted to the discussion of the results.

Materials and Methods
The butterfly samples used in the present work were obtained from the collections of the Hungarian Natural History Museum, Budapest, Hungary. All of the specimens investigated were males.

Optical Measurements
The wings of the Albulina metallica males were measured using an optical spectrophotometer (Avantes AvaSpec-HS1024x122TEC, Apeldoorn, The Netherlands). For illumination, an Avantes AvaLight-DH-S-BAL balanced light source was used. The illumination and detection angles were nearly perpendicular to the wing surface, and a slight modification (<5 • ) was used to maximize the reflected intensity.
The vapor sensing measurements were conducted by combining a vapor mixing setup with a home-built vapor-sensing cell and the spectrophotometer [31]. The butterfly wings were mounted inside the aluminum cell, which was covered with a quartz window to provide UV transmittance. The vapor concentration was set using computer-controlled mass flow controllers that were connected to the artificial air source and the bubblers containing the test liquids. During the vapor sensing measurements, we changed the concentration of the test vapor while monitoring the spectral variations in time. We applied concentrations in the 5-50% range, with 5% steps. The details of the procedure were published recently [27].
The wavelength shift of the main visible peak was 4.5 (6.5) nm for the dorsal (ventral) side for the maximal 50% vapor concentration, and the shift was a linear function of the vapor concentration in the above-mentioned concentration range.
When examining vapor-dependent optical spectra, it is practical to introduce relative reflectance spectra. This makes it easier to visualize the small spectral changes caused by the vapors. In principle, a relative spectrum can be obtained by calculating the ratio of the spectrum with a certain vapor pressure and of the spectrum with zero vapor pressure. In practical measurements, however, relative spectra are obtained directly from the spectrometer, by letting the spectrometer store the zero vapor pressure spectrum as a reference and measuring vapor-dependent spectra relative to that [31].

Electron Microscopy
Scanning electron microscopy (SEM) images were taken using an LEO 1540 XB (Carl Zeiss AG, Jena, Germany) microscope on wing pieces attached with conductive tape without any preparation. The TEM samples were prepared by incorporating the wing pieces in plastic blocks, followed by ultramicrotome sectioning, which resulted in 70 nm-thick slices. The samples were examined in TECNAI 10 TEM.

Finite Element Calculation
Reflectance spectra were calculated by numerically solving the three-dimensional Maxwell equations by the finite element method (FEM) [29]. We applied Floquet boundary conditions in the lateral X and Y directions. The FEM calculation was performed by the Wave Optics Module [32] of the Comsol Multiphysics [33] software on a uniform wavelength grid ranging from 200-800 nm, in 2-nm steps. The refractive index of chitin and ethanol was n chitin = 1.56 + 0.033i and n ethanol = 1.36 in the calculation. The small imaginary value was introduced to account for the pigment content of the scales of Albulina metallica [27]. The calculation was done for a normal backscattering arrangement, when both the illumination and the detection are normal to the sample surface. Figure 1 shows the top view SEM and cross-sectional TEM micrographs of a ventral wing scale of an Albulina metallica male. The dorsal side SEM and TEM images (not shown) were very similar to their ventral counterparts; only the characteristic sizes were somewhat smaller. A detailed comparison of the dorsal and ventral micrographs can be seen in [13,28]. The smaller length scale on the dorsal side corresponds to the color of the two wing surfaces: the dorsal wing surface of the Albulina metallica butterfly is blue, while the ventral wing surface is gold-green. No apparent order is seen in the arrangement of the small air voids on the top-view SEM image (Figure 1a), but as shown in [28], the two-dimensional radial distribution function (2D RDF) of the center of the air voids showed a marked first neighbor shell on both the dorsal and ventral sides, confirming the presence of a radial short-range order and the lack of angular order (a weak angular order could be seen on the ventral side). The cross-sectional TEM images (Figure 1b) showed, however, a layered structure; there were 3-5 pairs of layers where layer "A" is a perforated chitin sheet containing air voids and layer "B" is apparently solid chitin. The scale shown in Figure 1b has four "B" layers and three "A" layers; see Figure 2 for notation. This means that the 3D photonic nanoarchitecture of Albulina metallica is a 1D + 2D structure, where the vertical structure is a finite multilayer, i.e., a 1D photonic crystal (a "1D photonic cluster", meaning a periodic arrangement of a few layers), and the lateral structure is a quasi-random arrangement of air voids possessing a radial short-range order. According to our measurements on the SEM and TEM images [28], the vertical layer periodicity p ab = d a + d b and the average lateral distance of the air voids had a similar value on both sides, and it was 200 nm and 260 nm on the dorsal and ventral sides, respectively. d a (d b ) denotes the thickness of the "A" ("B") layers, respectively.

SEM and TEM Images
The top-view SEM image ( Figure 1a) shows that the topmost chitin layer ("B" layer) was porous, as we could "see into" the air voids of the topmost "A" layer (the pores are the darker spots on the SEM image). The size of the holes changed randomly; in the case of the larger holes, we could see through the "A" layer and see the second "B" layer, which was also porous. In the case of the TEM micrograph (Figure 1b), we could also see discontinuities in the chitin layers (the "B" layers) separating the perforated layers (the "A" layers), and also when the chitin layer seemed to be continuous, its gray level was not constant. Note that the TEM image is a cross-sectional image, but the slice is not infinitely thin; the TEM image is a kind of projection image through a 70 nm-thick slice, and the 70-nm thickness is in the same range as the size of the smaller pores. Hence, if this slice incorporated a part of an air void, the corresponding pixels became a lighter gray in the TEM image. These findings show that there are vertical channels present in this structural element, i.e., those connecting air voids in adjacent "A" layers through a "B" layer. Similar SEM and TEM images were seen on the scale structures of several polyommatine butterflies [34]. The TEM image also revealed the presence of lateral channels (those connecting adjacent spheres within the same "A" layer) because some of the air voids could be seen to be connected.
In the next subsection, we more carefully analyze the TEM image and thus build an appropriate model structure of the photonic nanoarchitecture filling the scales of the Albulina metallica butterfly.  [13]); the dark holes correspond to the air voids. (b) Cross-sectional TEM micrograph of a cover scale (from [28]). The light gray regions correspond to the air voids. The five triangular downward protrusions are the cross-sections of the ridges. The scale is displayed facing down, i.e., the ridges are pointing downwards to symbolize that it is a ventral scale. Note the edge of another scale in the upper right corner of the micrograph, which is a ground scale. The SEM and TEM images do not show the same scales.

Direct Space Averaging Method
Previously, we developed a direct space averaging (DSA) algorithm [28] in order to study the order-disorder effects in butterfly wing scales. The method was based on averaging the local environments of the repetitive units. When applied to the top-view SEM images of different butterflies, it was able to differentiate structures with a long range-, medium range-, and short-range order. In the case of the SEM images, the DSA analysis involved three consecutive steps: (i) find the middle of the holes; (ii) cut out a rectangular sub-image centered on each hole; and (iii) calculate an averaged grayscale image from these sub-images.
In this section, we further developed the DSA method in order to study the cross-sectional TEM images. As already noted in Section 3.1, the gray level of the "B" layers was not homogeneous in the TEM image, which indicates the presence of vertical air channels crossing the chitin layers. We wanted to determine whether there was a correlation between the position of the vertical air channels in the "B" layers and the position of the air voids in the "A" layer. In our original DSA algorithm, we simply cut out sub-images centered on each hole and averaged them. Now, we intended to do the same for the air voids of the middle "A" layer of the TEM image (A 2 layer in Figure 2a), but we had to compensate for two additional factors: (i) the layers were not perfectly flat; and (ii) the first neighbor distances of the adjacent air voids varied. Now, we intended to "transform out" both of these factors because we wanted to study not only the middle "A" layer itself, but also the relative position of the air spaces between the different layers, which include: (i) air voids in the "A" layers, (ii) horizontal channels in the "A" layers between the air voids, and (iii) vertical channels in the "B" layers. Hence, the enhanced algorithm of the DSA method is as follows: • FOR each hole in the middle "A" layer (A 2 layer) -FIND the center of the hole: r hole and the center of its left and right neighbors: r le f t , r right .
-CALCULATE the length and angle of the line section connecting the left and right neighbors, -SHIFT the image with r shi f t = − r hole , so that the present void is shifted to the center of the image. -ROTATE the image with −α, so that the line section connecting the left and right neighbors are horizontal. -SCALE (enlarge or shrink) the image so that the length d (distance of left and right neighbors) are the same on each image.
• AVERAGE gray values of all shifted images (pixel-by-pixel).
Thus, we cut out sub-images centered on the holes, rotated and scaled them, and then averaged them. The technical details of the procedure are the same as those in [28]. The centers of the holes are shown by red dots in Figure 2a. The leftmost and rightmost dots were excluded from the averaging because each void needed to have a left and right neighbor in the algorithm.

Results of TEM DSA Analysis
The resulting DSA image is shown in Figure 2b and its linecuts in Figure 2c. The DSA image is oriented in the same way as the original TEM micrograph (Figure 2a), the wing surface facing down, so this orientation symbolizes that it is a ventral scale. Its most prominent features are the horizontally-aligned three large white spots in the middle. The central white spot is the average central air void, and the other two white spots are its left and right neighbors. The two light gray horizontal bands on the bottom and top of the image are the average of the two "A" layers (A 1 and A 3 ; cf. Figure 2b) adjacent to the middle "A" layer (A 2 ). The gray level of these bands was nearly homogeneous, though some structure was present. There was noise due to the small number of averaged sub-images. While in the case of typical SEM images, there were several hundred holes, we only had 36 air voids in the TEM image in Figure 2a. This is because we need special, smaller magnification, but high resolution TEM images in order to produce a better statistics. Such measurements are currently underway. A weak periodic gray level variation was seen in the upper band (A 3 ), and we highlighted the maxima in Figure 2b with the dotted circles. These circles show that there was a weak correlation between the air void positions in the middle and upper "A" layers (A 2 and A 3 ), and the air voids in the upper layer tended to be located at sites above the walls separating the voids in the middle "A" layer (A 2 ).
The red and blue horizontal lines in Figure 2b show the positions of the two averaged "B" layers (B 2 and B 3 ), adjacent to the middle "A" layer (A 2 ). The red and blue linecuts, shown in Figure 2c, were taken along these lines. These linecuts show a marked quasi-periodicity, which indicates that there were air channels in the "B" layers and that the position of these air channels in the "B" layers was not random relative to the air voids in the "A" layers, but there was some correlation between the position of the air spaces in the adjacent "A" and "B" layers.
These findings led us to build the model structure; see Section 3.4 for details.

Building the 3D Model of the Scale
The TEM image analysis shown in the previous section (Section 3.3) indicates that the wing scales of the Albulina metallica butterfly had a 1D + 2D structure, meaning that the photonic nanoarchitecture is composed of A-B layer pairs, and in fact, both of the layers themselves have a porous structure; only the amount of porosity is different in the "A" and "B" layers. Based on the fact that most of the air voids in the "A" layers were completely white in the TEM images, we can conclude that the diameter of these air voids was larger than the TEM slice width (70 nm). Previously, we analyzed TEM slices cut in different directions [28] where the cutting planes were always perpendicular to the scale surface and rotated around a line perpendicular to the scale surface. That analysis showed that the air voids were spherical and had an isotropic arrangement inside the "A" layers. TEM DSA analysis (Section 3.3) indicated the presence of horizontal and vertical channels connecting these spherical airspaces. The average diameter of those channels was smaller than the diameter of the air voids, because the channels appeared gray in the TEM and TEM DSA images (cf. Figures 1 and 2b).
The model structure (Figure 3) utilized in this paper was constructed based on the findings summarized in the previous paragraph. First, we measured the average sizes and distances of the air voids in the "A" layers, as well as the d a and d b thicknesses of the "A" and "B" layers in the SEM and TEM images, then we inserted horizontal and vertical air channels. The diameter of these air channels was determined by calculating the main visible reflectance spectral peak wavelength with the algorithm given in [35] and fitting its value to the measured peak position by varying the diameter of the air channels as a parameter. This calculation was based on the so-called perforation factor. The perforation factor P is the ratio of the total volume of the air voids in the structure and the total volume of the structure, P = V air /V total = V air /(V air + V chitin ). With this, the effective refractive index is n e f f = Pn air + (1 − P)n chitin . We performed this analysis for the "A" and "B" layers separately and thus obtained the effective refractive indices of the two layers, n a and n b . The main visible reflectance spectral peak wavelength [7] is given as λ main = 2p e f f , where p e f f = n a d a + n b d b is the vertical effective layer periodicity.
Our geometry (Figure 3) was similar to the one we constructed earlier [36], where the "A" layers contained air voids and chitin "B" layers separated the "A" layers, but this time, we also inserted horizontal and vertical air channels. All air spaces were supposed to have a cuboid shape, though the real air voids and channels are spherical and helical, i.e., bounded by curved surfaces. This approximation makes the FEM spectrum calculations less demanding. We performed, however, a few test calculations with spherical air voids and tube-shaped channels constructed so that their volume was equal to our cuboid air voids and channels. The optical spectrum was practically unchanged in the visible range; however, considerable changes were seen in the UV range. Since in this paper, our main focus was on the visible part of the spectrum, the cuboid approximation was sufficient. It is not an easy task to demonstrate the 3D structure, because the air voids are inside the chitin. Hence, Figure 3 only displays the air spaces, rendered semi-transparent, without the chitin. The sizes of the model structure are given in the vertical and horizontal cross-sections shown in Figure 4 and in Table 1, where d a and d b are the thickness of the "A" and "B" layers; p z = d a + d b is the vertical periodicity (height of one story); p xy is the lateral periodicity; w a and w b are the thicknesses of the horizontal and vertical channels; and c is the length of the horizontal channel. The length of the vertical channel was identical to the "B" layer thickness, d b . The {v x , v y , v z } dimensions of the air voids could be easily derived as v x = v y = v xy = p xy − c, v z = d a .
Another major approximation we applied here was that our model structure was perfectly ordered, even though the real biological structure is always disordered. We address the validity and the ramifications of this approximation in Section 4, but the main conclusion is that for a normal backscattered configuration, the spectrum was unchanged in the first order when we introduced a lateral randomness into the positions of the air voids.

Three-Dimensional Finite Element Calculation of the Vapor-Dependent Reflectance Spectrum of a Scale
Reflectance spectra for the normal backscattered configuration were calculated by numerically solving the three-dimensional Maxwell equations by the finite element method; see Section 2.3 for details. The resulting reflectance curve is shown in Figure 5 by a solid line.
When the butterfly wing is subjected to vapors, liquid is condensed in the air voids. By increasing the vapor concentration, the amount of the condensed liquid also increases. The amount of condensed liquid also depends on the size of the air spaces, as the smaller air spaces are filled first. Therefore, according to the investigation of the SEM and TEM images in Sections 3.1 and 3.2, the air spaces in the "B" layers are smaller than the air spaces in the "A" layers, so we assumed that with an increasing concentration, the "B" layers are filled first before the "A" layers. Hence, if the filling factor (ratio of liquid filled volume to the total air volume) is f A for the "A" layer and f B for the "B" layer, then we can introduce, with a linear transformation, the constants c B and c AB , where f A = c AB and f B = c AB + c B . Here, c B corresponds to the process when the air spaces are filled in the "B" layer only, and c AB corresponds to the process when the air spaces are filled in both the "A" and "B" layers. The third theoretical possibility, when only the "A" layers are filled and the "B" layers are empty, does not have a physical significance because the voids in the "A" layers are larger than those in the "B" layers; hence, the voids in the "B" layers are filled first. When the wing is exposed to vapor, the liquid does not condense evenly in all air voids [37]; for a given concentration, when the smaller voids are already filled, the larger ones may be still empty, or only partially filled. Given that for small liquid concentrations, the average lateral distance of the liquid filled voids is larger than the coherence length [38] of white light, the reflectance components from the air-filled and from the liquid-filled voids are added incoherently. Hence, in our model of the vapor-dependent spectra, we calculated the reflectance of the partially liquid-filled nanoarchitecture as a weighted average of the spectra calculated for completely air-filled and completely liquid-filled configurations. First, we calculated the reflectance by assuming that all of the air spaces were filled with air, which is the R 0 (λ) spectrum (solid curve in Figure 5). Then, we repeated the calculation for the B filling process, which is the R B (λ) curve (dashed curve in Figure 5), and for the AB filling process, which is the R AB (λ) curve (dotted curve in Figure 5).
We also accounted for the swelling effect of the chitin matrix [30] in our calculation. Assuming an isotropic enlargement of the photonic nanoarchitecture by a factor of s, the swelling shifts the spectrum to larger wavelengths.
Summing up, the reflectance spectrum is given by:

Reflectance Spectrum of the Whole Wing
So far, the calculation provided the normal reflectance spectrum of one scale. In order to approximate the reflectance spectrum of a whole wing, two other important factors have to be accounted for: (i) the effect of the cover and ground scales and (ii) the effect of the random angular distribution of the scales.
(i) Lycaenid butterflies almost always contain absorbing pigments, mainly melanin [39]. In most cases, the butterfly wing contains several layers of scales [40], and the melanin is mainly concentrated in the ground scales [5]. Moreover, melanin can be eumelanin or pheomelanin, which have different absorption spectra [41]. Hence, in order to create a realistic optical model of the butterfly wing, we had to insert the wavelength-dependent absorption of the melanin [41][42][43] into the model. In our approximation, this was done by applying a linear approximation. Since we have no direct experimental data about the kind and distribution of melanin in this photonic nanoarchitecture, we relied on the experimental reflectance spectra. At the zeroth level, we added a small constant imaginary part to the refractive index of chitin for the scales containing the photonic nanoarchitecture. Assuming a wavelength-independent absorption is an oversimplification; hence, at the first level of approximation, we added a wavelength-dependent linear background to the reflectance spectrum calculated in Section 3.5. This calculation is the reverse of that applied [34] when analyzing the experimental whole wing reflectance spectra. The intercept and slope of the linear background line was obtained by fitting a linear function to the small-and large-wavelength parts of the experimental whole wing reflectance spectra.
(ii) The random angular distribution of the scales was accounted for by applying a Gaussian broadening to the spectrum. Figure 6 shows the effects of the three simulation parameters (s, c B , and c AB ) on the reflectance spectra. In order to reveal the effect of the parameters, we selected large parameter values, which do not occur in practical butterfly spectra. The insets show the spectra with the broadening already applied, but without the background correction (the background line is also shown in the insets). The left and right sub-figures show the direct and relative reflectance spectra, respectively.
• Figure 6a,b shows the effect of the s parameter, the swelling. The isotropic enlargement of the structure causes a red shift, R(λ) → R(sλ). The intensity of the peaks is unchanged (cf. Figure 6a inset), but after adding the background line, the intensity of the larger wavelength peaks also became larger (cf. Figure 6a). The red shift caused the appearance of a pair of negative and positive peaks in the relative spectra (cf. Figure 6b) on the left and right sides of the main peak wavelength. The amplitude of the negative peak was somewhat smaller, due to the background correction.
• Figure 6c,d shows the effect of the c B parameter. The intensity of the reflectance increased with an increase in c B because the ethanol filling increased the effective index of refraction of the "B" layer, and this increased the refractive index contrast of the "A" and "B" layers. The wavelength of the peak also increased, because an increase of the index of refraction increased the optical thickness of the structure. The interplay of these effects resulted in the left side of the reflectance curve being practically unchanged, but there was an increase in reflectance at the peak region and on the right side. This caused the appearance of a large positive peak in the relative reflectance spectrum (cf. Figure 6d) at the right side of the main peak wavelength.
• Figure 6e,f shows the effect of the c AB parameter. The main effect was the decrease of the reflectance, caused by the decrease of the refractive index contrast between the two layers. This in turn caused a decrease in the value of the relative reflectance (cf. Figure 6f) at the main peak wavelength with increasing c.  Figure 7a,b shows the measured dorsal and ventral relative reflectance spectra for increasing vapor concentration. Since the vapors only caused small changes in the spectra, especially for small concentrations, only the relative spectra are displayed. The dorsal and ventral spectra, apart from a red shift caused by the larger length scale of the ventral photonic nanoarchitecture (cf. Table 1) were similar: there was a negative (positive) peak at the left, bluer (right, redder) side of the main visual spectrum peak. The amplitude of the positive peak (at the redder side) was larger than the negative peak (that at the bluer side), however only in the case of the ventral spectrum. This indicates that the change of the spectrum was larger at the right, redder side of the main peak for the ventral side.

Comparison of the Measured and Simulated Reflectance Spectra
Next, we reproduced these measured relative reflectance spectra by fitting the three free parameters, s, c B , and c AB of the theoretical Equation (1) derived in Section 3.5, and we also applied the linear background (effect of the ground scales) and broadening (effect of the disorder) (Section 3.6). The fitting was performed by a least-squares fit by minimizing the integral of the square of the residuals: where "x" is the ethanol vapor concentration; λ min = 200 nm; and λ max = 800 nm. The residuals are linear in c B and c AB , but nonlinear in "s". For simplicity, we performed the fitting procedure only for the largest, 50% ethanol vapor concentration (x = 0.5) and obtained the values of the fitting parameters for the smaller ethanol concentrations by assuming that they depended linearly on "x". Figure 7c,d shows the calculated spectra for increasing "x" parameter values. It is perfectly possible to perform the fitting for each ethanol concentration separately, but as seen in Figure 7, the linear approximation sufficiently fits the measurements for these small concentrations. • The fitted parameters for the dorsal side were: c B = 0.08x, c AB = 0.08x, and s = 0.016x, where "x" is the ethanol vapor concentration, which means that c B = 0.04, c AB = 0.04, and s = 0.008 for the largest measured vapor concentration x = 0.5. By calculating the f A = c AB and f B = c AB + c B filling factors of the "A" and "B" layers (cf. Section 3.5), we obtained that for a x = 0.5 vapor concentration, there was a f B = 8% liquid ethanol filling of the air voids in the "B" layer and f A = 4% filling in the "A" layer, besides the 0.8% swelling.

•
The fitted parameters for the ventral side were: c B = 0.4x, c AB = 0, and s = 0.012x, which means that c B = 0.2, c AB = 0, and s = 0.006 for the largest measured concentration x = 0.5. Thus, we obtained that for a x = 0.5 vapor concentration, there was f B = 20% liquid ethanol filling of the air voids in the "B" layer and f A = 0% filling in the "A" layer, besides the 0.6% swelling.

Discussion
As we learned from the examination of the SEM and TEM images of the Albulina metallica butterfly wings (cf. Sections 3.1 and 3.3 and [13,28]), the photonic nanoarchitecture filling the wing scales had a 1D + 2D structure, i.e., it was composed of pairs of nearly flat chitin layers, and the layers themselves contained a quasi-random arrangement of air voids, while adjacent chitin layers were separated by layers where the principal component was air, with a small fraction of chitin (cf. Figure 1). Next, we examined the effect of the lateral disorder on the normal backscattered spectrum. Previously, we had performed goniometric spectrum calculations [28] for ordered and disordered lattices of scattering centers by utilizing the first Born approximation [44,45]. The results showed that in this approximation, the normal backscattered spectrum did not change by introducing a lateral disorder. The disorder only affected the angular dependence of the spectrum. This enabled us to model the photonic nanoarchitecture as an ordered model because we only analyzed the normal backscattered spectra in this paper.
Our model had three open parameters, c AB , c B , and s (the swelling factor), and the filling factors of the "A" and "B" layers depended on the c parameters linearly, f A = c AB and f B = c AB + c B . When the butterfly wing was exposed to a vapor of a specific substance with a specific vapor concentration and temperature, the detailed physical and chemical [19] interaction processes between the vapor and the photonic nanoarchitecture determined the specific values of the three parameters. As shown in [31], these physical and chemical processes are selective to the kind of the vapor; thus, the determination of the values of the c B (x), c AB (x), and s(x) functions (where x is the vapor concentration) helps to understand these microscopic processes. The different swelling factors associated with different substances [30] have an important contribution to the chemical selectivity. In the case of the experiments with ethanol for the Albulina metallica butterfly, the swelling was similar on the dorsal and ventral sides, but we saw a marked difference in the filling processes: the amount of the condensed liquid was larger at the ventral side, but it was confined only to the "B" layers, while both layers were filled on the dorsal side, though the filling of the "B" layers were somewhat larger. The negligible filling of the ventral "A" layers can be attributed to the large size of the air voids. The larger filling of the ventral "B" layers implies that the air channels in the "B" layers may be composed of several smaller diameter channels. The larger measured main peak shift on the ventral side (6.5 nm, as compared to 4.5 nm at the dorsal side) for the x = 0.5 vapor concentration can be easily explained by the larger vertical effective layer periodicity p e f f of the ventral scale structure. Indeed, the ratios of the peak shifted to the geometric thickness, 4.5 nm/200 nm = 0.0225 at the dorsal side and 6.5 nm/260 nm = 0.025 at the ventral side had a similar value, which shows that the larger peak shift on the ventral side was due to the larger size of the ventral structure.
A simplified model of the butterfly wing can be composed of the following structural elements: (i) cover scales (on the dorsal side); (ii) ground scales (on the dorsal side); (iii) wing membrane; (iv) ground scales (on the ventral side); and (v) cover scales (on the ventral side). The thickness of these objects are in the 1 µm range, and there are air gaps separating them. The distance between the above-mentioned structural elements is also in the µm range; hence, we did not expect coherent scattering between them, as coherent scattering may occur only inside the structural elements. Such incoherent scattering color calculations were performed for the case of vertebrate color patches [46], which are 3D structures that often contain multiple pigment types and structural features. In lycaenid butterfly wings, most of the melanin [47] is concentrated in the ground scales; therefore, considerable absorption only occurs there. In this paper, we modeled the melanin absorption by adding a linear background to the FEM calculated spectrum.
In the present paper, we only measured the normal backscattered spectrum and focused on the visible part of the spectrum. The fine details of the structure (e.g., the precise shapes and dimensions of the air voids) can, in principle, be revealed by measuring and analyzing angle-dependent spectra in the UV range. Such studies are currently underway.

Conclusions
Chemically-selective vapor sensors with an optical readout can be realized with the photonic nanoarchitectures occurring in butterfly wing scales [17,26] and by bioinspired sensors [48][49][50]. The Albulina metallica butterfly is exceptional because the wing scales of both its dorsal (blue) and ventral (gold-green) sides have a similar [28] 1D + 2D photonic nanoarchitecture.
We performed ethanol vapor-dependent relative reflectance spectrum measurements [27] on both the dorsal and ventral sides and analyzed the change of reflectance due to the vapor in the visible range. As we have recently shown [27] by principal component analysis (PCA), the structural colors of the two sides have different vapor sensing properties, and both the intensity and the wavelength can change. The change (increase) of the wavelength is larger on the ventral side, which can be explained by the larger air voids on the ventral side. The change of the intensity is different on the two sides because the tendency of the refractive index ratio is different.
In order to simulate these vapor-dependent relative reflectance spectra, we built a model, which was composed of four "A"-"B" layer pairs (corresponding to the number of layers observed in the TEM images). The "A" layers contain large air voids, which are interconnected by air channels, horizontal channels in the "A" layers and vertical channels crossing the "B" layers; these vertical channels are the air spaces in the "B" layers. When exposed to vapor, liquid is concentrated into the air voids, and the amount of concentrated liquid depends on the size and shape of the voids, as well as on the accessibility of the pores and the partial pressure. Therefore, the average liquid content of the "A" and "B" layers will be different for a given vapor concentration due to the different sizes of the air voids in the two layers. Next, we calculated the reflectance spectrum by exactly solving the 3D Maxwell equations with the help of a parametric FEM calculation. We found that the swelling of the nanoarchitecture was similar on the two wing sides, but the liquid filling of the "A" and "B" layers was different on the dorsal and ventral sides. This can be attributed to the different sizes and shapes of the air voids on the two sides.
It is worth pointing out that, to our knowledge, the qualitative agreement between the measured and computed response curves, when using lycaenid butterfly wings as chemical sensors, was achieved for the first time. The same structural model was able to reproduce the behavior of two distinct structures, one of blue and one of gold-green color, occurring on the wings of the same butterfly species. This is an indication that this type of model may be suitable for describing the behavior of numerous other butterfly scales of lycaenid butterflies possessing the so-called "pepper-pot"-type nanoarchitectures. Funding: This work was supported by the Hungarian OTKA K 115724 and OTKA K 111741. G.I.M. and G.P. wish to thank the Hungarian Academy of Sciences and the Belgian FNRS for financial support.
Acknowledgments: Helpful discussions with D.G. Stavenga are gratefully acknowledged.

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

Abbreviations
The following abbreviations are used in this manuscript: SEM Scanning electron microscopy TEM Transmission electron microscopy DSA Direct space algorithm FEM Finite element method UV Ultraviolet 1D One-dimensional 2D Two-dimensional 3D Three-dimensional PCA Principal component analysis