# Thin Film Solar Cell¶

Thin-film silicon solar cells rely on scattering structures diffracting incoming light into larger angles thus trapping the light by total internal reflection in the high-index absorber material. Their optimization is a topic of current research. In general the cell is modeled by applying periodic boundary conditions even if the scattering structures are statistical in nature as showcased in this example. This modelling choice accounts for light scattered out of the sides of the unit cell while possibly introducing periodification artifacts. This example relies heavily on the scripting interface describe in Matlab Interface.

This example features a model solar cell with a glass half-space from which the light is incident and a mirror boundary condition on top. A layer with the properties of a transparent conductive oxide (TCO) and an absorber layer mimicking silicon complete the layered solar cell geometry shown in the following.

Geometry Definition and Meshing for a Simple, Random Textured Thin-Film Solar Cell

The geometry is defined in the template file `data_analysis/layout.simple_setup.jcmt`

. It relies on the input defined in the script `data_analysis/run_wavelength_scan.m`

. The following section of `data_analysis/run_wavelength_scan.m`

loads the random interface defined as x-y-pairs in `data_analysis/interface.txt`

and sets most of the parameters required, such at the thicknesses of the TCO and absorber layers.

```
%% Setup solar cell geometry
% load interface and shift to zero integral to keep material volume
% constant
interface = load('interface.txt'); % simple list of x y pairs in keys.uol scaling
keys.pitch = diff(interface([1 end],1));
shift = trapz(interface(:,1),interface(:,2))/diff((interface([1 end],1)));
interface(:,2) = interface(:,2)-shift;
interface(:,1) = interface(:,1) - interface(1,1) - keys.pitch/2;
keys.interface = interface.';
keys.thicknesses = [ 600 ;... % tco thickness
300 ]; % silicon thickness
keys.h_offset = max(abs(interface(:,2)))+100;
```

These parameters are further processed in the scripting section at the beginning of `data_analysis/layout.simple_setup.jcmt`

. Here the appropriate key-value pairs are set up to be used throughout the template.

The cell geometry is built from three polygons: The computational domain, the TCO domain and the silicon domain. The computational domain is defined as a rectangular polygon.

```
# Computational Domain
Polygon {
Name = "ComputationalDomain"
Points = %(currPoints)e
DomainId = 1
Priority = -1
Boundary {
Number = 1
Class = Transparent
ExteriorDomainId=1
}
Boundary {
Number = 2
Class = Periodic
}
Boundary {
Number = 3
Class = DomainBoundary
BoundaryId=1
}
Boundary {
Number = 4
Class = Periodic
}
```

For this polygon `DomainId = 1`

and `Priority = ComputationalDomain`

are used. The specific value of `Priority = ComputationalDomain`

specifies that domain has a smaller priority than all overlapping objects. However, all parts of objects outside of the boundaries of this object are cut off / clipped. Furthermore, all boundary conditions specified on this boundary are inherited by overlapping objects with a higher priority. In this case we assign a `Class = Transparent`

boundary to the first `Boundary`

which are numbered counter clock-wise. The second and fourth boundary segments are set to `Class = Periodic`

. At the top of the model, a perfect electric conductor is assumed, i.e. a perfect loss-less mirror. This is modeled by a `Class = Domain`

boundary condition which is addressed by the specific `BoundaryId = 1`

. The further modelling of this boundary condition is explained in the section Boundary Conditions.

The remaining two polygons are built through a loop in the scripting block in `data_analysis/layout.simple_setup.jcmt`

. The `Name`

, `Priority`

, and `DomainId`

are simple to determine. The TCO layer, like the computational domain, is constructed as a simple quadrilateral whereas the silicon layer includes the textured interface which requires a little more work to account for the thickness offsets of the previous layers.

