1. Introduction

It is known for a long time that colors in nature are essentially produced in two different ways. Pigments contain molecules which absorb and/or transmit light to produce a color [1]. Structural colors, on the other hand, originate from nano- and microstructures that interfere with light [2]. Many interference and/or scatter effects have been discovered in the last century and eventually they led to the development of technical surfaces with fascinating optical features [3–8].

Understanding the actual optical principle of a structural color effect under review can be puzzling, especially, if the observed structure is more complex than a simple layered structure [9]. Consequently, most studies apply optical simulations to reveal the physical principle behind the observed optical effect. As the phenomenon to be uncovered is typically unexplored at the beginning of such studies, the best theory to describe it properly is unknown, too. Therefore, many simulation projects start with a simple approach like the Maxwell-Garnett theory [10] and combine it with the transfer matrix method [11]. Actually, this multi-layer technique works well in many practical cases (see, e.g., Refs. [12–16]) but not all, and important effects might be missed.

Consequently, advanced optical simulation methods are crucial for the analysis of lattice-like arrangements of nanoparticles because they might feature optical effects as observed in photonic crystals [17]. Such nanostructures are frequently found in nature and typically they consist of nanoparticles or air-voids in a surrounding matrix. Several birds, for example, have colorful feathers which are known to produce bright reflection through the lattice-like arrangement of nanoparticles [18–22]. The male peacock is a prominent case. It produces its colorful tail feathers with the help of 2D photonic crystals made of melanin rods connected by keratin, which are incorporated into the barbule [15,23].

Some reptiles utilize the principle of photonic crystals, too [24–27]. The most prominent one is, most likely, the male panther chameleon (Furcifer pardalis) which is famous for its ability to change its color. Teyssier and collaborators [26] showed that this famous dynamic color variation is caused by the change of distance between guanine crystals which are immersed in a cytoplasm matrix in the upper layer of the chameleon’s skin. In this way, it can change its color from blue to red by increasing the distance between the crystals. Another example are Phelsuma lizards [27]. These geckos feature skin colors which are generated via different physical origins. For example, the white or blue colors are caused by ordered arrangements of guanine nanocrystals in iridophores in the epidermis. Depending on the actual structure of the crystal, different colors are generated. Green skin is produced, either by structural green or by the interaction of structural blue with yellow pigments. Other colors, such as dark brown, are due to melanophores observed in the upper layers of the dermis.

Here, we show that the utilization of a rigorous electromagnetic approach like the Korringa–Kohn–Rostoker (KKR) method is well-suited to simulate and model the reflectance of such biological structures. In previous work, the KKR method has been successfully applied to study the color response of birds’ feathers [20,21]. This approach is inspired by the works of Korringa [28] as well as Kohn and Rostoker [29] which were originally intended for calculations in solid state physics. However, more than two decades ago, Stefanou et al. [30,31] introduced the software MULTEM allowing the straightforward optical simulation of photonic crystals based on the KKR method. Later, various extensions were developed to include coated spheres [32,33], quasi-periodic structures [34], non-spherical scatterers [35], and hybrid systems comprising metallic spheres [36], among others.

In this work we use the KKR method and its numerical implementation (MULTEM) to investigate the electromagnetic response of the skin of the panther chameleon [26] and of Phelsuma lizards [27] by proposing a simplified structure that captures the most relevant geometric characteristics of the tissue. We calculate the total reflectance and the efficiencies of the reflected diffraction orders that can be viewed from directions other than the specular. The KKR results reveal some features in the ultra-violet (UV) and short visible wavelengths region that might be overseen by less rigorous methods.

2. KKR method

The KKR method is an efficient rigorous formalism for the calculation of the reflectance, transmittance, and band structure of 3D photonic crystals. In particular, it permits to compute the individual efficiencies of the diffraction orders (defined as the power carried by each order normalized to the incident power) that propagate in directions other than the specular.

The photonic crystal can have a finite (in the $z$-direction) slab consisting of a number of unit slices (unit supercells), embedded in a homogeneous medium characterized by a real dielectric function and a real magnetic permeability. The unit slice consists of a number of different components, each of which can be either a homogeneous plate, or a multilayer of identical spherical particles in a homogeneous host medium. The multilayer is a stack of identical layers, each of which consists, in general, of a number of different planes of spheres with the same 2D periodicity with a particular type of arrangement, such as triangular or square, for instance. The dielectric function of the spheres is, in general, a complex function of the frequency, and therefore, either dielectric or metallic spheres can be considered. There is no restriction either on the size of the nanoparticles, nor on the filling fraction of the system, as long as the spheres do not overlap.

