Abstract

Self-consistent drift–diffusion model has been widely employed to simulate the device performance of intermediate band solar cell (IBSC) under practical device configuration. However, one of the remained issues in the drift–diffusion modeled-based works is the difficulty to reach the IB carrier continuity through the self-consistent manner. In most of the previous reports the constraints were relaxed or just partially satisfied; which render the unreliable performance results and misguide the device design strategy. In this work, in order to solve this issue and to validate our results, we performed extensive simulations to fully disclose the significant effect of the IB continuity constraints by taking InAs/GaAs quantum dot-based IBSC as a model device using the semiconductor modules in COMSOL Multiphysics combined with the Fortran codes. We found that under rigorous satisfaction of IB continuity constraint, the band potential profiles for the IBSC with either doped or nondoped IB under various light illumination conditions are nearly identical to those under the dark conditions. Moreover, from the simulated current–voltage curve dependence on the light concentration ratio, we found the device performance based on drift–diffusion under rigorous IB continuity constraint showed similar tendency to the features simulated based on detailed balance principle except the much-lowered power conversion efficiency. Our work demonstrated here, serves as an accurate and reliable IBSC device design approach toward better IB material screening, efficiency improvement, optical management, and extended application in the emerging field such as the perovskite material-based IBSC.

1. Introduction

High-efficiency solar cells which exceed the Shockley–Queisser limit have been expected as a third-generation photovoltaic concept [1]. One of them is the intermediate band solar cell (IBSC) which utilizes the two-photon step absorption by inserting a narrow band of electron states in the forbidden semiconductor bandgap. Therefore, IBSC can absorb the energy smaller than the host material energy bandgap, which in turn can generate larger current than the single junction solar cells. As a result, IBSC has a great potential to overcome the efficiency of conventional single junction solar cells [2, 3]. There are several ways to fabricate IBSC, for instance, multistacked quantum dots [4, 5], highly mismatched alloy materials [6, 7], and impurity bands [8]. Among them, quantum dot is most widely studied as the building block to create IB in the bandgap. Various quantum dots superlattices have been proposed to fulfill the operation principle of the IBSC, particularly, self-assembled InAs quantum dot arrays fabricated by MBE have been most widely studied as a building candidate to create the IB toward an InAs/GaAs-based IBSC device [911].

In addition to these experimental attempts, numerical simulations for the IBSC have also been intensively carried out. At the first stage, the detailed balance model has been applied. It is useful to acquire the theoretical limit and the optimal bandgap for the material [1214]. However, it is hard for the detailed balance model to include or extract the detailed information regarding the actual device structure. Therefore, at the second stage, a combination of modified detailed balance model and drift–diffusion models are proposed to simulate the IBSC with concrete material information by taking into the consideration of the actual working conditions such as the photon filling effect, IB mobility effect, nonradiative recombination effect as well as light concentration effect [1517]. However, the numerical simulation in reference simplified the IBSC structure in order to calculate the ideal condition of IBSC, that is, i-layer was assumed being flat. Moreover, p-layer and n-layer were excluded. These issues were addressed at the third stage. Here, self-consistent drift–diffusion model for simulating conventional semiconductor simulation approach were proposed to simulate the IBSC by including the important carrier transport and recombination effect, which were either ignored or not adequately addressed during the simulation approaches proposed in the second stage. By combining Poisson equation and carrier continuity equation, this method enabled us to extract electron/hole density depend on the position and optimal device structure [18, 19]. However, in this approach, there still remained an unsolved issue regarding the IB carrier continuity constraints. In this work, we performed extensive simulation to fully disclose the significant effect of the IB continuity constraints on the IBSC device characteristics.

This paper is organized as follows: Section 2 introduces the drift–diffusion method for IBSCs and explains the theoretical procedure to ensure the rigorous satisfaction of IB continuity constraint. In Section 3, we present the results and extended discussion regarding the comparison of the IBSC simulation with and without satisfaction of IB continuity constraints. Finally, conclusions are drawn in Section 4.

2. Methods

2.1. Drift–Diffusion Modeling of IBSC

IBSC model based on drift–diffusion principle contains the Poisson equation, kinetic equations of charge carriers and the continuity equations for both electrons and holes in all the three bands: valence band (VB), conduction band (CB), and IB.

2.1.1. Poisson Equation

