Machine learning classification for field distributions of photonic modes

Machine learning techniques can reveal hidden structures in large amounts of data and have the potential to replace analytical scientific methods. Electromagnetic simulations of photonic nanostructures often produce data in significant amounts, particularly when three-dimensional field distributions are calculated. An optimisation task, aiming at increased light yield from emitters interacting with photonic nanostructures, enforces systematic analysis of these data. Here we present a method that combines finite element simulations and clustering for the identification of photonic modes with large local field energies and specific spatial properties. For illustration, we use an experimental–numerical data set of quantum dot fluorescence on a photonic crystal surface. The application of Gaussian mixture model-based clustering allows to reduce the electric field distributions to a minimal subset of prototypes and the identification of characteristic spatial mode profiles. The presented clustering method potentially enables systematic optimisation of nanostructures for biosensing, bioimaging, and photon upconversion applications. Machine learning techniques are increasingly expanding their capabilities of making predictions on data across a variety of fields. The authors present a machine learning based approach capable of classifying the three-dimensional spatial electromagnetic field distributions of photonic crystals.

Machine learning is a rapidly developing discipline, which uses statistical approaches to learn from data without explicitly rule-based programming.Driven by today's massive increase in data amounts, the related techniques are extended and improved at a fast pace 1 .
Machine learning is currently applied to all aspects of science, from health sciences and psychology 2,3 , to biology [4][5][6] , to environmental 7 and material sciences 8,9 .But also to matters of everyday-life, from online security, to finance and insurance.While supervised learning has led to breakthroughs in computer vision 10 and speech recognition 11 , unsupervised learning is expected to become far more important in the future 12 .The latter techniques, such as clustering [13][14][15] , allow for the recognition of patterns in unlabeled data and can therefore reveal a hidden structure.They have been successfully applied to e.g.anomaly detection 16,17 or genetics 18 .
In the field of nanophotonics, increasing computer power, storage space and data throughput, as well as improvements in modeling techniques, greatly accelerated all-numerical system design.For nanostructures the following typical optimization tasks are met: • Simple design: Scalar parameters → scalar output (e.g.lengths/refractive indices → reflectivity) • Inverse design: Multivariate parameters → scalar output (e.g.permittivity distribution → reflectivity) • Qualitative design: Scalar parameters → multivariate output (e.g.lengths/refractive indices → 3D field distribution) Simple design tasks are generally solved by simulating the system for many different parameter combinations (i.e.grid search), or by applying function minimization routines.
More sophisticated techniques such as the reduced basis method 19 for finite element method (FEM) simulations have successfully been applied to speed-up this optimization process for large parameter spaces.Inverse design tasks introduce a high-dimensional input parameter space, typically by allowing for arbitrary changes in the permittivity distribution (r) of the nanostructure.Machine learning techniques have successfully been applied for this purpose, mainly using genetic algorithms [20][21][22][23][24] .Simple and inverse tasks have in common that they possess a scalar measure of success, i.e. they can be seen as minimization problems.The machine learning approach in inverse design therefore belongs to the field of supervised learning (more specifically regression).The third design task introduced above substantially differs in the way that the system should be optimized for a high-dimensional output.Due to the inaccessibility of a scalar success metric, we denote this problem as qualitative design.This is for example the case if the 3-dimensional spatial distribution of the electromagnetic fields has to be taken into account.Usually, such problems are solved by appropriate visualizations.But since any change in the input parameters leads to a change in the high-dimensional output, the data amounts quickly become extremely large.We will demonstrate below that machine learning, or more specifically clustering, is able to overcome these issues by reducing the output dimensionality.
As indicated before, an example of qualitative design is to optimize a photonic nanostructure, e.g. a photonic crystal (PhC), for an appropriate spatial field distribution.This is of high relevance whenever an interaction of the field with a (potentially vague) particle distribution is present, e.g. for emitters on nanophotonic surfaces or emitters embedded into the nanostructure.PhC slabs exhibit a phenomenon called leaky modes: resonances that can be excited using external radiation [25][26][27][28][29] .Leaky modes have been used to improve various applications (e.g.light trapping in photovoltaic devices [30][31][32][33][34] , light-emitting diodes 35,36 ), but can also affect near-surface emitters, such as QDs, atoms or molecules.Especially in the life sciences, the applications range from PhC enhanced microscopy and single molecule detection to enhanced live cell imaging, DNA sequencing and gene expression analysis [37][38][39][40] .
In this study we present a powerful technique based on machine learning for the classification of 3D electromagnetic field distribution data.This method can be of avail in any case where large amounts of electromagnetic field (or energy) distribution data should be reduced to a minimal subset of typical distributions.We will refer to these as distribution prototypes.We directly apply the technique to a specific dataset of our previous publication on fluorescence enhancement of lead sulfide (PbS) quantum dots (QDs) on a silicon PhC slab surface 52 , however, without loss of generality.A similar setup was used in previous studies 19,53,54 .
The effect is sketched in Fig. 1(a), depicting emitters (black dots) that interact with a leaky mode of the PhC excited by an external laser source.The illumination conditions introduce four parameters: the laser wavelength λ, the laser polarization P (TE or TM), the polar angle θ with the plane normal, and the azimuthal angle φ used to define the high symmetry direction (Γ − M or Γ − K).The latter is also indicated in the scanning electron microscope image of the sample (without emitters) in Fig. 1(b).An example of an electric field distributions E(r) of a leaky mode is depicted in Fig. 1(a).As mentioned, the energy density of the electric field of the leaky modes, w lm (r), can be larger compared to the energy density of the incident plane wave, w pw , known as field energy enhancement (w lm (r)/w pw > 1).To study this effect in large parameter spaces we usually define the volume-integrated field energy enhancement where V sup is the volume of interest.In our case V sup is the superspace of the computational domain, as indicated by the yellow dashed line in Fig. 1(a).The energy density of the plane wave has no spatial dependence and is proportional to the amplitude of the electric field, E pw,0 , and the refractive index n of the surrounding medium, i.e.
In the figure, a uniform random distribution of emitters is shown as an example.But depending on the coating process, emitters might have a very specific spatial distribution in a real application, e.g. a monolayer attached to the surface, or a higher concentration inside the holes, or at the plateaus between the holes.Consequently, the spatial distribution of the energy density w lm (r) becomes a determining factor and, therefore, the integrated field energy enhancement E + is not sufficient to quantify the effect on the emitters.
An optimized design for an application as sketched in Fig. 1(a) can hence be achieved by identifying a mode which has (i) a large volume-integrated field energy enhancement E + and (ii) an appropriate spatial field energy density distribution overlapping the locations of the emitters, at the same time.Task (i) is a "simple design" task, as defined in the introduction, while (ii) is a "qualitative task", i.e. an optimization of a multivariate output.
The results of the fluorescence enhancement measurements of our previous study 52 are shown again in the upper row of Fig. 2, as well as the results for the volume-integrated field energy enhancement E + in the center row.The latter results solve task (i), as described above.Task (ii) potentially enforces to take into account 3D field distribution data of all combinations of the illumination condition parameters λ, P, θ and φ.If the number of considered wavelengths N λ and the number of angles N θ becomes large, it is no longer feasible to directly visualize all the 3-dimensional field distributions for all points in the λ-θ-maps shown in Fig. 2. It is hence necessary to reduce the amount of field distribution data in an appropriate way.One possibility to achieve this reduction is to pitch on specific wavelengths and incident angles for which the field distribution is evaluated, as it was done in the previous study 52 .This way, however, information is mainly gained at random, so that general trends might be overseen.
A more systematical approach is to cluster field distributions which are similar, and to therefore derive typical distributions (i.e."distribution fingerprints").It is known that a certain undisturbed photonic band in the leaky mode regime will not significantly change its symmetry properties when crossing the λ-θ space 29,55 , as will be explained in more detail below.As a result, the entirety of field distributions are composed of a finite set of patterns which are basically caused by the finite number of bands.This feature space can efficiently be partitioned into the typical patterns using machine learning clustering techniques.
In the following, we will first reconsider the experimental and numerical results of the previous study 52 , highlighting aspects which were left unexplained by the prior analysis technique.Afterwards, we introduce the clustering technique and apply it to systematically analyze the 3D field energy distribution properties.The distributions are classified by as-signing them to distribution prototypes, which are consulted as representative solutions to fully explain the effects observed in the experiment.We further consider a mathematical method based on silhouette coefficients 56 to assess the clustering result.Based on these analyses we will explain how the method enables to solve complex optimization tasks with high-dimensional output, as indicated in the introduction.The all-numerical technique is of relevance for the design of nanophotonic structures for any application in which emitters interact with the electromagnetic field, e.g.highly-sensitive biosensors 37 , quantum dot solar cells [57][58][59][60] , or up-conversion devices 61,62 .

