Complex photonic response reveals three-dimensional self-organization of structural coloured bacterial colonies

Vivid colours found in living organisms are often the result of scattering from hierarchical nanostructures, where the interplay between order and disorder in their packing defines visual appearance. In the case of Flavobacterium IR1, the complex arrangement of the cells in polycrystalline three-dimensional lattices is found to be a distinctive fingerprint of colony organization. By combining analytical analysis of the angle-resolved scattering response of in vivo bacterial colonies with numerical modelling, we show that we can assess the inter-cell distance and cell diameter with a resolution below 10 nm, far better than what can be achieved with conventional electron microscopy, suffering from preparation artefacts. Retrieving the role of disorder at different length scales from the salient features in the scattering response enables a precise understanding of the structural organization of the bacteria.


Introduction
Iridescent structural colours originate from the interference of light scattered from transparent materials organized at a nanoscale level. In contrast to dye/ pigmentation-based coloration originating from absorption, such structural colours allow living organisms to fully control their visual appearance both in terms of colour and angular distribution. In fact, periodic arrangements of nanostructures in one, two or three dimensions can cause very bright, angledependent structural colours [1], as expected for photonic crystals [2][3][4][5][6][7], while matt, homogeneous colour arises from short-range correlation of scatterers [8] or resonant scattering phenomena in disordered structures [9][10][11].
In the case of periodic structures, diffraction has been used to assess bacterial dimensions as long as 100 years ago [12][13][14] and structural colour in bacteria was noted in the literature as long ago as 1941 [15]. Recently, a number of flavobacterial strains have been reported to grow in colonies that display structural colour [16][17][18]. However, a detailed explanation on how the interplay between order and disorder in their packing dominates the optical appearance in this system remain unclear. Such knowledge would allow insights into bacterial colony behaviour and enables their use as living optical materials.
In this manuscript, we used the Flavobacterium strain IR1 (wild-type, WT) [18] as a model system to study how the 3D polycrystalline organization of rodshaped bacteria creates structural colour (see figure 1a; electronic supplementary material, figure S1). Through analytical and numerical models, we assessed the lattice length of the packing with a resolution below 10 nm and estimate the value of the refractive index of their environment providing insights into the material used in their extracellular matrix. Moreover, we were able to extract information about the crystalline organization of the cells, their orientation, as well as the packing fraction by correlating all observed spectral features to physical mechanisms. Inducing deviations from the perfect crystal lattice arrangement of the bacterial positions in numerical simulations and the analytical analysis allows us to learn about the degree of disorder present on multiple length scales in the bacterial samples.
The results obtained enable further studies of the selforganizational capacity of bacterial colonies, as a thorough understanding of the structure-function relation of bacterial colonies opens new methods of modification and use of these bacteria for a range of technological applications.

