1. Introduction

Many species of butterfly display vibrant iridescent coloration that is widely believed to be caused by optical interference. This type of coloration is called structural color (see [1] for an overview of the field). Structural color can have several interesting features: The hue often changes with viewing angle, the color is often very intense and highly saturated, and the color can be changed or even eliminated by immersing the insect in a liquid with an appropriate index of refraction.

Scientific interest in structural color of insects has persisted for well over a century. Newton, Rayleigh, and Michelson all contributed theories to the field. The major breakthroughs in understanding occurred with the invention of the electron microscope [2]. For the first time, it became possible to observe the subwavelength features that are the source of structural color. Talented biologists, notably Ghiradella, published extensive studies of the structural features [3, 4]. Her studies, among others, were used to create a classification scheme for insects based on the underlying optical principles of the color [5].

A typical butterfly exhibiting structural color is the male Morpho rhetenor, shown in Fig. 1. This species is a very intense blue that fades to a deep violet when viewed near grazing. The wings are covered with a large number of small scales that are easily visible using a low-power optical microscope. As shown in Fig. 1b, it is actually the scales that produce the color. The scales, in turn, are covered with a series of ridges that are less than $1\mathrm{\mu m}$ apart, which are shown in Fig. 1c. Each ridge is a treelike structure with 4–12 branches on both sides. The structure is made of cuticle with an approximate complex index of refraction $\tilde{n}=1.56+i0.06$ [6].

In some species of butterfly, there are two types of scales. The first type, known as a ground scale, is similar to the scales on the M. rhetenor. In addition, a second type of scale, known as a cover scale, is layered over the ground scales. These scales generally have fewer ridges than the ground scales and are semitransparent. They are thought to diffuse the scattering from the ground scales [3]. This work primarily concentrates on ground scales, although a simplified cover scale is considered.

In this paper we use one of the most powerful numerical techniques for solving electromagnetic and optical problems, the finite-difference time-domain (FDTD) method [7], to study the scattering of light by the ridge structure described above. We perform a fully three-dimensional analysis that makes use of the most recent refinements to the method, which include a new approach for handling periodic structures [8, 9]. In the past, the FDTD method has been used to study butterfly scales; however, restrictions not introduced here, such as a reduced number of dimensions (two), were employed [10, 11]. In this paper we also present a novel method for displaying the color associated with the three-dimensional scattering pattern for the scale. We vary several of the features of the scale and observe the changes in color that result by using this new method. Our goal is not to create an exact computational model for the M. rhetenor or any other butterfly; rather, it is to study the effect of varying geometrical parameters on a generic scale.

2. Computational Model

2A. Discretized Geometry

By examining electron micrographs such as Fig. 1c, a computational model for a generic butterfly scale can be developed. An example of such a model is shown in Fig. 2a. In this model, the base, which is the relatively featureless cuticle sheet at the bottom of Fig. 1c, is modeled as a dielectric slab. The upper structure, known as a ridge, is a set of dielectric slabs with different dimensions and orientations. First, there is a central stem that supports the structure. On either side of this stem are branchlike arms called lamellae. The lamellae are often tapered, with those closest to the base being the widest. Usually, the lamellae are slightly tilted in the longitudinal direction (y) with respect to the base; however, some species have a very high tilt, as is discussed below. In addition to the lamellae, there is sometimes a second set of layers nearly orthogonal to the first. These layers are smaller than the lamellae and are called microribs. In Fig. 2, the components of the ridge are shown in different colors for ease of identification. In the model, they all have the same electrical properties: electrical permittivity $ϵ=2.43{ϵ}_{\mathrm{o}}$, magnetic permeability $\mu ={\mu }_{\mathrm{o}}$, and electrical conductivity $\sigma =6.4\times {10}^3\mathrm{S}/\mathrm{m}$. The permittivity and conductivity were obtained from the aforementioned complex index of refraction, $\tilde{n}$, at the wavelength $\lambda =488\mathrm{nm}$.

The structure in Fig. 2a is doubly periodic. The periodic cell, indicated by the black frame, is repeated three times in each lateral direction in this figure. The periodic cell is examined in more detail in Fig. 3. In this figure, the width of the periodic cell, $x_p$, is the period of the ridges, and the length, $y_p$, is determined by the spacing between the lamellae, g, the thickness of the lamellae, h, and the tilt angle of the lamellae with respect to the base, α, according to