To solve the electromagnetic problem, the multiple scattering between spheres of each single layer is calculated first. Then, the scattered response of multiple layers is determined by using a transfer matrix approach, a procedure similar to the one used to calculate the reflection and transmission properties of stratified media with planar interfaces.

The details of the KKR method and of its numerical implementation MULTEM can be found in [30,31]. During the last years, we added several extensions to the original MULTEM program. For this study, we emphasize the new option to calculate the individual efficiencies of the diffracted orders. The code is available from the authors on a reasonable request.

Finally, it is important to mention that the computing times of the method are quite moderate. For example, we used a Linux machine with an Intel Core i3-6100U CPU @ 2.30 GHz with 12 Gb RAM, 4 cores and 8 threads. The calculations for 1001 wavelengths between 200 nm and 1200 nm need less than 5 minutes considering 37 reciprocal vectors. Since the results for each wavelength are independent, the calculations are straightforward to parallelize. However, the simulations become more elaborate when several curves have to be averaged and might take some hours.

3. Diffracted orders of nanocrystal lattices

As mentioned above, the KKR permits to compute the individual efficiencies of the diffraction orders. For shorter wavelengths, the number of diffraction channels increases, and therefore, it is interesting not only to estimate the total reflectance but also to calculate the power carried by the non-specular orders, which might significantly contribute to the optical response of the structure, mainly within the UV. This could have biological implications in how reptiles are seen by other animals, like birds for instance [37]. In this section, we use an analytical approach to calculate the condition for the appearance of diffraction orders and show when higher order channels open up.

When a regular arrangement is illuminated by a linearly polarized plane wave of wavelength $\lambda$ and wavevector $\vec {k}_{i}$ (see Fig. 1), diffraction orders other than the specular can arise [38]. This implies that, in general, the total reflectance (transmittance) includes the contribution of the efficiencies associated to all propagating orders. The number of propagating orders depends on the characteristics of the incident wave as well as on the geometrical parameters of the lattice.

 

Fig. 1. (A) Scheme of a generic system of a single layer of spheres indicating the incident wave vector $\vec {k}_{i}$ and the angles $\theta$ and $\phi$; $\theta$ is the angle between $\vec {k}_{i}$ and the $z$ axis and $\phi$ is the angle between the $x$ axis and the projection of $\vec {k}_{i}$ on the $x-y$ plane. The inset shows an example of a two-dimensional Bravais lattice with primitive vectors $\vec {a}_1$ and $\vec {a}_2$, of length $s$ and $b$, respectively, and $\gamma$ is the angle between them. (B) Scheme illustrating the diffraction orders.

Download Full Size | PDF

Since the complete structure is formed by a stack of identical Bravais lattices, in order to analyze how many reflected and/or transmitted orders propagate, it is sufficient to consider the arrangement of scatterers in a single plane. It is known that the primitive vectors of an arbitrary two-dimensional lattice can be written as $\vec {a_1}=s\hat {x}$ and $\vec {a_2}=b \cos \gamma \,\hat {x}+b \sin \gamma \, \hat {y}$, with $b$ being the length of $\vec {a_2}$, $s$ the spacing between first neighbors (lattice constant), and $\gamma$ the angle between $\vec {a_1}$ and $\vec {a_2}$. The two unit vectors $\hat {x}$ and $\hat {y}$ are orthogonal to each other and to the stacking direction. The reciprocal vectors of such a two-dimensional lattice are $\vec {b_1}=\frac {2\pi }{s}\left (\hat {x}-\frac {\cos \gamma }{\sin \gamma }\hat {y}\right )$ and $\vec {b_2}=\frac {2\pi }{b \sin \gamma } \hat {y}$, which satisfy the orthogonality condition $\vec {a_i} \cdot \vec {b_j}=2\pi \delta _{ij}$. Now, every vector $\vec {g}$ of the reciprocal lattice can be written as $\vec {g}=p\,\vec {b_1}+q\,\vec {b_2}$, where $p$ and $q$ are integers. Therefore, $\vec {g}$ can be expressed as

