Introduction

Scientists and engineers have long been fascinated by the dazzling array of optical properties present in nature1,2,3,4,5,6,7,8,9,10,11,12,13. While coloration is a tangible reminder of the wonders of microstructures in animals, microstructures may also play a significant role in thermoregulation. Many animals take advantage of structural thermoregulation to thrive in diverse habitat climates around the world. Saharan silver ants, for instance, use surface microstructures to stay cool under desert sunlight by simultaneously reflecting solar heat and maximizing heat emission from their bodies in the mid-infrared (mid-IR) wavelengths14. Meanwhile, polar bears stay warm in the cold Arctic due to their dense fur, which is not just a warm blanket but is also a radiative layer that scatters IR radiation and minimizes heat loss15. As in other animals, we observe similar optical and thermal phenomena in butterflies, where wing microstructures affect not only coloration16,17,18 but also assist in the thermoregulation of these insects19,20. The variation in optical properties that facilitates thermoregulation in butterflies is the result of a combination of pigmentation and chitin-based photonic structures17,21,22,23,24,25.

The focus of the majority of optical and thermal studies on butterflies has been on visible to near-IR wavelengths (0.4–2.5 µm wavelengths) where absorption of heat from the incident sunlight takes place. Munro et al. 26 examined the correlation of the visible to near-IR optical properties of butterfly wings to the butterfliesʼ habitat temperature and precipitation, and showed that near-IR wing reflectance, rather than UV–Vis reflectance, is correlated with temperature. Nonetheless, our understanding of the optical properties of butterfly wing microstructures in mid-IR wavelengths beyond the visible to near-IR remains limited.

While the visible to near-IR wavelengths are important for heat gain via absorption of incident solar irradiation, the mid-IR wavelengths of 7.5–14 µm are critical for heat loss. The mid-IR wavelengths correspond to the wavelength range where thermal emission takes place as governed by Planck’s law27,28,29. Moreover, this wavelength region also hosts the window of atmospheric transmission30, which can lead to substantial radiative cooling in the ambient environment to outer space at 3 K (–270 °C)31,32. Importantly, this part of the spectrum (mid-IR) has not yet been systematically examined in butterflies although the effect of mid-IR radiative cooling has been examined in engineering systems31,33,34 and in other biological taxa14,35,36 in recent years.

Butterflies, being cold-blooded, regulate their body temperature using strategies ranging from behavioral adaptations to reliance on habitat and environmental conditions22,37,38,39,40. Thermoregulation is crucial in butterflies because they must maintain a body temperature of 20–50 °C in order to survive, regardless of where their habitat is located41,42,43. To attain a stable body temperature in such geographically and climatically diverse locations as the tropics and the Arctic, there is a need for passive thermoregulation in butterfly wings. Specifically, butterflies heat up their wings (and hence, bodies) by absorption of solar heat in the ultraviolet (UV) to near-IR wavelengths. Simultaneously, in order to avoid overheating, they cool down by utilizing re-emission of heat to outer space from the mid-IR wavelengths19,39,44. To effectively regulate their wing (and body) temperatures, butterflies utilize photonic structures2,21,24,45,46 to not only control visible and near-IR properties but also invisible mid-IR properties16,22,47,48 for thermoregulation. These microstructures may consequently be an adaptation of the different butterflies to their varying climates. Consequently, in order to understand how the butterflies thrive in diverse habitats, it is vital to understand the correlation between the microstructure-dependent mid-IR optical properties of the butterfly wings and their respective habitat climates.

Emerging studies have aimed to address our lack of understanding of the mid-IR optical properties of butterfly wings19,20,44 with increasing evidence being presented for the role played by wing microstructures in butterfly thermoregulation. Not only do the microstructures aid in the re-emission of heat from butterfly wings, but a matrix of living cells on butterfly wings coordinates with the microstructures to maintain a viable temperature for the butterfly20. Meanwhile, other findings have shown that different butterfly species have unique wing microstructures that are critical to how the mid-IR optical properties vary across butterflies from diverse climates19. These works have led to a greater knowledge of thermoregulation mechanisms in butterflies. However, there remains much to be learned about the adaptation of butterflies across vastly differing habitats to achieve thermoregulation via wing microstructures.