Poisson equation for modeling the electrostatic field in semiconductor structures is given as follows:here is the electrostatic potential, is hole density in VB, and are electron densities in CB and IB, respectively; are elementary charge, dielectric constant, vacuum permittivity, ionized donor, and ionized accepter densities, respectively.

By introducing the effective density of states, quasi-Fermi level and the assumption of Boltzmann distribution, the carrier densities shown in Equation (1) can be further described as follows:where , and are the effective density of states for CB, VB, and IB, respectively, and are schematically illustrated in Figure 1(a). , and are the band edges of CB, VB, and IB, respectively, as shown schematically in Figure 1(a). The quasi-Fermi levels for CB, VB, and IB: , and are also sketched in Figure 1(a) for better understanding.

2.1.2. Kinetic Equation of Charge Carrier

The current density of electrons and holes and are given as follows:

In both equations, the first term on the right side is a component of current density caused by the Lorenz force, taking into account only the influence of electric field , and establishes the effective mobility of electrons and holes . The second term on the right side of Equations (3) and (4) represents the effect of carrier diffusion in the direction of the gradient of carrier concentration and . The mobilities of electrons and holes and in CB and VB are related to the diffusion coefficients and by the Einstein–Smoluchowski equation as shown in Equations (3) and (4), in which is the Boltzmann constant and is temperature of the system which we have assumed here as room temperature 300 K during the whole simulations. Meanwhile, as shown in Figure 1(a), we assumed here that there is no current contribution from the IB, thus there is no formula for addressing .

2.1.3. General Continuity Equations

Under the steady state condition, namely, the time dependent carrier concentration , the general continuity equation used for modeling semiconductor can be defined as follows:

Here and in current continuity Equations (5) and (6) are the optical carrier generation and the recombination rates, where subscript stands for , which correspond to the transition between VB and IB, IB and CB as well as CB and VB. In order to account for the nonradiation recombination, the Shockley–Read–Hall (SRH) recombination rate is adopted in this work.

Next, we give explicit description regarding the calculation of the respective generation and recombination rate shown in Equations (5) and (6). For generation rates , if we assume , namely, the IB is located at an energy level which is closer to the conduction band edge than the valence band edge , the three generation rates correspond to the optical excitation among the VB, IB, and CB bands are given as follows:

Here , and represents the bandgap energy between CB and VB, the energy bandgap between CB and IB and the energy bandgap between IB and VB, respectively. , and are the optical absorption coefficients for the transition between CB–VB, CB–IB, and IB–VB, respectively. is incident photon flux density and is the position measured from the top emitter layer and directed downward to base layer as shown in Figure 1(b). Hyperparameters , and stand for the light concentration ratio, the solid angle subtended by the sun (1/46,050), the Planck’s constant, the speed of light, and the temperature of the sun (here we set as 5,800 K), respectively.

Similarly, the three recombination rates among the VB, IB, and CB are given as follows:

Recombination rates:

Here again under the assumption of Boltzmann distribution, the Equation (11) can be further approximated as follows:

Following the same treatment, the other two recombination rates can be given as follows:

Here , and are defined as recombination coefficient which have explicit forms as follows:

2.1.4. IB Continuity Constraint

After we have described all the simulation details regarding the IBSC modeling, we then move to tackle the issue of IB continuity constraint. As we have mentioned before, this is the key concept in this work because this constraint has not been adequately addressed since the first report about the drift–diffusion simulation of IBSC. We must mention here that the other IB dynamics such as carrier mobility in IB and the current contribution from IB which are important factors but are not addressed in the current work for simplicity and for better analyzing the effect of IB constraint conditions.

The governing equation for this constraint is given as below:

Equation (18) could be rearranged for easier interpretation as follows:

By refereeing to the sketch shown in Figure 1(a), this constraint states that the “net” generation rate between the VB and IB must be the same as the “net” generation rate between the IB and CB. In other words, there must have no surplus carriers in IB, namely, the carrier excited from the VB to IB will be reexcited completely to the CB without any loss. We must notice again that the key concept in this work is to discuss the theoretical IB continuity constraint in ideal circumstance, which is known to be greatly different from the practical application.

In this work, we handle the IB continuity constraint as follows to ensure the rigorous satisfaction of IB continuity constraints. By following the Equation (19), we rewrite it in the following form:

