Understanding Causalities in Organic Photovoltaics Device Degradation in a Machine‐Learning‐Driven High‐Throughput Platform

Organic solar cells (OSCs) now approach power conversion efficiencies of 20%. However, in order to enter mass markets, problems in upscaling and operational lifetime have to be solved, both concerning the connection between processing conditions and active layer morphology. Morphological studies supporting the development of structure–process–property relations are time‐consuming, complex, and expensive to undergo and for which statistics, needed to assess significance, are difficult to be collected. This work demonstrates that causal relationships between processing conditions, morphology, and stability can be obtained in a high‐throughput method by combining low‐cost automated experiments with data‐driven analysis methods. An automatic spectral modeling feeds parametrized absorption data into a feature selection technique that is combined with Gaussian process regression to quantify deterministic relationships linking morphological features and processing conditions with device functionality. The effect of the active layer thickness and the morphological order is further modeled by drift–diffusion simulations and returns valuable insight into the underlying mechanisms for improving device stability by tuning the microstructure morphology with versatile approaches. Predicting microstructural features as a function of processing parameters is decisive know‐how for the large‐scale production of OSCs.


Introduction
Power conversion efficiencies (PCEs), the feasibility of up-scaling manufacturing, and long-term stability of organic solar cells (OSCs) are major concerns in commercializing organic photovoltaics.Owing to the synergistic efforts on material design, morphology control, and interfacial engineering, etc., [1] device efficiencies of OSCs are getting closer to the commercialization benchmark PCE of 20%, [2] which showcases the great potential of this emerging photovoltaic technology for future applications.Commercialization of OSC products requires as well consistent high performance during operation for at least 10 years, which poses a significant challenge for OSC systems. [3]Central to performance and stability of OSCs is the morphology of blend films.It is known that the solidstate equilibrium morphology is dominated by the chemical structure, while the actual morphology is dominated by the processing conditions such as post-treatment, [4] processing additives, [5] and solvent choices, [6] as well as multicomponents. [7]Material and processing properties together control phenomena such as intrinsic aggregation behavior, nucleation, and crystal growth or phase separation as a function of miscibility. [8]With the rapid development of novel materials, numerous process parameters for morphological modification have to be explored in multidimensional optimization in order to achieve excellent device efficiency and stability, and thus a huge library of morphological features is generated.So far, traditional methods for morphological investigation of OSCs typically rely on sophisticated equipment, like grazing-incidence wide-angle X-ray scattering (GIWAXS), transient absorption spectroscopy (TAS), and transmission electron microscopy (TEM), etc. [7b,9] Electron microscopy or scattering methods even have been optimized to qualitatively but also quantitatively analyze the microstructure of complex composites over the last decade.Nevertheless, these sophisticated technologies require expensive equipment, experienced manual operators, and time.These limitations make it difficult to build up a large enough database with characteristic parameters over multiple material systems in order to analyze statistical trends, not to mention the difficulty in integrating the characterization routines into daily laboratory routines or even automatizing high-throughput research lines.
The advent of data-based research methods has opened a window of opportunity to train fast and cheap measurements towards providing equally valuable insight as the classical knowledge-based research methods.It is now the time to accompany the technological evolution of OSCs by establishing the fundamental understanding of causality between processing conditions and microstructure formation.We have recently demonstrated that Gaussian process regression (GPR) has the inherent capability to predict electrical performance and its degradation from features of the active layer morphology, [10] which we obtained from spectral decomposition of UV-vis spectra according to the Spano model. [11]Although parameterized for weak H aggregates, the method also works for small molecules with strong J aggregates and even for all-polymer systems. [12]n this work, we refine this technique with an embedded minimum redundancy-maximum relevance (mRMR) feature selection algorithm [13] (for details see Section E, Supporting Information) to identify the minimum feature set carrying the essential information to predict the target quantity.The causal connection between processing conditions, morphology, and long-term stability is established on the working horse polymer, PM6, and the non-fullerene acceptor (NFA), DTY6, with long side chains derived from the BTP family, by varying the processing solvents (halogenated chloroform and non-halogenated o-xylene) and the addition of a ternary component, [70]PCBM (full names are listed in the Supporting Information part).The PM6:DTY6 blend is attractive due to the fact that its device efficiency is not affected when switching the processing solvent from chloroform (CF) to the more eco-friendly o-xylene (o-XY).[6a] Experimentally obtained device stability shows that ternary devices exhibit a significant suppression of the fill factor (FF) loss, compared to those from binary devices, and o-XY-based devices exhibit a better operational stability property than those processed from CF.The morphological information of these samples is tracked with readily available UV-vis absorption spectra.We obtain morphological fingerprints specifically for electron donor and acceptor materials from automated spectral modeling. We ten combine a feature selection technique with Gaussian process regression to find deterministic relationships linking only the relevant morphological features to processing conditions and device functionality.This allows the distinction between direct and indirect influences on the stability and allows a quantitative assessment of each individual influence.Accordingly, our quick and elegant method unveils that the admixture of [70]PCBM increases the molecular ordering in the acceptor phases which makes the composite less prone to light-induced traps and is thus responsible for the enhanced operational stability.In addition, we model the effect of the active layer thickness and order by drift-diffusion simulations and obtain hints for the underlying mechanism for improved device stability.Our method provides a cheap and reliable strategy to efficiently identify and understand the underlying relationship between morphological stabilization and processing conditions.