Figures 3a, 3b, 3c show three variations on the generic model. The number of lamellae in all cases is N. In (a) the lamellae are all of the same width, a, and those on the left and right sides of the stem are aligned. In (b) some of the lamellae are tapered, so the maximum width is a and the minimum width is d. Finally, in (c) the lamellae on the left and right sides of the stem are offset vertically instead of aligned. Notice that there is also a dielectric layer above the ridges in (c). It is intended to be a highly simplified model for the aforementioned cover scale.

In the actual scales, there is also structure between the ridges and the base that supports the ridges. Our calculations show that this structure has little effect on the color of the scattered light, so it is omitted in this model.

2B. Finite-Difference Time-Domain Method

At this time, the FDTD method has been developed to a point at which many of the features are quite complicated. Our purpose here is not to cover the details of these features; rather, it is to give a simple explanation of the basic algorithm and an overview of the special features used for the current problem.

The FDTD method is used to directly solve Maxwell’s equations in the time domain [7]. In time domain, vector form, the equations are

The FDTD algorithm uses a procedure called “marching in time” to calculate future values of the electromagnetic field from current and past values. For example, the value of ${\mathcal{H}}_z$ at the future time $t=(n+1/2){\mathrm{\Delta }}_t$ is determined from the value of ${\mathcal{H}}_z$ at the past time $t=(n-1/2){\mathrm{\Delta }}_t$ and the values of ${\mathcal{E}}_x$ and ${\mathcal{E}}_y$ at the current time $t=n{\mathrm{\Delta }}_t$ using the update equation

Figure 4 is a schematic drawing for the FDTD computational domain used for our problem, one periodic cell. A field must be introduced into this cell before the algorithm in Eq. (5) can be used to advanced the solution in time. For this problem, the field is an incident plane wave, and it is produced by the “plane-wave injector” near the top of the cell. Here we use the recently developed “nearly perfect” plane-wave injector described in [13]. The injector surface is designed to produce the incident field below, toward the structure, and no field above. As a result, the injector surface naturally divides the domain into two regions: the region below in which the incident field and the scattered field (the field produced by the interaction of the incident field with the structure) exist, and the region above in which only the scattered field exists. These two regions are referred to as the “total field region” and “scattered field region,” respectively.

The incident field produced by the injector is a linearly polarized plane wave with an electric field of the form

Absorbing boundary conditions are applied to the top and bottom of the domain, as shown in Fig. 4. They prevent the reflection of outward propagating fields at these surfaces. For this work, we use the recently developed convolutional perfectly matched layer (CPML) absorbing boundary condition [14]. Periodic boundary conditions (PBCs) are applied at the other four faces of the domain (front, back, left, and right surfaces). These conditions make the calculation one in which the unit cell is effectively repeated ad infinitum in the lateral (x, y) directions. Here we use a recently developed technique for implementing PBCs that eliminates much of the previous complexity of modeling periodic structures in the FDTD method [8, 9].

2C. Modeling of Natural Light

We will consider butterfly scales exposed to natural light (daylight). Natural light is incoherent light that is unpolarized. Here it is modeled by making two separate calculations (1 and 2) in which the states of linear polarization of the incident plane waves, represented by Eq. (6), are orthogonal: ${\widehat{e}}_1·{\widehat{e}}_2=0$. All subsequent calculations based on these fields, such as reflection and transmission coefficients, involve the power, and they are obtained by averaging the appropriate powers for the two states.

3. Effect of Randomness in the Geometry

A periodic structure, such as the one shown Fig. 2, is a diffraction grating. When it is infinite in lateral extent, the scattered radiation is a discrete set of plane waves, the directions of which are determined by a grating equation. For a large grating (many periods), the scattered radiation consists of a number of very narrow lobes (grating lobes), with their directions again given by the grating equation. Observations of the scattered light from butterfly scales and wings do not show these narrow lobes. Using their measurements, Kinoshita and Yoshioka have determined that the absence of these lobes is probably due to irregularities in the structure, namely, random positioning of the ridges on a scale [15]. As a result of this randomness, the scattered fields from the ridges add incoherently, rather than coherently as in a perfect diffraction grating. The assumption of randomness, which we adopt, implies that the scattering from a whole scale can be determined from that of a single ridge. Below we describe numerical experiments that were performed to better understand this assumption and how to implement it in the FDTD method.

To reduce the computational burden (storage and run time), the numerical experiments were performed using a simplified version of the model shown in Fig. 2. The model was two dimensional; that is, it was invariant in the y direction ($\alpha =0$, $y_p\to \infty$) with a cross section similar to the one shown in Fig. 3b with the top four lamellae tapered, but without the microribs. The other dimensions were $N=8$, $c=60\mathrm{nm}$, $a=620\mathrm{nm}$, $d=140\mathrm{nm}$, $h=60\mathrm{nm}$, $g=140\mathrm{nm}$, and $x_p=750\mathrm{nm}$.