Differing environmental and climate conditions across habitats impose varying radiative demands on butterflies. Here our computational results for the optical properties of butterfly wings demonstrate the potential use of wing microstructures to regulate wing temperatures, enabling enhanced thermal performance and allowing for the survival of butterflies across the globe. The study of how nature adapts to extremes in climate is of vital importance to our own survival in such habitats. Recent advancements have been made in the fields of engineered materials and systems that make use of bio-inspired structures to achieve thermal control in various ambient climatic conditions49,50,51,52. By analyzing the thermoregulation strategies of butterflies around the world we might thus find better engineered solutions for our own lives.

Results

Structural control of spectral optical properties of butterfly wings

We investigated microstructure-driven optical properties by performing computations based on rigorous coupled-wave analysis (RCWA) and finite-difference time-domain (FDTD) methods. Both methods of analysis solve for Maxwell’s equations, with choices being dependent upon time and computational complexity requirements. We primarily used RCWA for parametric analysis and FDTD for validation32,49,52.

The computational models considered here have been used extensively in the optical physics literature53,54,55,56,57, and have also been previously used to analyze the experimentally determined mid-IR optical properties of butterfly wing microstructures19,20. The difference between measured emissivity and simulations does not exceed 0.1 in units of absolute emissivity, and our RCWA and FDTD simulations produce comparable results19, which validates their use in the current evaluation. Our studies make use of existing microstructure dimension data using methodology established in prior findings validated experimentally (Table S1, Figure S1)19,21,58,59,60,61,62,63,64 to newly evaluate the wing optical properties. The dimensions are used as input parameters for the RCWA computation, along with the refractive index and extinction coefficients of the materials involved (for instance, chitin65,66,67,68). We consider that the microstructures for all butterfly species are made of chitin, with the same refractive index values across the mid-IR wavelength range65,66,67,68. We find that in the mid-IR wavelengths, the butterfly wings attained emissivity values around 0.24–0.60 (Fig. 1). While butterflies are known to have regions of varying visible coloration across different regions of their wings, previous work has shown that regions of varying visible color have highly similar (or nearly identical) mid-IR emissivity values19. The likely reason for this is that the visible coloration and mid-IR emissivity arise from phenomena observed in two different wavelength ranges16,19,20,44,47. While visible coloration occurs in the 300–700 nm wavelength range, the mid-IR wavelength range extends from 7.5 µm to 14 µm wavelengths16,19,20. The interaction of light with photonic structures (such as the butterfly wing mesh) is largely scale-dependent; light of a given wavelength interacts strongly with structures of similar dimension scale. Consequently, the structural causes of coloration and mid-IR emissivity tend to be different in dimensional scale19,20,44,47. While structural coloration arises from diffraction and constructive interference caused by structures in the nanometer-scale, mid-IR optical properties are influenced by structures 10–100 times larger19,47. Diffraction and interference, primary causes of coloration and changes in emissivity, are thus highly wavelength selective responses of light interacting with photonic structures19. The varying values of emissivity observed in the mid-IR wavelengths could result in varying thermal performance of the butterfly wings under different environmental thermal conditions31,32,69,70,71. The findings indicate a possibility for radiative thermoregulation in butterfly wings by structural modulation of spectral emissivity.

Figure 1
figure 1

Butterfly habitat annual average air temperature and butterfly wing emissivity (a) Geographical mapping of the sampled butterfly species from around the world, generated using open source QGIS software v.3.20.3 (www.qgis.org). The coloration of the map indicates the annual average ambient air temperature90. The text adjacent to the data points indicates the habitat-related Köppen-Geiger climate classification72. Here, A is a tropical climate, B is arid, C is temperate, and D is continental, with W being a desert, w indicating a dry winter, s a dry summer, f meaning no dry season, a, b, and c pointing to hot, warm, and cold summers, and h indicating an overall hot climate year-round. (b) Correlation of average mid-IR emissivity of the wings of butterflies with annual average habitat air temperature. The computed mid-IR emissivity was averaged over 7.5–14 µm for each butterfly wing based on reported microstructure data in the literature19,21,58,59,60,61,62,63,64. The colored letters in the figure correlate to the individual emissivity data in (c). The error bars represent the standard error for each data point, with the temperature data being taken from 3–8 different locations for each species. A linear fit to the data gives a correlation of εmid-IR = 0.0103Tair – 2.5867 (the dotted line represents the fit, with the gray band representing the 90% interval from the fit). The corresponding coefficient of correlation is + 0.86 and the coefficient of determination is + 0.74. (c) Computational emissivity predictions for the wing structures of various butterfly species from around the world. The plot depicts the emissivity in the mid-IR wavelengths, up to 14 µm. The optical property values were computed based on structural dimensions from annoying electron microscopy (SEM) images (sample size ranging from 10–100 unique wing locations) in existing literature19,21,58,59,60,61,62,63,64. Complete genus and species names are given in Fig. 4.