Here the +(–) sign corresponds to a transmitted (reflected) wave, respectively, while $k$ is the magnitude of the incident wave vector and $\vec {k}_{\parallel }$ is its projection along the incidence plane [30]. In spherical coordinates $\vec {k}_{\parallel }$ can be written as

As already pointed out in Ref. [30], only wave vectors $\vec {K}^{\pm }_{\vec {g}}$ with real $z$ component represent propagating diffraction orders. An inspection of Eq. (2) shows that this condition is satisfied for $k^2 \geq |\vec {g}+\vec {k}_{\parallel }|^2$. Using Eqs. (1) and (3), this inequality can be used to calculate the cut-off wavelength in vacuum ${\lambda _{\textrm{cut}}}$ for a respective diffraction order $(p,q)$ from

4. Numerical simulations

Now, we present our numerical simulations based on the KKR method and investigate the optical response of the skin of two reptile species which are well known for their respective structural coloration. As mentioned in the introduction, the male panther chameleon (Furcifer pardalis) is famous for its ability to change its color [26]. Structural colors are also observed in the tropical geckos of the genus Phelsuma [27]. In both cases, the electromagnetic response originates from structures composed of guanine nanocrystals (refraction index $n_g= 1.83$) immersed in a cytoplasm matrix ($n_c= 1.33$). We examined these two well-documented cases and applied the KKR approach to calculate the optical response of the reptiles’ skin.

4.1 Male panther chameleon

To study the mechanisms of structural color generation in the skin of a typical panther chameleon, we consider the structure schematized in Fig. 2, which is composed of a stack of identical parallel layers of regularly distributed spheres of radius $r$ that represent the guanine nanocrystals ($n_g = 1.83$) immersed in a cytoplasm matrix ($n_c = 1.33$). The whole structure is surrounded by air.

 

Fig. 2. (A) Artistic view of a chameleon on a tree branch. The lower inset shows two transmission electron microscopy (TEM) images of the guanine crystal in the epidermis of the skin of a male panther chameleon in the artificially relaxed and excited state (scale bar 200 nm) [26]. Male panther chameleons are able to change their color in the respective part of the skin from bluish in the relaxed state to red in the excited state. (B) Model of the chameleon’s skin consisting of guanine spheres immersed in a cytoplasm matrix. The spheres have radius $r$ and the distance between nearest spheres is $s$. The distance between adjacent layers of spheres is $Z$ and $d$ indicates the distance between the cytoplasm-air interface and the centers of the spheres of the first and bottom layers. The drawing of the chameleon was obtained from [39] https://freesvg.org/chameleon-129790. The TEM images are reproduced from Fig. 2a in Ref. [26] (Creative Commons CC BY license).

Download Full Size | PDF

In the particular case of the chameleon’s structure, each layer is a 2D triangular lattice, i.e., $b=s$ and $\gamma =60^\circ$. Therefore, the distance between adjacent layers is $Z=(\sqrt {3}/2)\,s$. The layers are stacked in an ABC sequence so that for an infinite number of layers ($N \rightarrow \infty$) the structure becomes an FCC (photonic crystal) with stacking direction (111). We fixed the distance from the air-cytoplasm interface to the centers of the spheres of the top and bottom layer to half of the particle spacing, i.e., $d = 0.5 s$. In this way, we guarantee that the guanine spheres are completely immersed in the cytoplasm. We assumed values of 170 nm < $s$ < 340 nm and 46 nm < $r$ < 82 nm for the spacing and radius of the guanine spheres. With this choice the value range of $s$ is slightly larger than that reported for averaged values in Ref. [26]. Taking into account the thickness of 50 $\mu$m of the ordered iridophores, the number of layers can be estimated by dividing this value by the distance between the layers ($Z = \sqrt {3}/2 \,s$). For 170 nm < $s$ < 340 nm this gives the range 169.81 < $N$ < 339.61 for the number of layers. In the following, we chose $N=130$ for our simulations, taking into account that for $N \geq 130$ the spectral position and width of the main peak in the visible, which are the main characteristics that determine the observed color, do not change. However, it is important to mention that $N$ is not a convergence parameter for the KKR approach. As it increases, the total thickness of the structure increases, and this produces changes both in the periodicity of the Fabry-Perot oscillations at both sides of the main peak and in the fine details of the spectra within the UV region.