Results and Discussion
Figure 1A shows the chemical structures of the polymer donor PM6, the NFA DTY6, and the additive [70]PCBM component, as well as the two processing solvents (CF and o-XY) applied in this work.Device performance and device stability were explored based on an inverted configuration of ITO/tin oxide (SnO 2 )/active layer/molybdenum trioxide (MoO 3 )/Ag (Figure 1B).The binary PM6:DTY6-based OSCs from o-XY exhibit a promising PCE of 13.85%, which is comparable to that of the binary device from CF (13.95%).Furthermore, the PCEs of the ternary devices from o-XY with the addition of 20 wt.% [70]PCBM are slightly improved, mainly due to the increase of the opencircuit voltage (V oc ) by ≈20 mV (Figure 1C).In addition, all the photocurrent densities attained from external quantum efficiency (EQE) spectra agree with the short-circuit current density (Jsc) values obtained from the J-V characteristics (Figure 1D).The corresponding device parameters are summarized in Table S1 (Supporting Information).
Different processing conditions can regulate the bulkheterojunction (BHJ) composites into different structures and morphologies, [14] which are closely related to the device lifetime.Herein, to explore the influence of processing conditions on the operational stability, the device stability of CF-processed and o-XY-processed devices was monitored under illumination and thermal conditions in a home-built characterization setup in a dry nitrogen (N 2 ) atmosphere with oxygen and water level below 0.5 ppm, mimicking perfectly encapsulated devices.Figure 1E-I depicts the electrical performance evolution of four kinds of OSCs under ≈40 °C and continuous 1-sun-equivalent illumination provided by white light-emitting diodes (LEDs) (Figure S1A, Supporting Information).PM6:DTY6-based OSCs from CF and o-XY exhibit a different photostability.Specifically, the two binary devices show similar degradation of Voc and FF but exhibit a different trend for J sc degradation.Accordingly, the maintained PCE is improved from 62.0% of the initial PCE in CF to 66.9% of the initial PCE in o-XY after aging for 1550 h under 1-sun-equivalent illumination.
The significant photostability losses in the binary cells can be mitigated by introducing a small amount of a guest component, [70]PCBM, into the PM6:DTY6 binary blend.In line with literature on similar material systems, [15] both ternary OSCs processed from CF and o-XY show clearly enhanced operational stability, which originates from a remarkably improved FF stability, along with slightly enhanced Jsc and Voc stability.Moreover, the FF stability property of the o-XY-processed ternary OSCs outperforms that of CF-processed ternary OSCs.The stability data reported here were collected from more than seven devices (Figures S2-S5, Supporting Information).In addition to the stability studies at a moderate temperature, we further explored thermal stabilities by continuously aging at elevated temperatures of 65 and 85 °C for over 170 h.As presented in Figure S6 (Support-ing Information), the degradation trends for the various composites are independent of the harsh aging conditions.Clearly, degradation is accelerated at higher temperatures, but the relative stability of ternary devices is again higher compared to the binary ones.When increasing the aging temperature to 65 or 85 °C, the overall device parameters of o-XY-processed ternary OSCs are well preserved (the V oc drops due to its temperaturedependent nature). [16]Likewise, the stability properties of all devices tested under room temperature and metal halide lamps (MHLs) (Figure S7, Supporting Information) demonstrate a In order to find causal relationships between processing conditions, active layer morphology, and stability, we developed an integrated machine-learning workflow combining experimentation and model-based simulations, based on the extraction of a set of morphological features from spectral modeling of UVvis absorption spectra, shown in previous work. [10]As the trained spectral model considers both ordered and amorphous phases for both donor and acceptor individually, more than a dozen of morphological features can be obtained, some of which are in part redundant (see Section F, Supporting Information).Therefore, we embed a minimum redundancy maximum relevance (mRMR) feature extraction method in the modeling of the predictions by Gaussian process regression (GPR).A schematic of the workflow is shown in Figure 2. In the following, we first predict the FF loss, quantified by r FF,200hrs , the remaining FF after 200 h divided by the initial FF, from the morphological features.This will show us the most important morphological aspect of controlling FF loss.We focus on FF loss because NFAs are known to expose this kind of performance instability. [17]Then we will show cross-correlations between the main morphological feature and other aspects of morphology, which will give us hints for underlying mechanisms.Next, we show how these features can be controlled by processing conditions.Finally, we can use the obtained morphological trends to test scenarios for the mechanisms of degradation.
Figure 3A shows the result of the mRMR-GPR analysis of the processing-morphology-stability interdependencies in PM6:DTY6:[70]PCBM devices.The arrows point from the predictor to the target, with the contribution to the overall explanation of variance (coefficient of determination R 2 ) annotated.It is most instructive to read this graph structure from right to left.The FF loss (that is, r FF,200hrs ), is explained by one single predictor, namely the apparent Huang-Rhys factor h A , which is related to energetic disorder in the DTY6 phase.The assignment to disorder is based on the exciton theory and is unambiguous because there is no spectral congestion in the respective wavelength range.This predictor alone explains 58% of the total variance in r FF,200hrs .All other morphological predictors are correlated with h A and do not provide a significant additional explanation of variance.This also means that 42% of the variance in r FF,200hrs is not explained by any of the available features.The unexplained variance may either be due to uncertainty in the measurement of FF loss or due to further yet undiscovered predictors for FF stability for which UV-vis spectroscopy provides no extractable features.The distinction of these phenomena will require a larger and more varied dataset and the inclusion of further spectroscopic probes, e.g., photoluminescence spectroscopy. [18]alking now from right to left in Figure 3A, we consider h A as a target and predict it from the other morphological features by mRMR-GPR.We find three morphological predictors, together explaining 99% of the variance in h A , namely the total absorption of DTY6 excitons (A tot ), the energetic disorder in DTY6 (w A ), and the DTY6 exciton energy (c A ).The latter two correlations are to be expected from the Spano model and constitute an additional confirmation of its validity in our samples.A correlation between A tot and h A points to the often-observed fact that thicker films are more disordered. [10,19]The same is true for the correlation of D tot (the total absorption in the PM6 donor phase) with c A .Although care has been taken to keep the thickness of the samples in the same range, some fluctuations cannot be avoided.We will show below that by using appropriate surrogate functions, we can obtain quantitative trends with associated uncertainties from imperfect experiments.
Walking further to the left in Figure 3A, next we consider A tot , w A , c A , and D tot as targets and use mRMR-GPR to predict them from our applied processing conditions, namely the amount of [70]PCBM (X 70PCBM ), and the solvents, here represented by its -BP at a pressure of 1000 hPa (T B ). [18] We find that both A tot and D tot are mainly controlled by the solvent, only a small part being due to the [70]PCBM additive.X AmA is mainly controlled by admixture of [70]PCBM.Finally, the c A is mainly controlled by the solvent, simply because we already know that c A depends on thickness. [10]n Figure 3B, we show the surrogate function r GPR FF,200hrs = f (h A ) as blue solid line, showing that it reproduces the trend in the experimental data points, with different colors indicating different processing conditions.The dark blue shaded area gives the 95% confidence interval for the surrogate function, showing substan-tial uncertainty in the slope, especially for high values of h A where only a few data points are present.However, the confidence interval, while large, remains at negative slope, quantitatively confirming, on the basis of the available evidence, that energetic order in the ordered phase (low h A values) yields a higher FF stability.
Since Figure 3A has shown that thickness (via A tot and D tot ) influences h A , and that there is also a correlation between X 70PCBM and A tot , we need to distinguish how much of the observed FF stabilization effect is caused directly by the presence of [70]PCBM, and how much is indirectly caused by a different thickness.Usually, this is accomplished by "conditioning" (fixing) one of the parameters in question.As mentioned above, in material science this can be difficult to accomplish.However, even if such a dedicated experiment is unavailable, we can still use the GPRderived surrogate function to virtually condition (fix) one of the predictors (for details, see Section G, Supporting Information).In Figure 3C, we show the 2D surrogate function r GPR FF,200hrs = f (X 70PCBM , A tot ) (colored hypersurface), showing that even for constant A tot , there is a clear effect of X 70PCBM on r FF,200hrs .(For the confidence interval of the surrogate function see Figure S14C, Supporting Information).We thus confirm that the presence of [70]PCBM stabilizes the FF, despite the thickness variations in our samples.In Figure 3D, we further perform the same analysis for the solvent effect, showing the 2D surrogate function r GPR FF,200hrs = f (T B , A tot ).We find that at constant A tot there is very little variation of r FF,200hrs ; consideration of the confidence intervals in Figure S15C (Supporting Information) shows that it is insignificant.Thus, we learn that the contributions of the solvent to the FF loss in our dataset occur only indirectly due to the different film thicknesses, and no significant direct solvent effect can be identified.Finally, Figure 3E shows that both, admixture of [70]PCBM and the high-BP solvent o-XY, lead to a red shift of the DTY6 exciton, which is indicative of larger delocalization and thus higher order in the ordered phase of the acceptor.In summary, from Figure 3, we can explain the effect of admixture of [70]PCBM on the FF stabilization by an improved order in the acceptor phase, while we found that the solvent influences stability mainly indirectly, via its influence on the thickness of the active layer in our case.This shows the potential of our method yielding trends and their significance already from a limited dataset.This means that in the chosen deposition method, the solvent is neutral to both device performance and device stability, strongly suggesting the use of o-XY over CF for environmental and regulatory issues.Finally, we want to rationalize the results, obtained from an essential data-driven approach (mRMRembedded GPR), by comparing them with drift-diffusion simulations, a knowledge-based approach, which allows us to compare the expected FF loss from certain degradation scenarios to the experimental results.For the procedure see the right part of Figure 2. To accomplish a convergence with available experimental evidence, we will need to stipulate the model scenarios such that they yield a functional relationship between FF loss and a set of quantities for which we at least have qualitative predictors.
From Figure 3, we have evidenced that both the film thickness and the order in the acceptor phase influence FF loss.By calibrating some of our samples, we obtained the absolute active layer thickness for all our samples from their total absorption, by applying Lambert-Beer's Law (for details, see Section H, Supporting Information).According to the Spano model, the morphological feature h A scales with disorder in the acceptor phase; therefore, we use the inverse quantity 1/h A which should scale with order therein.Thus, we can display the FF loss as a function of active layer thickness and order, see Figure 4A.Although the number of available experimental data points (colored symbols) is limited, they are sufficiently independently distributed such that a GPR (colored hypersurface) finds clear trends for the FF loss with respect to both predictors.
In order to replicate these trends by our model, we assume that FF loss is caused by bulk traps, which is justified by Figure S10.We further assume a monotonous trend between the order in the acceptor phase and the extraction mobility for electrons μ n , because there is evidence in the literature that domain boundaries-and thus disorder-limit charge carrier mobility. [20]e distinguish between two different types of defects.On the one hand, one class of defects can cause structural disorder that dominantly impacts mobility but do not cause additional recombination.On the other hand, the second class of defects is causing traps that directly lead to recombination-these we classify as bulk traps.In the simplest possible scenario, we assume that the order in the acceptor layer has no influence on the increase of bulk traps.Under this assumption, we can model the FF loss during degradation by a tenfold increase in the concentration of bulk traps for all samples.In this scenario, the role of orderand thus extraction mobility -is just to allow for fast carrier ex- traction outperforming trapping.Since the active layer thickness controls the extraction field under the built-in voltage, and thus the drift velocity and the traversing time of the carriers, this simple scenario actually ascribes to both parameters the role of extraction accelerators, providing resilience against an increased trap concentration.From this model, the FF loss as function of bulk trap density tr B , active layer thickness d, and extraction mobility μ n , can be simulated by drift-diffusion (DD) simulations assuming a quasi-1D virtual semiconductor. [21]We have successfully modeled some of the current-voltage curves of the devices of the present dataset in their pristine state to obtain a reasonable set of system parameters (for details see Section J, Supporting Information).Using these parameters as starting point, we created a simulated dataset of 100 devices, randomly varying L between 75 and 225 nm, and log 10 (μ n /m 2 V −1 s −1 ) between −8.5 and −6.5.As Figure 3A suggests that mainly the acceptor phase is affected by processing conditions, we kept the hole mobility constant in the simplest scenario; however, we also studied scenarios where both hole and electron mobilities were varied (see Figure S18D, Supporting Information).For each {L,μ n } combination, we ran a DD simulation with high and another one with low density of bulk traps, and divided the resulting FF values to obtain the quantity FF tr(bulk)=2×10 19 ∕FF tr(bulk)=2×10 18 , which we can identify with the experimental quantity r FF,200hrs .
The result of this random variation is shown in Figure 4B.Qualitatively, the DD simulation shows the same trends as the experimental dataset, compared with Figure 4A.However, a closer inspection reveals quantitative differences.First, the experimental data show a much weaker thickness dependence than predicted by the DD simulation, compared with the colorbars from Figures 4A,B, respectively.As shown in Figure S18B (Supporting Information), this ratio depends on the assumption of the final trap densities after degradation.Second, the DD simulations show a saturation tendency for the FF loss for higher extraction mobilities (contour lines not equidistant along the μ n axis) while for the surrogate function in Figure 4A, no such saturation trend is found.This suggests that even for an order parameter as high as 1/h A = 5, the corresponding μ n values are still far lower than those in Figure 4B yielding a saturation trend at comparable thicknesses.This in turn suggests that a further FF stabilization is achievable by further increasing the order in the acceptor phase.For the simulation of more complex scenarios (trap density is dependent on extraction mobility, both hole and electron mobilities changing), we refer to Figures S18C,D (Supporting Information).In summary, the DD simulations show that an effect of molecular order on the FF stability can be expected also in other D:A systems provided that the DD parameters are similar to the ones assumed in the calculation of Figure 4D.This however also means that a general design rule cannot be deduced, as there may be D:A blends which due to a different electron configuration are not prone to trap formation, thus limiting the effect of disorder on the FF stability.