Full size image

Correlation of spectral optical properties with environmental conditions

To examine the relationship between the optical properties of butterfly wings and the butterflies’ habitat climatic conditions we focus on the spectral average emissivity values in an understudied wavelength range of significance to thermoregulation. Specifically, the mid-IR spectrum (from 7.5 µm to 14 µm wavelengths, the atmospheric transmission window) is expected to offer radiative cooling of the butterflies by re-emission of heat from the butterfly wings to outer space31,32. We thus compared the average mid-IR emissivity of the butterfly wings to the annual average air temperature of the corresponding habitat range of the butterflies.

Correlation of mid-IR emissivity with annual average air temperature

We first used the Köppen-Geiger climate classifications72,73 to take into consideration the monthly average, maximum, and minimum values of air temperature and precipitation across all locations of the world (Fig. 1a). The use of monthly average, maximum, and minimum values reduces uncertainties that arise by simply considering annual average temperature due to the fluctuation in the values recorded temporally and spatially.

We then looked at the individual habitat air temperature data for each of the butterflies in more detail. In general, the air temperature varied with the latitude of measurement and the altitude of the location73. For instance, the northern plains in India experience a different air temperature (corresponding to a different climate classification) compared to the elevated regions of the Tibetan plateau, although both locations are along the same latitude. In order to avoid spatial variations in the air temperature data, the reported values were taken as the average of the reported data across at least 3 different weather stations within the geographic range of each of the butterflies74,75. The temporal variations in data were avoided by considering the annual average ambient air temperature for the analysis.

We then compared the mid-IR emissivity of the butterfly wings with the corresponding annual average air temperature values of their respective habitats (Fig. 1b). The mid-IR emissivity values are averaged for each butterfly within the wavelength range of 7.5–14 µm. We found a correlation between mid-IR emissivity and annual average air temperature, in which linear regression yielded a coefficient of correlation (R) of + 0.86. The corresponding coefficient of determination (R2) is + 0.74, indicating that 74% of the changes in the mid-IR emissivity can be explained by changes in the air temperature.

The mean value of the mid-IR emissivity for the entire dataset was 0.40, with a standard deviation of 0.13, with all values ranging between a minimum and maximum of 0.21 and 0.60 respectively (Fig. 1c). The butterfly Celastrina echo has an average mid-IR emissivity of 0.21, with an average habitat air temperature of less than 280 K (7 °C) annually (climate classification: Dfc)72,73,74,75. However, Heliconius sara has an average mid-IR emissivity of 0.60 and its annual average air temperature exceeds 290 K (17 °C) (climate classification: Af)72,73,74,75. We hypothesize that the mid-IR emissivity offers cooling by re-emission of heat from the butterfly wings across varying habitat climates19,44. Higher values of mid-IR emissivity result in greater heat loss from the wings due to re-emission of heat to outer space (at a temperature of 3 K or –270 °C). The results show a correlation between the butterflies’ habitat annual average ambient air temperature and the mid-IR emissivity of the butterfly wings.

We also analyzed the possible correlation between mid-IR emissivity and precipitation, another factor considered in the Köppen-Geiger climate classifications. The analysis showed wide variations in precipitation globally, with no evident direct link to either latitude or altitude74,75. Butterfly wings have previously been extensively studied for their hydrophobic characteristics76,77,78. However, thus far there is no established evidence for thermal adaptation to precipitation. We analyzed the annual average precipitation across the various habitats and plotted them on a relative scale from low to high (light blue to deep red) (Fig. 1b). Based on the analysis, there was no direct link between mid-IR emissivity and precipitation. While studies document the overall hydrophobic nature of butterfly wings, the dominant geometric parameters controlling hydrophobicity could differ from those controlling spectral emissivity76,77,78.