To analyze the electromagnetic response and study the influence of the relevant parameters on the reflectance, we consider that the system is illuminated by a linearly polarized plane wave. To simulate natural light illumination conditions, the total reflectance, the specular reflectance and the efficiencies of the diffracted orders were obtained by averaging the corresponding results for transverse electric and transverse magnetic linearly polarized incident light.

4.1.1 Normal incidence

First, we study the influence of the radius and the spacing between spheres on the specular reflectance at normal incidence ($\theta =0^\circ$). In Fig. 3(A) we plot the specular reflectance as a function of the wavelength $\lambda$ for the spacing values $s = 170$, 255, and 340 nm for a fixed sphere radius $r=64$ nm. The calculated spectra exhibit a response typical for photonic crystals. The wide peaks correspond to photonic band gaps. It can be noticed that the respective main peaks are very sensitive to spacing variations. As the spacing $s$ increases, the peak redshifts and becomes narrower, see, e.g., the wide gap for wavelengths close to 450 nm for $s = 170$ nm and the much smaller gap for $s= 340$ nm for wavelengths around 800 nm. It is important to note that in addition to the main peak observed in the visible, an intensified reflectance is found in the near-UV region, which also redshifts as the spacing increases. The highly oscillating behavior observed at the right of each main peak, can be associated with Fabry-Perot resonances of the whole structure.

 

Fig. 3. (A) Specular reflectance at normal incidence as a function of the wavelength $\lambda$ for different spacing values $s$ for a sphere radius $r=64$ nm. (B) Color map of the specular reflectance as a function of $\lambda$ and the spacing $s$ for a fixed radius $r=64$ nm. The color coding ranges from blue (low) to red (high). The bright vertical red stripe in the middle of the map represents the movement of the gap which shifts from about 400 nm to 800 nm as the spacing $s$ increases from 170 nm to 340 nm. The left vertical red stripe corresponds to the intensified UV zone that redshifts, too. (C) Color map of the specular reflectance as a function of $\lambda$ and the radius $r$ for a fixed spacing $s=255$ nm. In this case the width of the gap (red stripe) close to 600 nm increases with $r$. Its position, however, does not change significantly. The UV intensified zone (left vertical red stripe) does not alter, too.

Download Full Size | PDF

The above described shift of gap position and width becomes clearer if we plot the optical response as color maps. In Figs. 3(B)-(C) we display the simulated specular reflectance at normal incidence as a function of the wavelength $\lambda$, for varying $s$ ($r=64$ nm) (B) and $r$ ($s=255$ nm) (C). In these maps, the ranges of spacing $s$ and radius $r$ were chosen according to experimental values reported in [26]. As observed in Fig. 3(B), the reflectance peak redshifts from 400 to 800 nm as the spacing $s$ increases from 170 to 340 nm. The intensified zone in the UV region, on the other hand, moves from 200 to 400 nm.

It is now interesting to compare this result with the color map in Fig. 3(C). Here, the spacing is fixed while the radius of the spheres is increased. In this case, the main peak becomes broader as the radius increases from 46 to 82 nm, but no significant shift is observed. The same observation holds for the intensified UV region. The sharp peak stays at 300 nm, although the radius is nearly doubled. These calculations demonstrate that the color response of these photonic crystals depends mainly on the spacing between the spheres while their radius is less important.

4.1.2 Different angles of incidence

In the following we explore the dependency of the reflectance on the incidence conditions. Figure 4(A) displays the specular reflectance for a structure with radius $r=64$ nm, spacing $s=340$ nm, and $\theta = 0^\circ$. Plotting three values of $\theta = 0^\circ, 30^\circ$, and $40^\circ$ we observe that the main peaks of the specular reflectance blueshift as $\theta$ increases.

 

Fig. 4. Reflectance spectra calculated for a regular structure with radius $r=64$ nm and $s=340$ nm. (A) Specular reflectance obtained for a fixed value of $\phi=0^\circ$, for three different incidence angles: $\theta =0^\circ, 30^\circ$, and $40^\circ$. (B) Reflectance calculated for $\theta =40^\circ$ and $\phi=0^\circ$. The total and specular reflectance and the efficiencies of the diffracted orders $(-1,0)$ and $(-1,-1)$ are shown for comparison. The cut-off wavelength $\lambda _{\textrm{cut}}=443$ nm is indicated as a vertical blue dashed line. (C) Zoom into the plot shown in (B) for $200\; {\rm nm} < \lambda < 443$ nm, i.e., the range in which the non-specular orders $(-1,0)$ and $(-1,-1)$ propagate.

