Enhancing quantum efficiency of thin-film silicon solar cells by Pareto optimality

We present a composite design methodology for the simulation and optimization of the solar cell performance. Our method is based on the synergy of different computational techniques and it is especially designed for the thin-film cell technology. In particular, we aim to efficiently simulate light trapping and plasmonic effects to enhance the light harvesting of the cell. The methodology is based on the sequential application of a hierarchy of approaches: (a) full Maxwell simulations are applied to derive the photon’s scattering probability in systems presenting textured interfaces; (b) calibrated Photonic Monte Carlo is used in junction with the scattering matrices method to evaluate coherent and scattered photon absorption in the full cell architectures; (c) the results of these advanced optical simulations are used as the pair generation terms in model implemented in an effective Technology Computer Aided Design tool for the derivation of the cell performance; (d) the models are investigated by qualitative and quantitativesensitivity analysis algorithms, to evaluate the importance of the design parameters considered on the models output and to get a first order descriptions of the objective space; (e) sensitivity analysis results are used to guide and simplify the optimization of the model achieved through both Single Objective Optimization (in order to fully maximize devices efficiency) and Multi Objective Optimization (in order to balance efficiency and cost); (f) Local, Global and “Glocal” robustness of optimal solutions found by the optimization algorithms are statistically evaluated; (g) data-based Identifiability Analysis is used to study the relationship between parameters. The results obtained show a noteworthy improvement with respect to the quantum efficiency of the reference cell demonstrating that the methodology presented is suitable for effective optimization of solar cell devices.


Introduction
In the last decades a lot of effort has been devoted to find optimal solar cells design, creating a very vast research field [1][2][3]. Three main kinds of solar cell devices had captured most of the attention: (I) photovoltaic (PV), which convert energy transported in light directly into electrical energy using devices based on semiconductor materials (this device will be analysed in this works), (II) thermophotovoltaic (TPV) [4], which convert heat into electricity by radiating photons that are then converted into electron-hole pairs via a photovoltaic medium, (III) nanophotonic thermophotovoltaic, promising devices recently developed in [5] which combines the best features of PV and TPV. However, optimization and analysis of solar cells have revealed to be a hard task, mainly because of the high number of (microscopic and conflicting) parameters to be set and the difficulties on finding a computationally-competitive accurate model.
In this work we present a combination of numerical algorithms especially designed for thin-film cell analysis and optimizations. Precisely, our approach is composed of four distinct steps: Sensitivity Analysis (SA) of model parameters, in order to study the effect that each parameter has on the model; Single Objective optimization (SOO) and Multi Objective Optimization (MOO) applied to the most sensitive parameters (we call these "SA + SOO" and "SA + MOO"), to find designs maximizing the efficiency of the cell along with trade-offs, balancing the efficiency of the cell and its production costs; robustness analysis, local, global and glocal, of the best designs found on the optimization cycle, to fully evaluate a device resistance against perturbation (e.g. usage, production imprecisions etc.); identifiability analysis, used to investigate the model for parameters functional relationships. We demonstrate the capabilities of our approach for the optimization and analysis of Tandem Thin-Film Silicon devices, using different Transparent Conductive Oxide (TCO) and Back Reflector (BR) materials. We show that we are able to obtain remarkable efficiency-cost improvements, up to 6.71% with respect to the reference cell considered.
The choice of tandem thin-film silicon solar cell in this study, is justified by the recent amount of work that has been spent in optimizing the efficiency of these promising devices. In fact the second generation of thin-film cells has recently arose as a valuable alternative to the more expensive devices, previously built, made of thick polycrystalline silicon wafers. However, the reduced thickness of the absorbing layer, as it will be demonstrated in the next sections, inevitably leads to a reduction in the device efficiency. Hence, several light trapping techniques have been proposed and introduced, during the years, in thin-film cells in order to balance such an effect. Examples are the introduction of a Back Reflector metals layer and, mainly, the chemical introduction of random nano-texturing at the interfaces separating the cell layers [6,7]. The former has been exhaustively investigated in this work. Indeed, a part of our study concerns the explorations of different doping dosage (four dosage level have been considered), together with two level of back reflector smoothness. We found that the best materials to be used, compared to the ones here analysed, are the ZnO for the T C O layer, and smooth Ag for the back reflector layer. However different doping dosage of ZnO should be applied to obtain different trade-offs balancing efficiency and cost of the cell.
The presence of random nano-texturing, whereas, lead to the computational intractability of this devices using Maxwell solvers [8], thus a lot of effort has been recently spent in the development of accurate, however, computational tractable models. The one we have implemented is based on a well balanced combination and generalization of the different approaches presented in [9,10].
Finally we perform a comparison among the results obtained by using the methodological approach here presented and by applying OptIA [11,12] , an immune system based optimization algorithm for single and multi objective optimization.
The remainder of the paper is organized as it follows. In the second section we explain the details of the algorithms composing our methodology. The optical model implemented is discussed along with the experimental results obtained using different materials in the third section. The simulation results obtained with our methodology are presented and commented in the fourth section. Finally, in the last section we present our conclusions and possible remarks for future works.