Effects of seasonal variations in air temperature

In similar fashion, we compared the butterfly wing mid-IR emissivity with summer and winter average habitat air temperatures74,75. In addition to the seasonal analysis, we also analyzed the air temperatures of the months during which each of the butterflies are found in abundance in their habitats19,21,58,59,60,61,62,63,64. The results (Fig. 2) indicate correlations of the mid-IR emissivity of their wings to the summer (Fig. 2a), the winter (Fig. 2b), and abundant months’ (Fig. 2c) average air temperature. The coefficient of correlation for the summer data was + 0.90, the coefficient for the winter data was + 0.82, indicating that 81% and 67% of the variation in the mid-IR emissivity, respectively, can be explained by variation in summer and winter average air temperatures. It was interesting to observe that the coefficient of correlation was higher for summertime, during which many butterflies are active. The coefficient of correlation for the butterflies’ abundant months was even higher, at + 0.95, indicating that 90% of the changes in mid-IR emissivity can be explained by variation in air temperature. We deduce that there is a strong correlation between the abundant months’ habitat air temperature and the butterfly wing mid-IR optical properties.

Figure 2
figure 2

Comparison of the butterflies’ average wing mid-IR emissivity in relation to their habitat average air temperature during summer, winter, and most abundant months. (a) Average mid-IR emissivity of butterfly wings in relation to the summer average air temperature of their habitats (June–September in the northern hemisphere and December–March in the southern hemisphere). The linear correlation corresponds to εmid-IR = 0.0116Tair – 3.0354, with a coefficient of correlation of + 0.90. (b) Average mid-IR emissivity of butterfly wings in relation to the winter average air temperature of their habitats (December–March in the northern hemisphere and June–September in the southern hemisphere). The linear correlation corresponds to εmid-IR = 0.0077Tair – 1.7892, with a coefficient of correlation of + 0.82. (c) Average mid-IR emissivity of butterfly wings in relation to the average air temperature of their habitats during the butterflies’ most abundant months (Supplementary Table S1)74,75. The linear correlation corresponds to εmid-IR = 0.0114Tair – 2.9091, with a coefficient of correlation of + 0.95. In all the plots, the gray bands represent the 90% interval from the fit. The plots use the same methods of calculation and the same sources of information as Fig. 1. The error bars depict the standard error for each data point, with the temperature data ample from 3–8 different locations for each species. The colored letters in the figures correlate to the individual species' emissivity data in Fig. 1c.

Full size image

Apart from seasonal variations, we also compare the effects of diurnal changes in air temperature (Figure S2). Though the results indicate noticeable correlations between daytime and nighttime air temperatures and mid-IR emissivity for the butterfly wings, it must be noted that butterflies are primarily active during daytime and spend the nights dormant. The nighttime air temperatures are also majorly influenced by daytime temperatures, and not much of significance may be deduced from such comparisons.

Correlation of mid-IR emissivity with other climatic factors

In addition to analyzing the correlation between the optical properties of butterfly wings and the air temperature of their habitats, we also compared the mid-IR emissivity values for the butterfly wings with their respective habitat annual average wind speed, annual precipitation, and altitude (as altitude can be related to the air pressure) values74,75 (Fig. 3). The annual average wind speed values for the various habitats ranged between 2.4–6.4 ms−1. We observed (Fig. 3a) no significant correlation between the annual average wind speed and the mid-IR emissivity of the butterflies, with the coefficient of determination being 0.13 and a coefficient of correlation of -0.37. Similarly, the annual precipitation values ranged from 200–2500 mm. There was no noticeable correlation of the mid-IR emissivity to the annual precipitation (Fig. 3b), with a coefficient of determination of 0.01 and a coefficient of correlation of + 0.09. Finally, we considered the habitat altitude, which may be related to the air pressure. The butterflies inhabit locations ranging from 0–2000 m above sea level. Their wing mid-IR emissivity values do not show a correlation to the habitat altitude (Fig. 3c), with a coefficient of determination of 0.03 and a coefficient of correlation of 0.16. The lack of correlation confirms that wind speed, precipitation, and altitude do not play dominant roles in radiative heat transfer phenomena of butterfly wings in the process of re-emission of heat from their wings. We thus conclude that habitat air temperature extensively influences butterfly thermoregulation among the various climatic factors examined.

