Microstructure design using graphs

Thin films with tailored microstructures are an emerging class of materials with applications such as battery electrodes, organic electronics, and biosensors. Such thin film devices typically exhibit a multi-phase microstructure that is confined, and show large anisotropy. Current approaches to microstructure design focus on optimizing bulk properties, by tuning features that are statistically averaged over a representative volume. Here, we report a tool for morphogenesis posed as a graph-based optimization problem that evolves microstructures recognizing confinement and anisotropy constraints. We illustrate the approach by designing optimized morphologies for photovoltaic applications, and evolve an initial morphology into an optimized morphology exhibiting substantially improved short circuit current (68% improvement over a conventional bulk-heterojunction morphology). We show optimized morphologies across a range of thicknesses exhibiting self-similar behavior. Results suggest that thicker films (250 nm) can be used to harvest more incident energy. Our graph based morphogenesis is broadly applicable to microstructure-sensitive design of batteries, biosensors and related applications. Representing microstructures with graphs, and mapping their properties mathematically, can help improve material properties. Morphogenesis is the process of optimizing performance in devices; representing the candidate structure with an appropriate mathematical framework in order to probe their properties and then mapping the structure to a property. Current approaches rely on computationally heavy methods for both stages, but now a team from Iowa State University and University at Buffalo simplify the process by representing the structures with labelled, weighted, undirected graphs. On this basis, one can create a “surrogate” model through generic physics graph descriptors (e.g. path lengths, domain sizes) and weighting functions for the particular property of interest. This approach reveals new designs for improved organic solar cells, but could expanded to other devices.


INTRODUCTION
Microstructure design, or morphogenesis, is the decision making process of determining material distribution to optimize performance. 1,2 Since the microstructure critically affects transport of mass, 3 charge 4 as well as reaction rate, 5 there has been a sustained focus on microstructure-sensitive design. The design goal here is to identify tailored microstructures that result in maximization of desired properties. Examples include designing the porous structure of membranes to enhance the filtering process for water reclamation, 6 designing the electrode microstructure in batteries to improve energy transport, 7 and controlling the microstructure in thin film organic electronics to improve energy harvesting 8 or sensing 9 capabilities.
Morphogenesis can broadly be decomposed into two distinct stages: the representation stage and the mapping stage. The first stage (representation stage) is to carefully construct a mathematical framework for representing and, most importantly, synthetically generating microstructures (by 'synthetically' we mean generated in silico). The second stage is to build the mathematical construct that maps the microstructure to a property. The first stage allows sampling candidate microstructures and perturbing them, while the second stage enables evaluating the performance of a microstructure. These two stages are then wrapped into an optimization framework to iteratively identify a microstructure with extreme performance. Current approaches to morphogenesis are limited in two ways: (a) the use of computationally complex methods for microstructure representation 10,11 -field based (phase field, level set, random fields), or feature based approaches-which make reconstruction or perturbing microstructures non-trivial, 12 and (b) the necessity of solving complex, non-linear partial differential equations based mappers 7,13 for evaluating the property of a microstructure. These limitations bottleneck rapid exploration of the search space, especially for geometrically constrained (thin films) and anisotropic systems where the number of constraints make exploration quite nontrivial.
Here, we show how morphogenesis can be substantially simplified and generalized by using concepts from graph theory. Labeled, weighted, undirected graphs can be used to represent microstructures. They provide a fast, accurate and extensible way to represent, perturb and generate microstructures. Furthermore, graph measures (or queries) can be used to construct viable surrogates for the microstructure-property mappers, which enable fast and efficient property evaluation. This is especially useful when considering microstructure-sensitive processes (like transport and reactions) that can be naturally represented in terms of pathways, connected components and interfaces, which are foundational concepts in graph theory. Additionally, expert intuition, available knowledge and data trends can be naturally fused into these surrogate graph metrics allowing easy integration of domain knowledge with data.
Finally, a graph based approach merges both stages of morphogenesis (microstructure representation and structureproperty mapping) in a computationally elegant and efficient way. We show that one can frame morphogenesis as a graph optimization problem and illustrate the approach by designing optimized morphologies for photovoltaic applications. 14 We demonstrate how this fast approach allows us to evolve an initial morphology into an optimized morphology exhibiting substantially improved short circuit current. We show optimized morphologies for thicker films (∼250 nm) than currently prevalent, which allows harvesting more incident energy. Our graph based morphogenesis is broadly applicable for microstructure-sensitive design for batteries, biosensors and filtration.

