S Le Roy

M Q Hoang

Bipolar charge transport model to simulate the impact of nanometric scale processes on the space charge behaviour in polyethylene

Keywords: charge transport model, space charge, roughness, physico-chemical structure, injection, polymeric insulation, interface

Predicting the electric field distribution in polymers used as electrical insulating materials remains the Holy grail, as the presence of charges disturbs the Laplacian electric field. Charges arising from the electrodes are one of the dominant mechanism of charge generation, particularly in polyethylene-based materials. Hence, nanometric scale processes at play at the interface have a non negligible impact on charge injection. In the present study, a bipolar charge transport model developed in 2D is used to simulate the impact of several nanometric scale processes, such as the variation of the barrier height linked to the chemical structure of the material at the interface, as well as surface roughness. Simulation results as regards net charge density, current, but also recombination rate, will be compared to a case where no specific supplementary hypothesis is set at the electrodes. At last, simulations have been performed for a combination of roughness and barrier height variation along the electrodes.

Introduction

Challenges in the field of polymer electrical insulations under DC stress are linked to the prediction of the electric field distribution, whatever the complexity of the system under study. Indeed, these polymer materials, used as electrical insulations, can be found in various systems such as cables for energy transmission for terrestrial [1] or aeronautic [2] applications, but also in microelectronic or power electronic devices. In these systems, when the applied constraints are low enough, macroscopic models [3][4][5] can be used to calculate the electric field distribution. However, when the applied constraints (electrical, thermal, …) are high, then space charges are present in the material bulk and disturb the Laplacien electric field. Macroscopic models fail to reproduce the space charge behaviour, as charge generation (i.e. injection, ionization), transport are not accounted for. To overcome these limitations, a more microscopic approach has been developed since two decades [6][7][8][9], based on bipolar charge transport models (BCT), also called mesoscopic models. These models are fluid or electro-hydrodynamic models, and account for charge generation, transport and accumulation, mostly in polyethylene based materials. One of the main issues with this microscopic approach remains to correctly describe charge generation at an interface, as it remains the dominant process providing charges inside the material. Researches on the subject have already shown that nanometric scale processes [10,[START_REF] Roy | Proc. IEEE Int. Conference on Dielectrics (ICD)[END_REF], as well as geometry (i.e. smoothness) of the surface [START_REF] Doedens | [END_REF]13] are key parameters in the interface description. The present paper focusses on modelling the space charge, current density and recombination rate behaviour when accounting for surface roughness, and physico-chemical variation of the energy levels at the interface, using a 2D bipolar charge transport model already published in [13]. At last, a coupling of these two hypotheses will be performed.

Model Description