```
# Layer Stack
<?
for layer=2:3
keys.mat_name=sprintf('Material_%d',layer);
keys.mat_id = layer;
keys.slc = .5*min(max(real(keys.vacuum_wavelength/keys.uol/keys.n(layer)),100),real(keys.vacuum_wavelength/keys.uol/keys.n(layer)));
if layer>2
dummy = keys.interface;
dummy(2,:)=dummy(2,:)+offsets(layer);
dummy(:,end+(1:2)) = [xLimits(2) xLimits(1); ymax ymax];
keys.currPoints = dummy(:);
else
keys.currPoints = [ xLimits(1) offsets(layer) xLimits(2) offsets(layer) xLimits(2) ymax xLimits(1) ymax];
end
?>
Polygon {
Name = "%(mat_name)s"
Points = %(currPoints)e
DomainId = %(mat_id)d
Priority = %(mat_id)d
MeshOptions{
MaximumSideLength = %(slc)2e
}
```

Note that we have adapted the mesh options for all domains to respect the wavelength as well as the material properties. This is done to resolve the wave phenomena accurately enough for the derived post processing quantities to remain meaningful.

Boundary Conditions

The cell modeling involves three types of boundary conditions available for electromagnetic modeling in `JCMsuite`

. Here a `Transparent`

boundary is used at the bottom modelling an infinite glass half space. The cell is set up to be a unit cell of an infinitely periodic array in X direction. At the top of the model, a perfect electric conductor is assumed, i.e. a perfect loss-less mirror. In the `Layout2D`

definition above these boundary conditions have already been assigned to the geometry. In case of the `Transparent`

and `Periodic`

setting this does not require any extra work. To model the perfect electric conductor correctly, further definitions are required. These are set in the `data_analysis/boundary_conditions.jcm`

and the `data_analysis/sources.jcmt`

files.

In `boundary_conditions.jcm`

a `BoundaryCondition`

is defined. On the specified `BoundaryId`

the electromagnetic field is set to include only a `TangentialElectric`

component.

```
BoundaryCondition {
BoundaryId = 1
Electromagnetic = TangentialElectric
}
```

In `sources.jcmt`

the source definition must then include a value for the `ElectricFieldStrength`

at the boundary. This is set to `[0.0 0.0 0.0]`

.

```
Source {
BoundaryId = 1
ElectricFieldStrength = [0.0 0.0 0.0]
}
```

Note

In case multiple `SourceBags`

are used, sources at the boundaries must be defined in every `SourceBag`

.

Source Definition

The illumination in this example is modeled as a z-polarized plane wave incident from the lower half space. This is considered a good enough approximation to solar illumination on earth at the length scales involved. As the illumination wavelength is parameterized, the template contains a scripting section which sets the appropriate keys for the k-vector `K`

and the `Amplitude`

. Note that in this a more simple definition of the `PlaneWave`

source is also possible. For details see Plane Wave

The special setting of a perfect electric conductor in this example requires additional care of the `Source`

definition. As described in the section Boundary Conditions an additional `Source`

on the boundary must be defined. Furthermore, this step necessitates the declaration of the source domain identified by `DomainId = 1`

which for problems without `DomainBoundaries`

the solver is able to determine automatically.

```
<?
keys.k_vector = 2*pi/keys.vacuum_wavelength*real(keys.n(1))*[0 1 0];
keys.U_s = [0 0 1];
?>
SourceBag {
Source {
DomainId = 1
ElectricFieldStrength {
PlaneWave {
K = %(k_vector)e
Amplitude = %(U_s)e
}
}
}
Source {
BoundaryId = 1
ElectricFieldStrength = [0.0 0.0 0.0]
}
}
```

Material Definition

The material properties have crucial influence on the results of the simulation. There are many different material databases available on the internet and in the literature. In this example we employ the freely available data sets from http://www.refractiveindex.info/ for the silicon material data. For glass we assume a constant refractive index of and for the TCO we assume . These are initialized in the `run_wavelength_scan.m`

script.

```
% refractive indices
keys.n(1) = 1.52; % glass
keys.n(2) = 2.1+.01*1i; % tco
keys.n(3) = 1; % silicon
```

To account for the dispersive nature of the material parameters, the material data is interpolated for the corresponding wavelength in the `run_wavelength_scan.m`

script. It is imported into Matlab??

```
% load material data for Si from file
Si = importdata('Si.txt',' ', 13); Si = Si.data;
Si(:,1) = Si(:,1)*1e-9; [foo,I] = unique(Si(:,1)); Si = Si(I,:);
```