Figure 3
figure 3

Comparison of average mid-IR emissivity of the butterfly wings to various climatic factors. (a) Average mid-IR emissivity of butterfly wings with respect to the annual average wind speed in their habitats74,75. The linear fit represents an equation of εmid-IR = −0.0381wind speed + 0.5376, with a coefficient of determination of 0.13 and a coefficient of correlation of -0.37. (b) Average mid-IR emissivity of butterfly wings with respect to the annual average precipitation in their habitats74,75. The linear fit represents an equation of εmid-IR = −1.0 × 10−5 precipitation + 0.3849, with a coefficient of determination of 0.01 and a coefficient of correlation of + 0.09. (c) Average mid-IR emissivity of butterfly wings with respect to the altitude in their habitats74,75. The linear fit represents an equation of εmid-IR = 3.0 × 10–5 altitude + 0.3852, with a coefficient of determination of 0.03 and a coefficient of correlation of + 0.16. The plots show a lack of correlation between habitat wind speed, precipitation, or altitude, and the optical properties of the butterfly wings. The plots use the same methods of calculation and the same sources of information as Fig. 1. The error bars depict the standard error for each data point, with the climatic data sampled from 3–8 different locations for each species. The colored letters in the figures correlate to the individual species' emissivity data in Fig. 1c.

Full size image

Phylogenetic analysis

While a direct comparison of the mid-IR optical properties of butterfly wings to the butterflies' respective habitat air temperatures shows a positive correlation, it is important to examine whether or not this correlation still holds when the phylogenetic relatedness of the butterflies is taken into account. In order to ascertain that these traits evolved as adaptations to the butterflies’ habitat climatic conditions rather than as a result of phylogenetic inertia, or the tendency of related organisms to have related traits, we performed a phylogenetic independent contrasts analysis79. We utilized the methodology of Felsenstein’s contrasts79,80,81 by comparing the phylogeny-corrected mid-IR emissivity with the average annual air temperature.

First, we made use of a phylogenetic tree available from existing literature82. Several of the species whose optical properties we examined are present on the tree. However, where phylogenetic data were not available for a specific species, we substituted our studied species for closely related or sister species on that tree. The substitutions were made with species proximity taken into account—for example, sister species were substituted, failing which, substitute species were taken within the same tribe. The following substitutions on the tree were made for those species with no directly available phylogenetic data: Papilio machaon in place of Papilio rumanzovia, Troides helena in place of Troides rhadamantus, Archaeoprepona demophoon in place of Prepona dexamenus, Heliconius sara and Heliconius doris in place of Heliconius melpomene, Euploeia mulciber in place of Danaus plexippus, Celastrina argiolus and Celastrina echo in place of Hemiargus ceraunus, Danis danis in place of Lepidochrysops patricia, Hypochrysops delicia in place of Lucia limbaria, Chrysozephyrus brillantinus in place of Artopoetes pryeri, and Arhopala japonica in place of Arhopala metamuta.

A trimmed phylogenetic tree (Fig. 4) was then imported into Mesquite83 for independent contrasts analysis using the PDAP:PDTree package. For the species under consideration, we mapped their mid-IR emissivity and the annual average air temperature of their habitat on the trimmed tree. We then ran a Felsenstein’s contrasts analysis on the pair of characters, in order to deduce the phylogeny-corrected correlations between the mid-IR emissivity and air temperature. The use of the method was first validated by checking for the correlation (or lack of) between the contrasts of the characteristics and the standard deviation for the data (Figure S3). Once the use of the method was validated, we then analyzed the phylogeny-corrected correlations between the optical properties and the habitat climatic conditions. In order to account for Heliconius melpomene and Hemiargus ceraunus being substituted for 2 species each, we performed 4 sets of analyses for each possible permutation in order to obtain a range of possible results to account for the substitutions.

Figure 4
figure 4

Phylogenetic tree used for the independent contrasts analysis for the species under consideration, based on an existing molecular phylogeny in the literature82. Where phylogenetic data were not directly available for the species in our study, closely related or sister species on the original phylogeny were substituted for Felsenstein’s contrasts analysis, as indicated in subsection 3 of our results. The numbers on the branches and the scale bar represent the branch length values in the original phylogeny. Stars indicate substituted species. References for species' SEM data are given as superscripts next to the species names.