Conclusion
We have studied the effect of processing conditions (solvent choice, and the admixture of [70]PCBM) on the electrical performance degradation of PM6:DTY6 based OSCs, with the goal to understand the influence of microstructure on device efficiency and stability.We obtained morphological features from UV-vis absorption spectra in an automated machine-learning pipeline.Using Gaussian process regression (GPR) embedded with minimum redundancy maximum relevance feature extraction (mRMR), we have been able to distinguish correlation from causation, individually quantifying the effects of blend film thickness, solvent, and [70]PCBM additive on the degradation-induced FF loss.We found that [70]PCBM increases FF stability by reducing energetic dispersion (improving on-chain order) in the ordered phase of the DTY6.On the other hand, the change of solvent (abstracting from the effect of thickness) did not lead to significant changes in FF stability in this case.Still, this work demonstrates that halogenated solvents can be easily replaced by "greener" solvents, achieving environmental protection without compromising device performance and stability.Furthermore, the effect of active layer thickness and the morphology order on stability is modeled by using drift-diffusion simulation, offering hints of the underlying degradation mechanisms.Using driftdiffusion simulations, we have shown that we can qualitatively reproduce the observed trends by assuming that molecular order and active layer thickness both accelerate charge extraction thus lending resilience to degradation-induced trap formation.Thus, we have shown that by our machine-learning-driven workflow, hints for possible mechanisms for electrical performance degradation can be obtained in a high-throughput workflow involving only inexpensive and fast experiments.The method is quite generally applicable to different target properties (such as J sc loss or transmissive properties), limited only by the choice of the optical probes, in the sense that phase specific morphologies can only be detected if the corresponding optical probes do not show complete spectral overlap.