First, a random array was simulated using $n=50$ ridges. The vertical positions of the ridges, the displacements ${\mathrm{\Delta }}_h$ in Fig. 5a, were varied randomly with a triangular probability density, such that $-275\mathrm{nm}\le {\mathrm{\Delta }}_h\le 275\mathrm{nm}$. The far-zone or Fraunhofer scattered field ($\lim kr\to \infty$) of the array was determined using a form of Huygens’ principle, which states that the far-zone scattered field can be determined from an integral involving the near-zone scattered field on a closed surface surrounding the array, such as the box marked “transform surface” in Fig. 5a. In the literature for the FDTD method, this procedure is called a near-to-far-field transformation. Now for a large array, the right and left sides of the box make relatively small contributions to the integral and can be ignored in the calculation. In addition, the total field below the array is approximately zero; all of the energy in the incident field is scattered or absorbed before reaching this surface. Thus the scattered field on the bottom surface of the box is approximately the negative of the incident field, which is unaffected by the randomness in the array. This coherent field on the bottom of the box produces a beam of radiation below the box (approximately the negative of the incident field) and little radiation above the box. From these observations, we conclude that it is mainly the near-zone scattered field on the top surface of the box that determines the far-zone scattered field above the box (the backscattered field). In the calculations that follow, only this surface will be used.

To obtain a statistical representation, the far-zone scattered field from the array of 50 randomly placed ridges was calculated 100 times using different random displacements for each case. The ensemble average of the 100 cases was then computed. This average can be thought of as what one sees when viewing a section of a wing containing many scales. The graphs in Fig. 6a are the results of those calculations, which are far-zone scattering patterns that show the time-average scattered power per unit area as a function of the direction θ in the x–z plane. Notice that the scattering is largest at the wavelengths that represent blue–violet light. The small “noiselike” variations along each curve are due to the random aspect of the calculation.

Next, we show how we can obtain results similar to those for the random array, shown in Fig. 5a, using the periodic FDTD method. For determining the far-zone scattered field above the array, the surface shown in Fig. 5b is equivalent to the top surface of the box in Fig. 5a. We ignore the small vertical portions of this surface; that is, we only consider the horizontal surfaces marked $S_1,S_2,...,S_n$ in our calculation. Because of the random displacement of the ridges, the scattering from these surfaces is incoherent, and in a statistical sense, it is equivalent to n times the scattering from one of the surfaces. Thus, knowledge of the near-zone scattered field on the surface immediately above one ridge can be used to obtain the far-zone scattered field from the random array of ridges. We approximate this field by the field from the periodic FDTD calculation based on the unit cell shown in Fig. 4, the field on the “transform surface” in the scattered field region. We note that because of the periodic boundary conditions in this calculation, the field within the unit cell includes the influence of all of the other ridges in the array.

The results based on the periodic FDTD calculation are shown in Fig. 6b, and they are to be compared to those for the random array in Fig. 6a. The scattering patterns for the two methods are very similar (the former is a smoothed version of the latter), with the shapes of the curves and relative values for the different wavelengths being essentially the same.

The power reflection coefficient is a measure of the total power scattered per ridge of the scale. It is obtained by integrating the normal component of the Poynting vector over the transform surface, then dividing by the total power incident on this surface. We note that, for the random array, the transform surface is the area above the 50 ridges that make up one element in the ensemble average, whereas for the periodic element, the transform surface is the area above the ridge (unit cell). The power reflection coefficients for the two methods are compared in Fig. 6c. As for the scattering patterns, we see that the results from the two methods are very similar.

The results presented in this section show that FDTD calculations based on the single periodic cell can be used to estimate the scattering, in an average sense, from a scale containing a large number of irregularly displaced ridges.

4. Interpretation of Scattering as Observed Color

In Section 3, we described the use of the periodic FDTD method to calculate the scattered field from a scale with irregularly positioned ridges. In this section, we explain the procedure used to obtain the observed color of a butterfly model from the results of the FDTD simulations. The goal is to calculate the observed color as a function of observer location when the structure is illuminated by an unpolarized light source located infinitely far away.

Referring to the flowchart in Fig. 7, the first step is to specify the incident field. Two simulations are performed using plane-wave excitations with orthogonal states of polarization. The scattered field is recorded on the transform surface above the structure. Although the computational method is a time- domain method, frequency-domain quantities are desired. Therefore the scattered field is converted to the frequency domain using an on-the-fly discrete Fourier transform (DFT) [7]. The frequencies used correspond to 81 wavelengths spaced in $5\mathrm{nm}$ increments between $\lambda =380\mathrm{nm}$ and $\lambda =780\mathrm{nm}$. This set of wavelengths covers the visible spectrum.