A 2D bipolar charge transport (BCT) model, already published in the literature [13,[START_REF] Hoang | Proc. IEEE Int. Conf. on Dielectrics (ICD)[END_REF], has been used in the present study. The sample considered is a low density polyethylene (LDPE) of thickness 150 µm, and a width of 45 µm. It is considered as homogeneous in the third dimension. An applied electric field of 40 kV/mm is set on the sample at room temperature (T=25°C). This corresponds to an applied voltage of 6 kV, which is applied at the top electrode, while the bottom electrode is grounded. The model accounts for electronic species (electrons and holes), generated by injection at each electrode, depending on the sign of the applied electric field. Transport is also accounted for using a field dependent mobility (hopping), as well as trapping into a single level of deep traps for each kind of carrier, from which they can detrap using a thermally activated coefficient. The model also features recombination of charges of opposite sign, with mobility dependent parameters, of the Langevin type. A schematic representation of the model hypotheses is proposed on Figure 1.

The equations to solve are of the form:

𝜕𝑛 𝑎 𝜕𝑡 + ∇. (𝑛 𝑎 µ 𝑎 𝐸 -𝐷 𝑑𝑖𝑓𝑓-𝑎 ∇𝑛 𝑎 ) = 𝑠 𝑎 (1)
∇. (𝜀 0 𝜀 𝑟 𝐸) = (+𝑞𝑛 ℎµ + 𝑞𝑛 ℎ𝑡 -𝑞𝑛 𝑒µ -𝑞𝑛 𝑒𝑡 ) = 𝜌 (2) Where na is the charge density (1/m 3 ), a refers to the charge carrier, i.e. electron (e) or hole (h), mobile (µ) or trapped (t). q is the elementary charge (C), t is the time, µa the mobility (m 2 /V/s), for each carrier. E is the electric field (V/m), ε0 the vacuum permittivity, εr the relative permittivity of the material (2.3 for LDPE), and ρ is the net charge density (C/m 3 ). All variables in equations (1-2) and in the following are space and time dependent, even if these variables have been omitted, for sake of simplicity.

The mobility for each type of carrier is of the hopping type, function of the electric field and the temperature, and is of the form:

𝜇 𝑎 = 2𝜆𝜈 𝐸 𝑒𝑥𝑝 ( -𝑞𝑤 𝑎 𝑘 𝐵 𝑇 ) 𝑠𝑖𝑛ℎ ( 𝑞𝜆𝐸 2𝑘 𝐵 𝑇 ) (3)
λ is the hopping distance (m), ν the attempt to escape frequency (s -1 ), we,h are the hopping barrier height (eV), for electrons and holes, kB is the Boltzmann constant, and T the temperature.

Ddiff-a refers to the diffusion coefficient, and follows the Einstein relation of the form:

𝐷 𝑑𝑖𝑓𝑓-𝑎 = 𝑘 𝐵 𝑇 𝑞 µ 𝑎 (4)
sa are the source terms, reflecting all the physical processes not linked to transport. An example of source term for mobile electrons is given in the following:

𝑠 𝑒µ = -𝐵 𝑒 𝑛 𝑒µ (1 - 𝑛 𝑒𝑡 𝑁 0𝑒𝑡 ) + 𝜈𝑛 𝑒𝑡 𝑒𝑥𝑝 ( -𝑒𝑤 𝑡𝑟𝑒 𝑘 𝐵 𝑇
) -𝑞𝑆 1_𝐿 𝑛 𝑒µ 𝑛 ℎ𝑡 -𝑞𝑆 3_𝐿 𝑛 𝑒µ 𝑛 ℎµ (5) Here, Be is the trapping coefficient (s -1 ), N0et the trap density for electrons, wtre is the detrapping barrier height for electrons (eV). S1_L and S3_L are the Langevin recombination coefficients, of the form:

𝑆 1_𝐿 = 𝜇 𝑒 𝜀 0 𝜀 𝑟 ; 𝑆 3_𝐿 = 𝜇 𝑒 +𝜇 ℎ 𝜀 0 𝜀 𝑟 (6)
The same kind of equation would hold for trapped electrons and mobile and trapped holes, and can be found in [7].

Charge generation is only due to injection at each electrode, function of the sign of the local electric field EX, and follows the Schottky law:

𝑗 𝑆𝑐ℎ𝑜𝑡𝑡 (𝑋) = 𝐴𝑇 2 𝑒𝑥𝑝 ( -𝑞𝑤 𝑎𝑖 𝑘 𝐵 𝑇 ) [𝑒𝑥𝑝 ( 𝑞 𝑘 𝐵 𝑇 √ 𝑞𝐸 𝑋 4𝜋𝜀 0 𝜀 𝑟 )](7)
A is the Richardson's constant, wai is the injection barrier height, for electrons or holes (eV), and the coordinate X refers to the anode or to the cathode. No extraction barriers are accounted for at the electrodes, so the extraction for holes at the cathode and for electrons at the anode are of the form:

𝑗 𝑎 = 𝑛 𝑎 µ 𝑎 𝐸 -𝐷 𝑑𝑖𝑓𝑓-𝑎 ∇𝑛 𝑎 (8) The total external current density is calculated by integrating over the dielectric surface the following Maxwell equation:

𝐽 𝑡𝑜𝑡 = 𝑗 𝑐𝑜𝑛𝑑 + 𝜀 0 𝜀 𝑟 𝜕𝐸 𝜕𝑡 (9)
Where jcond refers to the conduction current density, being the sum of the conduction current densities of each mobile carrier. The conduction current integrates the current arising from charge injection.

The model has been applied to a LDPE sample, for which optimized parameters have been published in the literature for a one-dimension model [7], and are listed in Table I. These parameters have been kept untouched in the present paper, as our goal is mainly to observe the impact of interface processes. Most equations are field and temperature dependent, so simulations could be performed for different electric fields and temperatures, or even for field and/or temperature gradients. This is however not the purpose of the present paper, so simulations are performed for a single applied electric field (40 kV/mm) and at room temperature. The present model, developed with COMSOL Multiphysics® uses the Transport of Diluted Species (TDS) module to solve the continuity equation for each kind of carrier [13]. This module already includes stabilization solutions to prevent oscillations. The Poisson equation module is used to couple these equations to the Poisson's equation. Backward Differentiation Formula solver is used for the time integration (maximum order 2, and minimum order 1). A zero flux is set as boundary conditions on the left and right sides. Before voltage application, the material is free of charges.

Results

Reference case

Simulations have first been performed for a reference case study where no specific hypothesis is set on charge generation. Only the Schottky law (i.e. Equation ( 7)) is accounted for charge generation. Parameters are the ones presented in Table 1. As it is difficult to compare the 2D charge density simulation results to experimental data, the charge density has been integrated over the full length of the sample (x-axis), and is then plotted as a surface graph as a function of sample depth and time. The aim is to be directly comparable with experimental data measured with the Pulsed Electro Acoustic (PEA) method as an example. In the present case however, the influence charges, sum of the capacitive charges and image charges, are not represented in order to see the charge behaviour at the electrodes. Figure 2a) presents the simulated net charge density as a function of space and time, for the reference case study. Positive charges are mainly observed on Figure 2a), arising from the anode, and travelling through the bulk with time. This is due to a lower injection barrier height for holes compared to the one for electrons. Negative charges are also observed at the cathode, but in a smaller amount. The bulk remains mainly positive during polarization, with only slight changes in the net charge density with time. This behaviour is typically what is measured experimentally in a LDPE under 40 kV/mm, at room temperature [7][8][9]. The simulated current density (Figure 2b) as a function of time is the integration over the sample thickness and over the electrode length of equation ( 9). It has been simulated for a polarization time of 100000s instead of 10000s, in order to observe a possible steady state behaviour. The current density in Figure 2b is decreasing with time for times below 4000s. The current is then slightly increasing, and achieves a steadystate value after around 40000s. At short times, the current density is mainly driven by the injection current density, 1.

