Theory of chirped photonic crystals in biological broadband reflectors

One-dimensional photonic crystals with slowly varying, i.e."chirped", lattice period are responsible for broadband light reflectance in many diverse biological contexts, ranging from the shiny coatings of various beetles to the eyes of certain butterflies. We present a quantum scattering analogy for light reflection from these adiabatically chirped photonic crystals (ACPCs) and apply a WKB-type approximation to obtain a closed-form expression for the reflectance. From this expression we infer several design principles, including a differential equation for the chirp pattern required to elicit a given reflectance spectrum and the minimal number of bilayers required to exceed a desired reflectance threshold. Comparison of the number of bilayers found in ACPCs throughout nature and our predicted minimal required number also gives a quantitative measure of the optimality of chirped biological reflectors. Together these results elucidate the design principles of chirped reflectors in nature and their possible application to future optical technologies.

Many of the brilliant colors throughout nature are produced via the interaction of light and microscopic structures, a phenomenon known as structural coloration [1][2][3]. In particular, certain species of beetles, butterflies, birds, spiders, and fish take advantage of structural color phenomena to achieve their striking appearance [4][5][6][7]. In many of these cases color arises even in the absence of pigments through the layered application of materials with differing refractive indices, structurally arranged so as to yield the desired color via light interference.
Perhaps the simplest such structure is a onedimensional photonic crystal, i.e. a stack of dielectric slabs with alternating high and low refractive indices whose periodicity leads to a band structure analogous to the electronic band structure of semiconductors [8]. The implication of this analogy is the existence of photonic stopbands in which incident light of certain wavelengths is highly reflected due to its inability to propagate within the crystal [9]. By preferentially reflecting a narrow range of wavelengths, photonic crystals can reliably produce a characteristic color when illuminated with white light.
A more sophisticated approach involves gradually changing the period of the crystal. These "adiabatically chirped" photonic crystals (ACPCs) act as a locally periodic crystal for a range of wavelengths and can therefore reflect a broad swath of colors strongly, allowing for reflection spectra more complex than monochromatic color, e.g. the shiny shells of certain beetles [10][11][12] and eye glow of some butterflies [13][14][15]. The striking colors that can arise from ACPCs are shown in the case of gold and silver beetles in Fig. 1. An ACPC in the eye of a deep sea crustacean is shown in Fig. 2, where the chirped structure is signaled by gradual variations in the layer thickness.
Replication of natural photonic structures has already yielded marketable benefit in the manufacture of paints, * calebqcook@gmail.com † arielamir@seas.harvard.edu fabrics, displays, and other technologies [5,17]. ACPCs in particular are well-suited for applications requiring complex, broadband reflection spectra, such as light manipulation in white LEDs [17], dispersion control of ultrashort laser pulses [18], and tailored optical behavior in fiber gratings [19]. Better understanding the design of biological photonic structures such as ACPCs may therefore provide inspiration for future technologies.
There are many well-known tools for calculating the reflectance of dielectric multilayers, including transfer matrix methods and recursion techniques [2,9,21]. Although numerically exact, these methods are iterative, and when there are many layers involved, as in most biological applications, carrying out these computations does not solve the inverse problem of designing reflectors with tailored optical behavior. In this Letter, we instead derive an approximate, but closed-form, expression for the reflectance of ACPCs by matching Bloch wave expressions for the electric field across stopband regions via a discrete Wentzel-Kramers-Brillouin (WKB) method. Analysis of Top: TEM image of an ACPC in the eye of the crustacean C. refulgens. Stained regions show layers, each on the order of 100 nm thick, of differing refractive index. Image taken from Ref. [20]. Bottom: Our model ACPC. Adiabatic chirp requires that the sequence of optical thicknesses m varies slowly. The forward-and backward-propagating fields E ± m are also shown, with an absorbing layer at the back enforcing this closed-form expression gives insight into several of the underlying the design principles of ACPCs.
Matuschek et. al. [18] apply a related WKB approximation to study ACPCs for ultrafast laser application. However, their work makes several simplifying assumptions, such as ignoring reflections from the back of the multilayer and thereby assuming near perfect reflection.
By instead performing asymptotic matching in full, we obtain a reflectance expression that gives better agreement with numerically exact results. Our method also follows a more elementary approach and does not invoke a needlessly cumbersome coupled-mode formalism [22].
Consider a stack of 2N dielectric slabs of refractive indices n L < n H arranged in an alternating fashion. We assume that within the m-th {n L , n H } bilayer, the slabs have fixed optical thickness ratio Ξ = ( L ) m / ( H ) m . The sequence of optical thicknesses m ≡ ( H ) m = ( L ) m /Ξ then encodes the chirp of the multilayer, which we assume to be monotonic and slowly-varying (compared to the wavelength of light) with bilayer index m. We also assume zero back-reflectance at the end of the ACPC, which models a final inhomogeneous absorptive layer, common to biological multilayers, whose light-diffusing character serves to strongly limit the amount of coherent radiation directed back into the structure [11,12]. A depiction of the model is shown in Fig. 2.
In the transfer matrix method [21], one decomposes the electric field at each bilayer interface m into forward-and backward-propagating components The fields E m at bilayer interfaces are related to one another via a transfer matrix product For incident light of wavelength λ, the m-th bilayer transfer matrix F m (λ) = F (x) of our model is [23] F where x = x m (λ) = (1 + Ξ) m /(λ/2), and ρ = (n H − n L )/(n H + n L ) is the Fresnel reflection coefficient, which we assume to be small as is the case in nearly all biological ACPCs (see Table I).
When x is approximately integer, the transfer matrix recursion (1) has evanescent electric field solutions, corresponding to high reflectance of incident light. Wavelengths λ and bilayers m for which x m (λ) is approximately integer therefore define the stopband of the ACPC. Conversely, when x lies between integers Eq. (1) has oscillatory solutions, corresponding to low reflectance and hence defining the passbands. The bilayer indices m t (λ) which separate stopband/passband regions define the turning points. To avoid reflections between multiple stopbands and simplify our analysis, we assume just a single stopband m ∈ [m 1 (λ) , m 2 (λ)] for all considered wavelengths λ. We also assume an x ≈ 1 stopband, since these require minimal layer thickness and are therefore most prevalent in biological systems.
In the case of no chirp m = , the vector recursion (1) has exact, Bloch wave solutions E ± m = exp [−m ln µ ± ] w ± where w ± m = w ± are the transfer matrix eigenvectors with corresponding eigenvalues µ ± m = µ ± . A natural generalization of these solutions are the modulated Bloch points. In particular, adiabatic chirp implies that near a given bilayer m, the fields E ± m approximately satisfy the recurrence which is formally a discrete Schrödinger equation in a potential determined by Tr(F m ) (see Section II of the Supplementary Material for details). We therefore apply a discrete WKB method [25,26] to approximately solve Eq. (4) near the turning points m 1,2 and match the mBloch wave (3) coefficients C ± across the ACPC stopband. We sketch this procedure in Fig. 3. The ratio of fields E − 0 /E + 0 in the mBloch solution (3) is fixed by enforcing the absorber condition E − N = 0 and applying the discrete WKB matching formulae, which are characterized by the tunneling exponent is the chirp rate at the center of stopband. In approximating Eq. (5), we have used the quadratic expansion of Tr(F m ) shown in Fig. 3.
The reflectance R = |E − 0 /E + 0 | 2 of ACPCs calculated via this discrete WKB matching procedure is then given to lowest order in ρ by analogous to the form found in quantum scattering [27]. This calculation is detailed in Section I of the Supplementary Material, where the next-order WKB reflectance is also given in Eq. (S18). A comparison between the numerically exact reflection spectrum of an ACPC and our corresponding WKB results (6) and (S18) is shown in Fig. 4. As Fig. 4 demonstrates, the WKB expression R(λ) (6) in fact captures a moving average of the exact reflectance spectrum. Dimensional analysis shows that q(λ) in Eq. (5), when written in terms of the wavelengths λ min = 2 (1 + Ξ) 0 and λ max = 2 (1 + Ξ) N over which an ACPC reflects strongly, is independent of the optical thickness ratio Ξ. The tunneling exponent γ(λ) (5) and hence the average reflectance R(λ) (6) is therefore maximal for Ξ = 1, i.e. when the optical thicknesses are equal in each bilayer. For periodic m = reflectors, the above fact is well-known, and such multilayers have been dubbed "ideal" [4]. Our results therefore generalize the claim that ideal bilayers elicit maximal reflectance to the case of periodicity-breaking adiabatic chirp.
Our results can also be used to solve the following design  (8). The percent difference between N and N γ∼1 exp/lin is given in the relevant column. All wavelengths λmin, max given in nanometers. The optical thickness ratio Ξ in each case is estimated via analysis of TEM images given in the corresponding reference (see Section IV of the Supplementary Material). FIG. 5. An ACPC with a designed average reflectance spectrum. The desired spectrum (red) is put into the chirp differential equation (7) with nH = 1.54 and nL = 1.34, which is solved numerically to give the chirp function m plotted in the inset (green). The resulting exact reflectance (blue) of the computed chirp m is then found numerically via transfer matrices and plotted against the input spectrum.
problem: given a desired average reflectance spectrum R(λ), what chirp function m is required to elicit R(λ)? Solving Eqs. (5) and (6) for d /dm reveals that this question can be answered by simply solving the differential equation The chirp differential equation (7) offers a novel method of designing multilayer reflectors of tailored optical behavior, and we exhibit such a design spectrum in 5. Moreover, Eq. (7) shows that achieving a flat spectrum R(λ) = R requires exponential chirp of the form m = 0 ( N / 0 ) m/N . Constant reflectance and hence exponential chirp are of special interest since reflectance above some threshold value is minimally achieved with a flat spectrum.
By computing q(λ) in Eq. (5) for specified chirp functions, we can in fact calculate the minimal number of bilayers N required to exceed some reflectance R over a given range of wavelengths. For comparison, we carry out this calculation for both exponential m = 0 ( N / 0 ) m/N and linear m = 0 + ( N − 0 )(m/N ) chirps, which gives the minimal bilayer numbers (8) is minimized in either case for Ξ = 1 and (ii) N exp ≤ N lin , consistent with our earlier conclusions. We also emphasize that the design principles contained in Eqs. (7) and (8) are not inferable via traditional, numerical (e.g. transfer matrix) approaches.
Comparison between the number of bilayers N found in several biological ACPCs and our predicted minimal number N exp/lin (8) gives a quantitative measure of the optimality of ACPCs in nature. In order to make a concrete comparison, we fix an arbitrary high reflectance benchmark of γ ∼ 1, i.e. R ∼ 86%. This comparison is shown in Table I for several ACPCs in nature and reveals that biological ACPCs are generally less optimal than exponentially chirped reflectors but generally more optimal than their linearly chirped counterparts.
Exceptions to the above conclusion are the butterfly tapeta highlighted in Table I, which have many more bilayers than expected per our analysis. This is consistent, however, with the fact that optical components such as tapeta likely require near perfect reflection, corresponding to γ > 1 and a greater number of bilayers in Eq. (8). Indeed, this explanation is supported by the fact that the reported reflection spectra of D. plexippus [15] and V. cardui [30] tapeta are essentially unity R ≈ 1 over the range of wavelengths [λ min , λ max ].
In summary, our closed-form results allow for tailoring ACPCs of desired reflectance spectra and estimating the optimality of chirped reflectors in nature. More generally, the mathematical methods developed herein are applicable to coupled wave propagation in slowly-varying media and therefore also relevant to the design of mechanical and acoustic metamaterials.
We thank Dvir Gur, Pete Vukusic, Gary Bernard, and Mark Arildsen for helpful discussions. C.Q.C. acknowledges support from the Harvard College Research Pro-gram and the Perimeter Institute. Research at PI is supported by the Government of Canada through Indus-try Canada and by the Province of Ontario through the Ministry of Research and Innovation. A.A. thanks the support of the Alfred P. Sloan Foundation.