After the simulations are completed, the far-zone scattered field is calculated using the previously described frequency-domain, near-to-far-field transform (Huygens’ principle). From the far-zone field, the radiant intensity, I, of the scattered light is found. In addition, the radiant flux of the incident light, $I^i$, is calculated. Note that both quantities are functions of the wavelength, but the scattered radiant intensity is also a function of the direction. The reflection factor, $R(\lambda ,\theta ,\varphi )$, is calculated as the ratio of the scattered radiant intensity to the incident radiant flux [16].

At each observation direction, the reflection factor is converted into tristimulus values $(X,Y,Z)$, for example,

The observed color is plotted as a three-dimensional pattern; the procedure is shown schematically in Fig. 8. The radial vector from the origin O is proportional to the luminance factor Y. The small patch of surface (angular size $\mathrm{\Delta }\theta$, $\mathrm{\Delta }\varphi$) at the end of this vector is the color that would be seen by an observer viewing the structure from that direction $(\theta ,\varphi )$. By examining such patterns, one can see how the luminance and color of the scattered light changes with the direction of observation.

5. Results

5A. Models for Morpho Butterflies

In this subsection, we examine the scattering from structures that are similar to a Morpho rhetenor scale (once again, the computational model is inspired by the specific insect, but it is not meant to be an exact model). The dimensions of the structure are based on the observations described in [19] and are summarized in the first column of Table 1.

Figure 9 shows the observed color for normally incident light (${\theta }^i=0°$) as a function of the location of the observer $(\theta ,\varphi )$. Each subplot corresponds to a change in the structure. Notice the small insets in Figs. 9a, 9d; they are simplified models that show the orientation of the ridge. In this figure, and in all similar figures that follow, logarithmic scaling is used; that is, the length of the radial vector is proportional to ${\log }_{10}(Y/Y_{\max })$, and the range displayed is two decades. Here, $Y_{\max }$ is the maximum value for each model for a given genus, and it is given for each subplot in the caption of the figure. Using these values, we can compare the relative luminance values for the subplots.

Figure 9a is for the simplest model; it is the one in Fig. 3a with a base and without microribs. The eight lamellae on the left and right sides of the stem are aligned, and all have the same width. Notice that the color is a highly saturated blue and is concentrated in a small angular sector, a beam. This beam is tilted at approximately the angle $\theta =2\alpha =20°$. That is, the tilt of the lamellae with respect to the plane of the wing (base) causes a rotation in the beam that is approximately twice the tilt angle. This rotation is the same as one would expect for specular reflection from a tilted dielectric slab.

To provide a quantitative measure for the shape of the beam in Fig. 9a, the scattered power is examined for two planes. Figure 10a is for the plane of the ridge (y–z plane), while Fig. 10b is for a plane perpendicular to the y–z plane and passing through the peak of the beam. Curves are given for four wave lengths in the visible range. The aforementioned tilt of the beam to $\theta \approx 20°$ is evident in Fig. 10a. For all wavelengths, the width of the beam is seen to be significantly greater in the perpendicular plane, Fig. 10b. This behavior is consistent with measured results for Morpho butterflies [6, 20].

The power reflection coefficient versus wavelength for this case is shown in Fig. 11a. There is a definite peak at approximately $\lambda \approx 470\mathrm{nm}$, which is in the blue range of the visible spectrum. The reflection coefficient is similar to the measurements of Berthier ([1], Fig. 8.20).

For Fig. 9b, the structure is identical to that for Fig. 9a, except the top four lamellae are tapered in width, as in Fig. 3b. The pattern of the scattered light is similar for both cases. However, as the power reflection coefficient in Fig. 11a shows, the intensity of the light is lower with the taper. This is also evident from the luminance scale factors for Figs. 9a, 9b: $Y_{\max }=1.0$ versus $Y_{\max }=0.30$.

For Fig. 9c, the model is changed once again by adding a cover scale (dielectric slab) above the ridge. Recall that the cover scale is prevalent in Morpho-type butterflies but not the M. rhetenor [6]. Two beams can now be identified in the pattern, a beam due to specular reflection from the ridge ($\theta =2\alpha =20°$), which is similar in shape and hue to the one in Fig. 9b, and a beam due to specular reflection from the diffusing layer ($\theta =0°$), which is of different hue and less saturated. The combination produces a pattern that is broader than in Fig. 9b, which is consistent with the concept of a diffusing layer. The decrease in saturation is consistent with observations of Morpho butterflies with cover scales, such as M. didius, for which the wings are less glossy than the M. rhetenor, and the observed color is whitish-blue [21].