Building blocks and literature review
In this section we introduce the problem addressed in this paper, the algorithms we used and the main reasons of our choices. The general form of a constrained global optimization problem is: . . , f n (x)) for some scalar functions f i (x), i = 1, 2, . . . , n (assuming that F is to be minimized is not restrictive), k is the domain space dimension, n is the number of objective functions (if n = 1 we talk about Single Objective Optimization (SOO), otherwise Multi-Objective Optimization (MOO)), and T is a subset of R k . A point x = (x 1 , . . . , x k ) ∈ R k is said to be a feasible (respectively unfeasible) point with respect to problem (1), if x ∈ T (respectively x / ∈ T ). The presence of constraints in an optimization problem may greatly increase its complexity [13].
In Fig. 1 we summarise the composite design methodology presented in this work. Initially the model is investigated by means of Sensitivity Analysis, that evaluate the magnitude of the effect that each parameter has on the solar cell efficiency, and possibly remove from the model input parameters that have only weak influence on the latter. The (reduced size) model is hence fed into black-box Single and Multi Objective Optimization solvers, respectively for full optimization of optical efficiency of the solar cell, and for trade-off analysis of efficiency Fig. 1 Flow chart of the composite design methodology applied for the analysis and optimization of the tandem thin-film silicon solar cells. Relative computational costs are also provided for each step and production costs. Finally optimal and Pareto-optimal designs are evaluated by means of Robustness Analysis, and the Identifiability Analysis is used to invert the relationship between model parameters and Pareto optimality.

Single objective Optimization
Almost all non-trivial engineering designing problems concern the optimization of an objective function characterized by a complex domain-codomain behaviour, for which traditionally technique cannot be applied [14][15][16]. Clearly there is currently no algorithm, which performs better than the others in every possible engineering problem [17], justifying the fact that new algorithms are continuously designed. In this work we use the Multilevel Coordinate Search (MCS) algorithm for solving the single objective optimization problem associated to thin-film silicon solar cell efficiency. However, in some particular cases MCS is strongly competitive and outperforms some state-of-art algorithms [18]. We here give just a quick explanation of the working principle of the method, refer to [19] for details.
The MCS algorithm is a combination of heuristics and calculus approximations; it is composed of two main sub-routines: a global routine and a local one. The global part is obtained via a branch-based algorithm. The domain is repeatedly split, along a coordinate direction, into, usually, three sub-boxes; each one of them containing a base point (a promising point for the optimization problem) and an opposite point, used to determine the sub-box. The splitting process is iterated up to an user-defined-number of times. The local search is whereas obtained via quadratic approximations of the function in the neighbourhood of some optimal points previously selected. Interpolating points used in the quadratic models are selected using a combination of triple and coordinate search.