Full size image

The results (Figure S4) depict a correlation between the contrasts for the mid-IR emissivity and annual average air temperature. For the correlation, Felsenstein’s contrasts correlation yielded an R2 value ranging between + 0.68 to + 0.71. The findings thus show that 68–71% of the variation in the mid-IR emissivity may be explained by variation in the air temperature. This result suggests that the optical properties of the butterflies’ wings evolved in response to the climatic conditions that they inhabit.

Discussion

Our findings reveal one way in which butterflies from diverse geographic and climatic regions have thrived in their environment. Previous works have observed that butterflies make use of radiative heat transfer via their wing surfaces (especially in the mid-IR wavelengths), which aids in maintaining a wing temperature within a habitable temperature range19,38,67. However, while earlier studies considered individual butterfly species within the context of their own habitat climates or microclimates, our present study attempted to evaluate global trends in butterfly wing mid-IR optical properties in relation to climatic conditions worldwide. Our findings suggest that the mid-IR emissivity values of butterfly wings are correlated with the habitat air temperatures for the butterflies. That is to say, butterflies from regions with higher annual average air temperatures possess higher average mid-IR emissivity values than butterflies from regions with lower annual average air temperatures.

Our analyses also evaluated whether or not the microstructures responsible for the varying mid-IR emissivity values are the result of evolutionary processes leading to the adaptation of these butterflies to their habitat climatic conditions. A Felsenstein’s independent contrasts analysis verified the presence of a link between the mid-IR optical properties of the butterflies and the corresponding annual average air temperature in their habitats. By controlling for the effects of phylogenetic inertia on this correlation, we infer that wing microstructures and their consequent mid-IR emissivity values are at least in part the result of natural selection-driven adaptation to ambient air temperature.

Across the habitat climate conditions examined, butterfly wing mid-infrared emissivity appears to aid in butterfly thermoregulation by enhancing or inhibiting heat loss. We have attempted evaluations of solar spectrum reflectivity (Figures S5, S6) of the butterfly wings, and observe a weak correlation of around 16% with annual average solar irradiation (Fig. S7). Further studies of the microstructure-dependent optical and thermal properties of the butterflies such as the current one may facilitate the development of bio-inspired optical metamaterials and photonic structures that aid in radiative thermal management31,70,84,85,86.

Butterflies have adapted to the varying thermal demands put on them by making use of varying microstructures. Similarly, we may be able to overcome inherent material-dependent limitations on the optical properties and thermal behavior of various materials we work with by making use of surface microstructures to cater to varying optical and thermal needs. For example, the microstructure-dependent surface emissivity may be optimized for radiative heating or cooling by controlling the optical properties across various regions of the electromagnetic spectrum. An ideal radiative cooling surface, where the emissivity is selectively maximized in the mid-infrared wavelength region32,69 will advance current breakthroughs in thermal management across a wide array of applications—from space systems to wearable devices, and solar panels to building thermal management via optical coatings87. Similarly, an ideal radiative heating surface, where the emissivity is selectively maximized in the solar spectrum wavelengths will improve upon current solutions for just as wide an array of applications—ranging from thermophotovoltaics88 to personal thermal management via light-weight insulated heating fabrics for Arctic and Antarctic expeditions and space explorations89.

The results reported here could thus lay the ground for novel thermal management solutions and ultimately connect us closer to nature. Learning from the adaptations of butterflies to their habitat climates, we could invent and improve upon means for our own survival across extreme climates around the world.

Conclusions

This work provides a new understanding of the relationship between the microstructure-dependent mid-infrared emissivity of butterfly wings to the butterflies’ habitat climates. The emissivity of butterfly wings in the mid-infrared wavelengths of 7.5–14 µm was computed using rigorous coupled-wave analysis and finite-difference time-domain methods based on existing microstructure microscopic images, a method previously validated experimentally19. The mid-infrared emissivity was then correlated to the butterflies’ habitat climate, via comparisons with the annual average air temperature, precipitation level, and wind speed. Apart from a comparison between mid-IR emissivity and the annual average air temperature, our analysis also examined the effects of diurnal and seasonal changes in climate. With butterflies being found in abundance in selective months of the year, and most being preferentially active during daytime, we observed that the correlation of these butterflies’ mid-IR emissivity values with the average air temperature for the time of year during which they are active was especially strong. We found a 90% correlation of the mid-infrared emissivity with the average air temperature for the months in which the butterflies are found in abundance.