This formula can be further expanded as follows by using Equations (13) and (14), notice we have omitted the position variable for easier reading.

Equation (21) can be further reshaped as follows:

By multiplying at both sides of Equation (22), we arrive at the following form for Equation (22):

This is actually a quadratic equation in terms of , thus we could easily obtain its analytical solution as follows:

Here, the negative solution is removed due to the positive feature of the exponential function. Equation (24) is extremely important since it provides an explicit form for simulating the quasi-fermi level of IB as follows by assuming:

We thus have

And finally, the can analytically calculated as follows:

By inserting shown in the Equation (27) into the Equation (2) for calculating the carrier concentration in IB, the can be calculated accordingly in the following form:

It is worth noting that can be analytically calculated in terms of , in a manner similar to Equation (27):

Here is given as follows:

Finally, we have a corresponding expression for the carrier concentration in terms of . Both methods are valid for deriving , but here we have chosen to use the form given in Equation (28).

2.1.5. IB Continuity Constraint Violation Design

In this work, we intentionally designed the violated IB continuity constraint in order to compare the difference between the cases where the IB continuity constraint is satisfied and the IB continuity constraint is not satisfied. The easiest way to realize violated IB continuity constraint condition is probably not to include the IB continuity constraint in the simulation. However, the computational model we are using does not allow us to remove the IB constraint. Thus, by following Equation (20), we express the violated IB continuity constraint as follows:

Here, and are violation coefficients of adjustment so that the “net” generation rate between VB and IB does not equal to the “net” generation rate between IB and CB. In this work, we determined the respective values to be and , where is created as an artificial condition for unsatisfied IB current constraint. Expanding the Equation (31) in the same way as in the Equation (20), the quasi-fermi level of IB is obtained as follows:

Similar to Equation (28), by inserting shown in Equation (32) into Equation (2), the carrier concentration in the violated IB continuity constraint can be calculated. Finally, using derived from Equation (32), the calculation converges through the self-consistent manner with the IB continuity constraint not being satisfied.

2.2. Device Structure and Material Parameters

After we have described the detailed theory employed in this work, we move further to provide the information regarding the device structure and material parameter. We adopted a device geometry with the IB located in the region. Information regarding the thickness and doping concentration for emitter layer, i-regime and base layer are shown in Figure 1(b). In order to have a practical contribution to the experimental study of IBSC, the most widely investigated InAs/GaAs quantum dot-based IBSC was used as a test bed except the IB level was not theoretically simulated but manually adjusted to the ideal (detailed balance) energy level for generating maximum power with the . Due to the special treatment of IB, there is no explicit structure information about the InAs/GaAs quantum dot array. These parameters along with the absorption coefficients profile listed below are required to calculate the generation and recombination rates shown in Equations (7)–(14):

Here, the overlap of the absorption coefficients among various transitions is ignored. Meanwhile, the Heaviside profile used here was chosen for its simplicity. More complicated absorption profiles can be found in our previous work [20]. Moreover, light management, specifically the suppression of reflectivity on the device surface, is crucial in the practical applications [21, 22]. However, for the sake of simplicity, we have assumed an ideal optical condition where the reflectivity from the device surface is set to zero. Two doping conditions in IB were considered in this work where the IB region is either nondoped or the IB region is predoped with a donor concentration of and all dopants are assumed to be completely ionized. The effective density of states for the three bands: VB, IB, and CB for the GaAs are set as follows: , which are required to calculate the three types of carrier concentration . The band diagrams and conversion efficiency were simulated under three different light concentration condition with . The initial values of the parameters required for conducting self-consistent drift–diffusion-based IBSC simulation are summarized in Table 1 for better reference. All simulation models are constructed using the semiconductor modules in COMSOL Multiphysics and Fortran 90 codes.

3. Results and Discussion

3.1. Band Potential Profile Comparison