Download Full Size | PDF

In Fig. 4(B) we show the specular and total reflectance for the same structure, for $\theta =40^\circ$ and $\phi =0^\circ$. As can be observed, the total reflectance coincides with the specular for wavelengths larger than $\lambda _{\textrm{cut}}=443$ nm (calculated using Eq. (4) and indicated with a vertical blue dashed line). For $\lambda < 443$ nm other diffracted orders start propagating. It can be observed that not only the specular order (magenta) contributes to the total reflectance (black), but also other diffraction orders such as the (-1,0) (red) and the (-1,-1) (cyan) give significant efficiency. It is important to mention that, for clarity, the efficiencies of other reflected orders activated for $\lambda < 443$ nm are not shown. However, in order to allow for a better comparison of the different spectra, in (C) we show a zoom of Fig. 4(B) for wavelengths at the left of $\lambda _{\textrm{cut}}$, where non-specular orders are activated, i.e., for $200\; {\rm nm} < \lambda < 443$ nm.

Here, it is important to mention that all peaks shift to the left for smaller spacing values and $\lambda _{\textrm{cut}}$ decreases, too. For $s = 255$ nm, for example, the cut-off wavelength is $\lambda _{\textrm{cut}} = 332$ nm and non-specular orders appear in the UV (see Fig. S1 in Supplement 1). Although the UV range is hardly visible for humans, it is well-known that several reptiles, birds, and insects are able to detect also this range of the electromagnetic spectrum and therefore, this characteristic could have biological implications in how chameleons are seen by other animals, like birds for instance [37]. It is worth noting that the reflectance spectra for different values of $\phi$ do not exhibit significant differences in the visible region (see Fig. S2 in Supplement 1).

We are interested in investigating the ranges of values of $\theta$ and $\phi$ that produce non-specular propagating diffraction orders within the visible and UV regions, which can be perceived by several animals, i.e., $\lambda _{\textrm{cut}} > 200$ nm. We, therefore, utilize Eq. (4) to calculate $\lambda _{\textrm{cut}}$ for different propagating orders. As an example, we plot in Fig. 5(A) a color scale map of the cut-off wavelength $\lambda _{\textrm{cut}}$ for the order $(p,q)=(-1,0)$ as a function of $\theta$ and $\phi$ for a structure with $s=340$ nm. When $50^\circ < \phi < 250^\circ$, $\lambda _{\textrm{cut}}$ decreases as $\theta$ increases. For other values of $\phi$, $\lambda _{\textrm{cut}}$ increases with increasing $\theta$, so it is more likely to get significant non-specular contributions of this order to the total reflectance. The limit value $\lambda _{\textrm{cut}} = 200$ nm is indicated with a black line. In the region enclosed by this line, the order $(-1,0)$ is not activated for any value of $\lambda$ in the range of interest. It is important to mention that the color maps obtained for the other first orders (not shown) are translations in $\phi$ of the map shown in Fig. 5(A) due to the symmetry of the lattice. This implies that depending on the incident conditions, some of these orders may or may not propagate.

 

Fig. 5. (A) $\lambda _{\textrm{cut}}$ vs. $\theta$ and $\phi$ for the diffraction order $(p,q)=(-1,0)$ and for $s=340$ nm; (B) $\lambda _{\rm Mcut}$ vs. $\theta$ and $s$ for the same diffraction order.

Download Full Size | PDF

For a given structure, $\lambda _{\textrm{cut}}$ of a given order $(p,q)$ is a function of $s$, $\theta$ and $\phi$ (Eq. (4)). If we fix $s$ and $\theta$, it is possible to find a value of $\phi$ which produces the largest achievable value of $\lambda _{\textrm{cut}}$, $\lambda _{\rm Mcut}$. However, it can be demonstrated that this value of $\phi$, $\phi (\lambda _{\rm Mcut})$, does not depend on $\theta$, i.e, is a function of $s$ only. For the parameters considered in Fig. 5(A), $\phi (\lambda _{\rm Mcut}) = 330^\circ$, as can be appreciated in this figure.