RESULTS
The upper two rows of Fig. 2 1)).The E + maps exhibit clearly visible bands of strong field energy enhancement, which partly correspond to regions of high measured fluorescence F + .Some deviations are caused by a Q-factor mismatch between the spectral bandwidths of the leaky modes and the excitation laser source.However, a few features of the measured F + maps remained unexplained, for example: The declining band after the anticrossing point, which is visible in the corresponding E + -map, but missing in the experimental fluorescence enhancement F + .
Please consult Ref. 52 for further details of the comparison.

Introduction and justification of the clustering technique
The E + -maps given in the center row of Fig. 2 only provide information about the volumeintegrated field enhancement over a characteristic volume V sup , marked by the yellow dashed line in Fig. 1(a).Therefore, regions of high E + can be regarded as a necessary condition for fluorescence enhancement, but not as a sufficient one.A high E + without a corresponding fluorescence enhancement F + hence indicates a lack in the spatial overlap of the emitters with the regions of enhanced field energy density.
However, it is known that bands of the photonic crystal have well-behaved spatial properties when varying the k-vector between two high-symmetry points of the irreducible Brillouin zone 29,55 .More specifically, the modes belong to the same symmetry point group as the system seen from the point in k-space, i.e. they exhibit the same spatial symmetry.
Consequently, it is theoretically justified to expect that the spatial properties of the bands only change smoothly with θ.For a fixed high-symmetry direction, e.g.Γ − K, we only expect two types of solutions, which are 1.regions that correspond to leaky-mode bands, and 2. regions that are off any photonic band, and therefore corresponding to the continuum of radiation modes.
The regions of the radiation modes are expected to exhibit solutions that resemble plane waves, i.e. show oscillatory behavior in the the exterior domain.All things considered, only a small number of different spatial symmetry types is expected, which is of the order of the number of bands that cross the parameter scan window.
This is where machine learning comes into play.If we consider the specific 3D electric field distribution for a single illumination setting as a sample, and the electric field values of each point in the considered volume as features, then clustering techniques are able to subdivide the entirety of field distributions into a finite number of field distribution prototypes.This In each of the silhouette plots (i.e.columns), silhouette coefficients for each sample are plotted as a bar in x-direction with a length corresponding to its value.The samples are sorted by their silhouette coefficients, with smaller values being located at smaller y-positions; and grouped and color-coded using the same colors as in Fig. 2. Red dashed lines mark the silhouette scores.
approach is reasonable because the data range is expected to contain a finite number of typical field patterns, and each of the real field patterns can be identified with one of those prototypes.Moreover, these prototypes have a sufficient "uniqueness", e.g. they considerably differ in their symmetry properties.The methods section "Clustering of electric field data" gives a detailed description of how the clustering is performed.In a nutshell, for each illumination condition, i.e. a set of (P, φ, θ, λ), the electric field strength E = (E x , E y , E z ) is derived from a FEM simulation.It is sufficient to export the fields on symmetry planes to reduce the data volume, for which we use the xy, xz and yz planes marked in Fig. 1(c).The validity of this approach was tested using a comparison to full 3D exports using a smaller dataset.Note that, in contrast to the the volume-integrated field energy enhancement Sckit-learn 63 .From the clustering itself two characteristics can directly be gained: the classification, which labels each observation with a cluster index i, and the distribution prototypes, usually denoted as "cluster centers" in the general clustering literature.
The latter are the average of the electric field distributions (on the chosen planes) of all samples that belong to a specific cluster i.We note that we averaged the normalized input data, i.e. the exact data used for the clustering, to calculate the prototypes.The prototypes therefore represent the actual mathematical cluster centers, with the tradeoff that the absolute field amplitude information is lost, because the samples are normalized individually.Another possibility would be to average the unnormalized fields, so that the amplitude information would be conserved, with the tradeoff that the prototypes derived that way are not exactly the cluster centers.We settled for the normalized fields, as the amplitude information is essentially included in the E + maps.
As in most clustering techniques, the number of clusters must be specified in GMM clustering, so that the appropriateness of this choice has to be validated.This aspect will be covered shortly.

