Solar cell simulation with angled incidence


@aalam Dear aalam, we adopted your provided EQE scripts in our structures successfully and so far, it´s working nice. Thank you very much for your effective support. BTW, I am gonna ask another question. We are trying to simulate solar cells at different incident angles (0 to 90 degree) correctly. We did use BCs as (Periodic) in X & Y axis and (Metal & PML) in Z axis, and did a sweep in the range of 0 to 90 degree incident angles. We are happy with the results so far but the problem is that we are not sure about the accuracy of the simulation results. Because, the FDTD absorption curve didn´t match with the analytical curves except the case of incident angle at 0 degree. To verify our results, we did two trails (incident angle 60 degree and polarization angles:0 & 90) and in that cases we used BCs as “BFAST”. In the case of BFAST, the FDTD absorption curve did close match with the analytical one but the absorption curve in the silicone layer (Absorber layer) was totally messed up. So, still we couldn’t made proper investigation of the accuracy of our simulation results. I am hoping you understand the our problem. Please, could you help me regarding this issue? How could we can be ensured that our simulation results (values of Jsc) that we got from the sweep (Conditions:= incident angles 0 to 90 degree, interval:5 degree, BCs: Periodic) are correct. If you need any more information or one of simulation file I can give you. N.B: one more thing, BFAST also takes so much time to complete the simulation. I am sorry for long writing. I am looking forward your help! Thanks in advance!

Script commands to calculate EQE of solar cell as a function of wavelength

Hi @mj.mendes, Thanks for letting me know that the script is working for you. The periodic boundary condition is valid only for normal incidence (of plane wave sources). For angled incidence the fields will not be periodic (there will be a phase difference between each period). You should use Bloch BC instead for single frequency simulations, or BFAST for broadband simulations. This is why your simulation with periodic BC only matches analytic result at 0 degree whereas the result from using BFAST will match analytic result for all angles (BFAST not Bloch since you are doing broadband simulation).

Can you please elaborate on why you feel the absorption curve in silicon is messed up when at the beginning of the sentence you mentioned that the absorption curve is in agreement with analytic calculations?

If you can share a simulation file then that would be great. Also if you can share the data you are trying to match with (I am not sure if you have any measured data though) then that will be even better.