In Fig. 5(B) we plot $\lambda _{\rm Mcut}$ vs $\theta$ and $s$ for the order $(-1,0)$. Note that for fixed $\theta$ ($s$), as $s$ ($\theta$) increases, $\lambda _{\rm Mcut}$ also increases. In the figure, the curves for four values of $\lambda _{\rm Mcut}$ (200 nm, 300 nm, 400 nm, and 500 nm) are drawn as contour lines. Within the parameters’ region above each of these curves, diffraction orders other than the specular are activated for $\lambda < \lambda _{\rm Mcut}$. Therefore, in the region shaded with red tones, that is, the one with the largest values of $s$ and $\theta$, non specular orders propagate both in the visible and in the UV region. Contrarily, for $s<230$ nm and $\theta < 20^\circ$ (region below the curve $\lambda _{\rm Mcut}=200$ nm), only the specular order propagates for incident wavelengths within the spectral region of interest.

The region at the top right of Fig. 5(B) has the largest $\lambda _{\rm Mcut}=500$ nm. Therefore, as non-specular diffracted orders appear in this region, diffracted orders take place in the near UV and partially in the visible. It is important to remark that the efficiency of these orders cannot be deduced neither from the diagram nor from Eq. 4. In order to analyze the relative relevance of non-specular orders in the UV, it is necessary to calculate their respective efficiencies.

4.1.3 Color impression for humans

It is known from experimental studies that in the top layer of the chameleon skin, the guanine arrangements are more or less randomly oriented [26]. Therefore, one possibility to simulate the color observed by humans is to perform an average of the reflectance for different incidence conditions. In Fig. 6, we plot the average specular reflectance for five values of $s$. Taking into account the symmetry of the hexagonal lattice used in the simulations and considering typical illumination conditions that occur in a natural environment, each reflectance spectrum was obtained considering $0^\circ < \phi < 55^\circ$ and $0^\circ < \theta < 40^\circ$. To calculate the observed color from the human perspective, we obtain the three standard tristimulus values X, Y and Z [40]. In the bar at the right of Fig. 6 we show the calculated color for each value of $s$. As already reported by Teyssier et al. [26] the color strongly depends on the lattice parameter $s$. Our results show that, as the spacing $s$ increases, the resulting color impression for humans goes from blue over green to red as already reported in Ref. [26] (see the Fig. S3 chromaticity diagram in Supplement 1).

 

Fig. 6. Averaged specular reflectance for different spacing values $s$. The right bar shows the calculated color for each value of $s$ as it would appear to the human eye.

Download Full Size | PDF

4.2 Phelsuma lizards

As second example, we analyze the structural color produced in the skin of Phelsuma lizards. Again, arrangements of guanine nanocrystals within the iridophores located in the epidermis cause interesting color effects [27]. A first approach to investigate the response of this structure was already presented by Saenko et al. [27]. In this study, it was assumed that the ordered microstructure is composed of stacks of layers of disks of guanine in cytoplasm, as schematized in Fig. 7. Based on this assumption the structure was modelled as a multilayered arrangement of 32 alternating layers of cytoplasm and layers with an effective medium (cytoplasm and guanine). The effective permittivity ($\epsilon _{\rm eff}$) was obtained using the Maxwell-Garnett formula [10], with a filling fraction of $60$%, which yields a value of $\epsilon _{\rm eff}=2.63$. The homogeneous layers had the same thickness $d_c=97 \pm 12$ nm while the composite ones had a thickness $d_g=83 \pm 11$ nm. Mean values and standard deviations of these parameters were obtained from Ref. [27] from a statistical analysis of several microscopy images of a Phelsuma laticauda specimen. The whole structure is immersed in air. To calculate the electromagnetic response of the multilayer system, we employ the transfer matrix method [11] utilizing the parameters and procedures described in [27]. A scheme of this structure is shown in Fig. 8(A).

 

Fig. 7. (A) Sketch of a lizard. The lower inset shows a TEM image of the longitudinal section of the guanine crystals in the epidermis of the lizard’s skin (scale bar, 1 $\mu$m). (B) Microstructure model of lizard’s skin of Saenko et al. [27]. The structure is assumed to be composed of a regular arrangement of guanine disks ($n_g = 1.83$) immersed in cytoplasm ($n_c = 1.33$). The relevant geometrical parameters are indicated by black arrows. The drawing of the lizard was obtained from [41] https://freesvg.org/lizard-silhouette-118766. The TEM image in (A) is reproduced from Fig. S3c in Ref. [27] (Creative Commons Attribution License 4.0).

