SolarDesign: An online photovoltaic device simulation and design platform^††thanks: Project supported by the Scientific Research Project of China Three Gorges Corporation (Grant No. 202203092).

To calculate physical properties of photovoltaic materials, first-principles calculations are performed within the framework of density-functional theory (DFT) [10] using the plane-wave pseudopotential method as implemented in the Vienna Ab initio Simulation Package [11, 12]. The electron-ion interactions were described by using the projected augmented wave pseudopotentials [13]. We adopted the generalized gradient approximation formulated by Perdew, Burke, and Ernzerhof as exchange-correlation functional [14]. The Heyd-Scuseria-Ernzerhof functional was then used for the bandgap correction [15, 16]. Structure optimization was achieved using the conjugate gradient technique until energies converged to $10^{-4}$ eV. For organic semiconductor materials, we constructed structurally ordered molecular crystals, taking periodic boundary conditions into account during the calculations. The hybrid perovskite materials with specific alloy compositions were simulated by constructing supercells, and for the chemical compositions that cannot be directly simulated, we used Vegard’s law to fit the properties. To properly take into account the long-range van der Waals (vdW) interaction that is non-negligible for hybrid perovskites involving organic molecules, the optB86b-vdW functional was adopted [17]. The process of DFT-based calculations and analysis of the calculated results were carried out using our in-house developed codes, the Jilin artificial intelligence-aided material design integrated package [18, 19].

The calculated physical properties include: bandgap, work function, electron affinity, dielectric constant, electronic band structure, density of states, effective mass, exciton binding energy, optical absorption spectrum, complex refractive index, and optical transition intensity. Most of these properties can be obtained through electronic structure analysis. Optical properties can be derived from the frequency-dependent dielectric function [20], exciton binding energy is calculated using the hydrogenic model [21], and optical transition intensity is proportional to the square of the transition dipole moments.

To calculate optical (absorption) properties of photovoltaic devices, mode-matching methods [22] are applied. Electromagnetic fields at each layer of the multi-layer device structure are expanded as plane waves, and Maxwell’s equations are solved using the continuity conditions for the tangential components of the electric and magnetic fields. Furthermore, the generation rate [23] is calculated based on the electric field distribution in the active layer, the corresponding complex refractive index, and the incident light spectrum. Additionally, for textured surfaces, we use the $4n^{2}$ (Yablonovitch) absorption limit [24, 25] and set the texture factor ranging from 0 to 1, with 1 indicating the $4n^{2}$ absorption limit.

For general models, non-uniform spatial discretization (Scharfetter-Gummel scheme), semi-implicit time stepping, and Gummel iteration [26, 27] are employed to solve the drift-diffusion model, which includes the Poisson’s equation and the current continuity equations for electrons and holes.

For specialized models, different methods should be adopted for different types of solar cells. Regarding quantum tunneling, we employ the non-local band-to-band tunneling models and WKB (Wentzel-Kramers-Brillouin) approximation, seamlessly integrated with the drift-diffusion model [1]. Regarding exciton dynamics, we use the Langevin bimolecular recombination model and a mixed dissociation model, which combines local diffusion-dissociation (Onsager-Braun model) with direct delocalization [2]. Regarding ion dynamics, we apply the ion migration and accumulation model, involving continuity equations for positive and negative ion flows (similar to carrier current continuity equations, including drift and diffusion terms, but excluding generation and recombination terms), initial concentration conditions, and annihilation boundary conditions [3].

Based on detailed balance theory, we compute radiative recombination using blackbody radiation spectrum and decompose non-radiative recombination into bulk/surface Shockley-Read-Hall (SRH) recombination and Auger recombination. Meanwhile, Ohmic loss from electrodes, low-mobility transport layers, and interfaces between the transport and active layers is represented by series resistance; Ohmic loss from shunt current caused by defects and pinholes is represented by shunt resistance. All of the above lead to the development of a new circuit model for intelligent quantification of solar cell efficiency losses [28, 29, 30].

Additionally, the modified nodal analysis (MNA) and high-dimensional Newton’s method are used to solve circuit models for both single solar cell and large-scale module configurations.

The input parameters include device structure (thickness of each layer), incident light spectrum, and complex refractive indices of all materials. Among them, the complex refractive index can be obtained from the Material Library, and the light spectrum can be the Solar AM 1.5G spectrum or user-defined spectrum.

The output results include optical absorptance and reflectance, carrier (exciton) generation rate, and maximum short-circuit current.

The input parameters of the non-electrode layers include relative dielectric constant, electron affinity, bandgap, electron/hole mobility, minority carrier lifetime of bulk/surface SRH recombination, Auger recombination coefficient, effective density of states in the conduction/valence band, and donor/acceptor doping concentration. Moreover, users can choose to import the generation rate from the optical output results or input the averaged generation rate themselves. The input parameters of the electrode layers include work function and contact type (Schottky/Ohmic). To simulate the transient current density-voltage curve, it is necessary to provide the voltage scanning range. To simulate the steady-state operating point, a given operating voltage is required. Among them, the relative dielectric constant, electron affinity, bandgap, and effective density of states can be obtained from the Material Library.