I. DERIVATION OF WKB REFLECTANCE
In this section, we outline the calculation of the WKB expression for the reflectance of adiabatic chirped photonic crystals presented in the main text. We also give the WKB reflectance calculated to the next order in ρ.
As presented in the main text, the wavelengthdependent bilayer transfer matrix F m (λ), relating the forward-and backward-propagating electric fields via the vector recursion relation for our model takes the form The eigenvectors w ± m and corresponding eigenvalues µ ± m of the general transfer matrix F m are explicitly where the D ± m are eigenvector normalizations, which we will discuss in more detail shortly.
In the main text it is shown that linear combinations of modulated Bloch waves approximately solve the transfer matrix recursion (S2) in the limit of adiabatic chirp. One could proceed by seeking higher-order corrections to the modulated Bloch wave solution (S6), but we instead invoke energy conservation to find the m-dependence of the eigenvector normalizations D ± m and thereby extract the dominant effect of these corrections [31]. From Eqs. (S5) and (S6) we calculate the Poynting energy flux density which a priori defpends on the bilayer index m. However, taking the eigenvector normalizations D ± m to scale as removes the m-dependence of P m = P (S7) and therefore enforces constant energy flux density in the multilayer, as required by energy conservation. The modulated Bloch wave solution (S6) with the eigenvector normalizations (S8) is a generalization of the Bremmer method for deriv-ing WKB-type solutions [32,33]. As a turning point m → m 1,2 in the multilayer is approached, the transfer matrix eigenvalues (S5) µ + m → µ − m converge to one another and the required eigenvector normalizations D ± m (S8) become arbitrarily large. This implies that our assumption of slowly varying transfer matrix eigenvectors w ± m breaks down and our approximate solution (S6) no longer holds near the turning points m ∼ m 1,2 , even in the case of adiabatic chirp. We must therefore write general solutions of the form (S6) separately in the three regions of the multilayer separated by the turning points m 1 and m 2 , giving where we have introduced the complex coefficients A, B, C, D, F , and G, in addition to the trace a m (λ) ≡ Tr [F m (λ)] whose magnitude determines the complex structure of the transfer matrix eigenvalues µ ± m , so that Since our modulated Bloch wave solution (S9) breaks down near turning points m ∼ m 1,2 , we return to the exact vector recurrence (S2), which may be approximated in the limit of adiabatic chirp as in the vicinity of a given bilayer m. This implies that we can treat E ± m as decoupled fields each independently satisfying the recurrence relation The discrete WKB method has traditionally found application in approximate calculations of quantum mechanical energy spectra [26,31,34,37]. The discrete WKB method in the present context is detailed in Section III, with discrete WKB solutions to the difference equation (S12) given by Eq. (S36).
Crucially, the discrete WKB solutions (S36) share the same functional form as the modulated Bloch waves (S9) in the vicinity of a turning point m ∼ m 1,2 , since in that limit the transfer matrix eigenvector normalizations (S8) become and Eq. (S9) subsequently reduces to the discrete WKB solutions (S36). We therefore identify the complex coefficients A, . . . , G in the modulated Bloch wave solution (S9) with the analogous coefficients in the discrete WKB solution (S36). Applying an asymptotic matching procedure to the discrete WKB solution coefficients, and hence the modulated Bloch wave coefficients A, B, F , and G, then allows for patching of the modulated Bloch waves (S9) across the stopband region. Asymptotic matching of the discrete WKB solutions is performed in Section III, and the resulting connection formulae across an x m ≈ 1 (or more generally, a m < −2) stopband m ∈ [m 1 , m 2 ] are given by where we have introduced the variables The parameter γ is analogous to the tunneling exponent in the Gamow theory of alpha decay [27], in which the probability T of He 2+ emission from an atomic nucleus satisfies T ∼ e −2γ . In the present case, γ(λ) (S17) characterizes how deeply light of wavelength λ penetrates an a m < −2 stopband of a chirped photonic crystal. The terminating absorber condition E − N = 0 assumed in our model gives a linear relation between the coefficients F and G in the modulated Bloch wave solution (S9) and thus a linear relation between the coefficients A and B via the discrete WKB matching formulae (S14). This linear relation fixes A/B and hence the field ratio r = E − 0 /E + 0 through the modulated Bloch wave solution (S9). The reflectance R = |r| 2 of ACPCs calculated in this way is given by where we have again introduced the zeroth-order WKB reflectance (S19) as well as the phase angles