Multi-objective optimization
SOO may reveal to be restrictive for the majority of real world problems; often mathematical models present more than one objective function to be optimized, usually those being in contrast with each other (e.g., generally the more efficient a device is the more expensive it is). We thus introduce the well known definition of Pareto optimality, and the partial order deriving from it. A point x ∈ T of the domain space, is said to dominate a point y ∈ T , with respect to problem (1) The Pareto-Front is then defined as the total set of non-dominated feasible points.
Evolutionary algorithms (EA) are a class of meta-heuristic methods that operate on a set of candidate solutions and subsequently, using biologically-derived evolution techniques (e.g. mutation, crossover), modify it. Because of their intrinsic nature, this class of algorithms may be easily adapted to face MOO problems [20]. Contrarily, exact algorithms for most of real world problems, as thin-film solar cell optimization process, may reveal to be unsuitable due to a high dimensionality of the problem considered. Indeed, in order to avoid an increase of the computational complexity in the optimization process, we do not consider a robust MOO [21][22][23] as a core of our optimization but, instead, we tackle this problem using a post-processing robustness analysis techniques (see 2.2.2). Among a vast literature on MOO algorithms [24][25][26][27], in this work we have used the NSGA-II algorithm, which we will briefly explain in the following.
NSGA-II [28] is a multi-objective evolutionary algorithm based on non-domination, elitism and crowding distances. At each iteration the algorithm assigns two attributes to each point (candidate solution) of the current population, namely rank and crowding distance, as it follows. The fast non-dominated sort algorithm implemented in NSGA-II, sorts the points of the current population P in different fronts F 1 , F 2 . . . , F s , such that F 1 contains all non-dominated points of the current population, F 2 contains the non-dominated points of the current population minus those ones which are in F 1 , and so on. We then say that a point x has rank p if x ∈ F p . So the rank is the measure of a point quality (the less the better). The crowding distance of x, estimates how isolated x is in the domain space. Hence the crowding distance of a point x is high if F(x) lies in unexplored regions of the objective space (thus the higher the better). The points are thus sorted according to the order induced by rank and crowding distance. The crowding distance of x, is an estimation of the density of points near x. It is calculated by selecting, for each objective, the two closest points on each side of x, calculating their distance and finally summing them over all the objective functions. Namely it is defined Front-wise as it follows: let x j , j = 1, . . . , l be all the points in the pth front F p . Let f i , i = 1, . . . , n, be the ith objective function; assume that the x j s are indexed accordingly to their ith objective function values, define: Finally the crowding distance of x is defined as: c dist (x) = n i=1 m i (x).

Sensitivity analysis
Sensitivity analysis (SA) may be loosely defined as the study of how output uncertainty of a model may be justified in terms of uncertainty in the model input parameters [29]. We use SA both to give an accurate mathematical description of the model parameters considered in the optimization problems and to reduce the domain space dimension, in order to obtain a computational reduction. In order to accomplish this task we implement two different SA algorithms, namely the Morris algorithm and the Sobol's indexes based approach. Methods as the Morris one are usually referred to in literature as screening methods or qualitative methods, because they tend to give a qualitative description of the parameters influence rather than a quantitative [30,31], indeed ranking the factors on an interval scale (i.e. it is a one-stepat-a-time method). The Sobol's indexes approach is whereas a Variance-based quantitative SA method, which indeed gives quantitative information about the output variance portion that may be explained by each parameters or even by groups of parameters [29]. The cost of this information is, of course, a significantly higher computational burden. It follows a review of the functioning principles of both methods. The basic idea behind the Morris method is that, given the model function and the ith parameter of the domain, low values of ∂ F ∂ x i may be interpreted as low influence of the ith parameter on the output. Hence the method consists of random sampling of the domain T (assumed to be a k-dimensional equally spaced grid) and approximation of the effect that every parameter has. Namely, given x = (x 1 , . . . , x k ) on the domain and a Δ ∈ R, multiple of the grid step, such that x +e i Δ = (x 1 , x 2 , . . . , x i +Δ, . . . , x k ) still belongs to the domain, where e i is a vector of zeros but with a unit as its ith component, the algorithm computes the elementary effect: As x spans the domain we obtain the set: The ith parameter overall influence is thus estimated considering the mean μ * i and standard deviation σ i of the absolute values of EE i elements [32]. Specifically, high values of mean stands for high influence of the ith parameter on the output, high values of standard deviation stands for strongly non-linear effects and/or correlation with other parameters.
The computation of Sobol's indexes [33] is based on a high dimensional model representation of the objective function F (known as ANOVA decomposition). It can be proved [34] that: where V is the total variance. The right-side terms of this equation represent the contribution of each and group of parameters on the total variance. The first order sensitivity index for the ith parameter is thus defined as S i = V i V . It gives information on the contribution that each factor (taken alone) has to the output. Whereas the total effect Sobol's indexes S T i is defined as where V −i is the variance computed on all parameters except x i , which is supposed to be fixed to its true (nominal) value [31].
To efficiently approximate the Sobol's indexes we implemented a Monte Carlo approach [29]. The complete procedure is illustrated in Algorithm 1.
where k is the domain space dimension and N is the number of samples considered) be two input sample matrices randomly generated in the domain space and let . , x k ) be the model output (i.e. in our case a n-dimensional vector), the firstorder indexes S i and the total effects S T i can be respectively computed as: where y A , y B and y C i are the evaluation of objective function on the matrices of data A, B and on C i , and f 0 is the mean, that can be either estimated by using sample A or B, defined as (e.g. using sample A) Algorithm 1 Sobol's Indexes Computation [29] procedure EstimateIndex(N)