and later interpolated for the corresponding wavelength.

```
keys.n(3) = interp1(Si(:,1), Si(:,2), keys.vacuum_wavelength, 'linear') ... % n
+ 1i*interp1(Si(:,1), Si(:,3), keys.vacuum_wavelength, 'linear'); % k
```

Note

A linear interpolation might not be a good choice, e.g. for coarse data sets with exponentially decaying imaginary parts of the refractive index. Here, a linear interpolation might overestimate absorption.

The template `materials.jcmt`

contains a very simple loop over all materials in which a `Name`

is derived from the `Id`

and the `RelPermittivity`

is computed from the refractive indices set in the script `run_wavelength_scan.m`

.

```
<?
for mat = 1:length(keys.n)
keys.permittivity = (keys.n(mat))^2;
keys.mat_id=mat;
keys.mat_name = sprintf('Material_%d', mat);
?>
Material {
Name = "%(mat_name)s"
DomainId = %(mat_id)e
RelPermittivity = %(permittivity)e
RelPermeability = 1.0
}
<?
end
?>
```

Project Setting

The `project.jcmpt`

files contains the solver settings required to solve the previously set up model. It is solved as a time-harmonic electromagnetic scattering problem.

```
Project {
InfoLevel = 3
Electromagnetics {
TimeHarmonic {
Scattering {
Accuracy {
FiniteElementDegree = %(fe_degree)i
Precision = %(precision)e
Refinement {
MaxNumberSteps = %(n_refinement_steps)i
PreRefinements = %(n_prerefinement_steps)i
}
}
}
}
}
}
```

To steer `Accuracy`

several settings are made in `run_wavelength_scan.m`

. A `Precision`

of 0.01 and a `FiniteElementDegree = 2`

is used. One adaptive refinement step is taken and assures a reasonably accurate simulation result.

```
% solver settings
keys.n_refinement_steps = 1;
keys.n_prerefinement_steps = 0;
keys.precision = 1e-2;
keys.fe_degree = 2;
```

Post Processing

The post processing steps required to compute absorption and reflection of the solar cell model are defined in `project.jcmp`

. The absorption can be derived from either the `DensityIntegration`

PostProcess or the `FluxIntegration`

PostProcess. The `DensityIntegration`

PostProcess integrates the `ElectromagneticFieldAbsorption`

quantity. This gives rise to the power loss within the absorber layer due to an imaginary part of the permittivity. As the glass block is assumed to be transparent, it would be sufficient to restrict the computation to `DomainIds = [2,3]`

.

```
PostProcess {
DensityIntegration {
FieldBagPath = "./project_results/fieldbag.jcm"
OutputFileName = "./project_results/field_energy.jcm"
OutputQuantity = ElectromagneticFieldAbsorption
}
}
```

Alternatively, the absorption can also be computed by looking at the net fluxes of the `ElectromagneticFieldEnergyFlux`

across domain boundaries. This also allows to compute the reflected flux leaving the computational domain.

```
PostProcess {
FluxIntegration {
InterfaceType = DomainInterfaces
FieldBagPath = "./project_results/fieldbag.jcm"
OutputFileName = "./project_results/energy_flux.jcm"
OutputQuantity = ElectromagneticFieldEnergyFlux
}
}
```

The Fourier Transform of the outgoing electromagnetic field allows for a more detailed account of the direction of the scattered light. It is also included in the `project.jcmpt`

file.

```
PostProcess {
FourierTransform {
FieldBagPath = "./project_results/fieldbag.jcm"
OutputFileName = "./project_results/propagating_fourier_coefficients_down.jcm"
NormalDirection = -Y
}
}
```

Simulation Run

The script `data_analysis/run_wavelength_scan.m`

sets up the simulation and executes the wavelength scan. Excerpts have already been discussed in the sections above. The central part of the script is the following loop over the different wavelengths.