which is of the Schottky type, function of the electric field at the electrode. This injection current is high just after voltage application, and then decreases, as the electric field at the electrode decreases due to trapped homocharges at each electrode. At longer time, the current density is driven by the bulk current, i.e. the transport, trapping and detrapping. The steady state is reached once there is an equilibrium between the injection current and the bulk current. The total recombination rate, which is related to the electroluminescence signal, is presented on Figure 2c as a function of time. It is the sum of all the interactions between mobile and trapped charges of opposite polarity. As an example, the recombination rate between mobile electrons and trapped holes can be written as:

𝑅 1 (𝑥, 𝑡) = 𝑆 1_𝐿 (𝑥, 𝑡)𝑛 𝑒µ (𝑥, 𝑡)𝑛 ℎ𝑡 (𝑥, 𝑡) (10) Which is function of the Langevin recombination coefficient, and of the density of mobile electrons and trapped holes. The total recombination rate is null for times below 10-20s, as there is no charge prior to voltage application. It is then always increasing, and reaches a steady-state value after 40000s. The time for which the steady state is achieved for the total recombination rate is the same as the one previously observed for the total current density. The recombination coefficients are function of mobility, but do not vary to a large extent in space and time. The main variations arise from the charge density. As charges are injected and quickly trapped at the vicinity of each electrode (i.e. it forms homocharges), the charge density for one type of charge varies along the y axis [7]. However, for one type of charge, and a specific position along the y axis, the charge density is almost always increasing with increasing time. As an example, the amount of positive trapped charges for y= 50µm is of the order of 0.6 C/m 3 at t=100s, and it increases to 6 C/m 3 at t=10000s. This behaviour is what explains the always growing recombination rate. Compared to the total recombination rate simulated with a constant mobility and constant recombination coefficients [7], which increases to a maximal value and then decreases, the total recombination rate in the present study is always increasing with polarization time. The factors affecting the shape are linked to the variation of the mobility, the recombination coefficients, and the absence of initial charges prior to polarization. However, the total recombination rate for the reference case study is always lower than 10 -3 C/m 3 /s, which is the value associated with the noise level of electroluminescence measurements [7]. Although the present simulated results differ slightly from the ones published for a constant mobility and constant recombination coefficients [7], the main behaviour for space charge and current are respected compared to the experimental data.