Classification maps
The classification can be visualized by assigning each point (θ j , λ k ) to a different color that corresponds to its label i.The procedure of determining the number of clusters will be explained in the next subsection.
When comparing the classification maps to the E + maps above, a striking accordance can be observed.The narrow bands of high field enhancement in the E + maps correspond to narrow areas at the same positions in the classification maps.Note that the E + maps and the classification maps are based on very different data sets: the former are derived from a spatial integration over the electric field energy density distribution w lm (r) in the superspace volume V sub only (Eq. ( 1)), while the latter uses electric field patterns E(r) on planes that include the PhC and glass domains.When observing the regions off the leaky-mode bands, i.e. the domains of the radiation modes, it is seen that these regions are multiply subdivided in some cases; e.g.Γ − K, TM: bottom left.In contrast, other parts are homogenous over large ranges, such as Γ − M , TE: top right.
Another detail of these plots are the different levels of saturation used for each point, obtained by alpha-blending with a black background.This additional layer of information illustrates the representation quality of the local solution by the assigned cluster, as determined using so-called silhouette coefficients 56 .The silhouette coefficients provide a way to assess the initial choice of the number of clusters, and how well the samples lie in their respective clusters, at the same time.The silhouette coefficient rates how well a sample fits into its own cluster.If it is far away from all other clusters and very close to the cluster center (i.e.prototype), the sample gets a positive rating.If the distances to a different cluster and its own cluster are comparable, it is rated with values close to zero.Finally, if it is much closer to a different cluster, a negative rating is assigned.See the methods section "Solution quality rating using silhouette coefficients" for a severe definition.
In all cases, we observe that the saturation decreases at the border of two clusters.This is expected, as silhouette scores close to 0 indicate a sample which is in fact close to the border of the neighboring cluster.It is apparent from this phenomenon, and important to stress here, that the clustering technique is a tool.The field distribution data is not categorical, so we expect superposed solutions which are badly represented by "pure" modes.That said, these intermediate parts are small, as it is seen from the saturation distributions, so that the clustering is still a valid and effective approach.

Silhouette analysis and the number of clusters
Before we investigate the field distribution prototypes, the quality of the clustering itself is evaluated using a mathematical analysis in the following.This can be done using the silhouette coefficients, using a scheme known as silhouette analysis 56 .Figure 3 depicts socalled silhouette plots for each combination using the same column order as in Fig. 2. In each of the plots, the silhouette coefficients for each sample are plotted as a bar in x-direction with a length corresponding to its value (negative values point into the −x-direction).The samples are sorted by their silhouette coefficients, with smaller values being located at smaller y-positions.In addition, the samples are grouped for each cluster k and color-coded using the same colors as in Fig. 2 (lower row).The red dashed lines mark the average of all silhouette coefficients, which is a measure for the absolute quality of the representation denoted as silhouette score.The results are the typical "sails" or "shark fins".The width of each fin in the silhouette plots is proportional to the area of the correspondingly labelled points in the classification maps.
Considering the distribution of the silhouette coefficients, fins which are not too sharp are observed, i.e. having broad plateaus of high silhouette coefficients.There is only a minimum number of values with negative coefficients.Both arguments together give a validation for the fact that the number of clusters is not underestimated: negative values would occur if there were to few clusters, leaving back samples which do not fit in one of the classes (s ∼ −1).Too many clusters could be identified by a large fluctuation in the fin widths.
But this does not fully apply here, as the areas occupied by the bands and the residual parts are unequal.Therefore, equally broad classes are not expected.A slightly too large number of clusters can be seen as unproblematic, because it would basically subdivide the radiation mode regions further, which are of limited relevance for the interpretation.Another point that suggests a good representation is that there are few clusters with below average silhouette scores.Using this reasoning, the optimum number of clusters was determined for each case by comparing the silhouette plots for different number of clusters (results omitted).