Figure 2 shows the simulated band diagram for the proposed IBSC under thermal equilibrium condition, namely, the dark condition using the device structure and material parameters described in previous section. Figure 2(a) represents the IBSC without IB doping simulated under condition where the IB constraint is satisfied. The comparative results for the same IBSC condition but with violated IB constraint are shown in Figure 2(b). As shown in the two figures, the band diagrams show typical feature and are completely identical to each other. The complete consistence is mainly due to the following two reasons. First, under dark condition, all the generation rate described in Equations (7)–(9) become zero. Second, under the thermal equilibrium condition, the three quasi-Fermi levels , and are completely overlapped with each other, which causes the three recombination rates described in Equations (12)–(14) to become zero. Due to these two factors, the IB continuity constraint defined in Equation (18) is constantly satisfied, thus there is no difference between these two figures. The band potential profiles of the IBSC with doped IB are shown in Figure 2(c) for the satisfied IB constraint and in Figure 2(d) for violated IB constraint. Different from the type potential profile for IBSC with nondoped IB, IBSC with doped IB showed typical potential profile and the flattened IB profile are solely due to the increased electron-carrier density at the IB region. Meanwhile, the explanation carried out for the potential profile equivalence between the IBSC with nondoped IB shown in Figures 2(a) and 2(b) could also be applied to the band potential profile of the IBSC with doped IB, which are identical for both the satisfied IB continuity constraint shown in Figure 2(c) and violated IB continuity constraint shown in Figures 2(d).

Dramatic differences are identified for the IBSC with satisfied IB constraint and violated IB constraint when the IBSC device is under light illumination. Since the system under nonequilibrium state, three quasi-Fermi levels , and start to split and finally end up at different energy potentials. Figure 3 shows the short-circuit band potential profile under 1-sun illumination for the nondoped IBSC with satisfied IB constraint in Figure 3(a) and violated IB constraint in Figure 3(b). It can be clearly seen here that unlike the thermal equilibrium dark state, the two IBSCs under 1-sun illumination showed large discrepancy between each other. For the case where the IB continuity is satisfied, the band potential profile maintained the similar feature as the dark state shown in Figure 2(a). In contrast, for the case where the IB continuity is violated, the band potential changed from the feature to a featured profile. The change in band profile can be further understood through the carrier concentration distribution shown in Figures 3(c) and 3(d). Compared to the satisfied IB constraint, significant increase of the hole concentration and IB electron concentration are found in the original region of the device, which cause the disappearance of the depletion region. The dramatic band potential profile difference between the satisfied IB constraint and violated IB constraint can be further understood by monitoring the two numerical simulation deviation errors, which are defined here by reformulating Equation (31) as follows:

For the case of satisfied IB constraint, we found that the deviation simulation error , as shown in lower figure in Figure 3(e), is at the magnitude below the order of which is consistent with our preset simulation error tolerance. However, for the violated IB constraint as shown in the lower figure in Figure 3(f) the simulation deviation error remains at the magnitude with the order of for most of the IB region. By plotting in the same figure of , it can be clearly seen that deviation error is completely negligible when compared to the deviation error . Meanwhile, from the upper graphs shown in Figures 3(e) and 3(f), it is found that greatly increased recombination rate in the is the main cause for the dramatic difference between and . Moreover, as shown in the Equation (14) for calculating the recombination rate , the increase of is mainly originated from the increase of the quasi-Fermi level splitting between and . By considering the Fermi level is maximized at around the energy level of IB, the increase of can only be achieved through lowering the , which inevitably cause the increase of the hole carrier concentration , as shown in Figure 3(d).

Similar to the nondoped IB simulation results shown in Figure 3, dramatic difference was also found for short-circuit band potential profile under 1-sun illumination for the doped IB simulation results and was summarized in Figure 4. The band potential profile for satisfied IB constrain maintained almost the same feature as the one under dark condition shown in Figure 2(c). The band flat profile shifted toward lower energy due to the fact that the electron concentration is about two orders higher than the hole carrier concentration as shown in Figure 4(c). We also found in here that the electron in the IB is depleted greatly at the interface with emitter layer, which induces strong electric field while damping the electric field in original device without doping (see Figures 3(a) and 3(b) for better understanding). We note that the band potential profile shown in Figure 4(b) for violated IB constraint holds a complicated nonlinear feature modified from the original linearized feature. We attribute this to the high electron doping concentration at the IB region as shown in Figure 4(d), where the IB carrier concentration reaches the maximum across the whole IB region. The depleted electrons from IB to emitter layer and base layer are much less than the one shown in Figure 4(c), which causes the reduction of the electric field near the IB edge and reappearance of the electric field originated from the original feature in the device without doping. The violated IB constraint results in the depletion region across the whole IB region, namely, an electric field in the region.