Unfortunately this delay in simulation is the price you pay for accuracy in the result. The BFAST simulation will take more and more time than simulation with regular plane wave source as the angle of the place wave increases. Please take a look at this KB page to learn about different simulation tips about simulation with BFAST (


Dear @aalam , I do apology for not well explained of our problem. For your kind understanding, I attached one of our simulation file along with few snapshots of the results of this simulation in a pdf file. I also enclosed a script file that we use to see and extract the results.
However, the main problem that we faced while using BFAST is that “Solar Generation” was not well aligned with the Absorber Layer (Silicon) and the data that created was diverged. As a result, Solar Generation was providing us weird Jsc value. You can easily see from the snapshots of the results.
Furthermore, I have read BFAST page from your given web-link but still I couldn’t figure out very well how could we can run a sweep in a larger scale (i.e: incident angles 0 to 90 degree with having 5 degree interval) in BFAST mode relatively more fast. Because of I already did use such values (dt multiplier: 0.5 to 1, auto shutoff 1e-06 to 1e-05 etc.) which are mentioned in the page but still it does take so much time to complete only one single simulation. Due to this fact, I am little concerned about time consuming of running a sweep in a BFAST mode as well as memory also can be another issue.
I really appreciate your effective support. Thank you so much! I am looking forward your valuable suggestions.
Script for analysis.lsf (2.0 KB)
Flat 300nm aSi_With ARC.fsp (358.6 KB)
Oblique Angle Study.pdf (216.0 KB)


Hi @mj.mendes, Just wanted to let you know that I am looking into your files. It looks like there is some issue with the convergence coming (I think) from the PML. I myself am not an expert in FDTD simulation with BFAST so I am consulting with my colleague @fgomez. We will hopefully get back to you by tomorrow. Thanks for waiting.


Hi @mj.mendes,

I discussed your simulation with @aalam and we found the issue behind the divergence in the simulation with BFAST. To simplify the simulation I made it quasi-2D by using only two cells along the y direction. Even though it seems like you want to include some structures afterwards, which won’t make this reduction possible, the simplified simulation illustrates the reason for the divergence.

If you run Flat 300nm aSi_With ARC - diverged.fsp (381.5 KB), the simulation will diverge and the log file will report some problems with the material fit:

Numerical discretization of material 'AZO_nk' might be unstable! Max. imaginary part of epsilon is 0.00101149!
Numerical discretization of material 'a-Si Fraun' might be unstable! Max. imaginary part of epsilon is 0.016145!
Numerical discretization of material 'Ag (Silver) - Palik (0-2um) Copy 1' might be unstable! Max. imaginary part of epsilon is 0.00856091!

These warnings tell you that the material ended up with some artificial gain (negative imaginary part of the permittivity) when fitting the material numerically for the BFAST calculation. Even though the original fit you see in the Material explorer does not show gain, the BFAST fit is slightly different because it depends on the time step dt. Therefore, the BFAST fit can end up with some small artificial gain that can cause an instability.

When you reduce dt by decreasing the BFAST dt multiplier, the artificial gain becomes smaller. For example, in Flat 300nm aSi_With ARC_converged.fsp (361.1 KB) I reduced the BFAST dt multiplier from the default 0.5 to 0.2 and the simulation converged.

To see more clearly the origin of the problem we can plot the material data, the FDTD fit from the Material Explorer and the fit used by BFAST. For the latter we can use the script function getnumericalpermittivity. You can simply run the script test_permittivity.lsf (1.2 KB) for the two simulations (which have different BFAST dt multiplier) to get the results shown below for the imaginary part of the permittivity near lambda=1100nm. Note how for smaller dt, the artificial gain decreases:

An alternative way to fix the problem is modifying the fit parameters until you find some settings where the BFAST permittivity does not have any artificial gain. You don’t need to run any simulations, as you can simply replot the BFAST permittivity using the method just described and check that there is no negative imaginary part. You can also add a small positive imaginary part to the original permittivity data to force the fit to move away from zero. This number can be close to the tolerance of the permittivity measurement, so that you do not increase the absorption in an unphysical way.

I have a couple of additional suggestions for your simulation:

  • Use a 2D simulation (or even a 1D simulation) to check that the planar structure is being simulated correctly. It is not only necessary to check that it is converging, you also want to check the accuracy of the results. For this purpose you can compare with results from the script function stackrt.
  • Because of the pronounced angle of incidence this simulation will be challenging once you move to 3D and add more structures. As explained here, the simulation time using BFAST increases with the angle of incidence. Furthermore, if you need to decrease the BFAST dt multiplier, the execution time will be longer. If this time becomes prohibitively large, you might want to consider sweeping over single-frequency simulations.


Dear @fgomez , Dear @aalam, Thank you so much for putting your lots of effort! We are thinking that using BFAST would be very complicated for all of our structures. So we are trying use bloch and run a 2D sweep to get a correct results. I will inform if we face any problem regarding 2D sweep in bloch mode. Have a nice time both of you!


Dear @aalam , Hoping, you are having a nice time. I need your some help regarding to the understanding of the results of 2D sweep in bloch mode. For your kind understanding, here I attached a pdf file in where I described the results and problems of two FDTD simulations. I am also uploaded those FDTD simulation files [one main (reference) file, and other one´s simulation file of 2D sweep]. If you need any more information, please let me know. Thank you so much for your effort! main simulation_Array TiO2 domes on 300nm aSi.fsp (356.2 KB)
bloch simulation_Array TiO2 domes on 300nm aSi.fsp (388.1 KB)
2D sweep with bloch FDTD simulation_MJM.pdf (575.9 KB)


Dear @mj.mendes,

Sorry for the wait. The reason why you are seeing different results is because the source bandwidth is different:

main_simulation: wavelength start = 0.4um and wavelenght stop = 1.1um
bloch_simulation: wavelength start = wavelength stop = 0.75um

Although in both cases the center wavelength is the same, the mesh will be different because when the “auto nonuniform” setting is used the mesh is calculated according to the smallest wavelength in the simulation bandwidth (which by default is the same as the source bandwidth). Therefore, the mesh in main_simulation is finer than in bloch_simulation, which probably means the results from main_simulation are more accurate.

If you set the source ranges to be the same in both simulations you get the same results because then the mesh is identical. If you want to have the same mesh without modifying the source settings you can use the advance option “set simulation bandwidth” in bloch_simulation:

This is recommended when running sweeps over single-frequency simulations as explained in this post:

One final suggestion regarding the PML: you should use the “steep angle” profile for the stretched coordinate PML as explained here:


Dear @fgomez , I am sorry for being late. Thank you so much for your extreme hard work that you did put for my simulations. Anyway, I couldn’t properly implemented your suggestion yet. At the moment, I am doing some other stuffs but I will again start angular simulations study in a month. May I am going to request you again for your some helps regarding this issue in coming days. Until then stay happy!