In general, we noted that the mid-IR emissivity values of butterfly wings from warmer climates (with higher annual average air temperatures) are higher than those of butterflies from cooler climates (with lower annual average air temperature), which could be a habitat-dependent-thermoregulatory adaptation of the butterflies. To rule out the effects of phylogenetic inertia and to examine the role of evolution in the adaptations of these butterflies, we performed a Felsenstein’s independent contrasts analysis. The phylogeny-corrected data suggested that up to 71% of the variation in the mid-IR optical properties of the butterflies could be directly explained due to the corresponding annual average air temperature in their habitats.

Our findings suggest that butterfly wing mid-infrared emissivity plays a critical role in butterfly wing thermoregulation and that the ability of butterflies to adapt to climate change will be limited by the speed with which wing microstructures can be modified evolutionarily.

Materials and methods

Computation of optical properties

This work computes the spectral optical properties of butterfly wing microstructures using the rigorous couple wave analysis (RCWA) and finite-difference time-domain (FDTD)57 methods32,49,52. The RCWA method considers topographical variations and handles rigorous solutions of Maxwell’s equations27,28,29. Direct results from RCWA yield scattering matrices in the forward and reverse directions, from which reflectivity (ρ) and transmissivity (τ) are computed. The emissivity is assumed identical to the absorptivity by Kirchhoff’s law27,28,29 and is computed from the reflectivity and transmissivity (α = 1−ρ−τ). The FDTD method discretizes the samples and solves for space- and time-variant Maxwell’s equation for each unit cell. RCWA is semi-analytical and treats the waves and fields as sets of gratings, hence making it very effective for mesh-like structures. On the other hand, FDTD is superior in modeling curved surfaces and spherical shapes in full three-dimension. While FDTD is time-consuming, RCWA is more efficient towards optimizing the designs by variation of geometrical parameters. Here we used FDTD and RCWA as means for validating each other’s results.

In the work presented here, the butterfly wing microstructures are modeled as periodic arrays of mesh-like structures, whose structural parameters and dimensions are listed in Supplementary Section 1, Table S1. The spectral permittivity and complex refractive index of the chitin65,66,67,91,92 for the wings are input optical parameters for the spectral computations. Both the structural parameters (ridge periodicity (a), cross-link periodicity (b), ridge width (c), cross-link width (d), and height (e), detailed in Table S1 of the Supplementary Information) and optical parameters (spectral permittivity and complex refractive index of chitin) are inputted into custom built MATLAB (version 2020b) code in the case of RCWA, and commercially available Lumerical software (version 2020a) in the case of FDTD. The models are first described as mathematical arrays constituting the different materials that comprise the microstructures, and then the optical properties of the materials in each structural region are specified, which is comparable to building an outline first, followed by filling it in with material.

The simulation models are enclosed by boundaries that conform to the structural limits of the models themselves on four sides (along the X- and Y-axes)—the X–Y dimensions being defined by the ridge periodicity (a) and the cross-link periodicity (b). The remaining two sides, the top and bottom boundaries, extending for 10 times the model height (e), above and below (along the Z-axis) respectively. The boundary conditions used are as follows: perfectly matched layer boundary conditions at the top and bottom boundaries to effectively negate any erroneous reflection of light back into the model, and on the 4 sides of the model, periodic boundary conditions to emulate periodic repetitions of the mesh-like microstructures. The model was illuminated by a plane wave light source, with the wavelength ranging from 0.2 µm (200 nm) to 20 µm, although we primarily focused on the mid-IR wavelengths of 7.5–14 µm for our analyses in the work above. The models are discretized into meshes of 20 nm by 20 nm squares.

Once built, running the simulation yields two matrices as results, a forward propagation matrix and a backwards propagation matrix, respectively signifying the transmissivity and reflectivity of the structure. We then infer the emissivity as unity minus the sums of transmissivity and reflectivity27,28. The modeling code and corresponding optical properties can be accessed from the Materials and Correspondence section.