Analysis of mesh density in electrical simulations of heterojunction silicon solar cells Jošt Balent, Janez Krč and Marko Topič University of Ljubljana, Faculty of Electrical Engineering, Tržaška 25, 1000 Ljubljana E-mail:jost.balent@fe. uni-lj.si Abstract. Solar cells convert solar energy into electricity. Numerical modelling has been proven as an efficient tool in the process of optimisation of solar cell performance. In this contribution we focus on electrical modelling of a silicon heterojunction solar cell, employing the Sentaurus TCAD simulation tool. In particular we studied the impact of an adaptive mesh density on accuracy of the J-V characteristics and the main output parameters of the solar cell. It is shown that, for a broad range of densities (134 - 15962 discrete elements), parameters vary below 3 % for the analysed solar cell structure. Relatively low variation is a consequence of the adaptive character of the mesh -denser mesh in the vicinity of interfaces and in the regions of high gradient changes of the internal quantities, such as defect regions. 1 Introduction Energy demand is on the rise as population soars above 7.6 billion. Meeting this demand by fossil or nuclear fuels proved to be unsustainable. Photovoltaics (PV) presents a promising renewable energy solution to supply the world with clean and low-cost electrical energy as it uses solar cells to convert sunlight directly into electrical energy. Among various PV technologies [1] crystalline silicon (c-Si) wafer-based solar cells are dominating the market today (>90 %) [2], with commercial PV modules reaching conversion efficiency in the range of 17-20 %, whereas the research cell record efficiency is currently 26.7 % [3]. This record-braking Si cell design has a structure of a so called heterojunction (HJ-Si) device. The term hetero refers to the use of different energy band-gap materials. In particular, HJ-Si combines the advantages of Si thin-film and wafer-based technology. The HJ-Si solar cells are also the subject of our numerical investigation in this paper. State-of-the-art HJ-Si devices are quite complex, both optically and electrically, as they are comprised of ultra-thin (few nanometres) and thick layers (hundreds of micrometres), made of semiconductor, oxide and metal materials. Various textures and contacting schemes are applied to obtain superior performance of the device [1]. Numerical modelling proved to be an indispensable tool in the R&D cycle as it offers an in-depth analysis of such complex devices which consequently leads to cost reduction and improvement of efficiency. In this paper we are going to use both optical (SunShine [4]) and electrical (Sentaurus TCAD [5]) numerical simulations, with the emphasis on the latter, to simulate optoelectronic properties of the HJ-Si solar cell. In particular we will show the importance of meshing density when doing finite element method analysis (FEM) on such layered structures. This analysis is a key step towards reliable simulations of realistic devices and further optimization of their operational performances. 2 Structure and principle of operation of a HJ-Si device The HJ-Si cell (Figure 1) employs a lightly doped p or n type crystalline silicon (in our case n-type) as an absorber with a thickness of around 100 - 200 ^m. Both sides of the absorber are capped with a thin (1-5 nm) passivation layer (typically intrinsic hydrogenated amorphous silicon, i-a-Si:H or other alternative material) in order to passivate c-Si surface defects - fill the free, dangling bonds of Si. Passivation layers are usually covered with a highly doped p-type (one side) and n-type (other side) a-Si:H material (10-50 nm thick) to assure selectivity of electron/hole extraction and to supply them to the front (top) and rear (back) contact. o-+ o- Figure 1: Structure of a typical HJ-Si solar cell. A highly doped n-type transparent conductive oxide (TCO) is used as a front contact in combination with a metal grid (the grid is not shown in the figure). Indium tin oxide (ITO) is usually employed in HJ-Si devices as the TCO. Also at the rear side a TCO layer is often used in combination with Al or Ag metal contacts to ensure good optical and electrical characteristics of the rear part of the device. A cell with contacts solely at the rear has also been developed in the form of an Inter-digitated-back contact (IBC) scheme (not shown here). The device employs some form of light trapping by Ligbt def. reg. 2 nni 1.12 eV n-c-Si 150jim 1.12 eV? ITO 70 Dm p-a-Si:H 10 nm 1.75 eV i a Si:H 5 nm 1.75 eV 2 mu 1.12 «V i a Si:H 5 nm 1.74 eV n-a-Si:H 10 nm 1.80 eV ITO 100 nm Al 500 nm ERK'2018, Portorož, 182-327 324 interface micro-texturing and anti-reflection coating (ARC) to maximize the absorption of light in c-Si and minimize the reflection loss (in Figure 1 the texturing of the interfaces is not shown, but it was considered in optical simulation of the device). Optical simulations were performed on the entirety of the schematic shown in Figure 1, whereas electrical simulations omitted the ITO layers at this stage and assumed ideal ohmic contacts instead. 3 Numerical models In electrical simulations with the Sentaurus TCAD we used the drift-diffusion based model to solve the transport equations of free electrons and holes in the HJ-Si device. The densities of electrons and holes were calculated with Fermi-Dirac statistics in both low and highly doped regions. Mobility and diffusivity were modelled by the Einstein relation. We also included the Masetti model for doping dependence of mobility [6]. Band-gap narrowing effects were considered negligible and were not included in the simulations. We incorporated Auger, Shockley-Read-Hall (SRH) and radiative recombination mechanisms. Passivation quality at front and rear surfaces of the c-Si absorber was modelled by the 2 nm thick defective regions (def. reg. in Figure 1). In these two regions amphoteric Gaussian-like trap distributions were introduced, similar to those used in modelling of the a-Si:H layers, which included exponential tail-states and Gaussian-like dangling bonds. Amphoteric nature of dangling-bonds was modelled as donor and acceptor Gaussian functions whose peak trap densities were apart by correlation energy. Quantum tunnelling was included at layers thinner than 10 nm. The model used is based on Wentzel-Kramers-Brillouin (WKB) approximation and takes into account only intra-band tunnelling. Doping ionization was assumed to be 100 %. All material parameters were taken from our previous publication [7]. The generation rate profile across the structure, due to absorbed sunlight, was calculated by the SunShine simulator. The performed electrical simulations of the solar cell were two-dimensional (2-D), although the structure of the device would allow for a 1-D simulation (device planar dimensions are much larger than its thickness and all parameters are also uniform in the planar directions). However, a 2-D model was used instead since it will be required in our further work, where the device structure will be changed (e.g. different lateral contacting schemes). Here, a 2-D model with appropriately small planar dimensions was applied. In finite element method (FEM) linear interpolation between the calculated points (nodes of the discrete triangles) was applied, which was in our previous studies indicated to be an effective solution compared to higher order interpolations. In the analysis presented here the attention was paid on variation of mesh density - variation of number of discrete elements in 2-D FEM simulations with the Sentaurus TCAD tool. In case of a 2-D simulation these elements are represented by triangles. There is a trade-off between validity and accuracy of simulations on one hand and computational time and processing power of computing on the other. Many simulations have to be performed during the optimization process of a device, therefore it is very important to find optimal settings for reliable, accurate and time-efficient simulations. 4 Simulation results Figures of merit of a solar cell performance are the current density - voltage (J-V) characteristic and the corresponding parameters under full sun (1000 W/m2) illumination: the short-circuit current density (Jsc), the open-circuit voltage (Voc), the maximum power-point (MPP) current density and voltage (Jmpp,, Vmpp), the fill-Factor (FF, see definition in Figure 2) and the conversion efficiency (n) To demonstrate these characteristics we show in Figure 2 the J-V curve obtained by the optoelectrical simulations of the HJ-Si solar cell where 5654 discrete triangular elements were used in electrical simulations. As shown later, this number of elements used in an adaptive mesh is recognized as a proper discretization setting in simulations of the device under scope. V [mV] Figure 2: Simulated J-V curve and parameters of the HJ-Si cell. 5654 of discrete elements were used in the simulation. The Jsc is the maximum current density a cell can produce under given illumination (short-circuit condition, load with zero resistance). The Voc is the maximum voltage obtained (under open circuit condition-no load). Fill-factor tells us what percentage of the product of Jsc and Voc we can utilize when applying optimal load to obtain the maximum power point (MPP) condition. Conversion efficiency, as a final parameter, defines how much of the incoming solar energy is being converted into electricity. All parameter values shown in Figure 2 are given for the simulated solar cell reaching a conversion efficiency of 23.8 % (the cell was not optimized). 183 4.1 Analysis of mesh density in electrical simulations In FEM analysis mesh density can play a crucial role when it comes to validity and accuracy of simulations. It is prudent to use denser mesh near material interfaces or regions with high gradients of internal quantities, such as defect densities near interfaces. The Sentaurus TCAD tool allows for a definition of an adaptive mesh according to various conditions (doping density gradients, material or region interfaces etc.). Meshing starts by defining maximum and minimum triangle element dimensions (x-lateral and y-vertical axis for 2-D). In our case minimum and maximum values for the x axis were both equal to 1 ^m. Depth of the device was assigned to the y axis and element dimensions were restricted to the interval between 1 pm-10 ^m. We implemented an adaptive mesh for all material interfaces by defining a starting triangle element height-yo and an enlargement factor-k when going away from the interface or high gradient region. The adaptive mesh algorithm uses a material interface as its starting point where it applies triangles with the height yo (see Figure 3). Next step repeats the procedure by multiplying yo with k and forming larger triangles. The height y remains constant when it reaches the maximum value defined at the beginning of this section. Interface edge Figure 3: First two steps of the adaptive mesh algorithm. In our analysis 20 different mesh densities were included with the total number of discrete elements ranging from 134 to 15962 obtained by varying the factor k from 1.005 to 200 (yo =1 pm and was kept constant) An example of low, medium and higher schematic mesh densities can be seen in Figure 4. No. of elements: p, i, def. reg. 132 3568 5654 n-c-Si def. reg., i, n Figure 4: Examples of different mesh densities. Results of our analysis are shown in Figure 5 (a-d). Simulated values of the parameters Jsc, Voc, FF and n are presented as a function of number of mesh elements. The span of No. of elements is relatively broad. We can see that the results of all parameters vary the most for the smallest No. of elements. a) > E b) c) 39.7 39.6 39.5 39.4 39.3 39.2 39.1 39.0 38.9 — ■■ » • • 1.5 % • • • 100 250 500 1000 2500 500010000 25000 No. of elements ---------------- 0.3 % • 100 250 500 1000 2500 500010000 25000 No. of elements 83.0 82.8 82.6 -82.4 82.2 -82.0 81.8 -81.6 81.4 -»--------------- 1.2 % • -----■»■■ il« «mti 100 250 500 1000 2500 500010000 25000 No. of elements d) Figure 5: Parameters of the solar cell vs No. of elements. 747.5 747.0 746.5 746.0 745.5 745.0 0 184 1.5 1.0 0.5 - a o.o -4 -1.0 --1.5 - 0.000 0.005 0.010 0.015 50.000 60.000 70.000 80.000 90.000 150.010 150.015 150.020 150.025 150.030 Depth [i m] Figure 6: Energy-band diagram of a reference cell at MPP with 5654 elements. 2 -0.5 - -2.0 The values of JSc (Figure 5a) and n (Figure 5d) settle down in two regions (250 - 500 and 550 - 1100) before they finally reach their constant value after 5654 No. of elements. The values of Voc vary slowly after the first two points up until 2500 No. elements when the value changes to the final and constant value at 5654. Values for FF are quite constant right after the first two points due to it being a kind of composed parameter (in the definition products of J and V are present) and the variations of the components can cancel each other out. Maximum relative variation was observed for n and was 2.7 %, followed by Jsc with 1.5 %, FF with 1.2 % and Voc with 0.3 %. We have to point out that these values are really low (low sensitivity to the mesh density) due to the utilization of an adaptive mesh. Mesh density as low as 134 elements can be used when doing preliminary research, but for an accurate and still time efficient simulation we recommend at least 5654 elements (triangular element dimensions range from 1.15e-6 ^m2 -1.17 ^m2). We also have to take special care if the gradients of internal parameter distribution change significantly during optimization - in this case the mesh has to be changed accordingly during the optimization process. 4.2 A view into internal quantities After our simulations pass the test of validity and accuracy we can use them to get an inside view in the device. Distributions of internal quantities, which often cannot be accessed by measurements, can be determined by simulations. This presents a key advantage when we do further optimisation of the devices. These internal quantities are for instance the energy band diagram, free charge carrier densities, recombination rates, electric field and others. Here we show, as an example, only the energy band diagram (Figure 6) calculated with 5654 elements and under MPP. The result of the 2-D simulation is presented as a function of the vertical direction (depth) only, belonging to a vertical cross-section taken at the middle position of the 2-D structure. Results corresponding to cross-sections taken at different positions differ only in the range of numerical errors. From the diagram we can see, for example, how conduction and valence bands change throughout the front and rear heterojunctions. We can also see the splitting of quasi-fermi levels which is the voltage we have at the terminals (Vmpp). The shape of the bands at material interfaces (i-a-Si regions in our case) can tell us the level of carrier selectivity (energy barrier for electrons in p,i-a-Si regions at the front and an energy barrier for holes in i,n-a-Si at the rear). We can also see regions where there are high gradients of internal quantities which allows us to redefine the mesh at this critical points to insure accurate simulations (defect region interfaces in our case) 5 Conclusions In numerical simulation of solar cells, mesh type and density play an important role trading off accuracy and validity for time efficiency and computational power. We showed, that by using an adaptive mesh in electrical simulations of a HJ-Si solar cell, main output parameters vary with mesh density (134-15962 No. of elements) up to 2.7 % of maximum relative difference. The values of all parameters become constant when we use at least 5654 elements in the mesh. When an adaptive mesh is used, we can do preliminary optimization with as few as 134 points for a structure shown in Figure 1. To perform time efficient and accurate simulations, an adaptive mesh of at least 5654 points are to be used. Literature [1] A. Luque and S. Hegedus, Eds., Handbook of photovoltaic science and engineering, 2nd ed. Chichester, West Sussex, U.K: Wiley, 2011. [2] F. ISE, 'Photovoltaics report', vol. 2016, p. 45. [3] K. Yamamoto et al., 'Record-breaking efficiency back-contact heterojunction crystalline si solar cell and module', 33rd European Photovoltaic Solar Energy Conference and Exhibition, vol. 2017, p. 4. [4] J. Krč, S. Franc, and T. Marko, 'One-dimensional semi-coherent optical model for thin film solar cells with rough interfaces', in Informacije MIDEM 32, Ljubljana, vol. 2002. [5] Sentaurus Device, https://www.synopsys.com/silico n/tcad/device-simulation/sentaurus-device.html [6] G. Masetti, S. Maurizio, and S. Sandro, 'Modeling of carrier mobility against carrier concentration in arsenic-, phosphorus-, and boron-Doped Silicon', IEEE Transactions on Electron Devices, no. 7, p. 6, Jul. 1983. [7] J. Balent, J. Krč, A. Čampa, F. Smole, and M. Topič, 'Electrical simulations of hetero-junction silicon solar cells', in 53rd International Conference on Microelectronics, Devices and Materials & the Workshop on Materials for Energy Conversion and their Applications: Electrocalorics and Thermoelectrics, Ljubljana, Slovenia, 2017. 185