II. TIGHT-BINDING ANALOGY
In this section, we make an analogy between the electric fields E ± m in an adiabatically chirped photonic crystal and the zero energy eigenmodes of a one-dimensional chain, described quantum mechanically by a tight-binding Hamiltonian. Naively applying the continuous WKB method to the quantum chain yields results consistent with those of the discrete WKB analysis near turning points and offers insight into the functional form of E ± m and our matching formulae. By taking a continuum limit of the quantum chain tight-binding Hamiltonian, we also make concrete the analogy between the reflection of light from a photonic stopband and quantum scattering from a potential barrier.
We note that the recurrence relation (S12) is precisely the equation one would obtain by seeking the zero energy eigenmodes E m of a Hamiltonian H given by We can view the matrix H as a tight-binding Hamiltonian of a finite, one-dimensional quantum chain with nearestneighbor hopping elements [38]. Let us first begin with the simple case that a m = a is constant. In this case, we can write the exact eigenmode solution as a discrete plane wave E m = e ±iκm with wavevector κ. Substitution into our recurrence relation (S12) then tells us that which implies that our wavevector κ is given by (S24) The general eigenmode solution E m in the constant a m = a case is then just a linear superposition of the discrete plane waves of wavevectors ±κ. When a ∈ [−2, 2] and κ is real, the zero energy eigenmode components E m are oscillatory in the bilayer/chain site index m. Conversely, when |a| > 2 and κ is imaginary, the components E m are the sum of growing and decaying exponentials in the index m.
In the more general case that a m varies slowly with bilayer/chain site index m, we can write the continuous WKB solution to the above eigenmode problem as [27] Substituting this expression into our recursion relation gives, to lowest order in the variation of κ m with respect to m, the same consistency equation as before which implies (variable) wavevector solutions (S28) Again, the magnitude of a m tells us whether the the eigenmode components E m are oscillatory or exponential in the index m.
Comparing the continuum WKB result (S26) with the form of E m obtained via the discrete WKB method (S36), we see that the two differ only by the m-dependence of the passband magnitude |E m | of the electric field/zero energy eigenmode. In particular, in passband regions |a m | < 2. However, these two forms are consistent to lowest order in (a m ± 2), i.e. near turning points, which is precisely the regime in which the difference equation (S12) is applied in Section III to match the modulated Bloch waves (S9). This explains the strong similarity (exact agreement) between the standard WKB connection formulae [27] and the discrete WKB matching formulae across an odd (even) stopband and given in Section III. We can gain further intuition from this tight-binding analogy by rewriting the tight-binding Hamiltonian H in the continuum limit as where m → z represents the continuous coordinate parametrizing the depth in the multilayer and/or length along the quantum chain, and we have introduced the potential Examining the continuum quantum chain potential V (z) reveals why one obtains evanescent solutions when a (z) > 2; in this case the potential (S31) is positive, and therefore the zero energy modes we seek are classically forbidden by energy conservation. In quantum mechanics, however, such classically forbidden solutions to HE (z) = 0 are allowed, but are exponentially damped in z due to quantum tunneling [27].
An identical analogy can be made for odd a (z) < −2 stopbands by transforming the recurrence (S12) in E m to the recurrence in the transformed variable E m = (−1) m E m , where have defined a(z) = −a(z). Repeating the same analysis as before, we arrive at the odd stopband potential which is plotted in the main text. We therefore conclude that stopbands in adiabatically chirped photonic crystals are analogous to the classically forbidden regions on a one-dimensional tight-binding quantum chain, with light scattering from the former corresponding to quantum tunneling through the latter.