Field distribution prototypes
As the second essential outcome, the clustering procedure yields the field distribution prototypes.As the input data for the GMM algorithm has been electric field values on three planes, namely xy, xz and yz, the prototype data is available on these planes as well.
The prototypes for all clusters of each combination, and on all three planes, are depicted in Fig. have strongly localized energy distributions (see panel 1, column 6 in Fig. 4).In contrast, clusters that belong to radiation modes have energy distributions that increase away from the PhC, e.g.cluster 2 of the Γ − K, TE case (dark blue, see panel 2, column 2).Further comparisons will be considered below.

Putting the pieces together
To explain the measured fluorescence enhancement effects shown in the upper maps of Fig. 2, it is necessary to combine all information gained from the numerical analysis.This is, the volume-integrated field energy enhancement maps (E + , Fig. 2, center row), the classification maps (Fig. 2, lower row) and the prototype maps (Fig. 4).A guide on how the different aspects of the results can be connected to yield a complete interpretation may read as follows: 1. Select a feature in the volume-integrated field energy enhancement (E + ) maps.For these features the simulation suggests a possible excitation enhancement effect.
2. Check whether there is an according feature in the experimental fluorescence enhancement (F + ) maps.
3. Afterwards, observe the corresponding region in the classification maps and determine the cluster label from the color using the color bar.
4. Using this label or color, locate the related column in the prototype map that belongs to the direction/polarization combination (the prototype maps are ordered according to the columns in Fig. 2, from left to right).Check if the field energy distributions on the three planes can explain the observed fluorescence enhancement (this may necessitate to take into account all cases, because the QD distribution is unknown).
For ease of comprehension, we will analyze the results for selected cases in order of increasing As the QDs in this experiment are distributed inside the holes and particularly in an about 100 nm -300 nm thick film on top of the structure, they overlap with the leaky mode volume very well.When observing other columns of the prototype map, there are other interesting patterns which could potentially increase the emission of the QDs, for instance in columns 1 and 5.These two modes gather their energy in the center of the hole and at the flanks, respectively.However, when looking at the classification maps again, these correspond to the red and yellow regions, which are outside the measurement window.Another important point is given by the patterns of solutions that are related to radiation modes.These would be expected at regions off any band, e.g. in the dark blue and cyan regions with labels 2 and 7. Returning to the prototype maps, these modes in fact have energy distributions that increase with the distance from the PhC surface, but less dominant in case of the the prototype of cluster 2.Moreover, these modes do not fulfill the necessary condition 2 of the guide given above, i.e. they do not exhibit an integrated field energy enhancement (E + ).
Γ − K, TM.We observe a steep band of high fluorescence and a large field energy enhancement (condition 2) at a similar position.The clustering approach found the band as well, labeled as 6 (pinkish).Column 6 of the prototype map shows that this band localizes its energy at the center of the hole (slightly above the PhC in z-direction) but also at the From the prototype map it is observed that this band has a node along the x-direction and concentrates its energy at the flanks and the plateaus in y-direction.The energy distribution is therefore comparable to the one of the orange band in the Γ−M , TM case, which is clearly seen in the experimentally measured fluorescence as well.The shallower mentioned band undergoes a transition from the green cluster (label 3) to the red cluster (label 1) in the classification maps.For the green parts the energy is strongly localized at the flanks, while