RESULTS AND DISCUSSION
Graph-based microstructure design framework Standard approaches for computationally predicting the performance of microstructure sensitive systems involve interrogating the microstructure of interest using computationally intensive frameworks. Incorporating these strategies into any inverse design problem becomes computationally intractable. We address this challenge in two stages. In the first stage, we develop a physicsbased 'graph surrogate model' of the performance of the device. The 'surrogate model' extracts the microstructure-aware performance measures and is built (a) using the underlying equivalence between a multi phase morphology and a labeled, weighted, undirected, graph, (b) using concepts and algorithms from graph theory, (c) to seamlessly allow incorporating domain expert knowledge and intuition. In the second stage, this simple and fast surrogate model is linked with a parallel probabilistic optimization algorithm that enables the inverse design of morphologies to maximize certain performance criteria. The choice of the probabilistic approach again allows seamless incorporation of domain knowledge in terms of the prior distribution (see Fig. 1).
The key idea here is to represent microstructure as a graph. A graph-based representation allows us to leverage well studied algorithms from graph theory along with their optimized implementations. We note that graphs have been used in many material science models ranging from molecular structure 15,16 to grain boundaries 17,18 and porous microstructure. 19 Discretized microstructure is represented as a labeled, weighted, undirected graph G = (V,E,W,L) (Fig. 1a). The graph captures all imperative information inherent to a digitized morphology. Each pixel (or voxel) in the microstructure becomes a graph vertex, with V being the set of all vertices in the microstructure. An associated labeling function, L, assigns a label to each vertex in V. For instance, considering a two-phase microstructure, L could associate a color (red) to represent one phase, and another color (blue) to represent the second phase. Thus L acts on each vertex in V. It should be noted that the set of vertex labels can be easily extended to account for multi phase system or to represent multiple characteristics (like electrical, thermal and magnetic properties, or crystalline vs amorphous properties) of each vertex. Each vertex is connected to its neighboring vertices through edges. The set of all edges in the graph is denoted as E. Finally, the weight function, W : E ! R þ , assigns a non-negative real weight to each edge in E. For instance, the weight function could simply encode the physical distance between voxels. Alternate weight functions could encode hopping frequencies of charge/mass/momentum transfer between vertices or other transport characteristics. This G = (V,E,W,L) based representation provides a lot of flexibility in representing a variety of microstructures and associated spatially varying properties. On top of the versatility of this graph representation of a microstructure, we gain the ability to use very efficient graph algorithms to extensively characterize the graph (morphology). In general, most transport and reaction characteristics that depend on the underlying morphology can be recast as graph measures. 19,20 Given graph-based representation of the microstructure, a surrogate model of material properties is constructed directly on the graph-based representation. Traditionally, there have been a spectrum of approaches for constructing a mapping from a given microstructure to its associated device performance indicator. On one end of the (complexity) spectrum are the full-physics based approaches that use detailed microstructure-aware propertystructure simulators as the mapping function. These typically come at a very high computational cost precluding optimization except on large computing resources. On the other end of the spectrum are the approaches that use a small set of salient microstructural features from the microstructural data that is then mapped to the performance indicator. These features encode domain knowledge, are intuition driven and/or a result of carefully designed hypothesis driven structure-property studies. A classic example of the latter approach is the Hall-Petch relationship, where the average grain size (salient feature) governs the yield strength of the material. 21 We emphasize that representing a digitized morphology as an equivalent graph allows one to construct a 'graph surrogate model' (based on graph metrics) that can span this full spectrum of structure-property mappers.
A surrogate model is an approximate (but fast) mapping from the microstructure to the output, usually constructed with the aid of statistical and/or data-driven approaches. 22 If a small set of fullphysics structure-property simulation results are available, the surrogate model is created using this limited set of results. The surrogate model construction (or 'training') can be done by extracting a large set of graph descriptors/measures from the graph and identifying the best set of features that create a good regression to the full-physics predictions. 23 Representing the full microstructure as a graph helps avoid any a priori feature extraction, and builds the surrogate model on the same graphbased representation. Additionally, intuition driven or domain   Fig. 1 Elements of graph-based microstructure design framework. a The concept of representing morphology as an equivalent graph allows construction of fast graph based surrogates to full-physics based device performance indicators; b Scheme of population based incremental learning: An initial (or prior) distribution is sampled to generate candidate graphs (i.e., microstructures), which are rapidly evaluated and ranked. This ranked information is used to update the distribution. Iterating this process evolves the distribution towards graphs exhibiting better performance. Such an approach naturally allows incorporation of expert information into the prior knowledge based features can be easily represented as graph metrics and incorporated into a surrogate mapper. 20 Examples include evaluating the connectivity in porous media as graph connectivity metrics and graph paths, domain sizes in ceramics and alloys as graph connected components, among several other examples that originated from network-based computational models for characterizing transport properties. [24][25][26] To keep the approach generic and allow diverse applications to a variety of problem specific traits, we separate construction of the surrogate model into two steps: defining a set of physics relevant graph descriptors that are generic (i.e., application agnostic), and weighting functions that encapsulate the specifics of the property of interest (i.e., application dependent). The generic graph descriptors include traits such as path lengths (from specific locations, say, boundaries to all the internal nodes or from one phase to another), centrality measures (of, say, reactants or initiators), connectivity measures (of, say, phases in the microstructure with the boundary or to each other), domain sizes, among others. 27,28 Next, weighting functions are constructed to capture the detailed trends (say, kinetics aspects, or material specific aspects) of the transport. For example, effect of Brownian motion (or random walks) can be encoded by defining a weighting function that associates smaller weights with larger distances that an ion/charge/molecule has to travel on the graph. This is integrated with the graph path length measures to construct a viable surrogate for the property of interest. This generality provides substantial flexibility to creatively construct reliable surrogate mappers.
Probabilistic graph-based optimization The graph based surrogate model provides an efficient map from microstructure to performance. We link this surrogate model with a probabilistic optimization algorithm. We specifically choose a probabilistic approach for the following reasons: (a) the outcome of the optimization is a distribution, i.e., a microstructure class (rather than specific microstructure), (b) such probabilistic approaches perform global searches in a gradient free manner and are embarrassingly parallelizable, (c) it is straightforward to encode domain knowledge as a prior distribution, and (d) the optimization is performed directly on graphs, hence eliminating any conversion steps.
In this work, we use a probabilistic optimizer called populationbased incremental learning (PbIL). 29 The basic steps of the algorithm are outlined in Fig. 1 (b). PbIL is an optimization algorithm that estimates the explicit distribution of the optimal solution rather than searching for a specific solution, as in typical meta-heuristics (e.g., genetic algorithms, 30 hill-climbing 31 or simulated annealing 31 ). The multi variate probability distribution is stored as a probability vector P of vertex labeling. In particular, with each vertex in the graph (variable) we associate the probability of a given phase (or component, or more generally the label, L) appearing at this location in the morphology. The probability vector P is updated in each iteration of the optimization based on the feedback from the cost functional. The optimization starts with a given probability vector of vertex labeling (step S1 in Fig. 1(b)), which is either uniform or reflects domain expert knowledge (see Fig. 3). In each iteration, n morphology realizations (i.e., instances of a colored graph) are sampled from the current probability vector P (step S2 in Fig. 1(b)). For each realization, j, the graph surrogate model is deployed to evaluate the performance indicator (cost functional), f(j) (step S3 in Fig. 1(b)). The n b best samples n b ( n ð Þare identified (step S4 in Fig. 1(b)). These best samples are used to calculate the update vector, P u , of the probability distribution (step S5 in Fig. 1(b)). Specifically, the update vector P u represents the maximum likelihood estimation of vertex labeling from the n b best samples. An iteration ends by updating the probability vector P ¼ P Á ð1 À l r Þ þ P u Á l r , where l r is the learning rate (step S6 in Fig. 1(b)). Intuitively, the update reinforces features present in the morphology, and dampens those missing, at the constant rate prescribed by l r . The algorithm iterates until standard termination criteria are met (iteration limits, improvement bounds). The resulting method has a small memory footprint with only one vector of vertex labels being stored and updated at each iteration. The method is also embarrassingly parallel, since multiple samples drawn from the evolving distribution vector P can be evaluated simultaneously. The cost functional f expresses a single objective to optimize. However, we note that in practical settings, f can encapsulate multiple criteria, as we demonstrate in the next section.
Application-specific exploration We illustrate the power of this approach by optimizing the morphology of organic solar cells (OSC). OSCs have the potential for widespread usage due to their low cost-per-watt, mechanical flexibility, and ease of manufacturability. 14 Their widespread use and commercialization, however, are bottlenecked primarily by relatively low solar efficiencies. Efficiency can be improved either by changing the material chemistry (molecular engineering), microstructure tailoring (microstructure engineering) or architecture (device engineering). It has been argued that most advantages from tuning material energy levels (i.e., molecular engineering) have already been achieved, 32 leaving microstructure design as a promising avenue for efficiency increase. 33,34 Traditional approaches to microstructure engineering in OSC have relied on systematic trial-and-error experiments to generate diverse microstructures and their device performance (which is resource intensive), or via using computationally intensive simulators to virtually interrogate the candidate morphologies to identify features that an optimal morphology should exhibit. 35,36 It is generally believed that an interdigitated columnar morphology (Fig. 3) is an optimal morphology. This is because the interdigitated morphology exhibits a graceful combination of two competing requirements-having a large donor-acceptor interface to enable efficient exciton dissociation, while still providing straight, fully connected domains to enable efficient charge transfer and collection. We substantially improve upon this current state-of-the-art using our graph based morphogenesis (Fig. 3) and show that a whole class of morphologies exist that not only exhibit a much higher photovoltaic performance, but which may also be easier to fabricate.
Following the template for graph based morphogenesis outlined in subsection 'Graph-based microstructure design framework', we represent a OSC microstructure-consisting of two phases, electron-donating and electron-accepting-as an labeled, weighted, undirected graph. 20 We then utilize a finite set of full physics simulations to identify three physically meaningful graph descriptors that together produce an exceptional surrogate for device performance (specifically, the short circuit current) 23 (see Supplementary Information for more details). Figure 2 illustrates the quality of the graph surrogate model. Our surrogate for the short circuit current produced by a morphology is given as: F ¼ f abs f diss f out , where these three descriptors are: (i) f abs : weighted fraction of photoactive material. This accounts for that fraction of the microstructure that absorbs incident light and produces excitons. This also accounts for decaying light intensity as light propagates through a semitransparent media (light intensity decays with distance from the top electrode); where n D is the total number of electron donor vertices, d H is the distance to the top surface, w abs is the weighting function that encodes the probability of light absorption. We consider three scenarios here: w abs = 1 homogeneous light absorption, w abs ¼ expðÀd H =λÞ and w abs ¼ P P p¼1 a p ðd H Þ p is the polynomial of order p fitted to the results of Maxwell equation for varying thickness of the device and reflection from the bottom surface accounted for (details provided in Supplementary Information).
(ii) f diss : weighted fraction of donor with weighting based on the distance to the donor-acceptor interface. This feature accounts for the ability of the microstructure to successfully dissociate the produced excitons. Excitons generated in the donor phase diffuse (i.e., via Brownian motion) and are dissociated only if they reach a boundary. This graph based metric accounts for the average (weighted) distance from any donor vertex to the nearest donor-acceptor interface; where d I is the shortest distance from any donor vertex to the nearest acceptor, w diss ¼ expðÀd I =L d Þ is the weighting function that encodes the probability of exciton reaching the interface, and L d is the exitonic diffusion length. Note how domain knowledge (the exciton dissociation length) is gracefully incorporated into the surrogate model. (iii) f out : fraction of complementary paths to the electrodes. This graph feature encodes the availability (and length) of charge transport pathways for the produced electrons and holes to move through microstructure and reach the electrodes where they are collected to form useful current: where n C D is the total number of donor vertices that are connected to anode, n C A is the total number of acceptor vertices connected to cathode. Both set of vertices are identified using the connected component algorithm, 27 which is a standard graph algorithm. For each vertex the shortest path towards the respective electrode is determined along with the tortuosity. Tortuosity is defined as the ratio between the shortest path via the appropriate domain (donor or acceptor) versus the shortest straight path without any obstacles. A tortuosity of one corresponds to a straight rising path, where the shortest path is equivalent to the straight path without any obstacles. The tortuosity is calculated for each vertex and subsequently weighted by the function w tort ¼ expðÀ4ðτ À 1ÞÞ. The weighting function penalizes the curvy pathways where more recombination can occur and reinforce the straight paths.
To demonstrate the robustness of the surrogate model, Fig. 2 depicts the correlation between graph-based descriptors (surrogate model) and performance metrics from exitonic drift diffusion model (physics-based model). 37 We plot three metrics, full performance evaluation and two major contributions to the The morphogenesis from an initial guess towards the final optimized structure. (Bottom) Effect of initial probability field on the morphogenesis outcome. Three initial probability fields: bilayer, uniform and interdigitated initial field evolve towards consistently similar morphologies with statistically similar properties and fractal dimensions performance, namely exciton dissociation and charge transport. The correlation coefficient for three metrics are 96.65%, 96.95% and 58.81%, respectively. In each subplot, we plot data for three resolutions of the adaptive exploration: 20 × 20, 40 × 40, and 80 × 80. One data point corresponds to one iteration of the optimization. This analysis was performed on data collected during optimization. In particular, for each iteration of the optimization, we saved the best performing intermediate candidate morphology and run the physics-based simulation. Our correlation studies show that the surrogate model is well correlated with the physics-based model at all stages of the optimization. Figure 3 (top panel) illustrates results for the evolution of the morphology from an initially uniform random distribution. We consider a device thickness of 100 nm, which is typical for OSC devices. We assumed a homogeneous light absorption scenario (i.e., a non-reflecting bottom surface). The PbIL algorithm identifies a dendritic type structure, which is uncommon in OSC's. Interestingly, this structure exhibits a bicontinuous interpenetrated network of well balanced donor and acceptor domains with domain size comparable to the exciton diffusion length (which is set as L d = 10 nm). This dendritic structure exhibits a graceful combination of two competing requirements-having a large donor-acceptor interface to enable efficient exciton dissociation, while still providing straight, fully connected domains to enable efficient charge transfer and collection. In fact, the currently believed optimal morphology-an interdigitated columnar morphology-is though to be optimal for precisely this reason.
We virtually interrogate the dendritic microstructure using a morphology-aware physics-based device simulator. 37 The device simulator calculates the short circuit current density, J sc by solving the associated exciton-drift-diffusion equations (see Supplementary Information). The predicted J sc for the dendritic structure is 8.9 mA cm −2 , which is substantially larger than the J sc = 5.28 mA cm −2 exhibited by the interdigitated columnar morphology. As anticipated, this is due to enhanced exciton dissociation with minimal increase in recombination. More specifically, the increased donor-acceptor interfacial area results in enhanced exciton dissociation, while the charge transport pathways remain relatively non-tortuous, resulting in minimal performance degradation due to recombination.
Incorporating domain knowledge into morphogenesis We next account for domain knowledge about features of an optimal microstructure into the graph based optimization. Instead of starting from a random initial distribution, we repeat the morphogenesis by initializing the material distribution with the interdigitated columnar morphology (Fig. 3). The ensuing optimized microstructure improves upon the initial guess, and exhibits a J sc = 9 mA cm −2 . Interestingly, the optimization only changes the material distribution very slightly (see movie M1 in Supplementary Information) by creating the secondary dendritic arms that increase the donor-acceptor interfacial area. We also used a bad initial guess by using a bilayer morphology (which exhibits a J sc = 0.8 mA cm −2 ) as the initial distribution. Here, the optimization produces both primary and secondary dendrites (see movie M2 in Supplementary Information) as the optimization proceeds (Fig. 3).
Universal scaling Searching for a distribution rather than a specific solution seems to make the framework agnostic to initial distributions, with all three optima exhibiting similar J sc . We qualitatively check this by rerunning the optimization five times, with different random seeds, for each of the three initial distributions considered (random, bilayer, interdigitated). All 15 morphologies exhibited very similar J sc with a 3% variability across the optima. While visually dissimilar, all exhibited similar statistical properties, especially a fractal dimension around 1.5. This suggests that a whole class of morphologies could potentially be optimal. Moreover, we found the fractal dimension to be close to that of quadratic cross (1.49 38 ). This is consistent with dendritic structure where perpendicular side branches mimic quadratic cross topology. Although morphologies we report in this paper have columnar dendritic structure, it is interesting to see consistent fractal dimension with quadratic cross topology.
Effect of changing device architecture We next explore variations in the optimal morphology when light reflecting boundary conditions are imposed at the bottom electrode. As we emphasized earlier, accounting for this new aspect of physical phenomena is trivial: we simply augment the graph surrogate with a weighting term that encodes a thickness dependent light intensity (see Supplementary Information for Fig. 4 Effect of changing device architecture: light absorption scenario and device thickness. Two device thicknesses and three light reflection scenarios evolve towards consistently similar morphologies. Optimized morphologies demonstrate superior properties for a range of device configurations details). Figure 4 depicts the optimized microstructure for this scenario. Notice that the absorption of the reflected light favors more donor phase in the middle of the active layer. However, the morphology still exhibits a dendritic structure with nearly the same fractal dimensions.
Optimal morphology for thicker devices We finally explore variations in the optimized microstructure with increasing device thickness. This is a pressing problem in commercializing OSCs as increasing the device thickness can potentially enable absorbing (and converting) more incident radiation, but usually results in poor performance due to the concomitant increased recombination resulting from increased charge transport pathways. Figure 4 compares the optimized microstructures for thin (100 nm) and thicker devices (250 nm). The thicker device has 28.9% more incident radiation that is absorbed. Interestingly, most of this additional incident energy is converted to collected current, with the short circuit current improving from J sc = 8.05 mA cm −2 for the 100 nm device to J sc = 10.09 mA cm −2 for thicker devices. Notice that the morphology for the thicker device exhibits a similar dendritic structure with similar fractal dimensions (∼1.5). However, we notice two types of dendrites within the thicker devices. Dendrites with only first order side branches are promoted close to the top, while higher order branching is promoted in the bottom of the device. The fractal dimension of this structure is marginally larger than that exhibited by the thinner device, due to this complexity. The heterogeneity of the optimized structure demonstrates the robustness of the surrogate model in the search of structures with spatially varying complexity.
In summary, microstructure-sensitive design has become an indispensable part of the materials discovery paradigm. However, the identification of optimal microstructures (morphogenesis) has traditionally been hindered by twin bottlenecks, namely nontrivial approaches for microstructure perturbations (i.e., searching the microstructure space), and complex approaches for microstructure-to-property mappings. Here, we illustrate how both these issues can be efficiently resolved via a graph-based strategy. Treating microstructures as graphs allows efficient, modular, and extensible representation which results in simple approaches to perturb and explore the microstructure space. Simultaneously, graph based measures (that can be rapidly, and efficiently computed) make excellent surrogates of microstructure-to-property mapping. This graph based representation and graph-to-property surrogate mapper can be integrated with a probabilistic optimization strategy to efficiently identify optimal microstructures. A graph based approach furthermore allows natural incorporation of domain knowledge into the design process. We illustrate a proof-of-concept application by designing the optimal structure for OSC. This resulted in the identification of a new class of microstructures that exhibit better performance than the currently hypothesized optimal microstructure. We anticipate and look forward to widespread use of this microstructure design strategy by the materials community.

METHODS
For the design framework, an in-house C++ code is used that implements probabilistic graph-based optimization. In each iteration of the optimization, n = 4000 samples are generated from the current distribution. Each realization is evaluated according to the cost function, and n b = 100 best microstructures are chosen to compute the update vector P u . The learning rate is set as 0.1. The search terminates when no improvement in cost is recorded over 100 consecutive iterations. The cost function is defined through the surrogate model evaluated directly on the graph as discussed in the text. The model is constructed to approximate a short circuit current of the device. Specifically, the full model consists of the exitonic drift diffusion model and is used to calculate the short circuit current J sc . More details can be found in 23 and Supplementary Information. To determine statistical similarity, all microstructures are characterized with fractal dimension. In particular, the Hausdorff (box-counting) fractal dimension has been computed for binary image using MATLAB function hausDim. More details are given in Supplementary Information.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.