III. DERIVATION OF DISCRETE WKB CONNECTION FORMULAE
In this section, we first derive the discrete WKB connection formulae across an even, i.e. a m > 2, stopband and then apply a transformation to the recurrence (S12) to obtain the analogous connection formulae across an odd, a m < −2, stopband. These results connect scattering states E ± m across a classically forbidden region and augment the discrete WKB matching procedure used by Braun [34] to derive the quantization rules for bound states.
Since both E + m and E − m (approximately) satisfy the same recursion relation near a turning point m ∼ m 1,2 , it suffices to consider just E + m . Approximate discrete WKB expressions for the field E + m can be obtained by assuming a solution to the above recurrence relation of the form where the unknown functions Φ j m are labeled according to their rate of change with respect to the index m, so thaṫ Φ j m = O m −j ,Φ j m = O m −j−1 , and so on. Substituting the discrete WKB ansatz (S35) into the recurrence relation (S34) and solving for the unknown functions Φ j m up to order O m −1 yields [26,34]: Our goal is to match these solutions across the stopband, i.e. find a linear relationship between the coefficients A, B and F , G. Performing this matching procedure will require analysis of approximate solutions to original recurrence (S12) in the vicinity of each turning turning point, which we obtain by linearizing a m near m 1 and m 2 .