The validity of the simulation results for the satisfied IB constraint and violated IB constraint could be further verified by monitoring the constrain deviation error and shown, respectively, in the Figures 4(e) and 4(f). The is found within the order of , which is compatible with our simulation error tolerance. On the other hand, the violate IB constrain deviation error is at the range above the order of . By plotting in the same figure with , it can be clearly seen here that is completely negligible. We have also performed the simulation for the IBSCs under higher light illumination intensities as and . Since the band potential profile maintain the similar feature as the one for except a slightly increased quasi-Femi levels, the plots are not presented here due to the limit of space. However, the device performance dependence on light illumination intensity will be given in subsequent section for better comparison.

3.2. Net IB Carrier Generation Rates Comparison

We define the net IB generation rates as the difference between the generation rate and recombination rate via IB as follows:

Figure 5 shows the net-generation rates IB as a function of position for IBSCs with doped and nondoped IB under various light concentration ratios: and . Figure 5(a) shows the simulation results with the satisfied IB constraint while Figure 5(b) presents the one with violated IB constraint based on the Equation (31). All the results have been normalized by the light concentration ratio . As shown in Figure 5(a), the net generation rate is completely identical for all the simulation under either different doping concentration in IB or the different light concentration ratio. These results again verified the validity for our simulation results where a IB continuity constraint is rigorously guaranteed. To the best of our knowledge, this is the first work of the same kind which succeeded in a complete realization of IB continuity. On the other hand, as shown in Figure 5(b), the results simulated with the violated IB continuity constraint display totally different tendency from Figure 5(a). The low between in IB region has been attributed to the low IB concentration and high IB concentration between and the edge of base layer is due to the increased IB concentration as shown in Figure 2(b).

3.3. Electric Field, Current, and Power Device Distribution

In this section, we provided a brief description of the electric field, current, and power distribution within the device. The electric field distribution, simulated under a 1-sun short-circuit condition for the IBSC with satisfied IB constraints, showed a significantly higher value for the doped IB compared to the nondoped IB, as illustrated in Figure 6(a). This observation agrees with the band profiles depicted in Figures 3(a) and 4(a). This correlation can be easily validated by taking the first derivative of the results displayed in these two figures. Conversely, for the IBSC modeled with violated IB constraints, the electric field was much higher for the nondoped IB than for the doped IB. This trend is the exact opposite of what was observed with the satisfied IB constraints. Similarly, this electric field distribution aligns with the profiles shown in Figures 3(b) and 4(b). For the current and power distribution depicted in Figure 6(b)6(d), it is evident that the current distribution profile aligns with the power distribution profile, as illustrated by the solid and dotted lines in Figure 6(d) for simulations under satisfied IB constraints. On the other hand, with the violated IB constraints, a significant deviation is observed between the current and power distributions. Since the power distribution profile corresponds to the product of current and voltage, these results suggest that the simulation for the violated IB constraints might not have been accurately conducted when the device was under an external bias voltage. Such findings hint at potential poor device performance; for instance, a higher shunt resistance might be present in the device simulated under violated IB conditions compared to the one under satisfied IB constraints. We will address this further in the subsequent section.

3.4. Current–Voltage Characteristics and Light Concentration Dependence

At last, we show the current–voltage characteristics () curves for the IBSCs with nondoped and doped IB under light concentration ratio and . Notice that the current density is normalized by the light concentration ratio for fair comparison. Figure 7(a) shows the results simulated under the condition where the IB continuity constraint is rigorously satisfied. The open circuit voltage is found increased when the illumination concentration ratio increase, which agrees with the conventional solar cell device. It is also found here that the simulated under various conditions are quite close to each other. However, the small difference can be clearly identified through the close-up of the region near the region as shown in the inset of Figure 7(a). As a general tendency, we found here the for the IBSC with doped IB showed slightly higher value than the one without predoped IB. We attribute the increase to the large carrier concentration doped in IB. However, by carefully checking the value shown in Figure 5(a), we found the net-generation rate in IB for the doped IB is actually lower than the one with nondoped IB in most of the region of IB except in the narrow region near the emitter layer. A consistent feature in the curves with the one shown in Figure 7(a) is that again we found here the for the IBSC with doped IB showed higher value than the one without predoped IB.