Surface roughness

A certain roughness is always present at the surface of an insulating material, hence participating at the variation of the electric field distribution at the interface. It has already been published in [START_REF] Doedens | [END_REF]13] that concave protusions increase the electric field at the interface, while convex one will decrease the electric field. In the present paper, a Schottky law ( 7) is used to simulate charge generation at each electrode. This equation varies exponentially with the electric field, with a large variation of the current density for high electric fields in the range 10 6 -10 8 V/m. On the contrary, there is quasi no variation of the injection current density for electric fields in the range 0-10 6 V/m, due to the thermal term of the Schottky injection. Hence, concave defects will have more impact than convex ones. In the present case, only concave defects have been simulated. This was also motivated by the need to decrease as much as possible the calculation time, which is large is the present case. The shape of each protusion is elliptic, ranging from 200 to 800 nm, while the width of the ellipse ranges between 50 and 400nm. This roughness is within the range of what can be measured for a lab-made sample [START_REF] Taleb | Proc. IEEE Conf. on Electrical Insulation and Dielectric Phenomena (CEIDP)[END_REF], and what has also been measured for different samples at the industrial scale [START_REF] Doedens | [END_REF]. An example of surface roughness has already been presented in [13]. Here the surface roughness has been simulated at the top and at the bottom electrode. The roughness is not symmetric at the top and bottom electrode and does not cover the whole electrode for sake of calculation time. The surface roughness is located from 10 to 21 µm at the top electrode, while it is located between 7 to 28 µm at the bottom electrode, for an overall electrode length of 45 µm. Profilometry measurements have been performed on a LDPE sample, to measure the mean surface roughness, and the maximal 'height' [START_REF] Taleb | Proc. IEEE Conf. on Electrical Insulation and Dielectric Phenomena (CEIDP)[END_REF]. Other data related to the present configuration (i.e. number, shape, and location of the ellipses) are totally arbitrary. Figure 3 presents the 2D net charge density as a function of space for a time of 100s (Figure 3a) and for 10000s (Figure 3b). Charge injection is highly enhanced at the defect maximal height, i.e. where the electric field is maximum due to the particular geometry. The most injecting defects are the ellipses having a high ratio b/a, b being the ellipse vertical radius, i.e. the defect height, and a being the horizontal radius, i.e. the defect half-width. The maximal net charge density overcomes 200 C/m 3 , although the net charge density in Figures 3 is limited to ±50 C/m 3 . The injected charges penetrate then fast within the bulk, and after 100s, charges injected at the vertical of the ellipses having the higher b/a ratio have already crossed the sample and achieved the cathode (150 µm). The electric field is decreased next to each electrode soon after voltage application due to a large amount of injected charges. However, the electric field is enhanced in the volume ahead of the charge cloud, leading to an increased charge mobility, and hence an increased charge velocity (v=µE). After 10000s (Figure 3b) positive charge injection is also visible at the top electrode where no cavities exist. Negative charges are also injected at the bottom electrode where no cavities exist, but they are not visible due to the choice of range for the net charge density. These simulated results show how generation processes at the electrodes controls the bulk behaviour. This has however already been highlighted for such charge transport models using a parameters sensitivity analysis (Soboll indexes) [17]. Integrating these results along the x axis leads to Figure 4, i.e. a 2D cartography as what is measured during space charge measurement. At short times (lower than 500s), positive charges, injected from the anode, and transporting through the dielectric thickness, are observable on Figure 4. This is only due to the surface roughness, that leads to large amounts of injected charges on small regions. As observed previously on Figures 3a and3b, positive charges are dominant compared to negative ones, due to the choice of parameter value for injection barrier height. Negative charges are also observable on Figure 4, and where not observed in the case of a smooth rough electrodes at short times is the amount of charges next to the electrode (positive charges at the anode, and negative charges at the cathode), which decreases with time for a rough electrode. Charges, being injected in a large amount, decrease the electric field next to the electrode, leading to a decrease of the injection at each electrode. It is to note that a variation of one decade in the electric field in the range 10 6 -10 8 V/m implies a variation of more than one decade on the injection current density. Meanwhile, the electric field is enhanced inside the bulk, so the charge mobility is enhanced. This can be observed on Figure 5, which compares the net charge density, integrated over the x axis, as a function of the material thickness for a time of 500s, and for the different case studies presented in the paper. It is even noticeable at the cathode (left) that positive charges, transported from the anode, arrived at the opposite electrode, and may accumulate here as the electric field is decreased.