A. Even Stopband Connection Formulae
Suppose we have an a m > 2 stopband for m ∈ [m 1 , m 2 ] (where m 1 , m 2 are not necessarily integers), surrounded on either side by propagation bands |a m | < 2.

Matching Across Even m = m1 Turning Point
When linearized near m 1 , the recurrence (S12) takes the form for some α 1 > 0 determined by a m , or equivalently where we have defined ν ≡ x + (m − m 1 ) and x ≡ 2/α 1 . In terms of the variables ν and x, the solution to this recurrence is a linear combination of the Bessel functions E + m = aJ ν (x) + bY ν (x). Now, the asymptotic expansions of the Bessel functions J ν (x) and Y ν (x) for large argument x (consistent with small α 1 and hence slow chirp) and large order ν (consistent with m far from m 1 ) depend on the relative magnitude of x and ν. In each case, the asymptotic expansions are given by (see [39], 9.3.2 and 9.3.3): We can then express these Bessel function expansions in terms of the original variables m, m 1 , and α by noting ν 2 − x 2 = a 2 m − 4 /α 2 1 in the linearized region and applying the change of variable k → k + m 1 − x to obtain Combining these results, the Bessel function asymptotic expansions can be rewritten as so that the linear combination of Bessel functions E + m = aJ ν (x) + bY ν (x) therefore becomes  Using γ, we can then rewrite the discrete WKB solution (S36) as (S48) This reformulation is useful for matching near m = m 2 , since the discrete WKB solution in the stopband is now expressed as a function of right turning point m 2 . The recurrence (S12) linearized near m = m 2 takes the form or equivalently, where we have defined ν ≡ x−(m − m 2 ) and x ≡ 2/α 2 . In terms of the variables ν and x, the solution to this recurrence is a linear combination of the Bessel functions E + m = aJ ν (x) + bY ν (x), just as in the m = m 1 case. In particular, the mapping from the m = m 1 case is exact under the switch m ↔ m 1 (so as to obtain the modified m-dependence of ν) and subsequent replacements m 1 , α 1 → m 2 , α 2 . Applying this transformation to the linear combination of asymptotic Bessel functions (S45) gives These are the matching relations across the m = m 2 turning point.