```
for ii = 1:length(wavelengths)
keys.vacuum_wavelength = wavelengths(ii)*keys.uol;
keys.n(3) = interp1(Si(:,1), Si(:,2), keys.vacuum_wavelength, 'linear') ... % n
+ 1i*interp1(Si(:,1), Si(:,3), keys.vacuum_wavelength, 'linear'); % k
%% Solve
results = jcmwave_solve('--solve', 'project.jcmp', keys,'jcmt_pattern','simple_setup');
%% Absorbed energy
absorption(:,ii) = results{2}.ElectromagneticFieldAbsorption{1}/p_in;
%% Fourier Transform
e_field = cell2mat(results{3}.ElectricFieldStrength);
km = norm(results{3}.K(ceil(end/2),:));
angles = real(acosd(-results{3}.K(:,2)/km));
refl = sum(conj(e_field).*e_field,2).'*cosd(angles);
reflection(ii) = sum(refl,1)*.5/Z0*keys.n(1)*keys.pitch*keys.uol/p_in;
%% Flux outgoing
energy_flux = cell2mat(results{4}.ElectromagneticFieldEnergyFlux);
index = results{4}.DomainIdFirst == 2 & results{4}.DomainIdSecond == 3;
absorption_via_flux(ii) = real(energy_flux(index))/p_in;
index = results{4}.DomainIdFirst == 2 & results{4}.DomainIdSecond == 1;
reflection_via_flux(ii) = real(p_in+energy_flux(index))/p_in;
%% Statistics
total_time(ii) = results{1}.computational_costs.TotalTime(end);
unknowns(ii) = results{1}.computational_costs.Unknowns(end);
mem(ii) = results{1}.computational_costs.TotalMemory_GB(end);
end
```

Within this loop the variable `keys.vacuum_wavelength`

is set to the correct value and the interpolation of the material data is executed. Subsequently the project is solved by the appropriate call to `jcmwave_solve`

with the keys structure and the `jcmt_pattern`

as additional arguments.

The obtained results from the post processes are returned in the cell array `results`

in the order defined in `data_analysis/project.jcmpt`

. A scaling with of the electromagnetic field absorption in the absorbing layers gives rise to the absorbed fraction of incident power to the domain.

In similar ways, the output of the `FluxIntegration`

and `FourierTransform`

are processed and scaled to obtain the reflected and absorbed fractions of the incident power.

Finally the simulated results are plotted as implied EQE and 1-reflection over wavelength and as an area plot with absorptance in the difference cell layers as well as reflections. The `sampling_rate`

variable allows to refine the sampling rate of the `wavelength_interval`

.

The final section of the script executes the weighting with the AM1.5 solar spectrum to compute implied photo current densities. To compute this, the spectrum is first loaded from a *csv* file *data_analysis/ASTMG173.csv* which is a slightly reformatted version of the file *ASTMG173.xls* available at http://rredc.nrel.gov/solar/spectra/am1.5/. From the data in this file the photon count per second is computed and the `spectral_weight`

determined at the simulated wavelengths. Numerical integration and scaling with the elementary charge allows to display the simulated implied photo current densities for the silicon absorber as well as the parasitic losses in the TCO and the reflection losses.

```
%% Solar Spectrum
solar_spectrum_data = importdata('ASTMG173.csv',';',2); % Data is reformatted from ASTMG173.xls available at http://rredc.nrel.gov/solar/spectra/am1.5/
planck = 6.62606957*10^-34; % Planck's constant
energy = planck*c0/keys.uol./solar_spectrum_data.data(:,1); % photo energy
photon_number = solar_spectrum_data.data(:,3)./energy*1e-4; % number of photons in cm-2 nm^-1 s-1
spectrum_weight = interp1(solar_spectrum_data.data(:,1),photon_number,wavelengths);
eC = 1.602176565*10^-19; % elementary charge
currentDensity = eC*trapz(wavelengths,[absorption.*repmat(spectrum_weight,3,1)]')*1000;
currentDensity(3) = eC*trapz(wavelengths,reflection.*spectrum_weight)*1000;
fprintf('\n Implied photo current densities in (%dnm, %dnm):\n %g mAcm-2 (TCO)\n %g mAcm-2 (Si)\n %g mAcm-2(R)\n\n',round(wavelengths([1 end])), currentDensity)
```