In Figs. 9d, 9e, 9f, 11b, the previous study is repeated, but the lamellae on the left and right sides of the stem are offset from each other, as is the case for some Morpho butterflies. Except for this change, the models are identical to the ones presented above. As seen in Fig. 9d, the main beam due to specular reflection is now split into two, with a relative null in between, and it is less intense than for the case with aligned lamellae: $Y_{\max }=0.21$ versus $Y_{\max }=1.0$. A small beam in the backscatter direction ($\theta =0°$) is also present; this beam is caused by scattering from surfaces that are parallel to the x–y plane (for example, the base layer). The split in the beam is consistent with the observations of Vukusic et al. [6] and Berthier et al. [20] and the calculations of Plattner [10]. It is caused by the destructive interference for the scattering from a pair of adjacent lamellae, one on each side of the ridge. These lamellae are offset vertically by the distance $(h+g)/2=100\mathrm{nm}$, which is close to $\lambda /4$ for wavelengths in the blue— violet region of the spectrum. This causes a relative phase shift of about $180°(\lambda /2)$ for specular scattering from adjacent lamellae. Notice that the hue of the light is more violet than for the case with aligned lamellae, which suggests a shift to shorter wavelengths. This observation is confirmed by compar ing the power reflection coefficients in Figs. 11a, 11b. The reflection coefficients computed using this method are similar to earlier results computed using lamellar grating theory [22].

For Fig. 9e, the top four lamellae are tapered in width. The taper is seen to shift the hue slightly toward violet. The effect of the taper on the intensity is more dramatic and similar to that seen for aligned lamellae. It causes a significant decrease, as can be seen in the power reflection coefficient in Fig. 11b and in the scale factors for the luminance: $Y_{\max }=0.05$ versus $Y_{\max }=0.21$.

Finally, the effect of a cover scale (dielectric slab) is shown in Fig. 9f. The specular reflection (beam) from the cover scale is significantly larger than the scattering from the ridge below and produces a fairly unsaturated blue color. Comparing this result with the one for aligned lamellae in Fig. 9c, we can directly see the effect of offsetting the lamellae. The scattering from the cover scale and ridge are comparable when the lamellae are aligned; however, offsetting the lamellae significantly reduces the scattering from the ridge, making the scattering from the cover scale the dominant feature.

Notice that in all cases shown in Figs. 9, there is a change in hue with viewing position (the angles θ and ϕ). This is a clear indication of the iridescence observed for these butterflies.

Next, the effect of adding microribs to the Morpho-like models is examined. We recall that the microribs are the small layers that are perpendicular to the lamellae, as shown in Figs. 2, 3. The microribs are added to the models with aligned and offset lamellae with taper; results for these cases are presented in Figs. 9b, 9e. The results in Fig. 12 show that the addition of the microribs causes very little change in the shape of the patterns and only a small change in the hue for the Morpho-like models. The reason is straightforward: The lamellae are tilted at a slight angle $(10°)$ relative to the base; thus the microribs are nearly perpendicular to the base. Therefore specular scattering from the microribs is mainly reflected into the adjacent structure rather than back at the observer.

In summary, the calculations for Morpho-like scales have demonstrated several interesting characteristics of the structural color for these butterflies, some of which are present in the measurements of previous investigators. The calculations have shown the following:

5B. Troides magellanus Butterfly

Next, the scattering from a structure similar to the male Troides magellanus butterfly scale is studied. This insect, shown in Fig. 13a, is remarkable due to the interplay of structural and pigmentary color [23, 24]. When viewed at any angle other than near grazing, the lower wings of the butterfly are yellow. However, when the observer and the light source are collocated in a direction near grazing and pointing toward the outer margin of the wing, the observed color changes with angle to blue, as shown in Fig. 13b.

The parameters used in the model for this butterfly are given in the second column of Table 1; they are based on values given in [23, 24]. The key feature is that the tilt angle of the lamellae is large, $\alpha \approx 54°$. Recall that the tilt angle was only $10°$ for the Morpho-like structures. Due to the larger tilt, the period along the y direction, $y_p$, is shorter than for the Morpho models by about a factor of four. Therefore, to keep the length for coherence the same for both butterflies, four periods along the y direction are used for the calculations for this insect. The other parameters for the lamellae are similar to those for the Morpho-like models; therefore it is expected that the observed color will be similar under certain conditions. Specifically, when the direction of the incident light is near normal to the lamellae (approaching grazing to the wing), then the observed color should be blue.