Robustness analysis
In this section we introduce the concept of robust design and the three different indexes we have used to evaluate thin-film devices robustness. It is indeed clear that the quality of a solar cell cannot be evaluated only considering its Pareto-optimality with respect to the objective functions. We argue that an important factor that has to be taken into account, especially in such a design optimization problem (i.e. solar cell), is the robustness of the device. Intuitively, a device (more generally a point of the domain) is said to be robust with respect to problem (1) if small changes on the design parameters correspond to small changes on the objective function computed on the "perturbed design" [35][36][37]. In order to precise this intuitive definition of robustness we introduced three robustness indexes, which will reveal to be essential for the Decision Making.
The Local robustness (LR)and Global robustness (GR) indexes are based on the idea that perturbations, interpreted to be production errors or usage, may be both assumed to be Gaussian distributed. We set to zero the distribution mean, and the standard deviation of the distribution to 1% of the corresponding parameter value. The computation of the Local index is as it follows: for i = 1, . . . , k, the ith Local robustness index of x (L R i (x)) is computed via a Monte Carlo approach by perturbing N times the ith parameter of x (say x + e i Δ the perturbed point) using Gaussian noise (set as explained above), and then computing the fraction of perturbed points for which F(x + e i Δ), where e i is a k-vector of zeros but with a unit as its ith component, does not deviate more than 1% from F(x). In other words: Analogously the Global robustness index of x (G R(x)) is computed by perturbing N times all the parameters of x with Gaussian noise and calculating the fraction of perturbed points (say x + Δ every each perturbation) such that F(x + Δ) do not deviate more than 1% from F(x); in short: The Local robustness index is computed only considering perturbation along coordinate directions, it does work well as a screening of a point robustness, and quantifies the robustness of each parameter of the point. The Global robustness index quantifies, on the other hand, the point resistance against random perturbations. However these methods work with a fixed perturbation strength (σ ) that bounds the region in which perturbed points are explored.
The "Glocal" robustness index (proposed in [38]) is a hybrid method for the evaluation of a point robustness. Given an initial point, x ∈ T , it iteratively tries to find the biggest viable region (which is a region in which F(x + Δ) does not deviate too much from F(x)) in a neighbourhood of x; then a Monte-Carlo integration is used to approximate the volume of the Viable region, and finally the volume is normalized in order to get a robustness index which is independent on the number of parameters composing the model. Hence the "Glocal" index depends on the biggest "stable" region around x. An obvious drawback of this method is the tremendous computational burden. In order to reduce the computational effort in the Viable region search, the latter is guided by a Principal Component Analysis (PCA), which is a standard statistical tool for non-parametrically extracting relevant information from complex data-sets (see for example [39]).