Matching Across Entire Even Stopband
Eliminating the coefficients C and D in the matching relations (S46) and (S52) finally gives These are the discrete WKB connection formulae across an entire a m > 2 stopband and are identical to the connection formulae across a potential barrier in the standard, continuous WKB approximation [27].

B. Odd Stopband Connection Formulae
Now consider an a m < −2 stopband for m ∈ [m 1 , m 2 ] (where m 1 , m 2 are not necessarily integers), surrounded on either side by propagation bands |a m | < 2. We can transform the solution E + m in this odd stopband region, given again by Eq. (S36), to a solution in the region near an even stopband by taking which leaves the recursion relation (S12) unchanged. This transformation can be achieved by multiplying (S36) by (−1) m and using the identities cos −1 (a m /2) = π − cos −1 ( a m /2) and (−1) m e ±imπ = 1. Effecting this transformation yields where we have defined These are the discrete WKB connection formulae across an entire odd stopband and are precisely the relations contained in Eq. (S14). Interestingly, odd stopbands have no analogy in the continuous WKB approximation, and the odd stopband connection formulae (S14) are not identical to the even stopband/continuous WKB connection formulae (S53). In particular, the presence of additional turning point phases exp (±im 1,2 π) in the odd stopband connection formulae (S57) has important consequences in quantum mechanical applications [26].

IV. TEM IMAGE ANALYSIS
In the main text, we compare the bilayer number N found in several ACPCs in nature to our predicted minimal number of bilayers N γ∼1 exp/lin in the case of exponential and linear chirp. This comparison requires an estimate of the optical thickness ratio which we estimate for each biological ACPC via analysis of TEM images given in the relevant reference. In this section we detail the image analysis process, a sketch of which is shown in Fig. S1, using the elytra of the golden Chrysina aurigans beetle [28] as an example.
First, an anisotropic Gaussian filter is applied to the TEM image along the assumed axis of the ACPC, which smooths the image while highlighting the low-and highindex layers. Then a curvature flow filter is applied in order further smooth the image while preserving the edges between low-and high-index layers. Finally we binarize the image, where the binarization thresholds are determined locally in the image in order isolate the the lowand high-index layers.
After binarization, we take a thin (few pixel high) horizontal cross-section of the image, from which the layer thicknesses d L/H m and hence local optical thickness ratios Ξ m (S59) can be determined. The value of Ξ m for each bilayer m in the C. aurigans TEM image is plotted in Fig. S2. The fixed ratio Ξ is finally estimated by taking the average optical thickness ratio Ξ = µ Ξ over all the ACPC bilayers m.