The brilliancy and variety of structural colors found in nature has become a major scientific topic in recent years. Rapid-prototyping processes enable the fabrication of according structures, but the technical exploitation requires a profound understanding of structural features and material properties regarding the generation of reflected color. This paper presents an extensive simulation of the reflectance spectra of a simplified 2D Morpho butterfly wing model by utilizing the finite-difference time-domain method. The structural parameters are optimized for reflection in a given spectral range. A comparison to simpler models, such as a plane dielectric layer stack, provides an understanding of the origin of the reflection behavior. We find that the wavelength of the reflection maximum is mainly set by the lateral dimensions of the structures. Furthermore small variations of the vertical dimensions leave the spectral position of the reflectance wavelength unchanged, potentially reducing grating effects.
©2012 Optical Society of America
1. Introduction and motivation
A wide range of species such as butterflies [1–3], moths , peacocks , beetles  or fungi  utilize structural features to generate brilliant colors. Their unique intensity, gloss or angular color isotropy together with the absence of bleaching are major topics of scientific interest [8,9]. Some animals are even able to vary their color when exposed to stress by adapting the distance between structures or by varying the amount of liquid in gaps between them . The concept of structural colors is best explained in contrast to pigments. Pigmentary color is based on absorption or emission of light. In contrast, the spatial assembly of structural components plays the central role for structural color [10,11]. Furthermore, plasmonic effects can give rise to intense color effects . The practical use of structural colors could reach from textiles and decoration , solar cells  to sensors . Inspired by natural concepts and the large diversity of structural colors, the question arises how the spectral reflectance of structural colors can be tailored.
The well-established finite-difference time-domain (FDTD) method provides an excellent tool to analyze the interaction between electromagnetic waves and optical structures. Only a few papers such as Ref [16–18] connecting FDTD simulations with structural colors of Morpho butterfly wings exist.
This paper focuses on an extensive parameter study of a structure inspired by the wing scale structure of the Morpho rhetenor butterfly. The structure is idealized to a rectangular scalable 2D-model which includes all important geometrical parameters. As the reflection spectra are strongly connected to the shape of the structure, those parameters are optimized to achieve selective reflection in a given color range. The optimized spectrum is compared to the reflection spectrum of the Morpho rhetenor butterfly.
By comparison to simpler models such as a planar layer system the influence of various structural components on the spectrum was analyzed. Finally, the influence of the different parameters on the reflection spectrum and their potential relation to the angular color isotropy of Morpho butterflies is discussed. It was found that the wavelength of a narrow reflectance peak can be easily adjusted by a variation of geometry parameters defining the horizontal dimensions. A small variation of vertical dimensions leaves the wavelength of the reflectance peak unchanged. By randomizing these parameters this behavior can be used for enhancing the angular color isotropy of such structures reducing grating effects.
Based on the large amount of analyzed data the presented comprehensive parametric study is beyond a simple investigation of the structural features of the Morpho butterfly. The simulation was used as a tool for studying both, individual features of the intricate butterfly wing structure and the coaction of the combined structure with respect to the resulting spectral reflectance. It must be pointed out that this work was not only done for the feature sizes and shapes found in the real butterfly wing, but also for derived geometries taking into account fabrication issues of such structures. Therefore, a systematic variation of parameters determining such a derived geometry enabled the collection of spectral data in a database or look up table that correlates color and structure. This is a prerequisite and of highest importance regarding the fabrication of structures for a tailored spectral reflectance, e.g. a desired color impression. Due to technical limitations of lithographic processes, it is not possible to manufacture any arbitrary designed structure and compromises have to be made. Hence, we have put the focus of the simulations and the database on the comprehensive investigation of a class of structures derived from the Morpho butterfly that can on one side cover the whole visible spectral range and on the other side can be fabricated by e.g. 3D laser lithography.
2. Materials and methods
For the following simulations the commercial software tool FDTD-solutions  was used to simulate the reflection spectrum of a simplified structure similar to the “bookshelf-structure” of the wing scales of the Morpho rhetenor butterfly. The analysis and graphical representation of the simulation data was done with Matlab .
The basic 2D-model consists of a central rectangle (“pillar”) and attached rectangles on both sides (“shelves”). All important geometric parameters are varied within a certain range (Table 1 ). In horizontal direction the model is extended to infinity with periodic boundary conditions (Fig. 1 ). Periodic boundary conditions in combination with a 2D-model (extending the model structure infinitely perpendicular to the projection plane) reduce the simulation time significantly, making a wide range parameter screening possible. The boundary conditions in vertical direction are absorbing (perfectly matched layer, PML). The refractive index of the polymer, which will be used for a two-photon absorption rapid prototyping process in future experiments, shows only a weak dependence on the wavelength in the visual range (n = 1.54 at 300 nm and n = 1.52 at 700 nm) as determined from ellipsometric measurements (Vase Ellipsometer, J.A.Woollam Co. Inc.). Hence the simulations were carried out assuming a constant value of 1.53 + 0.005i, corresponding to the refractive index at 500 nm. For comparison, the refractive index of the cuticle of the Morpho butterfly is 1.56 + 0.06i. Furthermore, we assumed a SiO2-substrate with a refractive index taken from Palik . A collimated source is located above the structure and emits a spectrally broad (300 - 800 nm) plane transverse magnetic wave  with a Poynting vector (red arrow, Fig. 1) perpendicular to the substrate surface. The spectral reflectance R was determined for each wavelength by calculating the surface integral of the real part of the Poynting vector over the detector surface (Fig. 1) and normalizing it to the incident source power. This means that for periodic simulation regions with a detector extending through the simulation region the whole angular range is included in the reflectance. The accuracy of the FDTD-mesh was chosen to give a good tradeoff between computer memory requirements and simulation time while ensuring the convergence of the results. The source distance (300 nm) and the detector distance (400 nm, to the upper part of the pillar) were chosen to be close to the structure to avoid errors due to large simulation regions .
Initially all spatial parameters as defined in Fig. 1 were set to 150 nm and 7 shelves attached to each pillar were considered. An independent variation and optimization process was then performed iteratively for each of the parameters. The light spectrum was divided into six spectral ranges, corresponding to the colors violet, blue, green, yellow, orange and red. As detailed below, the first structural parameter was then varied between 0 and 1000 nm and for each color range the total spectral reflectance was calculated. The total reflectance of the color to be optimized was then divided by the total reflectance of the other color ranges. After that the maximum value of these ratios was determined and the corresponding structural parameter set for the next iteration step. Subsequently, the entire model parameters were varied one after the other finally yielding an optimized set of parameters (Table 1). As it was not possible to minimize the reflection in all other spectral ranges simultaneously, we focused on maximizing one color while minimizing its complementary color. A restriction to geometry values above 100 nm was chosen due to the limited resolution of the rapid-prototyping process. The optimized spectrum for the blue (λcenter = 450 nm, Δλ = ± 30 nm) color range was the basis for a further analysis of the system.
The focus in the second part of the paper lies on analyzing the origin of the reflectance maximum of the optimized spectrum. In order to reveal the mechanisms of color generation, the resulting structure was compared to simpler models by using a “reverse engineering process”. This was done by inducing symmetry breaks to a planar layer system and proceeding to more complex structures (Fig. 2 ) while simulating the corresponding reflection spectra. The results are plotted as 2D contour plots, where the x-axis shows the varied parameter, the y-axis displays the light wavelength and the color indicates the reflectance R.
In the third part we propose a different way to analyze the influence on the angular color isotropy of the Morpho rhetenor butterfly. The central topic is how and if it is possible to vary certain parameters without changing the reflectance behavior. Starting from the color-optimized spectrum, the spatial parameters (Fig. 1, Table 1) were screened separately between 0 and 1000 nm. Depending on the behavior of the reflectance maximum and their potential connection to the angular color isotropy the parameters were categorized into two groups.
We conclude with a far field simulation involving 100 randomized structures comparing them to the results of the previous chapter.
3. Results and discussion
3.1. Parameter optimization leading to a reflection in the blue color range
The rightmost column in Table 1 summarizes the optimized structural parameters for reflection in the blue color range (420 - 480 nm). Examining the reflection spectrum (Fig. 3(a) ) of such a structure, we find a sharp peak at around 445 nm, reflecting almost 80% of the incident light. Furthermore, a broader feature exists centered at approximately 700 nm. Confirming the optimization process, the spectrum reaches a minimum at the yellow spectral range (560 - 580 nm). The intensity ratio color / complementary color plays an important part for the intensity of the color perceived by the human eye . As a comparison the reflection spectrum of a real Morpho rhetenor butterfly (Fig. 3(c)) was measured for an angle of incidence of 0° to the surface normal of the plane of projection and various detector angles. The measurements were done with an Instrument Systems GON 360 Gonometer connected via a fiber optic cable to an Instrument Systems CAS 140 CT Compact Array Spectrometer. Inserting the structural parameters and the refractive index of the Morpho rhetenor butterfly  into our model gives a simulated reflection spectrum (Fig. 3(b)) which agrees well with the experimental results (Fig. 3(c)). The simulated spectrum shows a primary peak at 430 nm and also covers the shoulder of the peak close to 500 nm. Furthermore the reflection in the red color range is suppressed.
3.2. Investigation of the optimized reflectance spectrum by comparison to simpler models
In order to investigate the origin of the reflection peak in Fig. 3(a) a series of over 1200 FDTD simulations was done by means of a “reverse engineering process”. Starting from a simple dielectric layer stack the system was transformed into the final “bookshelf structure” structure (Fig. 2(e)). First, a continuous layer system consisting of 3 layers (of the same material which was used above) with gaps of air in between (Fig. 2(a)) was analyzed. The gap between the layers (“layer distance”) was chosen to be constant at 255 nm according to the optimized structural parameter “shelf distance”. Afterwards, the height of the layers was varied between 0 and 1000 nm and the simulated reflection spectra displayed in a contour-plot (Fig. 4(a) ). The contour plot shows the typical broad reflection pattern due to multilayer interference. As a confirmation of the FDTD simulations the spectra were compared to an analytical transfer-matrix-algorithm model and the results were found to be consistent. Comparing the contour plot of the planar layer stack with the optimized reflectance spectrum (Fig. 3(a)) we find that the reflection features are wider and located at different wavelengths. We can conclude that the optimized spectrum cannot only come from the alternating layers but there must be a dependence of the other structural components.
The white vertical lines highlighted in all of the contour plots correspond to the parameters in Table 1 and furthermore define the starting conditions for the following transformation step. In Fig. 4(a) for example the white line indicates the spectrum where the “layer height” has reached the optimized value for the parameter “shelve height” of 140 nm (Table 1). This value is then set as the starting value and the next transformation step is performed.
Cutting the planar system vertically into two halves or “shelves” (Fig. 2(b)) and thereby introducing a vertical offset transforms the 1D into a 2D photonic crystal system . The contour plot displays the spectral reflectance in dependence of the offset (Fig. 4(b)). At low offsets regular elliptic shaped patterns exist throughout the Fig. with a wavelength similar to the one of the undisturbed planar system. The closer the shelves get to the centered position (white line, shelf offset = (shelf height + shelf distance) / 2 = 197.5 nm) the more irregular the features get, especially at low wavelengths. At this transformation step we see the narrow reflection peak (Fig. 3(a)) for the first time. Before performing the next transformation, to get the correct number of five shelves, the lowermost shelf was removed from the simulation.
By increasing the distance to the boundaries of the simulation region (Fig. 2(c)) a parameter corresponding to “structure distance” was introduced. We find (Fig. 4(c)) that the narrow peak which was observed before is shifted to higher wavelengths while its intensity is reduced the further the structures are apart. With further increasing the structure distance a broader low-intensity reflection region appears at wavelengths between 600 - 700 nm. Simultaneously, the regular feature just below 500 nm fades away while keeping a constant wavelength. This behavior leads us to the assumption that the broad horizontal features in Fig. 4(c) exist due to the layer system and the diagonal features come from the vertical separation of the shelves corresponding to a photonic crystal structure. At distances of 200 nm, 600 nm and 900 nm three additional diagonal reflectance peaks exist which seem to be almost parallel to each other. In comparison to the results of a square lattice of dielectric columns in  a photonic crystal based reflection may be assumed. When reaching 300 nm - the optimized distance between the structures - the second peak has shifted to a wavelength of 320 nm, which is the starting condition for the next parameter simulation.
The separation of the shelves in the horizontal direction (Fig. 4(d)) shows a similar behavior. Diagonal “photonic-crystal-features” coexist together with the “planar interference-feature” at a wavelength of 460 nm, which continuously gets less distinctive the further the shelves are separated. At a distance of 155 nm, equal to the optimized parameter of “pillar width”, the first sharp reflectance peak has almost faded away and the second not yet appeared. The broader feature can still be found located between 600 and 700 nm. We find that the structure is at this point identical to the final structure, except that the refractive index of the central pillar is set to 1. Still no reflecting behavior identical to Fig. 3(a) can be found.
Increasing the height of the central pillar (Fig. 2(e)) lets the reflection peak slowly reappear (Fig. 4(e)) but shifted to a higher wavelength of 445 nm while the broader feature can now be found at wavelengths of 700 nm. The dashed line is now identical to the spectrum displayed in Fig. 3(a). Inspired by their behavior on a variation of “structure distance” we assume that both observed reflectance patterns are related to the photonic crystal structure.
A second variation of the vertical offset of the shelves, this time for the finalized structure (Fig. 2(f)), underlines the importance of the 2D photonic crystal structure. The reflectance at 445 nm disappears the closer the shelves are placed together, but a slight variation of the offset will not change the reflection spectrum (Fig. 4(f)).
3.3. Influence of structural parameters on the reflectance spectrum
As the influence of the different structural components was investigated before, the next part focuses on how certain parameters alter the reflection spectrum of the “bookshelf structure”. Each time starting with the optimized parameters as starting conditions, the six spatial variables in Table 1 were varied between 0 and 1000 nm. The effect of the different parameters on the reflection spectrum was again visualized by the use of contour plots. In all of the following plots the white dotted lines indicate the location of the optimized spectrum shown in Fig. 3(a). According to the behavior of the peak at 445 nm the variables can be basically classified into two groups, horizontal (x-dimension) and vertical (y-dimension) parameters.
Horizontal parameters show a similar linear wavelength-dependence as found before (Fig. 5(a) -5(c)). Hence, the location of the reflection peak can easily be tuned by adjusting the parameters “structure distance”, “pillar width” or “shelf width”. The intensity of the reflectance reaches a maximum at a value close to the optimized parameters. If the parameters are increased further, additional reflectance peaks appear having a lower total intensity. When observing the parameter variation of “shelf width” it occurs that two broader reflectance patterns intensify at widths of 300 and 800 nm.
Looking at vertical parameters like “shelf distance” and “shelf height” (Fig. 5(d)-5(e)), the wavelength of the peak is almost not sensitive to a variation of those variables but the intensity reaches a maximum approximately at the optimized parameters. Changing the parameter “shelf y0” does not significantly change the reflection spectrum.
According to Ref [10,17]. the randomization of the starting height (“shelf y0”, Fig. 5(f)) of the shelve structures on each pillar and the height of the pillar itself (“pillar height”, Fig. 4(e)) play an important part in canceling out the grating effects, improving color isotropy in a broad angular range. From our simulations we conclude that a slight variation of the distance and height of the shelves (Fig. 5(d), Fig. 5(e)) together with a randomization of the offset of the shelves (Fig. 4(f)) leads as well to an improvement of the color isotropy.
The narrow width (Fig. 3(a)) of the reflection peak has to be further investigated. The first time it appeared was when introducing a vertical shelf offset (Fig. 4(b)) to the planar layer system. In addition to that for the specific optimized parameter set the width of the peak also results from the small number of shelves. Increasing the parameter “shelf number” leads to a broadening of the reflectance region at 445 nm. Simultaneously, the pattern at 700 nm intensifies to reach reflectance values of approximately 0.8. Similar to a multilayer system with a large number of layers additional lower intensity side maxima appear. We conclude that the narrow shape of the peak is mainly influenced by the 2D photonic crystal structure and the low number of shelves.
3.4. Far field calculations for randomized structural parameters
To prove the influence of randomization on the angular distribution of the reflection spectrum, the electric far field (|E|^2) at a distance of 1 m and a wavelength of 445 nm of 100 “bookshelf structures” was calculated for several randomizations. The FDTD-model was adjusted in a way that for each separate tree the values of vertical parameters randomly differ from each other, given by the parameter Δr. We looked at three cases, a) no randomization, b) randomization Δr = ± 10 nm and c) randomization Δr = ± 25 nm for each of the parameters “shelf height”, “shelf distance” and “shelf y0”. Figure 6(b) and Fig. 6(c) clarify that the intensity at viewing angles θ other than the grating orders (at 0° and ± 35°) is significantly increased the larger the randomization gets.
4. Summary and outlook
The simplified structure based on the “bookshelf-structure” of the wing scales of the Morpho rhetenor butterfly was analyzed and optimized by means of the FDTD method. A spectral optimization led to a narrow reflection peak in the blue color range, reflecting 80% of the incoming light at a wavelength of 445 nm. To investigate the origin of the spectral reflectance a transformation of a simple planar layer system to the “bookshelf structure” was performed including a wide range parameter study of several parameters. This analysis led us to the assumption that the reflection features of the optimized spectrum mainly result from a 2D photonic crystal structure. A parameter study of the spatial parameters of the model clarified that the wavelength of the narrow reflection peak can easily be tuned by the variation of the parameters “structure distance”, “pillar width” and “shelf width”. A small randomization of the parameters “shelf distance”, “shelf height”, “shelf y0”, “pillar height” and “shelf offset” does not influence the reflectance wavelength and can improve the cancellation of grating effects. A far field simulation for randomized structures shows that at angles differing from the grating orders the electric field is intensified.
All of the simulations were done in respect to a rapid prototyping process utilizing two-photon absorption. The fabrication of such high-resolution structures with two-photon absorption requires profound knowledge of laser beam shaping to be able to reduce the voxel size polymerized by the laser focus. Computer-aided design software tools provide a perfect way to model randomized structures based on simulation results.
We gratefully acknowledge the financial support of the Austrian Research Promotion Agency (FFG) within the framework of project No. 824903 (MorphoColor).
References and links
1. A. L. Ingram and A. R. Parker, “A review of the diversity and evolution of photonic structures in butterflies, incorporating the work of John Huxley (The Natural History Museum, London from 1961 to 1990),” Philos. Trans. R. Soc. Lond. B Biol. Sci. 363(1502), 2465–2480 (2008). [CrossRef] [PubMed]
3. D. Pantelić, S. Curčić, S. Savić-Šević, A. Korać, A. Kovačević, B. Curčić, and B. Bokić, “High angular and spectral selectivity of purple emperor (Lepidoptera: Apatura iris and A. ilia) butterfly wings,” Opt. Express 19(7), 5817–5826 (2011). [CrossRef] [PubMed]
5. S. Yoshioka and S. Kinoshita, “Effect of macroscopic structure in iridescent color of the peacock feathers,” Forma 17, 169–181 (2002).
8. M. Kolle, P. M. Salgard-Cunha, M. R. J. Scherer, F. Huang, P. Vukusic, S. Mahajan, J. J. Baumberg, and U. Steiner, “Mimicking the colourful wing scale structure of the Papilio blumei butterfly,” Nat. Nanotechnol. 5(7), 511–515 (2010). [CrossRef] [PubMed]
9. J. D. Forster, H. Noh, S. F. Liew, V. Saranathan, C. F. Schreck, L. Yang, J.-G. Park, R. O. Prum, S. G. J. Mochrie, C. S. O’Hern, H. Cao, and E. R. Dufresne, “Biomimetic isotropic nanostructures for structural coloration,” Adv. Mater. (Deerfield Beach Fla.) 22(26-27), 2939–2944 (2010). [CrossRef] [PubMed]
10. S. Kinoshita, S. Yoshioka, and J. Miyazaki, “Physics of structural colors,” Rep. Prog. Phys. 71(7), 076401 (2008). [CrossRef]
11. A. R. Parker, “515 million years of structural colour,” J. Opt. A, Pure Appl. Opt. 2(6), R15–R28 (2000). [CrossRef]
12. J. Zhang, J.-Y. Ou, N. Papasimakis, Y. Chen, K. F. Macdonald, and N. I. Zheludev, “Continuous metal plasmonic frequency selective surfaces,” Opt. Express 19(23), 23279–23285 (2011). [CrossRef] [PubMed]
14. W. Zhang, D. Zhang, T. Fan, J. Gu, J. Ding, H. Wang, Q. Guo, and H. Ogawa, “Novel photoanode structure templated from butterfly wing scales,” Chem. Mater. 21(1), 33–40 (2009). [CrossRef]
15. X. Yang, Z. Peng, H. Zuo, T. Shi, and G. Liao, “Using hierarchy architecture of Morpho butterfly scales for chemical sensing: Experiment and modeling,” Sensor. Actuat. A-Phys 167, 367–373 (2011).
16. D. Zhu, S. Kinoshita, D. Cai, and J. B. Cole, “Investigation of structural colors in Morpho butterflies using the nonstandard-finite-difference time-domain method: Effects of alternately stacked shelves and ridge density,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 80(5), 051924 (2009). [CrossRef] [PubMed]
18. 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(2), 97–103 (2007). [CrossRef] [PubMed]
19. Lumerical, “FDTD Solutions,” http://www.lumerical.com/.
20. MathWorks, “Matlab,” http://www.mathworks.de/.
21. E. D. Palik, Handbook of Optical Constants (Academic Press, 1985).
22. J. B. Schneider, “Plane waves in FDTD simulations and a nearly perfect total-field/scattered-field boundary,” IEEE Trans. Antenn. Propag. 52(12), 3280–3287 (2004). [CrossRef]
23. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton University Press, 1995).