Mode B "�lank mode"
Mode C "hole mode" FIG. 5. Full-3D volume renderings of selected modes for the Γ−K, TE case.Semi-artistic ray tracing images depicting multiple periods of the photonic crystal as a grayish material.The upper row shows a topview of the full-3D E-field energy density, color-coded using a heat map (not comparable between the figures).The lower row shows a closer view and indicates the same random distribution of quantum dots (bright small spheres), emitting white light with an intensity proportional to the field energy density at their specific positions.The columns relate to three different modes of the Γ − K/TE case, denoted as modes A, B and C, as marked in Fig. 4. The modes are the actual solutions from the finite-element solver that have the smallest deviations from the assigned prototype (i.e.cluster center), determined using the silhouette coefficients.Incident angle θ and wavelength λ for each mode are given in the headings.The figures use real physical proportions.
the red one can be identified with a radiation mode.The energy is therefore less well confined to the surface in the red case, which is exactly seen in the fluorescence maps, where only at the location of a broad green region a fluorescence enhancement can be observed, but no enhancement is seen at the prosecution of the mode at higher incident angles.Γ − M , TE.The E + -maps cover two bands that show anticrossing at the long wavelength/steep angle end of the data window (top left).The lower band and most parts of the interaction zone are labeled as cluster 3 (green), while the upper band has label 7 (cyan).
The green band has a node in the yz plane, while the energy in x-direction is strongly localized at the flanks.It is therefore likely to be seen in the fluorescence enhancement, when comparing to the previous results.However, a fluorescence enhancement for the green band is mainly seen in the interaction zone.This can be explained when considering the energy distribution of the cyan band: this band does not have the node in the yz plane.In the interaction zone, the two bands basically overlap, so that the node is partially erased.Therefore, a stronger effect in the fluorescence enhancement is expected, just as it is observed.Interestingly, there is moreover a small effect at ≈ 1080 nm for large angles in the F + map.The clustering reveals another band with strong localization at the flanks: the orange region with label 4 (and the similar yellow one which couples more strongly to the radiation modes).It is likely that the measured effect is actually caused by this band, as it has a similar energy distribution as other bands which are clearly seen.
To give a clear idea of the 3D energy density distribution for three selected modes, and also to show how well the clusters match the actual physical fields, Fig. 5 shows full-3D renderings 64 .The images depict multiple periods of the photonic crystal as a grayish metallike material, without showing a superspace material.The upper row shows a topview of the volume-rendered electric field energy density color-coded using a heat map, which is not comparable between the figures.The lower row shows a closer view and indicates a random distribution of QDs as bright small spheres, emitting white light with an intensity proportional to the field energy density at their specific positions.The QD distribution is the same for all three images.The columns relate to three different modes of the Γ−K, TE case, denoted as A, B and C.They correspond to clusters 8, 3 and 6, as also marked in Fig. 4.
The modes are the actual solutions from the finite-element solver that have the smallest deviations from the assigned prototype (i.e.cluster center), determined using the silhouette coefficients.Incident angle θ and wavelength λ for each mode are given in the headings.
Note that these images have an illustrative character, but can be very helpful to imagine the actual physical situation.Modes A and B are the ones which have been discussed recently, and it is clearly seen that the former concentrates its energy at the plateaus, while the latter has high energy densities at the flanks of the holes.A third type is shown with mode C, which focusses the energy directly inside the holes.The illustrations in the lower row give a notion of how these modes activate different QDs, depending on their position.Only a small density of QDs is used for the images for purposes of visibility, and they are randomly distributed in a layer that fills the holes and extents 100 nm in z-direction.Mode A very efficiently excites QDs at the plateaus, just as expected, while modes B and C do the same at the flanks and inside the holes, respectively.Consequently, these renderings completely confirm the results of the clustering approach.