This hypothesis is confirmed by examining the plots in Figs. 14, 15. When the angle of incidence is well away from grazing, the scattering from the structure is very weak. This is the case in Fig. 14a, which is for ${\theta }^i=40°$. Notice that the scale factor for the luminance is $Y_{\max }=0.05$. The color is brown, not yellow as for the actual butterfly, because the yellow color is due to pigmentation, which is not included in the model. As the angle of incidence is increased slightly toward grazing, a lobe that is distinctly blue in color appears, as in Fig. 14b, which is for ${\theta }^i=50°$. This blue lobe is even more pronounced when ${\theta }^i=70°$, as in Fig. 14c. Notice that the scale factor for the luminance is now $Y_{\max }=1.0$. As for the Morpho, an intense blue appears when the light is nearly normally incident on the broad surfaces of the lamellae.

There are factors in addition to the direct scattering from the broad surfaces of the lamellae that contribute to the color shown in Fig. 14. As pointed out by Lawrence et al., the exposed ends of the highly tilted lamellae form a grating on the horizontal (x–y) plane [23]. For our dimensions, the period of this grating is $p_e=(g+h)/\sin \alpha =260\mathrm{nm}$, and there are two lobes in the scattering pattern. One is a specular lobe (zeroth order) that is fairly wavelength independent, and is probably responsible for the dark lobe on the right of all of the plots in Fig. 14. The other lobe (first order) occurs close to the backscattering direction ($\theta ={\theta }_i$) for wavelengths shorter than $\lambda =2p_e=520\mathrm{nm}$, and it contributes to the blue-violet color.

In Fig. 15, the power reflection coefficient is plotted versus wavelength for two angles of incidence, ${\theta }^i=40°$ and ${\theta }^i=70°$. These results clearly show the increase in scattering that occurs as the angle of incidence approaches grazing. They also clearly indicate the blue color for the structural scattering.

5C. Ancyluris meliboeus Butterfly

The final model examined is for the Ancyluris meliboeus butterfly. In contrast to the previous butterflies, it is the ventral side (bottom) of the wing that displays iridescence. Only small regions on this side are iridescent, and we restrict our discussion to one of these regions, that appears yellow. This butterfly has been studied by Vukusic et al. [25], and the dimensions in column three of Table 1 are based on their estimates for the ventral scales. The computational model is similar to the Morpho-like model with aligned lamellae discussed above. The most significant changes are that the lamellae are tilted by $30°$ relative to the plane of the wing, and the spacing between the lamellae is increased to $260\mathrm{nm}$. As for the previous case, multiple periods along the y direction (two in this case) are used to keep the length for coherence the same for all of the butterflies. The increased spacing between the lamellae causes the scattered light to shift toward longer wavelengths. As a result, the observed color, shown by the color pattern in Fig. 16a, is yellow instead of blue.

The observed color for the model with microribs added is shown in Fig. 16b. The shape of the beam, which is largest near the specular direction ($\theta =60°$), is not significantly changed, but the color is shifted toward longer wavelengths. In addition, a second lobe appears. This lobe is believed to be caused by the scattering from the microribs acting as a second set of lamellae. The shift in color with the addition of the microribs is also indicated by the power reflection coefficients in Fig. 17. The observed color is similar to that in photographs of a specimen taken by Vukusic et al. [25].

6. Conclusion

In this work, a set of detailed electromagnetic simulations for the structural color of butterfly wings was presented. The simulations were performed using the finite-difference time-domain method with periodic boundary conditions. The effect of the irregularity between the ridges of a scale was examined, and a method for using a periodic model to predict the response from an irregular structure was developed. The results from the simulations were interpreted using colorimetric theory in terms of human color vision. A novel method for presenting the three- dimensional color pattern was introduced. The techniques were used to study computational models that are similar to the structures of three butterflies: Morpho rhetenor, Troides magellanus, and Ancyluris meliboeus. The effect of variations in the structure of these models was explored.

This work was supported in part by the John Pippin Chair in Electromagnetics within the School of Electrical and Computer Engineering at the Georgia Institute of Technology.

Table 1. Dimensions for Butterfly Models

View Table

 

Fig. 1 Images of a male Morpho rhetenor butterfly. (a) Photograph of dorsal side (top) of insect. (b) Detail of wing showing scales as viewed under an optical microscope. (c) Cross section of scale showing individual ridges as viewed using a scanning electron microscope. SEM image courtesy S. Kinoshita [[26], p. 13, Fig. 12].

Download Full Size | PDF

 

Fig. 2 Computational model. (a) Portion of a scale showing three ridges. The volume within the black frame is the volume used for the periodic FDTD calculations. (b) Detail of a single ridge showing the individual FDTD cells. (c) Arrangement of the electromagnetic field components within a single FDTD cell.