Bacterial colonies acting as photonic systems
Flavobacterium IR1 was cultivated on Artificial Sea Water Black (ASWB, see Methods) agar in a Petri dish. As shown in figure 1a, a strong angle-dependent colour appearance is observed. For small illumination angles with respect to surface normal (θ i = 0°, 10°) the colony appears blue with a very low scattering intensity, while for larger angles (−45°, −60°) we can observe a strong green scattering response (see Methods for details and figure 1b for a schematic of the measurement geometry). In the growing edge of the colony, the colour is red-shifted, however, in the following we will consider only the centre of the colony as the growing edge continues to organize and reorganize as well as similar green appearances have been reported in earlier studies [16,[18][19][20].
Structural organization of bacteria. An overview of the structural organization of IR1 bacterial colonies is obtained with a combination of electron microscopy (EM) techniques (scanning EM (SEM), cryogenic-SEM and transmission EM(TEM)) as shown in figure 1c-f . The SEM image in figure  1c shows an edge of the colony allowing insights into bacterial organization. The positions of the bacteria are strongly correlated in space (figure 1d) and this correlation can be observed across all the colony (see also electronic supplementary material in [18]). Hexagonal packing is clearly recognizable in cryo-SEM cross-section imaging (figure 1e, see Methods for details) and TEM (figure 1f , see Methods for details). However, large-area imaging both via cryo-SEM images (electronic supplementary material, figure S2a and b) and TEM (electronic supplementary material, figure S2d) reveals that the bacteria are aligned in a hexagonal lattice locally, but that there is no long-range correlation in their orientation. In fact, each domain shows local order on a length scale which is comparable to the coherence length of incidence light. Therefore, when measured on a larger scale (approx. 100 μm [18]) all possible in-(surface)plane rotations of the lattice can be expected. This observation is in agreement with the structural colour appearance that shows no angular dependency with in-plane rotation (electronic supplementary material, figure S4 in [18]) and allows us to model the bacterial colony as domains of packed bacteria with local order and random orientation on a larger scale as visualized in electronic supplementary material, figure S2c.
However, while the EM imaging techniques provide a very good qualitative description of the organization of the cells, the precise measurement of the value of the lattice constant d (inter bacteria distance) from these images remains challenging due to sample preparation artefacts. In fact, TEM image analysis of the structure in figure 1f reveals a lattice constant (inter bacteria distance) of d = 179 ± 17 nm, which is less than half of what we measure from the cryo-SEM cross-section (figure 1e, d = 414 ± 29 nm). Such values differ also from the cryo-SEM top view images where the lattice constant is d = 374 ± 36, nm if measured by a line plot (see electronic supplementary material, figure S3) and is d = 396 ± 36 when estimated by autocorrelation image analysis as described by Johansen et al. [18], where the authors obtained a value of 357 nm from a SEM top view image. Interestingly, while the absolute values of d differ between the different methods, the ratio between the bacterial size a and lattice distance d remains similar for all methods within errors (a/d = 0.9 ± 0.05), revealing that the artefacts are mainly connected to shrinking caused by dehydration of the colony. Cryo-SEM imaging suffers less from dehydration artefacts but quantitative values of d are difficult to achieve deep inside the colony, as ice crystal growth becomes an increasing issue away from the sample surface [21]. The overall sample thickness is estimated from cryo-SEM cross sections (electronic supplementary material, figure S2e) to be 10-15 μm. Moreover, as the cell orientation varies within the colony, the measured d is often not measured in a cross-section image that is perpendicular to the orientation of the cells. This can be seen in electronic supplementary material, figure S2d, where the same sample is imaged at different locations in cross-section leading to different absolute lattice values and sizes. On the other hand, when the structure factor of fixed colonies is evaluated from TEM and SEM images (see electronic supplementary material, figure S4), good hexagonal packing is observed. Therefore, this technique provides a very good overview of the ordering of the sample especially deep inside the colony, where the cryo-EM methods fail. Goniometry as an advanced characterization tool. In order to evaluate the value of d without perturbing the colony we recorded angle-resolved reflectance spectra with an optical goniometer set-up. This technique has the advantage of probing a large area of living colonies without the need for any sample preparation. Figure 2a shows the schematics of the geometry used in our experiment for a scattering configuration where the illumination is fixed at one angle and the scattering from the sample is measured at all the other angles (see Methods for details of the set-up). Similarly, a specularreflection configuration can be measured (figure 3a), where the angle of illumination always equals the angle of observation. The specular reflection (θ in = −θ out ) is measured in the range θ in = −45°to 45°, while the scattering properties of the sample are collected for the light incidence angles θ i = 0°, −45°and −60°with respect to the sample normal. The scattering goniometry data shown in figure 2c-e together with the structural characterization by EM, will be the foundation for the subsequent analysis.

Analytical analysis allows structure determination
Figure 2c-e shows three scattering goniometer measurements for the light incidence angles of θ i = 0°, −45°and −60°of a sample cultured for 2 days in an incubator at 25°C, measured in the centre of the colony. For the investigation, a piece of agar (approx. 2 × 2 cm) with the colony was simply cut out from the Petri dish, attached to a glass slide using doublesided tape and then mounted on the goniometer set-up. In all measurements, a region of approximately 5°around the angle of incidence shows no signal caused by the limitation in the design of the set-up (the detector arm overlaps the source). In the following logic, we aim to uncover the structural origin of all observed features by assigning them to their physical mechanisms. Small changes in the optical appearance can be detected by goniometry [18] and if related to their structural change, this will reveal changes in bacterial organization.
As expected from the EM image analysis, we observe a scattering response that can be attributed to a polycrystalline two-dimensional (2D) structure with hexagonal packing in cross-section. In our configuration, only the scattering from bacterial domains aligned perpendicular to the goniometer detector plane will contribute to the recorded scattered data (see sketch in electronic supplementary material, figure S2c). In fact, the measured angle-resolved reflection for any domain and any illumination is consistent with conical diffraction (see figs 1-3 in [22]). In more detail, in figure 2c-e, we can recognize several distinct features in the optical response of the bacterial colony. In the following, we provide a phenomenological, but quantitative description of these features derived from closed-form expressions and use this to extract the relevant parameters determining the optical response of the colony. (c)  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200196 Lattice constant retrieval from the diffraction spot. Diffraction grating-like structures are formed by the bacterial colonies (see for example the cryo-SEM top view figure 1d). The lattice constant, corresponding to the inter-bacterial distance, can be directly obtained from the experimental goniometer data (figure 2). The angles of constructive interference from a diffraction grating can be expressed by the grating equation where m ∈ [0, − 1, + 1, − 2, + 2, …] is the order, λ is the wavelength of light, d is the period of the structure, θ i is the angle of incidence and θ m is the reflection angle for a given order. By matching this equation with the high-intensity reflection spots in the scattering measurements (see white lines in figure 2c-e) we find that the spot location is indeed governed by the grating equation. The equation estimates the periodicity to be 395 ± 5 nm. Note that this simple method allows us to achieve very high resolution combined with high reproducibility. An uncertainty of only ± 5 nm is obtained from fitting three different WT IR1 samples on day 2 as shown in electronic supplementary material, figure S6. For comparison, a variation of Δd = ± 20 nm is shown in electronic supplementary material, figure S6a as dashed white lines. This shows a remarkable degree of reproducibility for an in vivo measurement and confirms that goniometry can be far more accurate than TEM, SEM or cryo-SEM for obtaining averaged quantitative values of the packing. The possibility to achieve such a high accuracy allows us to map efficiently different life stages of the bacteria colony. As an example, a lattice constant of d = 395 ± 5 nm was obtained for day 2 (figure 2e). The same sample is measured at day 1 (figure 2f) where the colony shows a slightly more red shifted diffraction spot, corresponding to a lattice constant of d = 425 ± 5 nm.
Diffraction streak as an indication of domain tilts in the colony. The diffraction spots in figure 2 show a certain width which can be related to a tilt of the sample surface by ± 10°(dashed white lines in figure 2e). This tilt in the orientation of the domains can in fact be observed in figure 1a, electronic supplementary material, figure S3e and S3f on the sample surface. Note that the bacteria still grow mainly as a flat film on the agar plates but with a local deviation from that is expected to occur in bacterial colonies. Similarly, earlier reports show that buckling is a common phenomenon in bacterial colonies [23][24][25] and calculations by Kientz et al. [20] suggest such effect to be present in structurally coloured bacterial colonies.
Remarkably, in the proximity of the diffraction spot, we observe a significant spread of intensity (diffraction streak). This spread stretches over the whole angular range and crosses the specular reflection even for normal incidence (θ in = 0°) where the higher wavelength diffraction spot is out of the measured angular range. Such spread cannot be caused by tilts in the sample surface, variations in the lattice constant of the periodic arrangement of the bacteria or other types of disorder, such as orientational disorder within one crystalline unit cell. We expect this streak to originate from tilted domains within the colony. To confirm that tilted domains within the colony are responsible for the streak line of the diffraction peaks, we extracted the low wavelength scattering features in figure 2c (black line), d (red line) and e (blue line) and corrected them for the angle of incidence of the light. From geometrical considerations, it follows that a shift in input angle corresponds to a similar shift in output angle (electronic supplementary material, figure S5a). In figure 2b, we observe that the extracted scattering behaviour overlays on each other once corrected for the angle of tilt, suggesting that the streaks indeed originate from tilted domains. The observation of the spread of intensity around the diffraction spots over the whole angular range allows us to conclude that a large variety of tilts are present in the sample. Note that the diffraction spots still govern the signal intensity, meaning that most domains grow flat and parallel to the agar surface.
Refractive index retrieval from angle-dependent specular reflection. The optical data allows the extraction of another key parameter for photonic systems, its effective refractive index. The specular reflection peaks for a specific angle θ out are caused by the constructive interference of waves travelling in the medium at an angle given by Snell's Law u 2 ¼ arcsin (sin u in =n avg ), where θ in is the illumination angle and n avg is the volume average effective refractive index of the total material composite in the photonic crystal. Inside the sample, the optical path length changes with the z-projection of the wavevector, meaning that it is the cosine-projection of the distance the wave travels that determines the constructive interference criterion (electronic supplementary material, figure S5b). Therefore, the peak reflection wavelength λ p for normal incidence is shifted in angular space to a shifted peak wavelength λ s as follows: λ p can be extracted from the specular measurement (figure 3b) by reading out the wavelength where the specular line crosses 0°. Then equation (2.2) can be used to fit n avg to the specular reflection curve profile (red dashed lines in figure 3b). The fitting procedure determines the average refractive index to n avg = 1.4 ± 0.05. This is close to the reported values of n bac = 1.38 and n agar = 1.34 for another Flavobacterium grown on an agar plate [20]. For measurements on a dried sample, we obtain n avg ≈ 1.2 ± 0.05 (see blue line in figure 3c), which indeed confirms the drying out of the material compared to the live, hydrated state. We conclude that this method is a quick way to estimate or confirm the refractive index range of an unknown system, and can directly be used to investigate the composition of the extracellular matrix although more precise techniques [26,27] are needed to better resolve the refractive index precisely.

Numerical analysis of a crystalline bacterial arrangement
In the following, we provide a full description of the optical response from the bacterial colony by finite-element calculations. We further present a bandgap analysis of the photonic crystal structure to provide a more intuitive insight on the propagation of the electromagnetic modes inside the bacterial colony. These results uncover the physical origins of the optical response of the studied sample.
In figure 1c, we observed that the bacterial stack in tightly packed layers, forming a hexagonal arrangement in crosssection ( figure 1e,f ). In a first approximation, the optical response can be obtained by using a 1D transfer matrix approach [28,29] (see Methods) with a refractive index structure as pictured in figure 4a (right). In this case, we assumed the refractive index of the bacteria to be n bac = 1.38 and the surrounding matrix n env = 1.34 in agreement with previous reports [20]. The angular dependency of the multilayer system is studied by changing the angle of incidence and measuring the wavelength-resolved reflectance. The result shown in figure 4b reveals that this simplified simulation of the 2D photonic crystal structure already predicts the specular reflection response measured experimentally.
Bandgap analysis confirms 2D hexagonal photonic crystal. To further understand the origin of the specular reflection and diffraction peaks, a photonic bandgap analysis on a hexagonal 2D photonic crystal structure was performed using the MPB code [30] (figure 4b). In this calculation, the input parameters are the lattice constant (d = 395 nm) and the tabulated refractive indices. The result of this analysis shows similar behaviour for both transverse electric (TE) and transverse magnetic (TM) polarization. Therefore, in this section, we will only discuss the results for the TM modes in figure 4b and refer to electronic supplementary material, figure S7b for the TE modes. In figure 4b, note that the dashed lines correspond to asymmetric modes while the solid lines correspond to symmetric modes. The symmetry of the modes was identified in electronic supplementary material, figure S7a. Asymmetric modes can be disregarded when searching for bandgaps as they cannot be excited by an incident plane wave [31].
From the bandgap diagram, three partial band gaps in the visible light range are found manually. The two bandgaps in the ΓK-direction are at 309 nm and 471 nm (light blue background) and correspond to the direction where light is incident perpendicular to the colony surface. The gaps are both less than 1 nm (see figure 4b zoom on the right), owing to the low refractive index contrast, such that mode banding cannot be excluded. Despite their small spectral width these two band gaps can still lead to a noticeable signal in the experiments as they form the origin of the specular reflection lines. It is worth noticing that the visible light reflections in the experiments originate from band gaps between higher bands and not the first band, which would require roughly twice the bacterial cell size. This seems to be different from other natural photonic crystal systems [32][33][34].
Another bandgap is found in the ΓM-direction at 336 nm (magenta background) corresponding to a 30°(or −30°) tilt of the structure. Interestingly, the appearance of one of the band gaps in the UV part of the optical spectrum may provide a functional adaptation of the bacterial colonies in providing a UV protection mechanism. To further understand the optical signal corresponding to this gap we simulated the angular resolved reflectance of a multilayer tilted by 30°and −30°a pproximating the hexagonal structure as sketched in electronic supplementary material, figure S8a. Although the experimental case is different, the observed features can be related to two weak signals observed in the specular measurements seen as crossing lines of the high wavelength specular line (see arrows to rescaled specular measurement in electronic supplementary material, figure S8). In the experiments, we observe that the signal coming from tilted layers with respect to the sample normal (see electronic supplementary material, figure S8a) is weakened in intensity compared to the symmetric specular lines.
Finite-element method reveals photonic bandgap lines. In order to confirm the parameters established from the analytical analysis and prior literature, as well as to check our predictions from the bandgap analysis and multilayer simulations, we simulated the specular and diffracted reflection using the finite-element method (FEM). Combining the FEM with a transfer matrix formalism, we can resolve Maxwell's equations for the full system, including the air interface and the substrate interface. In this model, we only assumed a periodic, hexagonal royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200196 packed bacterial colony with a period of 395 nm as extracted from equation (2.1), and the refractive indices of n bac = 1.38, n agar = 1.34. We use an in-house finite-element solver [35,36] combined with a transfer matrix method [37] for efficient reflection calculations of the periodic layers of bacteria, meaning that only one layer of bacteria were simulated and the solution then repeated N times to obtain the overall scattering characteristics (figure 4c). This result shows excellent agreement with the specular measurements as shown with a rescaled intensity in electronic supplementary material, figure S8c. All specular lines are recovered by the FEM modelled hexagonal packed crystal structure. Given the model assumptions, the prediction of the reflection wavelengths is fairly accurate. The reflection bandwidths are much more confined, which is to be expected for a completely periodic model system. In contrast to all earlier introduced studies, the FEM calculations allow us to study the effect of the photonic band gaps on higher order diffraction. By plotting the first-order diffraction spectra for varying angles of incidence (see electronic supplementary material, figure S8d), we find that the high wavelength diffraction giving rise to the strong green appearance of the bacterial colonies only sets in for incident angles larger than θ in ≈ 20°.

Learning from disorder
In the bandgap analysis as well as the FEM study we assumed, the bacteria were forming a perfect hexagonal crystalline lattice from touching hard disks in cross-section. By contrast, in EM cross-sections (figure 1e,f), we observed that the bacteria resemble soft disks with a certain degree of packing and orientational disorder. To account for such effects we, therefore, perform finite-difference time-domain (FDTD) simulations [38,39].
Perfect lattice reveals crystalline orientation. In all FDTD calculations, we study the integrated total reflection for normal incidence of a plane wave. As shown in figure 4, the direction of normal incidence on the hexagonal structure corresponds to the KΓ direction. We first investigate the effect of the orientation of the hexagonal packing in respect to the incident light. The total reflectance of the perfect hexagonal lattice with lattice constant of d = 395 nm is simulated for two crystal orientations in figure 5a.
For incidence in the KΓ direction, two strong peaks at 336 nm, 473 nm and one weak peak at 313 nm are obtained corresponding to experimentally observed features while the MΓ incidence leads to a spectrum with five peaks in the visible range, in contradiction to the situation observed in the experiments.
In figure 5b, we compare the total integrated reflectance with incidence in the KΓ direction with the spectrum of the specular reflection and find that the peaks at 313 nm and 473 nm correspond to the specular reflection while the peak at 336 nm correspond to the first-order diffraction in excellent agreement with the observation in the 0°incidence measurement in figure 2. Note that in the experimental spectra (figures 2 and 3), the specular reflection peaks show higher intensities than royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200196 the diffraction peaks. This is not fully recovered in the numerical calculations as the boundary of the sample has not been taken into account and the intensities strongly depend on the domain thickness inside the colony assumed in the numerical calculations. Taking into account the sample boundary conditions, the domain thickness inside the bacterial colonies and the absorbing substrate in transmission is a challenging task, both to investigate experimentally and to calculate numerically, but would allow a quantitative comparison of intensities and give further insights into the bacterial assembly.
To conclude we related the features observed in the spectra to a specific orientation of the crystalline structure as pictured in figure 5a (right inset). The orientation found here is expected as it allows a better packing in a 2D geometry, forming a flat surface, consistent with surface energy minimization arguments.
Disorder analysis. To study the effect of disorder on the total integrated reflectance, we introduce lattice disorder (figure 5c) and crystalline orientation disorder (figure 5d) separately (see Methods for details). The average value of the lattice constant is assumed to be the one retrieved from the analytical analysis. We expect the rod-shaped bacteria to have a diameter smaller than this value to allow for small amounts of disorder. This claim is further strengthened by the cryo-SEM cross-section shown in figure 1e,f . In these images, we extracted averaged sizes at least 5% smaller than the distance. Therefore, in all FDTD simulations, we used a particle diameter of 375 nm while the particle distance (lattice constant) is kept at 395 nm.
In figure 5c, we observe that small amounts of disorder in the lattice distance (hardly visible in the simulated structure) already lead to a strong decrease in the peak intensity and a broadening of the peak width, similar to the reported literature for another 2D photonic crystal [40]. A similar effect is also observed for the disorder in angular alignment within the hexagonal crystal as introduced in figure 5d. Here, already small amounts of disorder can be observed in the real structure (electronic supplementary material, figure S9). Complete angular disorder (see red line and corresponding structure as well as isotropic angular space distribution in figure 5d) leads to a broad spectrum with no distinct peaks in the visible range. Despite the low refractive index contrast in our samples, a remarkably bright angle-dependent colour appearance is observed in the measurement. From this, we conclude that a high degree of order within each crystalline domain remains in the studied samples.
The FDTD method can further be applied to real structures imaged by TEM or cryo-SEM in cross-section. A TEM image is converted to a binary image (electronic supplementary material, Fig S4) and the total angular integrated reflectance spectrum is calculated. The result is shown in figure 5 as a black dotted line. A very broad peak around 510 nm is observed strongly deviating from the experimental results. This emphasizes that goniometry in comparison with royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200196 numerical tools allows more insights into the structural arrangement of living microorganisms than EM imaging methods alone.

Conclusion
In summary, we studied the interplay between order and disorder in the structural organization of iridescent bacterial colonies. By correlating angle-resolved measurements with simple closed-form expressions, we were able to deduct the inter-cellular spacing as well as estimate the average refractive index of cells and their matrix. The combination of analytical scattering analysis, multilayer simulations, bandgap analysis and finite-element calculations allowed us to recover the origin of all optical features observed in specular and scattering goniometry measurements. Moreover, the quantitative characterization of the broadening of the diffraction peak led us learn about the degree of disorder on domain alignment and enables further studies of bacteria communication under different growth conditions: as an example it might be useful in estimating the physical interactions of cells during adaptation of the colony structure in response to environmental conditions.
In conclusion, our in vivo, non-invasive method is particularly important for biofilm characterization as bacterial colonies and biofilms are known to be differentiated dynamic structures [41,42] and traditional methods rely often only on imaging studies confined to monolayer microcolonies [43] or static sectioning approaches. Additionally, our finding can be extended to any other 2D photonic tissue in other organisms as the response only depends on the size ratio between the scattering elements and the optical wavelength.

Sample isolation
Strain IR1 was isolated on ASW agar during a screening of estuarine sediment samples from the Neckarhaven region of Rotterdam harbour [18]. Strain IR1 is a yellow-pigmented Gram-negative bacterium culturable on ASW agar under aerobic conditions from 2 to 30°C with salinity matching the location where it was isolated (0.5-1.5). On ASWB plates, containing nigrosine as a contrasting agent, IR1 formed intensely structurally coloured green colonies with red/orange edges as seen in electronic supplementary material, figure S1.

Sample preparation
Bacteria were grown on ASWB agar plates, containing 17.5 g agar (Invitrogen), 10 g KCl, 1 g yeast extract (Sigma: Y1625), 5 g peptone (Sigma: 70 173) and 0.33 g nigrosin (Sigma: N4763) per litre. Plates were allowed to air dry for 30 min after which bacteria were added, and then air dried for another 15 min. Bacteria were inoculated on the plates by suspending them in a 1% w/v KCl solution, of which 5 μl was deposited on the plates. These plates were placed into an incubator at 25°C for 1/2 days before measuring.

Photographic images
A directed beam from a xenon lamp (HPX2000, Ocean Optics) was used as a light source. A white sheet was put in between the beam and the sample to simulate a diffusive, but directed illumination. All images where taken with a camera (Nikon D3200).

Scanning electron microscopy
The SEM studies (figure 1c) were performed on fixed colonies that resulted in the dried material maintaining structural colours described elsewhere (see electronic supplementary material in Johansen et al. [18]).

Cryogenic SEM
Cryo-SEM was performed on an FEI Verios 460 scanning electron microscope with a Quorum cryo-transfer system PP3010T. A piece of agar (2 × 2 cm) with bacteria was cut out and placed in a specimen holder with colloidal graphite suspension. The specimen holder plus sample was then plunge frozen in liquid ethane and transferred to a specimen preparation chamber cooled down to approximately −140°C. Samples for cross-sectional view were then freeze-fractured, all samples were sublimed at −90°C and sputter-coated with platinum. All samples were imaged at 2.00 kV acceleration voltage.

Transmission electron microscopy
Samples were embedded for TEM by chemical fixation. Small pieces of the bacteria colony on agar gel were cut out and fixed with glutaraldehyde (2 wt%) and formaldehyde (2 wt%) in 0.1 M sodium cacodylate buffer (pH 7.5) overnight at 4°C. During the fixation the bacteria colony had detached from the agar which was subsequently removed. Samples were then rinsed with the buffer solution and post-fixed in an OsO 4 solution buffered at pH 7.4 for 2 h at 4°C. The samples were then washed with deionized water and dehydrated through a graded ethanol series (30-100%) to dry acetonitrile. They were left in a mixture of 50% dry acetonitrile and 50% Quetol 651 epoxy resin overnight, followed by infiltration in Quetol resin for 2 days. Samples were placed in a silicon mould with Quetol resin and cured for 2 days at 65°C. Ultrathin sections were cut using an ultramicrotome Leica Ultracut E. with a 45°diamond knife (Diatome). Sections were placed on carboncoated copper grids and post-stained with 2 wt% uranyl acetate aqueous solution and Reynolds lead citrate solution. They were imaged in a FEI Tecnai G2 transmission electron microscope operated at 200 kV and equipped with an AMD CCD camera.

Goniometry
We used a xenon lamp (HPX2000, Ocean Optics) coupled through a fibre to a reflective collimator (RC08SMA-F01) as the light source in order to produce a collimated beam with a spot size of 5 mm. The same type of reflective collimator and fibre was used to couple the light to a spectrometer (AvaSpec-HS2048, Avantes) for detection. For normalization, a white diffuser was measured at 0°incidence with the detector at 5°. For further details, see [44].

Generation of two-dimensional structures for finite-difference time domain simulations
Structures with different degrees of disorder were generated following an inverse-design approach [39]: first, the desired number of hard particles were added to the system, then their positions were gradually changed in order to minimize the difference between the targeted structure factor and that of the ensemble of particles.

Finite difference time-domain simulations
The simulations were performed using LUMERICAL 8.18 (Lumerical Solutions Inc., Vancouver, BC, Canada), a commercial software based on the FDTD method. We used a fixed number of scatterers for a simulation volume of 5 × 5 μm with a filling fraction (ff ) of ff = 0.6. The value of ff royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200196 was decided to guarantee the same number of particles both in the disordered and in the perfect hexagonal packing. Periodic boundaries conditions in the lateral direction, i.e. perpendicular to the incoming beam, Y and perfect matching layer (PML) boundaries in the X direction were used in all the calculations.
A plane incident wave is sent into the sample and the angular integrated reflection (−90°to 90°) is measured. This corresponds to an angular integration of the 0°incidence measurement in scattering goniometry (figure 2c). For each disorder configuration, we averaged over seven independently generated samples. In these calculations, the effect of the air-sample interface was not taken into account.

Finite-element method
An in-house periodic 2D finite-element solver based on [37] and implemented following the recipe in [45] was used. It has been thoroughly tested in several publications [35,36]. In the following, a short explanation of some important aspects for these particular simulations are given, but the reader is referred to [37,45] for full details.
Due to the periodic boundary conditions, the electric and magnetic field at the boundaries can be resolved by infinite sums (as opposed to integrals for the general case) of exponential basis functions multiplied by a set of coefficients. These are referred to as reflection and transmission modes at the boundaries. In the FEM implementation, this sum of modes is truncated to a reasonable number, since the lowest orders are propagation modes and higher orders are evanescent modes, whose coefficients normally tend to zero. Due to our implementation, this means that the non-periodic interfaces of the simulation model have to be of a homogeneous material in order to be described by the basis given in [37,45].
The highest computational cost in solving the FEM system is to decompose the so-called stiffness matrix in the linear FE system. It is, therefore, extremely fast to re-use this decomposition in order to solve for other combinations of transmission and reflection coefficients. By doing this, we can obtain a complete (up to the order truncation) description of transmission and reflection for any possible combination of incident waves (on both sides) [37]. This is represented in a transfer matrix, and such matrices can be combined to represent combined geometries. As a side note, one has to remember to calculate incident angles based on Snell's law when doing this (not mentioned in [37]). Therefore, by dividing the bacteria model as seen in the sketch in figure 4c, we have an extremely fast way of calculating light interaction through just one period of bacteria, and then extending the results to 2N number of repeated bacterial cells in the depth direction. Compared to the bandgap analysis, we therefore include the effects at the interfaces (in particular the top interface, where the plane wave couples to the bacterial photonic crystal) and the effect of a finite height of the bacterial colony. Furthermore, the approach reduces computational time dramatically (compared to a full model) and makes it feasible to use a very fine resolution (Δλ = 2 nm)-and also to include different incident angles. All FEM simulations for this paper were performed on a laptop computer from 2015 with two physical Intel i7 3.1 GHz cores, taking less than 24 h. The approach has a few limitations. Firstly, the scattering matrices may become imprecise or unstable due to large differences in numerical values and since matrix inversions are needed [37]. Secondly, the interfaces have to be homogeneous (so we cannot make the bacterium smaller than the period). In practice, this means that some numeric values may be slightly off compared to full models, but we found the differences to be insignificant. The scattering matrix approach was verified by simulating the reflecting of N repetitions for one wavelength using a full model for low numbers of N, and then compared to the transfer matrix result. A deviation of less than 1% was generally observed.

Photonic bandgap calculations
Bandgap solving was performed using the open-source software MPB [30]. A 2D model was run for both TE and TM modes. The ctm file and post-processing scripts for reproducing the results can be found in the data repository accompanying the paper. Since a plane wave can only couple to symmetric modes [31], we needed to detect symmetric and anti-symmetric modes. This was done manually by plotting the field of the modes and search for symmetries (electronic supplementary material, figure S8).

One-dimensional transfer matrix method
The multilayer simulations were performed by using a 1D transfer matrix method python package [46]. It simulates light propagation in planar multilayer thin films, including the effects of multiple internal reflections and interference.
Data accessibility. All data needed to evaluate the conclusions in the paper are present in the paper and/or the electronic supplementary material and/or available at the University of Cambridge data repository (https://doi.org/10.17863/CAM.46768).