Identifiability analysis
The last step of the methodology we present is the Identifiability Analysis (IA). In the most general form (see [40] for further details on the definition), consider a system governed by a set of (possibly differential) equation. Let z be the state variable of the model, with z 0 the initial value for z. Let t be the time variable and x be the parameter of the model (i.e. the design parameter of the solar cell in our case). Identifiability analysis is applied to a model of the form where ϕ is a function of the state variable x, its derivatives, a k-dimensional parameters vector x, and time t. Equation (3) enforces the initial condition for the model considered. y, the observation variable, is a function of z, t and x. A parameter x i is said to be identifiable (with respect to the model considered) if there exist, for given values of z 0 and y , a unique solution for x i from equations (2)-(4); otherwise it is said to be (structural) non-identifiable. Non identifiability manifests itself through functional relations between the non-identifiable parameters, thus most of the IA algorithm analyse the model looking for those. In this work we use the MOTA (mean optimal transformation approach) algorithm, a databased Identifiability Analysis algorithm, which are a class of a-posteriori model-independent method, hence perfectly suited for our application. The MOTA [40] algorithm is based on consecutive iterations of the ACE (alternating conditional expectation) algorithm [41], which estimates optimal transformations that maximize the linear correlation between parameters.
Briefly by choosing a parameter as response, say x j , and the others as predictors, assuming that there exist transformations such that being Gaussian noise, the ACE algorithm estimates, via an iterative process, optimal transformations ψ and φ j such that: However the ACE algorithm itself cannot be used for IA: if, say, φ i (x i ) = 0 (i.e. x i is not correlated to x j ) then it may occur that φ i (x i ) = 0, because x i will be used by the algorithm to justify the noise in Eq. (5). The idea behind the MOTA algorithm is then easy to grasp: running the ACE algorithm more than once, if the ith transformation, φ i (x i ), remains "quite" stable (the authors introduced a precise index for that), then parameter x i is considered to be correlated to x j . Finally parameters that have been found correlated with each other are retained to be non-identifiable.
In the next section we will introduce the optical model we have used for the evaluation of the optical efficiency of tandem thin-film silicon solar cell. As conclusion of this methodology section we present in Algorithm 2 the pseudo-code of our methodology, in the MOO case.

Thin-film silicon solar cell
The structure of a Thin-Film Silicon Solar Cell used as a case of study is shown in Fig. 2. The cell structure is based on a two coupled p-i-n devices, one made by amorphous Si Namely instead of cutting thick wafers of polycrystalline Si and then using them to build the device, the layers are subsequently grown using the Plasma-assisted chemical vapour deposition technique [42]. This makes possible to design devices with a Si layer (more generally the absorber layer) thickness of the μm order, which is 2 order of magnitude less then previously used devices. However, clearly, less thickness means less absorption, so light trapping techniques have been used to improve the photon optical path in the layers. Examples of this are the introduction of a back reflectors layer in the device architecture and the nano-texturing in the cell layer interfaces. The latter is usually obtained using alkaline solutions, as KOH or NaOH, which corrode the silicon forming randomly positioned square pyramids. The depth of those can be accurately tuned controlling the temperature and the time length of the corrosion process (an overview of the recent advances may be found in [43]). The optimization of the cell geometry in term of light harvesting can be performed modifying meso-scale features (thickness of the different layer) and nano-scale features (nano-texturing of the interface). However, the meso-scale features cannot be optimized without considering the processing time (and the related cost), which is a monotonic function of the layer thickness (especially the μc − Si : i one [1]). It is worth noting that there are alternatives to random nano-texturing: the work in [44] demonstrates that periodic nano-texturing may strongly improves the efficiency of the solar device. Understanding the trade-off between the different effects in the cell design is, of course, a complex problem and it needs a suitable numerical approach. We now give a quick overview of the two methods [9,10] we have combined in order to compute devices quantum efficiency.

Electromagnetic fields computation
Let us assume that the electromagnetic radiation is normal to the device. First, we define: -I, R, T the amplitude of incident, reflected and transmitted electric fields.
σ i , i = 1, . . . , 11, as the roughness of the interface between the (i − 1)th and the ith layer (where the 0th layer is the external ambient); d i , i = 1, . . . , 11, the thickness of the ith layer; -N i = n i − ı k i the complex index of refraction of the ith layer; . . , 12, respectively the amplitude of the incident wave on the ith interface travelling in a forward direction, backward direction, left and right.
According to Fresnel's equations for nano-rough interfaces the complex reflection and transmitted coefficients (forward and backward) are given by The second factors in (7)-(8) are due to the roughness of the interfaces. Hence, the relationship between the incident, reflected and transmitted amplitude is where where i = 1, 2, . . . , 11, 12 and β i ≡ 2π n i λ − 2π k i λ ı. Having calculated the amplitudes, the respective intensity may be easily computed. Using (9) we can determine E R F i and E R B i and hence the resultant forward-going wave E F m (x) and the resultant backward-going wave E B m (x). Finally, we can obtain the absorption profile in the ith layer and the corresponding integrated absorption and sum them over all the wavelength to obtain the total absorption profile and relative integrated absorption. Thus, as explained in [10], the scattered light in the ith layer can be approximated by Refer to the original work for details on the computation of α.