Figure 6a presents the comparison of the total current density as a function of time for rough and smooth electrodes. In the case of rough electrodes, at short time, the current density is higher of more than one decade compared to the one for smooth electrodes. As the quantity of injected charges is higher (Figure 5), the calculated current density is also higher, particularly at short times where the injection current controls the total current. The current density in the case of rough electrodes decreases then and achieves a quasi-steady state value, also higher than the reference one. It is to note that the sharp decrease of the current density for rough electrodes is comparable at short times (<100s) to what is normally observed experimentally, being either attributed in the literature to an initial charge present inside the sample, or to polar processes. Two regions can be observed on the total recombination rate (Figure 6b) as a function of time, when a rough electrode is considered, with a threshold at around 1000s. The total recombination rate soon after voltage application is high (5 10 -4 C/m 3 /s) and slightly decreases with time. When a rough electrode is considered, positive charges penetrate faster inside the dielectric, reaching the cathode in less than 100s. Electrons, injected at the cathode, being mobile or trapped, can recombine with these positive charges. This allow having a high recombination rate particularly at the vicinity of the cathode. This recombination rate decreases with time as these fast injected positive charges can be extracted, or are 'consumed' by recombination. The second part of curve (for times above 1000s) is mainly due to the charges arising from the rest of the surface where no defects are present, with an increase of the total recombination rate, as in the reference case. However, there is a non-negligible effect of the defects on the overall electrode electric field. The total recombination rate is always higher that in the reference case. In the case of rough electrodes also the total recombination rate remains below the 'threshold' recombination rate (10 -3 C/m 3 /s) previously set as representing the experimental noise level.

Variation of the injection barrier height along the x axis