DISCUSSION
The aim of the numerical approach presented here is the systematic identification of suitable leaky modes of nanophotonic structures for interaction with near-surface emitters.For instance, (a) a monolayer of emitting species attached at the surface of the nanophotonic structure is expected to strongly interact with a leaky mode with shallow field distribution.
In contrast, an experiment with (b) emitters in a coating on top of the photonic nanostructure, the interaction with a shallow leaky mode will be rather small due to the limited spatial overlap of the mode volume with the emitting material.Here, a leaky mode with an energy density enhancement in a large volume outside the photonic nanostructure would be better suited.In a third scenario, (c), emitters fill the voids of the photonic nanostructure, for example if they are solved in a liquid solution and dropped onto the structure.In that case, leaky modes with strong field enhancement inside these voids are expected to cause the strongest effects.In the chosen dataset, the fluorescence enhancement experiment of PbS quantum dots on a silicon PhC slab in nanohole geometry, the distribution of QDs resembles a mixture of case (b) and (c).The clustering technique revealed that the modes which have the best spatial overlap with the QD distribution effectively cause the strongest fluorescence enhancement effects in the measurements.
In the previous study 52 we used a selection of a small number of points for which the field energy distributions was analyzed.The clustering technique confirmed the results that were achieved this way, but it also helps to explains complicated details, e.g. as caused by the superposition of two modes.Therefore, the clustering approach gave a much more coherent and detailed explication of the underlying physical phenomena.It emphasizes the interesting parts automatically and systematically, e.g. by revealing regions of rapidly changing field distributions through individual clusters, or through large deviations from the assigned prototype.Moreover, the clustering technique seems to be applicable to even more complicated cases, e.g. in windows with more bands for which an analysis using selected points is not reasonable any more.
The presented technique composed of (i) the field energy enhancement maps and (ii) the 3D electric field distribution clustering provides a versatile tool for the analysis and design of photonic nanostructures for applications that utilize near-field enhancement effects for increased emission.For any known distribution of near-surface emitters that should be affected by leaky modes, optimum values for all relevant parameters can in principle be determined.It is e.g.possible to define a wavelength range for the excitation of the emitters by considering their absorption properties, and to numerically calculate the field energy enhancement E + and field values in 3D for clustering (as provided in the center and lower row of Fig. 2, here).By choosing the mode with the largest spatial overlap of high field energy with the emitter distribution from the prototypes (as in Fig. 4), an optimum mode can systematically be determined.This process can moreover be repeated for possible geometrical parameters of the photonic nanostructure, e.g. the lattice constant, slab thickness or hole radius.Alternatively, if the geometrical parameters should be varied extensively, the technique could be applied for an initial set of geometrical parameters to select a potential mode and to reduce the wavelength and angle window.Successively, only the field energy enhancement E + may be calculated in the scan over the possible geometrical parameters to determine the absolute maximum of the enhancement.
The clustering technique is extremely flexible.It is not limited to uniformly sampled feature spaces as shown in our example application.It would also have been possible to choose arbitrary snapshot points in the θ-λ space, e.g. with a higher density in regions of high field energy enhancement E + .It is further not limited to the shown number of feature parameters, i.e. we could have added a variation of the hole diameter or other geometrical parameters as well.But the method is even more powerful, because the trained classifier can be used to classify field distributions that it has not "seen" yet, known as prediction.In contrast to the clustering itself, this is a computationally cheap process, and the classifier can even be persistently stored on disk for later use.To make these considerations more clear, it would have been possible to choose a smaller number of possibly non-uniformly sampled points in the θ-λ space for efficient clustering.The silhouette analysis can be used to make sure that the number of samples is sufficient to reach an appropriate clustering result.From this clustering the prototype field distributions can be derived and the classifier can be stored to disk.Afterwards, an e.g.uniform scan over θ, λ, and other parameters that are expected to not change the field distributions considerably, (e.g.hole diameter, slab thickness, refractive indices, . . . ) could be performed.The resulting new solutions could then be assigned to the prototypes using the classifier from disk with minimal computational effort.
Numerous applications could benefit from these optimization abilities.In the field of biosensing, photonic nanostructures have become an important platform for e.g.label-free biosensing or for the enhancement of the output of photon emitting tags used in the life sciences and in vitro diagnostics.A recent review article 37 shows that nanophotonic enhanced biosensors are yet extremely relevant, even commercially and potentially on industrial scale.Exploiting leaky modes with large Q-factors enables for narrow bandwidths (< 1 nm) and extremely high sensitivities, e.g. for detection of disease biomarkers in serum with concentrations of ∼ 1 pg ml −1 .The numerous applications that are described in the mentioned review article have in common that the nanophotonic structure is designed for a very specific mode, i.e. a specific illumination condition and a determinable distribution of the molecules/cells/virus particles in question.This is where the technique presented here could be utilized for a systematic optimization in the design process, and hence to further increase the sensitivities of related sensors.Photon upconversion 65,66 in biomedical imaging and solar energy is another application that could benefit from the discussed all-numerical design abilities.Recent publications 61,62 demonstrate upconversion using thin emitter layers, which as well could potentially be improved using specifically tailored nanophotonic structures.
In summary, we have developed a numerical method that allows to systematically optimize nanophotonic structures pertaining to the 3D field distribution and field energy enhancement of modes.The method uses a combination of FEM simulations and postprocessing using machine learning clustering.We showcased the modelling power of the method by explaining experimentally measured fluorescence enhancement of QDs on a photonic crystal slabs surface.The method yielded information that was not easily accessible using e.g. a visualization-based analysis for selected parameter combinations, and which allowed to fully explain the experimental results.Consequently, the presented technique could be of great avail for applications that utilize effects that depend on the spatial field distribution of nanophotonic modes, such as in the fields of biosensing 37,66 , quantum dot solar cells [57][58][59][60] , or up-conversion in solar energy 61,62,65 .

Clustering of electric field data
The clustering is executed on an input matrix X of shape N s × N f , where N s is the number of samples and N f the number of features.A sample is the solution for a specific set of input parameters, in our case incident angle θ and wavelength λ.The features, in the present case, are absolute values of the electric field components E j with j ∈ {x, y, z} for a number of points r i ∈ R 3 , i.e. of the form |E j (r i )|.Consequently, if the field is evaluated at N p points, these are N f = 3N p features.To avoid exporting the electric field on a full Cartesian grid in 3D, which would cause huge amounts of data when trying to achieve a reasonable resolution, data is only exported on the symmetry planes marked in Fig. 1(c), respectively.More symmetry planes could be used as well, but based on these three planes a reasonable classification can be reached, as tested using smaller data sets and comparing to a full 3D field output.A field pattern of a single simulation holds data for each of the 3 spatial directions, and for each component j of the electric field (altogether a 4D data set).As each sample X i must be a 1D row vector with observations of single scalar values x 0 , . . ., x N f −1 , it is necessary to flatten these data sets in always the same way, yielding "1D representations" of the fields.The data is moreover normalized by scaling each sample to unit norm individually.The field export is performed for each point in each map of Fig. 2, center row, so that the samples are unique simulations for a given direction/polarization combination, wavelength λ and incident angle θ.The number of samples for a single map is given by N s = N λ • N θ .To give an expression for the complete input matrix X we abbreviate , where the additional indices m = 0 . . .N θ and l = 0 . . .N λ have been introduced, and where the hat denotes the absolute value and normalization.The input matrix then reads For the wavelength and angle resolution values of 0.5 nm and 0.3 • have been used, respectively.For each clustering procedure the input matrix X had a size of N s × N f = 47 034 × 8616.This is a comparably large problem size, especially because the large feature dimensionality (N f ), so that the procedure took more than 10 hours on a hexa-core workstation with roughly 40 GB of memory consumption.