Figure 1 .
Figure 1.A) Chemical structures of materials and solvents and B) the device structure employed in this study.C) J-V characteristics and D) EQE spectra of devices based on binary PM6:DTY6 from CF and o-XY, as well as ternary PM6:DTY6:[70]PCBM from CF and o-XY, respectively.Photostability of devices based on binary PM6:DTY6 from CF and o-XY, as well as ternary PM6:DTY6:[70]PCBM from o-XY, respectively.E-H) Normalized recorded current-voltage curves with normalized initial short-circuit current density.I) Normalized photovoltaic parameters of devices based on three different active layers as a function of illumination time.Device stability data are collected over seven devices.Note that the V oc curves of CF-and o-XY-processed binary devices almost overlap.

Figure 2 .
Figure 2. Machine-learning workflow integrated with high-throughput experimentation and characterization and drift-diffusion simulations, with the goal to find quantitative structure-property relationships.

Figure 3 .
Figure 3. Gaussian process regression (GPR) with embedded feature extraction by minimum redundancy maximum relevance (mRMR).A) Graph display of mRMR-GPR analysis.The arrows point from predictor to target with contribution to explanation of variance annotated.B) Partial dependence plot to predict r FF,200hrs from the most relevant feature h A ; the colors refer to the different blends (see inset).C,D) GPR to predict r FF,200hrs from A tot and the [70]PCBM weight fraction X 70PCBM or the solvent BP, T B , respectively.E) GPR to predict c A from the [70]PCBM weight fraction X 70PCBM and the solvent BP, T B .

Figure 4 .
Figure 4. A) FF loss as a function of the active layer thickness, obtained from the total absorption (see Section H, Supporting Information) and the order parameter 1/h A .B) Drift-diffusion simulation of FF loss as a function of the extraction mobility μ n and the active layer thickness L. For other parameters see Section I (Supporting Information).Symbols: simulation; hypersurface: surrogate model by GPR.C) Bayesian Network representation of the probability distributions used to predict FF loss from experimentally available quantities (blue arrows) and of the deterministic relationships assumed in the drift-diffusion model for degradation (black arrows).