Download Full Size | PDF

 

Fig. 8. Direct comparison of the Maxwell-Garnett (MG) approach and the KKR method. (A) The MG approach of Saenko et al. [27] assumes a multilayered arrangement of alternating layers of cytoplasm and of an effective medium comprising cytoplasm and guanine. (B) The KKR method, on the other hand, considers alternating layers of cytoplasm and of guanine spheres immersed in cytoplasm. (C) Simulated reflectance obtained with the KKR method (solid orange line) and with the MG method (dashed blue line) for $d_c=97$ nm and $d_g=83$ nm. As observed, the MG method correctly predicts the band gap around 550 nm, but fails to calculate the secondary peak in the UV. (D) Averaged reflectance for both methods utilizing an averaging procedure described in the text.

Download Full Size | PDF

To get a more detailed insight into the optical response of this structure, we apply the KKR method. Each layer of disks is modeled as a layer of spheres whose diameter matches the height of the disks. We consider 32 alternating layers of cytoplasm and layers of spheres centered at the nodes of a triangular lattice with $s=89$ nm and a filling fraction $f\approx 60{\% }$ (Fig. 8 (B)). There is no lateral displacement between the spheres of different layers.

In Fig. 8 (C) we compare the reflectance obtained by using the two models considering normal incidence and the mean values of the thicknesses. Both curves exhibit a typical photonic crystal behavior with a bandgap between 500 and 550 nm. However, in addition to this peak, the KKR curve exhibits an enhancement in the UV region near $\lambda = 250$ nm, which is not present in the reflectance calculated by the MG approach. This peak is the result of multipolar interactions of higher order that occur at short wavelengths, which are not taken into account within the MG formulation [42]. This comparison evidences the importance of using a rigorous electromagnetic method which takes into account higher order multipolar interactions, such as the KKR.

As previously stated, the measured thicknesses of the layers present variations throughout the sample. For this reason, in order to make a more precise estimation of the reflectance which could be comparable to experimental measurements, we applied an averaging method in both the MG and the KKR approaches, to obtain the curves in Fig. 8 (D). We employed the procedure followed in Ref. [27] to perform the averages. First, we fixed $d_g = 83$ nm and obtained the reflectance as the weighted average for 100 equally spaced values of $d_c$ in the range 97 $\pm$ 24 nm. Second, the same procedure was repeated for the averages over $d_g$ in the range 83 $\pm$ 22 nm and considering $d_c=97$ nm. Finally, the reflectance was calculated as the arithmetic mean of both averaged reflectances. In Fig. 8 (D) we compare the averaged reflectances obtained with the MG and the KKR methods. In comparison with Fig. 8 (C), notice that the oscillations smooth out as a result of the averaging procedure. Both spectra in Fig. 8 (D) are practically identical for $\lambda > 400$ nm, with a main peak at $\approx 525$ nm. The KKR curve, however, features an additional secondary peak in the UV that is not predicted by the MG method. This additional peak is the result of multipolar interactions of higher order and evidences the need of rigorous electromagnetic methods to accurately simulate the optical response of nanostructured biological tissues. And indeed, there are significant additional peaks in the UV reported in Fig. 3(a) of Saenko et al. [27] . However, to perform a detailed comparison with their experimental results, additional characteristics of the structure (presence of pigments, irregularities in inclusions’ shape and distribution, curvature, etc.) have to be taken into account.

5. Conclusions

To conclude, we analyzed the optical response of reptiles’ skins that exhibit structural color by the example of the male panther chameleon and Phelsuma lizards. In both cases, we rigorously modelled the microstructure responsible of the color response as a quasi-regular arrangement of guanine spheres immersed in cytoplasm. By applying the KKR method, we analyzed their reflectance response for various cases. Since several animals have a larger visibility range than humans, we explicitly considered the visible and near UV regions. We investigated the dependency of the reflectance on the geometrical parameters and on the angle of incidence, and performed a detailed analysis of the generation of diffraction orders. The application of a method such as the KKR permitted us to predict important characteristics of the reflectance spectra, especially within the UV and short visible wavelengths region, that cannot be explained by approximate models like the often applied transfer matrix approach. Our study highlights the importance of employing rigorous electromagnetic methods to simulate the optical response of biological systems, especially for structures that act like photonic crystals.