The output results of transient simulation include current density-voltage curve and characteristic parameters (short-circuit current, open-circuit voltage, fill factor, power conversion efficiency). The output results of steady-state simulation include electron/hole current, electron/hole concentration, electric potential, non-radiative recombination, and energy band diagram.

The additional input parameters for perovskite materials include initial concentration of positive/negative ions, mobility of positive/negative ions, and voltage scanning rate. The transient output results include current density-voltage hysteresis curves (forward-reverse scan) and corresponding characteristic parameters. The steady-state output results include the net ion concentration distribution in addition to the general output.

The additional input parameters for organic materials include the delocalization ratio of excitons, exciton lifetime, exciton mobility, exciton radius, and bimolecular recombination coefficient. The steady-state output results include the exciton concentration distribution and exciton dissociation probability in addition to the general output.

The additional input parameters for tunneling layers include the electron/hole tunneling coefficient (relative effective mass).

Whether for traditional or new circuit models, the input parameters are the experimentally measured current density-voltage curves, and the output results are the numerically fitted current density-voltage curves and fitting errors (root mean square error and mean absolute error). Furthermore, the output results of traditional circuit model include five parameters of the model (short-circuit current, saturated dark current, ideality factor, series/shunt resistances). The output results of the new circuit model also include five parameters of the model (saturated dark currents of bulk SRH, surface SRH, and Auger recombination, series/shunt resistances) and the quantified proportion of device efficiency loss (bulk SRH recombination, surface SRH recombination, Auger recombination, series/shunt resistances).

Device structure of a tunnel oxide passivated contact (TOPCon) silicon solar cell is set as: Ag/SiNx [70 nm]-c-Si ($\rm{p}^{+}$) [300 nm]-c-Si (n) [200 $\rm\mu$m]-SiO2 [1 nm]-c-Si ($\rm{n}^{++}$) [30 nm]-Ag [1000 nm]. Compared to Silvaco, the calculation error of current density-voltage (J-V) characteristics by SolarDesign is lower than 1 % (See Fig. 1). Particularly, the computer time and memory cost are reduced by more than 10 times. Moreover, the short-circuit current (Jsc), open-circuit voltage (Voc), fill factor (FF), and power conversion efficiency (PCE) of the TOPCon cell are $42.69\,\,\mathrm{mA/cm^{2}}$, 0.71 V, 0.85, and 25.64%, respectively. The Shockley–Queisser limit with 33 % efficiency could be approached by both near-perfect Sunlight absorption above the bandgap of silicon and ignorable non-radiative recombination in the device. Additionally, the continuity of classical total current (that is the sum of the drift and diffusion currents) can be satisfied at all device regions except for the tunnel oxide layer, where quantum tunneling occurs (See Fig. 2).

Device structure of an organic solar cell is set as: ITO [70 nm]-ZnO [30 nm]-P3HT:PCBM [200 nm]-PEDOT [50 nm]-Au [100 nm]. Table 3 lists characteristic parameters of the organic cell with different exciton delocalization ratios. As illustrated in Fig. 3, higher Jsc is expected with a larger delocalization ratio of exciton. Moreover, due to reduced built-in potential, exciton dissociation probability decreases when the organic cell is operated from the short-circuit to open-circuit condition (see the inset of Fig. 3).

Device structure of a perovskite solar cell is set as: ITO [50 nm]-Ni$\rm{O}_{X}$ [20 nm]-MAPb$\rm{I}_{3}$ [500 nm]-$\rm{C}_{60}$ [50 nm]-Ag [120 nm]. The Jsc, Voc, FF, and PCE of the perovskite cell are $21.68\,\,\mathrm{mA/cm^{2}}$, 0.99 V, 0.81, and 17.38%, respectively. Compared to COMSOL, the calculation error of J-V characteristics by SolarDesign is lower than 1% (See Fig. 4). Particularly, the computer time is reduced by more than 100 times. Figure 5 shows energy band diagram of the perovskite cell under the open-circuit condition. The quasi-Fermi levels of electron and hole remain constant, which indicates zero currents. Furthermore, hysteresis behaviors are simulated under a series of voltage scanning rates as depicted in Fig. 6. The J-V characteristics are traced in forward (from short-circuit to open-circuit) and reverse (from open-circuit to short-circuit) scans. The hysteresis effect can be ignored at very fast or very slow scanning rate. The significant difference of the Jsc between the forward and reverse scans indicates dominated bulk recombination at the perovskite active layer. Additionally, the ion migration results by SolarDesign have been benchmarked against those by IonMonger [3, 4].

The efficiency losses of a perovskite solar cell incorporating star-shaped polymer [31] are quantified by a modified diode model [30] in the SolarDesign platform. The averaged fitting errors for the measured J-V characteristics are only $\sim 0.1$ $\rm{mA}/\rm{cm}^{2}$ for both the control and polymer incorporated cells, where hysteresis effects are ignorable and thus only forward scanning results are analyzed. Figure 7 shows that the incorporation of polymer reduces the bulk recombination and increases both the series and shunt resistances. The dominated surface recombination in the polymer incorporated cell results in a very high FF of 0.86.