Download Full Size | PDF

 

Fig. 3 Parameters used to describe the geometry of a periodic cell and the coordinates for a point P in space.

Download Full Size | PDF

 

Fig. 4 Schematic drawing showing elements within an FDTD periodic cell.

Download Full Size | PDF

 

Fig. 5 (a) Geometry for two-dimensional array of ridges with random heights. (b) Equivalent transform surface for obtaining the far-zone scattered field above the array.

Download Full Size | PDF

 

Fig. 6 Comparison of results from an ensemble average of randomized finite models and a periodic model. (a) Far-zone scattering pattern for the ensemble average of 100 simulations with 50 randomly spaced ridges per simulation. (b) Far-zone scattering pattern for the periodic model. (c) Comparison of the reflection coefficients for the ensemble average and the periodic model.

Download Full Size | PDF

 

Fig. 7 Flowchart illustrating the steps used to calculate the observed color using the periodic FDTD method. ℱ indicates the Fourier transform, and η is the intrinsic impedance of free space.

Download Full Size | PDF

 

Fig. 8 Illustration for the three-dimensional representation of the observed color as a function of the direction of observation.

Download Full Size | PDF

 

Fig. 9 Observed color as a function of direction for Morpho-like models. For all cases, there are eight lamellae and a base. (a) Aligned lamellae without taper, $Y_{\max }=1.0$. (b) Aligned lamellae with the top four lamellae tapered, $Y_{\max }=0.30$. (c) Aligned, tapered lamellae plus a cover scale (dielectric slab), $Y_{\max }=0.28$. (d)–(f) Same as (a)–(c), but with offset lamellae. (d) $Y_{\max }=0.21$. (e) $Y_{\max }=0.05$. (f) $Y_{\max }=0.22$.

Download Full Size | PDF

 

Fig. 10 Far-zone scattering patterns for the Morpho model shown in Fig. 9a. (a) Pattern in the y–z plane. (b) Pattern in a perpendicular plane that passes through the peak of the beam (angle measured from $\theta =20°$, $\varphi =90°$).

Download Full Size | PDF

 

Fig. 11 Power reflection coefficients (a) for structures with aligned lamellae, as in Figs. 9a, 9b, 9c, and (b) for structures with offset lamellae, as in Figs. 9d, 9e, 9f.

Download Full Size | PDF

 

Fig. 12 Observed color as a function of direction for Morpho-like models that include microribs. For both cases, there are eight lamellae with the top four lamellae tapered and a base: (a) aligned lamellae, $Y_{\max }=0.26$; (b) offset lamellae, $Y_{\max }=0.05$.

Download Full Size | PDF

 

Fig. 13 Photographs of the dorsal side of a male Troides magellanus butterfly. In both views, the light source is located near the observation point. (a) From most observation directions, the observed color of the lower wings is yellow. (b) When illuminated and viewed near grazing, the color of the lower wings changes abruptly, with a change in angle, to blue.

Download Full Size | PDF

 

Fig. 14 Observed color as a function of direction for Troides magellanus model with six aligned lamellae. The details for the model are the same for all cases, and only the direction of the incident light is changed: (a) ${\theta }^i=40°$, $Y_{\max }=0.05$; (b) ${\theta }^i=50°$, $Y_{\max }=0.14$; (c) ${\theta }_i=70°$, $Y_{\max }=1.0$.

Download Full Size | PDF

 

Fig. 15 Power reflection coefficient for Troides magellanus model. The rapid increase in the intensity of the iridescent color is shown by comparing the reflection coefficient for ${\theta }^i=40°$ and ${\theta }^i=70°$.

Download Full Size | PDF

 

Fig. 16 Observed color as a function of direction for Ancyluris meliboeus model with six aligned lamellae and a base: (a) without microribs, $Y_{\max }=1.0$; (b) with microribs, $Y_{\max }=0.48$.

Download Full Size | PDF

 

Fig. 17 Power reflection coefficient for Ancyluris meliboeus model. The addition of the microribs causes the peak in the reflection coefficient to shift to longer wavelengths.

Download Full Size | PDF

1. S. Berthier, Iridescences: The Physical Colors of Insects (Springer, 2007).

2. T. F. Anderson and A. G. Richards, “An electron microscope study of some structural colors of insects,” J. Appl. Phys. 13, 748–758 (1942). [CrossRef]  

3. H. Ghiradella, “Structure of iridescent lepidopteran scales: variations on several themes,” Ann. Entomol. Soc. Am. 77, 637–645 (1984).