In contrast to the results where the IB continuity constraint is strictly adhered to, the simulation under the violated IB constraints, as shown in Figure 7(b), reveals an unstable current–voltage profile. As highlighted in the conclusion of the previous section, depicted in Figure 6(d), this instability in the current–voltage relationship can be attributed to the improper conduction of the current under external bias when the IB constraints are not met. This is consistent with the elevated shunt resistance observed in Figure 7(b).

We speculate here that for the IBSC simulated under rigorous satisfaction of the IB continuity constraint, the device performance behaves similarly to the device simulated under the detailed balance simulation condition. Under the detailed balance principle, the normalized showed no dependence on the concentration ratio . Figure 7(c) shows the detailed simulation results regarding the device performance dependence on the light concentration ratio . It can be clearly seen here that the normalized short circuit currentof showed no or little dependence on in both the cases with doped IB and nondoped IB, which is fairly consistent with simulation results based on detailed balance principle-based simulation. The open circuit voltage showed logarithmic relation with the concentration ratio , which again aligns with detailed balance principle. The power conversion efficiencies which are ∼32% at 1 sun and ∼39% under 1,000 suns, showed also a nearly logarithmic relation with except some deviation appeared in the results for the device with doped IB. The slight deviation away from logarithmic relation is probably due to the poor fill factor for the device at low concentration and recovered due to the so-called photon filling effect at higher light concentration ratio. We must note here that although the device performance presents similar feature as the device operation under detailed balance principle, all the conversion efficiencies simulated under drift-–diffusion model are found dramatically lower than the maximum theoretical efficiency simulated under detailed balance principle. On the other hand, the simulated efficiency efficiencies which are ∼32% at 1 sun and ∼39% under 1,000 suns are still higher than the experimentally reported results of InAs/GaAs QD IBSC, indicating a lot of rooms to improve in the future experimental work [11, 23]. At last, it is worth mentioning that when using the stratified IB continuity in our simulations, the efficiency obtained is comparable to that achieved using a hybridized detailed balance model and drift–diffusion models, as shown in literature [20]. As noted in the introduction, the hybridized method is suitable for examining effects such as photon filling, IB mobility, nonradiative recombination, and light concentration. However, the numerical simulation in the hybridized approach simplified the IBSC structure to compute the ideal conditions of an IBSC, assuming a flat i-layer. Additionally, both the p-layer and n-layer were omitted from the device model, limiting its practical application in the practical IBSC device design.

4. Conclusion

In this work, we performed extensive simulation to fully disclose the significant effect of the IB continuity constraints on the IBSC device characteristics. The band potential profiles for the IBSC with doped and nondoped IB under the dark state without light illumination are identical to each other because all the generation and recombination rates in the IB are zeros and IB continuity constraint is automatically satisfied. However, we found dramatical difference in the band potential profile occurred between the satisfied IB constraint and violated IB constraint when the device simulation is performed under light illumination. For the case where the IB constraint is satisfied, we found the band potential profile maintained basically the same feature as the one under dark conditions for IBSCs with both the doped and nondoped IB. On the other hand, the simulation results with violated IB constraint showed inappropriate band profile with dramatic deviation from the dark condition. Our results based on satisfied IB constraint were further verified by the deviation simulation error lower than the order of , which is completely negligible when compared to the order of for the case of violated IB constraint. Meanwhile, we also found the normalized net generation rates showed completely different features between the case with satisfied IB constraint and violated IB constraint. From the simulated J–V curve and light concentration ratio dependence, we found that the device performance based on drift–diffusion under rigorous IB continuity constraint showed similar tendency to the features simulated based on detailed balance principle except the much-lowered power conversion efficiency. Our work demonstrated for the first time the appropriate band potential profile and the device performance for drift–diffusion based IBSC device simulation. All the results verified in this work will have significant impact on the proper device design of IBSC for better IB material candidate screening, efficiency improvement, optical management for both the III–V compound-based IBSC as well as the emerging state of the art perovskite solar cell, perovskite-based IBSC, and other state-of-the art silicon-based photovlotaics [2428].

Data Availability

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Conflicts of Interest

The authors declare that there is no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors gratefully acknowledge the funding from New Energy and Industrial Technology Development Organization (NEDO) Grant number: JPNP20015 and the Ministry of Economy, Trade and Industry (METI), Japan.