An easy way to account for the physico-chemical structure of the material at an interface is to make the injection barrier height value vary along the electrode surface, i.e. the x axis in the present case. This could account for the variation of energy states along the electrode, due to chain conformation, available energy levels for injection, etc. In this section, the value of the barrier height for electrons and holes has been calculated randomly using Matlab®. A small variation on the barrier height implies a large variation of the injection current density [START_REF] Roy | Proc. IEEE Int. Conf. on Properties and Applications of Dielectric Materials (ICPADM)[END_REF]. A minimal and a maximal value for barrier heights have been set., as being ±10% of the injection barrier height presented in Table 1. These barrier height values are then extrapolated along the x axis on the top electrode for holes and on the bottom electrode for electrons. This percentage value is arbitrary, and other values could have been set, as it is not possible in reality to have access to such informaiton. However, calculations already performed using the Schottky equation with different barrier heights [START_REF] Roy | Proc. IEEE Int. Conf. on Properties and Applications of Dielectric Materials (ICPADM)[END_REF] show that a too high barrier height value (>1.5 eV) leads to insignificant injection current densities, i.e. a low injected charge density. Barrier height values lower than 1.05 eV lead to non negligible impact in the calculations (problem of timestep calculation, sharp variation of the electric field at the electrode…). An example of barrier height distribution for electrons as a function of the x axis is given on Figure 7. Figure 8a presents the net charge density as a function of space and time for this case study. The global space charge behaviour is comparable to the one simulated for a constant injection barrier height, with a large injection of positive charges at the anode, while negative charge injection at the cathode is lower. The bulk remains then mainly positive until the end of the polarization time. However, the dynamic is faster in the present case, with a higher net charge density, for holes and for electrons. When the barrier height is decreased, there is an enhanced charge injection, while in other regions the injection barrier height is increased, i.e. the injected charge density is lower. The overall net charge density is however higher compared to the reference one (see Figure 2). In the present case, the variation of the injection barrier height is at the origin of the electric field distribution. A large amount of injected charges leads to a local decrease of the electric field. The barrier height variation is however small (less than 10%). The electric field next to the electrode does not achieve values observed for the case of rough electrodes, where high positive charge density were siulated (>200 C/m 3 ), so the decrease of the injected charge is not as pronounced as in the case of rough electrodes (see Figure 5).

A large number of simulations have been performed, with a)

b) Figure 6. Comparison of a) the current densities and b) the total recombination rates, as a function of time, for the case of smooth and rough electrodes. Parameters of Table 1. different barrier height distributions for electrons and holes at each electrode, in order to verify if the results were reproducible. Table 2 presents the mean values of the barrier heights for electrons and holes for some of the simulated cases that have been performed. Although these variations are small, this leads to a non-negligible impact on the simulated results, and even on the integrated net charge density as a function of time and space (see Figure 8a) to d)).This, among other processes, could explain the small differences that are always observed experimentally for the same kind of materials but different samples. Figure 9a presents the current density as a function of time, for the reference case and for the different simulations as presented in Table 2. The current density is always higher when injection barrier height varies compared to the reference case, but the shape is globally similar. This higher value of the current density is explained by a higher charge density in all case studies compared to the reference case (see Figures 8). It is to note that the current density value does not converge to the same steady state value after longer polarization times. Figure 9b presents the total recombination rate as a function of time, for the reference case and for the different simulations as presented in Table 2. The total recombination rate has the same shape as the reference one for all case studies, with higher values. These higher values are as previously explained by the higher charge density for all case studies compared to the reference case. The total recombination rates in the present case stay lower than the one simulated in the case of rough electrodes.

Coupling electrode roughness and injection barrier height variation along the x-axis

Both interface roughness and variation of the barrier height can take place at the same time, as the first one is linked to the 'geometry', while the second one is linked to the material structure at the surface. Simulations have then been performed with both interface hypotheses, the surface roughness being the same as the one presented before, while the variation of barrier height for electrons and holes is random, but the conditions (±10% of the optimized barrier height values) are the same as in section 3.3. The global space charge behaviour is comparable to the one when only roughness is accounted for (Figure 5). Figure 10 compares the net charge density as a function of sample depth for specific polarization times for the two case studies. There is only a small variation of the net charge density in the present case compared with the case of roughness only. The maximal difference between the two case studies is observable at short time, while the net charge density 2. Other parameters as in Table 1.

is almost the same after 3 hours under voltage. The authors thought it would be possible, either with different roughness, or with different interface processes such as what is presented in the present section, to reproduce specific charge packets behaviour observed particularly for polyethylene based materials under relatively high electric fields [START_REF] Hozumi | [END_REF], and often simulated using an arbitrary threshold in field on the charge injection. This is actually not the case, even for higher values of applied electric fields. However, other interface processes need to be accounted for.