Gaussian mixture model clustering
Simple clustering techniques, such as the k-means algorithm 67 , can be extremely robust, but also have their disadvantages.E.g. k-means assumes that the clusters are circular, i.e. representable by a (hyper-)sphere in feature space.The center of this sphere defines the cluster center (i.e.prototype), while the radius acts as a hard boundary used to decide which samples belong to the cluster.In contrast, the GMM 67,68 is a so-called soft method.That is, a "score" for each cluster is assigned to the samples, which account for the probability that the sample belongs to a specific cluster.In GMM clustering, the clusters are represented by Gaussian distributions of the dimensionality of the features space (i.e.N f ).
In general, a superposition of N multivariate Gaussian distributions of the form can be used to approximate almost any continuous density to arbitrary accuracy (this is intuitive with 1D Gaussians, which can fit almost any 1D signal if enough Gaussians are superimposed).Here, the N i (x) are multivariate Gaussian distributions of the form 67 for a D-dimensional vector x, the D-dimensional mean-vector µ, and the D × D covariance matrix Σ with determinant | Σ|.Equation ( 4) is called a Gaussian mixture, the N i (x) are called components of the mixture, and the c i are weight factors.Loosely speaking, the distribution of sample points is "fitted" using a set of high-dimensional Gaussians.A GMM can therefore represent much more complex data sets and can be seen as a generalization of the k-means algorithm for non-circular clusters.One can imagine that it would be straight-forward to fit the multivariate Gaussians to a data set for which the labels are known.With unlabeled data the case is more difficult, and enforces to take into account another step.In the literature, this problem is commonly denoted as to find out which (latent) component is "responsible" for a certain sample, -which is somehow a different way of asking to which cluster the sample belongs.But it underlines that the GMM clustering is a probabilistic approach, because it calculates the probability that the sample was generated by cluster i for all clusters.These probabilities, which are also called responsibilities, are simply the weight factors c i of Eq. ( 4).In the implementation that was utilized here, the cluster assignment is solved using a method known as expectation-maximization 69,70 .This algorithm starts with a random Gaussian mixture (i.e.random components), which is typically initialized using a prior application of k-means to improve the convergence.In the next step it determines for each sample the probability of being generated by each component of the mixture.Based on these probabilities, the parameters of the Gaussian distributions are fitted to give the best approximation of the data by maximizing their likelihood 67 .This process is executed iteratively and is guaranteed to converge to a local optimum.

Solution quality rating using silhouette coefficients
To give a definition of the silhouette coefficient, let X k i be a sample that was assigned to the cluster k and a(i) be the average dissimilarity of X k i to all other members X k j =i of this cluster.The measure for the dissimilarity is usually the Euclidian distance.Let d(i, m) be the average dissimilarity of X k i to all members of the cluster m = k and b(i) be the minimum of d(i, m) for these clusters, i.e.The cluster m for which this minimum is obtained is called the neighboring cluster of X i .If the number of clusters N k is > 1, we can define the silhouette coefficient s(i) for the sample

b(i) = min
From this definition it is seen that the silhouette coefficient s is in the range −1 ≤ s ≤

FIG. 1 .
FIG. 1. Overview of the nanostructure.(a) Light incident on a silicon photonic crystal (PhC, gray) on glass (cyan) excites a leaky mode that exhibits enhanced electromagnetic near-field energies in the superspace volume (marked by the yellow dashed line).Emitters (black dots) in the vicinity of the PhC surface interact with the local electric field distribution.(b) Scanning electron microscopy image of the PhC sample with denoted high-symmetry directions Γ − K and Γ − M .(c) A unit cell of the PhC system as used in the simulation.Yellow, green and red rectangles mark the planes used for the field export.