Scattered light computation
In order to evaluate the total absorption profile of scattered light we have used the Monte Carlo approach presented in [10]: a certain number of photons are traced until the photon absorption in a layer or the photon reflection to air takes place. Photon possible behaviours depend on his position within the cell: -A photon in a layer may be absorbed with a probability that depends on k i , d i and θ , the latter being the convex angle between the photon direction and the normal direction to the layer. -A photon incident with a rough interface may be (i) transmitted and scattered, (ii) reflected and scattered, (iii) refractively transmitted and specularly reflected; these events probabilities depend on θ i , θ i+1 (the refractive angles that can be calculated using Snell's law), n i and σ i .

Thin-film silicon solar cell problem
As previously explained the parameters we consider are: (i) the roughness of each interface: σ i , i = 1, . . . , 11; (ii) the thickness of the micro-crystalline silicon layer: d μc−Si . The multiobjective optimization problem can be summarized as: (10) in which Q e is the average overall quantum efficiency of the cell in the ideal charge collection conditions (100 % collection efficiency independently on the frequency) calculated as the average absorption in the two intrinsic μc − Si:i and a − Si:i layers, d μc−Si is the thickness of the intrinsic μc − Si layer and d i is the thickness of the ith reference cell material ( Table 1). The minimization of the micro crystalline layer thickness is justified, as already mentioned, by the high fabrication cost of this particular layer. Constraints in (10) are introduced to mitigate the scattering weight in case of film much thinner than the first ZnO layer roughness. In these cases the approximation made by the modelling approach used in this paper does not hold any more; an improved description can be obtained using the effective medium approximation [45] and it will be included in the next version of the algorithm. The result shows that the most important parameters considered are d μc−Si , as it could be easily expected, σ 2 and σ 4 , hence these three are the parameters to be considered during the analysis and the designing. Parameter as σ 5 , σ 8 , σ 9 have a weak effect on the model output

Results
In this section we present step-by-step all the results obtained by our methodology.

Sensitivity analysis
The results of the Morris analysis with a number of trials N = 10,000, plotted in Fig. 3, show a qualitative ranking of the most relevant factors (represented by high values of μ * ) are: (i) the thickness of the intrinsic layer, demonstrating what previously stated; (ii) the interface roughness above the ZnO-layer, which is the first textured materials encountered by an incident photon and shall, thus, have a relevant effect in terms of scattering events; (iii) and the one of the intrinsic amorphous silicon layer, which is the first absorbing layer encountered by an incident photon. It is worth noting that their respective values of σ are high too, i.e. the relationships between these three parameters and the model objectives are strongly non-linear. Parameters σ 5 , σ 8 , σ 9 have a weak effect on the model output (less than 1% of maximum effect). The Morris analysis results are confirmed, extended and quantified by the Sobol's indexes for the first order effects and the total effects, listed in Table 2 along with their relative errors (50% probability). Once again parameters as σ 3 , σ 5 , σ 6 , σ 8 and σ 9 are ranked as non-effective parameters (less than 1% of maximum S T i value) confirming the preliminary hypothesis of the Morris analysis qualitative ranking; whereas d μc−Si , σ 2 and σ 4 are the most influent parameters (the total effect of the other parameters is at least reduced of a 10 −1 factor). Notice however that results for σ 11 and σ 7 are different for the two SA here performed. The latter are ranked as insensitive by the Sobol analysis; whereas the Morris analysis ranks them as the fifth and sixth important parameters for the model. We here perform a conservative choice, selecting as in-sensitive only the parameter for which both analysis agree. We hence discard (i.e. fixing to some predefined values) the non-influent parameters. However in order to demonstrate the power of the "SA + Optimization" approach we ran the optimization algorithms twice: considering the full parametrized model, and the reduced one. We find that computational effort in the "SA + Optimization" is remarkably lower with respect to the full Confirming the preliminary Morris analysis parameters as σ 3 , σ 5 , σ 6 , σ 8 and σ 9 are ranked as unsensitive.
The low values of S T j − S j reveal that a large portion of the variance on the model output uncertainty may be explained in terms of first order effects. Negative signs in the table and S T j of d μc−Si > 1 are due to numerical errors (i.e. negative sensitivity indices close to zero stand for unimportant factors) parametrized model, however a little loss in the optimal point is found in the SOO, whereas a significant improvement is found in the MOO case.

Single objective optimization
The termination criterion we used in the single objective optimizations is either an upper bound on functions evaluations (50k 2 , k being the domain space dimension) or the stagnation of the algorithm. Notice that in our simulations the algorithm always terminates for the stagnation of the solution. The efficiency gain of the best design found, compared to the reference cell one, is definitely remarkable holding a relative improvement in Q e about 6.71%. The best total absorption profile obtained is plotted in Fig. 4. The improvement, compared to the reference cell, is of 5.88% (Quantum efficiency and thickness of the intrinsic layer are 0.604169 and 2210). However, notice that almost all parameter values are maximized (with respect to constraints in (10) = 1305). Note that the design minimizing the micro-crystalline intrinsic layer thickness (30% less than the reference cell) still holds a 1.27% absorption improvement and a relative higher Q e . Finally the results obtained trough the "SA + MCS" reach a 5.55% total absorption improvement. Therefore, using this type of optimizations we have a little loss on the performance with respect to the "full optimization" MCS, yet we want to underline that the computational burden of the "SA + SOO" optimizations is significantly reduced: 24.3% on average.

Multi objective optimization
The Pareto-Front approximation obtained applying NSGA-II to the multi-objective, fullparametrized, problem is plotted in Fig. 5a (total number of generations is 700, population size is 100). Q e of non-dominated designs, found by the algorithm, spans from ≈ 0.567 to ≈ 0.6, whereas d μc−Si spans from 1190 to ≈ 2210 (which is the maximum range allowed by the constrains). Thus we have obtained a wide spread Pareto-front approximation. We may notice that the designs obtained using MOO are definitely more robust than the previously obtained. The total absorption profile of the cell with the maximum value of quantum efficiency (Q e = 0.6000, d μc−Si = 2209) obtained in these runs is plotted, along with the reference cell one, in Fig. 6: the improvement, which mainly lies in the [700, 900] wavelength region, is of 5.083%. The best absorption improvement obtained in these runs is 5.17%.
The Pareto-front approximation obtained in the "SA + MOO" optimization (using NSGA-II set as above), along with the feasible designs explored, is plotted in Fig. 5b, whereas the Pareto Pareto-front approximations obtained in the "MOO" and in the "SA + MOO" approach are compared in Fig. 5c . The improvement, with respect to the "full-parameterized" Pareto-front approximation, is remarkable. Figure 5c clearly shows how the "SA + MOO" Pareto-front dominates the one obtained using the simple MOO approach. This significant improvement may be due to the particular way in which EAs explore the objective space; a reduced domain space seems to significantly enhance the search. This result demonstrates the consistency of our combined algorithmic approach, which in complex problem (such as the one associated to thin-film cell design) shall enhance the capability of a heuristic optimization approach. In particular, SA reveals to be a fundamental tool as a pre-processing phase to investigate parameter effects on the model. The convergence rate of the algorithm is graphically analysed in Fig. fig:convergence, where we plot the non-dominated fronts computed during the iterations of NSGA-II (Fig. 7).
The Pareto-orientation of the optimization allows us to perform a parametric analysis of the absorption profile of optimal trade-offs. Figure 8 shows the absorption profile of Pareto- optimal solutions as a function of the quantum efficiency. This plot demonstrates that most of the gain in the optical efficiency of Pareto-optimal cell is due to improved absorption capabilities in long wavelengths. Indeed, the effect is reflected by a shift in the figure, which is appreciable even for quantum efficiency variation of the order of 10 −2 .  Fig. 9 depicts the absorption profile associated to the closest-to-ideal design found by OptIA for this optimization problem. This has a Q e of: 0.5719 and a d μc−Si of 1190, i.e. a variation of − 2.29 and − 25.63% ,respectively, compared to the closest-to-ideal cell design found by NSGA-II.

Robustness analysis of notable designs
The Local and Global robustness indexes of the trade-offs found by NSGA-II are listed in Table 3. Even if the Local indexes are very high, the Global Robustness indexes do not even reach 70%. Thus, the designs analysed are very resistant to perturbations along coordinate directions, but become unrobust (fragile) if the perturbations are combined. This means that, at least locally for the points analysed, the objective functions of the model are strongly non-linear, depending mostly on parameters combinations, rather then disjoint effects due to single parameters. The only not maximized Local indexes are the ones associated to d μc−Si , σ 2 and σ 4 . This locally confirms the global results given by the SA. In Fig. 10 are highlighted the knee points of the Pareto front, along with their respective Global robustness indexes. The highest index found so far is the one of Trade-off 3, namely 67.37%; which makes it the most robust design analysed with respect to the optical model. Its "Glocal" index is 1.8035, which, compared to the total volume of the feasible region, indicates a relatively big Viable region, confirming and improving our hypothesis on the robustness of this device.

Identifiability analysis
In Table 4 we list the results of MOTA applied to 10 −5 non-dominated points. The important relationships found using these points are between p4 and p5 and between p7 and p9, plotted in Fig. 11a, b respectively. Notice that these hidden functional relationships cause complex behaviour of the model output with respect to the domain space, and thus making harder the task of an optimization algorithm. However, all of them do not show up in the "SA + MOO"     11 Scatter plots of the optimal transformations F found by the MOTA algorithm, using the 10 −5 non dominated points as input data set, for p4 and p5 and for p7 and p9 respectively. a p4 and p5. b p7 and p9 Optimization (at least one parameter involved in these relationships is considered fixed in the latter), and this may be one of the reason for the improved behaviour of the optimization algorithms.

TCO and back reflector materials
In the previous sections we have discussed about the results of our optimization and the analysis of the most common structure of tandem solar cells: highly doped ZnO used as TCO and smooth Ag used as a BR. Of course this is not the only possible choice, there is a huge amount of different possible materials (having price comparable with Si) and different doping dosage that may be used in solar cell design. In this section we will discuss the results that we have obtained exploring some alternatives: (i) varying the dope dosage of the ZnO layers in five conditions: non-doped, low doped (resistivity larger than 100 mΩ× cm), normally Each simulation is performed with a different ZnO, used for the TCO layer, doping dosage (optimal, normal, low, not) and Ag ,the back reflector used, roughness (smooth, rough) combination  smooth Ag combination, is the device analysed in Sect. 3.2. We applied NSGA-II algorithm setting the population size to 100 and the maximum number of generation to 200, in order to find which of these devices better perform. The eight Pareto-Front approximations obtained are showed in Fig. 12. The results clearly show that, for low efficiency, the best combination of material is given by normal doped ZnO and smooth Ag. For high efficiency the devices seems to perform in a very similar way. In Tables 5 and 6 we report comparison among the results obtained by using NSGA-II and OptIA in term of optimality of the closest-to-ideal trade off for the eighth optimization performed in this section (Table 7).

Conclusions
In this work we have given a complete analysis of a tandem thin-film Silicon solar cell model. The combination between the matrix method and the calibrated photonic Monte Carlo simulations allowed us to have an accurate fairly-quick optical model that could thus be analysed.
The results obtained in all the different optimizations performed, not only widely demonstrate the applicability of our methodology but also show the importance of every step: the understanding of a high dimensional model would have not been possible without a preliminary sensitivity analysis; MCS algorithm was used both to find extremely efficient designs, up to a 6.71% Quantum Efficiency improvement and a 5.88% gain on the total absorption profile, and fair efficient designs, up to a 3.70% gain on the total absorption profile, considering cost-savings bounds; NSGA-II was used to find Pareto-Optimal points of the problem: Maximize Quantum Efficiency, Minimize Intrinsic Layer Thickness; hence balancing Efficiency and cost-savings. An interesting fact about the points found in NSGA-II runs is the higher Global Robustness indexes they hold compared to SOO results, even exceeding 80%. Hence MCS results behave better in terms of efficiency, while NSGA-II results shows a better behaviour in terms of stability. Of great importance was the improved behaviour of NSGA-II algorithm when used in combination with the SA results. In fact "SA + MOO" NSGA-II optimization, performed significantly better then the full-parameterized optimization in terms of both optimality of the results obtained and reduced computational effort. Finally Identifiability Analysis gave us higher order information on the domain parameters and their behaviours with respect to the model outputs.
It was also demonstrated that our approach could be systematically applied on designing thin-film solar cell devices, ensuring remarkable improvements in terms of efficiency, cost savings and robustness. In particular, the robustness designs found in our optimization process might be useful for a decision maker and should be take in consideration in a design process. Future works will concern the exploration of alternative materials and solar cell types. An interesting direction of work lies in experimental validations of the optical efficiency of several Pareto-optimal designs obtained by our methodology, which could lead to new possible prototypes in collaboration with industrial partners.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.