Conclusions

A 2D bipolar charge transport model has been used to simulate the impact of different nanometric scale processes on macroscopic measurable variables such as the net charge density, the current density and the total recombination rate. Surface roughness has been simulated using concave ellipses with different ratio, while a variation of the barrier height value for injection has been taken into account to simulate the variation of the physico-chemical structure of low density polyethylene at each interface. Simulations show that accounting for surface roughness changes the macroscopic behaviour to a large extent. It allows observing a high current density at short times, which has always been attributed in previsous researches to polarization or initial charges already present in the material. Rough electrodes also allow observing a fast charge injection at the electrodes, a lowering of the injection at short time, and a fast penetration in the bulk, which is not observable otherwise. A variation on the barrier height for injection leads almost to the same charge density behaviour, appart from the lowering of the net charge density at each electrode at short times. The current density has however a different shape, consistent with the reference case. The current density value is higher of more than one decade compared to the reference case. Different random distributions of the barrier height value lead to slighly different net charge density cartography and current density behaviour. The same observation holds in the case of different surface roughness. This could totally explain the small differences always observable when making the same measurement on two different samples. Lastly, the total recombination rate has a total different behaviour as a function of time when rough electrodes are accounted for. Roughness seems to have a dominant impact compared to a variation of the barrier height for injection. Other nanometric scale processes, such as an interfacial zone, or contact charges, need to be also accounted for, in order to observe which processes are necessary in the future simulations. 1.

Figure 1 .

 1 Figure 1. Schematic representation of the two-levels transport model. Conduction is by free charge in the transport level, associated with a field dependent mobility, for each kind of carrier.

Figure 2 .

 2 Figure 2. a) Simulated net charge density as a function of time and space, b) current density as a function of polarization time,and c) total recombination rate as a function of polarization time, in the case where no specific hypothesis is set on the interfaces. Parameters of Table1.

Figure 3 .

 3 Simulated net charge density as a function of space in a LDPE electrode, with top and bottom rough electrodes, for a time of a)100s, and b) 10000s. Parameters of Table1.

Figure 4 .

 4 Figure 4. Simulated net charge density as a function of time and space for the case where a rough zone of around 15 µm is considered at each electrode. Parameters ofTable 1.

Figure 5 .

 5 Figure 5. Net charge density as a function of sample depth for t=500s, and several case studies presented in the paper. Parameters ofTable 1.

Figure 7 .

 7 Figure 7. Example of distribution of the barrier height values along the x axis for electrons at the cathode. Distribution calculated as being in the range ±10% of the barrier height value in Table1.

Figure 8 .

 8 Simulated net charge density as a function of time and space for a) case study 1, b) case study 2, c) case study 3, and d) case study 4. Parameters are those of

Figure 9 .

 9 a) Simulated current density as a function of time and b) simulated total recombination rate as a function of time for the reference case and for different case studies related to barrier height variation as presented in Table

Figure 10 .

 10 Figure10. Net charge density as a function of sample depth for different polarization times, for the case where only surface roughness is considered, and the case where surface roughness as well as a variation of the injection barrier height is accounted for at each electrode. Other parameters of Table1.

Table 1 .

 1 Parameters used for the simulations

	Symbol	value	units
	Trapping coefficients
	Be electrons	1. 10 -1	s -1
	Bh holes	2. 10 -1	s -1
	Hopping barrier
	we for electrons	0.66	eV
	wh for holes	0.6	eV
	Trap densities
	Noet for electrons
	Noht for holes

Table

Table 2 .

 2 1, and mean parameters for electrons and holes injection are those of Table2. Mean values of the injection barrier height for electrons and holes calculated for four cases study.

		Mean electron barrier	Mean hole barrier
		height at the bottom	height at the top
		electrode (eV)	electrode (eV)
	Case study 1	1.286	1.161
	Case study 2	1.288	1.172
	Case study 3	1.264	1.142
	Case study 4	1.279	1.160

Acknowledgements

This research is partly funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.02-2019.33.