FIG. 2 .
FIG. 2. Comparison of measured quantum dot fluorescence enhancement F + , simulated volume-integrated field energy enhancement E + , and corresponding classification maps.(Upper row) Measured fluorescence enhancement F + as a function of vacuum wavelength and incident angle θ of the laser source (logarithmic color scale; see Fig. 1(a) which indicates θ; see supplementary material of Ref. 52 for experimental setup).The columns correspond to the four combinations of sample orientation (Γ − M and Γ − K) and source polarization (TE and TM).(Center row) Simulated volume-integrated electric field energy enhancement E + for the same conditions as in the upper row.For a definition of the volume V sup see Fig. 1(a).The white lines mark the experimental data limits.(Lower row) Classification maps depicting the cluster assignments (labels) using different colors independently for each plot, and the respective silhouette coefficients using alpha-blending with a black background (color bar omitted; see Fig. 3 for value ranges).More saturated colors denote larger silhouette coefficients.(Note: Upper and center rows of the panel repeat the same results as already shown in Ref. 52 for a larger angle and wavelength range.) repeat the main findings of the prior fluorescence enhancement study 52 .The upper row shows the fluorescence enhancement (F + ) maps obtained by tilting the QD-coated PhC sample along the respective high-symmetry directions of the irreducible Brillouin zone (Γ−M or Γ−K, adjusted using φ), and by using transverse-electric (TE) or transverse-magnetic (TM) polarization P of the incident laser radiation.Each measured spectrum (for a single incident angle) was first integrated over the fluorescence peak from λ = 1200 nm to λ = 1700 nm and normalized to the measured incident laser power and the absorption profile of the QDs, yielding the fluorescence F .A minimum estimate for the fluorescence enhancement F + is obtained from dividing by the minimal value in each of the maps.The maps feature regions of enhanced fluorescence.The measured fluorescence enhancement is caused by increased energy densities of the fields at the emitter positions.Concerning task (i) of the introduction, the center row of Fig. 2 maps the electric field energy enhancement E + integrated over the simulated superspace volume V sup , which contains the QDs (see Eq. (

FIG. 3 .
FIG. 3. Silhouette analysis plots for the different direction/polarization combinations.In each of the silhouette plots (i.e.columns), silhouette coefficients for each sample are plotted as a bar in x-direction with a length corresponding to its value.The samples are sorted by their silhouette coefficients, with smaller values being located at smaller y-positions; and grouped and color-coded using the same colors as in Fig.2.Red dashed lines mark the silhouette scores.
. 1(a)), the fields for the clustering are also considered in the dielectric materials (silicon PhC and glass substrate, Fig. 1(c)).To account for the different cluster sizes (narrow bands) and unknown cluster shapes in the data set, the flexible Gaussian mixture model (GMM) clustering technique is used (see the methods section "Gaussian mixture model clustering" for details), implemented in the Python library Recall that the clustering is carried out individually for each combination of polarization P and azimuthal angle φ (=high symmetry direction Γ − M or Γ − K).Plotted in the same fashion as the E + maps of the center row of Fig. 2, we denote the resulting figures as classification maps.These classification maps are shown in the lower row of the figure.The color scale relates the colors to the labels and, hence, identify the corresponding cluster.Note that the classification maps cannot be compared among each other, although the same colors have been used.The clusterings for the Γ − K-cases used 8 clusters, while the Γ − M -cases only required 7 (i.e.there is no gray region in these maps).

4 .
For each direction/polarization combination a prototype map is shown in one of the four horizontal panels.Each panel consists of three rows for the xz, yz and xy planes, respectively (from top to bottom).Each of these rows has the same number of columns that accounts for the number of clusters, and each column has a colored edge in the top-most row that corresponds to the color used for that label in the classification maps shown in the lower row of Fig.2.The cluster label is further given in the title of the xz-row.Each distribution plot depicts the electric field energy distribution E 2 in the respective plane.The distribution plots further feature semitransparent markings for the glass superstrate (blue) and the silicon of the PhC (gray) in the case of xz and yz; and a white circle indicating the hole circumference in the case of xy.Recall that the color scales do not give absolute values, as the prototypes are based on normalized data and, therefore, cannot be compared with respect to their absolute amplitudes.For each prototype, the field energy plots on the three planes give a notion of the 3D field energy distribution.The solutions with the same label (color) in the classification maps of Fig.2all share this distribution type.Lower saturations quantify how much the individual solutions deviate from the prototype.Clusters that correspond to leaky mode bands with strong field enhancement, such as cluster 6 of the Γ − K, TM case (pinkish),

FIG. 4 .
FIG. 4. Prototype maps for the different direction/polarization combinations.For each direction/polarization combination (4 horizontal panels) a prototype (i.e.cluster center) map is shown, which consists of 3 rows for the xz, yz and xy planes, respectively (from top to bottom).Each of these rows has a number of columns that accounts for the number of clusters, and each column has a colored edge in the top-most row that corresponds to the color used for that label in the classification maps shown in the lower row of Fig. 2. The cluster label is further given in the title of the xz-row.Each distribution plot depicts the electric field energy distribution E 2 in the respective plane.The distribution plots feature semitransparent markings for the glass superstrate (blue) and the silicon of the PhC (gray) in the case of xz and yz; and a white circle indicating the hole circumference in the case of xy.(Color scales do not give absolute values, as the prototypes are based on normalized data.)

Γ
− M , TM.The experimental fluorescence enhancement (F + ) map features a single stripe of increased fluorescence with a high contrast.This stripe excellently corresponds to the single leaky-mode band causing a high volume-integrated field enhancement in the E + map.The classification map reveals this band accordingly with label 4 (orange), for which the field distributions are shown in column 4 of the prototype map.The xz and yz patterns show that the energy of this band is accumulated at the plateaus between the holes.
plateaus in the xz plane.It therefore potentially affects QDs relatively independently of whether they are gathered inside the hole or on the plateaus.Again, as the emitters are distributed inside the hole and in a bulk film on top of the structure in the experiment, the overlap with the leaky mode field distribution is good.Returning to the clustering, it is seen that a second, much broader band running from top left to bottom right is seen in the E + maps.The classification maps reveal that the field distribution of this band undergoes a change when crossing the pinkish band, from label 7 (cyan) to label 1 (red).This band is basically not seen in the fluorescence enhancement.Γ−K, TE.The TE cases both feature a more complicated band structure.In the Γ−K, TE case, there are two very clear bands that show anticrossing, a steeper band crossing the complete wavelength range from roughly 20 • to 40 • , and a shallower one coming from top left.The former is very clearly seen in the clustering by the gray region with label 8.

1 .
Values near 1 indicate that the sample is far away from the neighboring cluster and accordingly fits well into its own cluster.A value of 0 indicates that the sample is on or very close to the boundary between its own and the neighboring cluster, and negative values indicate that it might have been assigned to the wrong cluster.A sorted diagram of all silhouette coefficients can thus be used to visualize the representation quality of a clustering.In addition, the average silhouette coefficient for all samples -usually denoted as silhouette score -can be used to compare the representation quality for different clusterings, e.g. using different N k -values.It hence even provides a single numeric value for solution quality assessment.