4. H. Ghiradella, “Hairs, bristles, and scales,” in Microscopic Anatomy of Invertebrates, W. H. Frederick and L. Michael, eds. (Wiley-Liss, 1998), Vol. 11A, pp. 257–287.

5. P. Vukusic, J. R. Sambles, and H. Ghiradella, “Optical classification of microstructure in butterfly wing-scales,” Photonics Sci. News 6, 61–68 (2000).

6. P. Vukusic, J. R. Sambles, C. R. Lawrence, and R. J. Wooten, “Quantified interference and diffraction in single Morpho butterfly scales,” Proc. R. Soc. London Ser. B 266, 1403–1411 (1999). [CrossRef]  

7. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech, 2005).

8. R. T. Lee and G. S. Smith, “An alternative approach for implementing periodic boundary conditions in the FDTD method using multiple unit cells,” IEEE Trans. Antennas Propag. 54, 698–705 (2006). [CrossRef]  

9. R. T. Lee, “A novel technique for incorporating periodic boundaries into the FDTD method and the application to the study of structural color of insects,” Ph.D.thesis(Georgia Institute of Technology, 2009).

10. L. Plattner, “Optical properties of the scales of Morpho rhetenor butterflies: theoretical and experimental investigation of the backscattering of light in the visible spectrum,” J. R. Soc. Interface 1, 49–59 (2004). [CrossRef]  

11. S. Banerjee, J. B. Cole, and T. Yatagai, “Colour characterization of a Morpho butterfly wing-scale using a high accuracy nonstandard finite-difference time-domain method,” Micron 38, 97–103 (2007). [CrossRef]  

12. K. S. Yee, “Numerical solution of initial boundary value prob lems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]  

13. J. B. Schneider, “Plane waves in FDTD simulations and a nearly perfect total-field/scattered-field boundary,” IEEE Trans. Antennas Propag. 52, 3280–3287 (2004). [CrossRef]  

14. J. A. Roden and S. D. Gedney, “Convolutional PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media,” Microw. Opt. Technol. Lett. 27, 334–339 (2000). [CrossRef]  

15. S. Kinoshita and S. Yoshioka, “Structural colors in nature: The role of regularity and irregularity in the structure,” Chem. Phys. Chem. 6, 1442–1459 (2005). [CrossRef]   [PubMed]  

16. “Standard terminology of appearance,” ASTM E 284 08 (ASTM International, 2008).

17. R. S. Berns, Billmeyer and Saltzman’s Principles of Color Technology, 3rd ed. (Wiley, 2000).

18. “Standard practice for computing the colors of objects by using the CIE system,” ASTM E 308 06 (ASTM International, 2006).

19. R. A. Potyrailo, H. Ghiradella, A. Vertiatchikh, K. Dovidenko, J. R. Cournoyer, and E. Olson, “Morpho butterfly wing scales demonstrate highly selective vapour response,” Nat. Photon. 1, 123–128 (2007). [CrossRef]  

20. S. Berthier, E. Charron, and J. Boulenguez, “Morphological structure and optical properties of the wings of Morphidae,” Insect Sci. 13, 145–158 (2006). [CrossRef]  

21. S. Yoshioka and S. Kinoshita, “Wavelength-selective and anisotropic light-diffusing scale on the wing of the Morpho butterfly,” Proc. R. Soc. London Ser. B 271, 581–587 (2004). [CrossRef]  

22. B. Gralak, G. Tayeb, and S. Enoch, “Morpho butterflies wings color modeled with lamellar grating theory,” Opt. Express 9, 567–578 (2001). [CrossRef]   [PubMed]  

23. C. Lawrence, P. Vukusic, and R. Sambles, “Grazing-incidence iridescence from a butterfly wing,” Appl. Opt. 41, 437–441 (2002). [CrossRef]   [PubMed]  

24. J. P. Vigneron, K. Kertesz, Z. Vertesy, M. Rassart, V. Lousse, Z. Balint, and L. P. Biro, “Correlated diffraction and fluorescence in the backscattering iridescence of the male butterfly Troides magellanus (Papilionidae),” Phys. Rev. E 78, 021903 (2008). [CrossRef]  

25. P. Vukusic, J. R. Sambles, C. R. Lawrence, and R. J. Wootton, “Limited-view iridescence in the butterfly Ancyluris meliboeus,” Proc. R. Soc. London Ser. B , 269, 7–14 (2002). [CrossRef]  

26. S. Kinoshita, S. Yoshioka, and J. Miyazaki, “Physics of structural colors,” Rep. Prog. Phys. 71, 076401 (2008). [CrossRef]