# Study of finite volume number density fluctuation at RHIC energies in PNJL
model for the search of QCD Critical Point

A. Sarkar<EMAIL_ADDRESS>P. Deb<EMAIL_ADDRESS>R. Bose
<EMAIL_ADDRESS>*‡Department of Physics, Ramakrishna Mission Residential
College (Autonomous), Narendrapur, Kolkata - 700103, India  †Department of
Physics, Indian Institute of Technology Bombay, Powai, Mumbai- 400076, India

###### Abstract

The Quantum Chromo Dynamical (QCD) phase diagram can be explored in the heavy
ion collision experiments. The non-monotonic behavior of the conserved number
fluctuations are believed to be a signatures of the QCD phase transition and
the Critical End Point (CEP) as a function of the collision energy. Due to the
experimental limitations the elementary constituents produced, are evanescent,
it is very difficult to measure the observables directly. So, we need
theoretical approaches to estimate the values of different observables. The
Polyakov-Nambu-Jona-Lasinio (PNJL) models of QCD is such an effective model,
which posses the benefit of having characteristics similar to the observables.
We present the number density fluctuation at RHIC energies and quantitative
estimates of its contribution of critical fluctuations in PNJL model.

###### pacs:

12.38.AW, 12.38.Mh, 12.39.-x

##  I Introduction

In order to study the early universe, we must have the knowledge of the
matter, that existed at that time. Upon reproducing the matter with the exact
density and temperature, existing at that time, in laboratory conditions, the
characteristics of the early Universe can be studied. So far, the only
possible way to do that is, to collide two heavy atomic nuclei, accelerated to
energies of orders of a hundred GeV, inside a collider. The Relativistic Heavy
Ion Collider (RHIC) at the Brookhaven National Library (BNL) and The Large
Hadron Collider (LHC) at CERN are the two biggest colliders for this. Using
the results of head-on collisions, it is possible to model the density and
temperature that existed in the first instants of the formation of the
Universe. Over the past few decades, physicists around the globe have been
trying to recreate the matter, named quark-gluon plasma (QGP), by colliding
together nuclei of atoms, with enough high energy, to produce trillion degree
temperatures. At these high temperatures, matter stays in the QGP state for a
fraction of a second.

When two particle beams collide head-on inside the tunnel of colliders, some
of the matter changes state from hadronic state to QGP state. In this chiral
phase transition, a first order phase transition may exist. The possible
existence of a first order phase transition and the smooth cross-over implies
that the transition line will end at a point. This point is known as the QCD
critical point or critical end point (CEP) STAR ; STAR1 . The CEP at the end
of first order line is a critical point where the phase transition is of
second order. With respect to the phase diagram, the critical point is the
point where the boundary between phases disappear. To obtain the critical
point from a phase diagram, physicists observe the enhancement of fluctuations
of final state observables, for example, the number of pions emitted in a
collision or the ratio of the numbers of kaons and pions emitted, higher order
moment of conserved quantity, multiplicity, etc..

The QGP matter formed in the heavy ion collision experiments has a finite
volume which depends on the size of the colliding matter, center of mass
energy and collision centrality. In the different centrality measurements of
HBT radii, several efforts have been made to estimate the finite volume of the
system adamova . These results suggest that the volume increases with
centrality during the freeze out of the system. The effect of finite volume
system has been studied in many theoretical models such as non-interacting bag
model elze , chiral perturbation theory luscher ; gasser , Nambu–Jona-Lasinio
(NJL) model nambu ; kiriyama , linear sigma model braun1 ; braun2 , Polyakov
loop extended Nambu–Jona-Lasinio (PNJL) model fukushima ; ratti ; pisarski ;
fukushima1 ; hansen ; ciminale ; ghosh ; deb1 ; abhijit ;
bhattacharyya2010investigation ; deb2009mesonic  and by the first principle
study of pure gluon theory on space time lattices bazavov ; boyd ; engels ;
fodor ; allton ; forcrand ; aoki ; megias . In the NJL model the finite size
of a dense baryonic matter has been described by the induction of a charged
pion condensation phenomenon. Recently in PNJL model, it has been observed
that as the volume decreases, the critical temperature for the crossover
transition decreases. For lower volumes, CEP is shifted to a domain with a
higher chemical potential $\mu$ and lower temperature $T$ deb1 . All the
theoretical results indicate that the location of CEP strongly depends on the
size of the fireball of the heavy ion collision system.

In this present work, our aim is to find out the critical end point of the QCD
phase diagram. To do so we will use an effective model - Polyakov-loop
extended Nambu-Jona-Lasinio (PNJL) model. We will focus on the calculation of
quark number density and will try to find the critical region from the quark
number densities and compare the results with the experimental beam energy
scan results tlusty2018rhic ; odyniec2019beam . In order to study the quark
numbers at finite temperature and chemical potential, it is found that the
quark-number density is determined by the corresponding dressed quark
propagator only zong . The quark number densities are important thermodynamic
quantities because the quark-number susceptibilities at finite temperatures
and densities can be determined from the derivatives of the quark number
density with respect to the chemical potential. Since the quark number
susceptibilities are good signatures for predicting the critical end point,
the quark number densities should also provide a good prediction for the
location of CEP. Thus we can investigate the number densities in PNJL model in
the search for QCD critical point.

In this paper, first, we discuss PNJL model for infinite and finite volume
systems followed by the calculation of the quark number densities for light
quarks and strange quarks for different beam energies of the infinite volume
and finite volume system. From the nature of the quark number density, as a
function of the center of mass energy, we will try to analyse the location of
the critical point.

##  II PNJL Model

Effective models are very important in studying high energy physics. They are
constructed with Lagrangians of simple mathematical structure having
properties similar to that of QCD. PNJL model is such an effective model which
gives good results at finite temperature and chemical potential. We shall
consider the 2+1 flavor PNJL model with six quark interactions. In the PNJL
model the gluon dynamics is described by the chiral point couplings between
quarks (present in the NJL part) and a background gauge field representing
Polyakov Loop dynamics. The Polyakov line is represented as,

$L(\bar{x})={\cal P}{\rm exp}[i{\int_{0}}^{\beta}d\tau
A_{4}{({\bar{x}},\tau)}]$ (1)

where $A_{4}=iA_{0}$ is the temporal component of Eucledian gauge field
$(\bar{A},A_{4})$, $\beta=\frac{1}{T}$, and $\cal P$ denotes path ordering.
$L(\bar{x})$ transforms as a field with charge one under global Z(3) symmetry.
The Polyakov loop is then given by $\Phi=(Tr_{c}L)/N_{c}$, and its conjugate
by, ${\bar{\Phi}}=(Tr_{c}L^{\dagger})/N_{c}$.

The Lagrangian of the PNJL model with three flavours and $U_{A}(1)$ anomaly is
given by :

$\begin{split}\cal{L}&=\sum_{f=u,d,s}\overline{\psi}_{f}(i\gamma^{\mu}D_{\mu}-m_{f}+\gamma_{0}\mu)\psi_{f}+\frac{g_{S}}{2}\sum_{a=0}^{8}[(\overline{\psi}\lambda^{a}\psi)^{2}\\\
&+(\overline{\psi}i\gamma_{5}\lambda^{a}\psi)^{2}]-g_{D}\\{det[\overline{\psi}_{f}\frac{(1+\gamma_{5})}{2}\psi_{f^{\prime}}]+\\\
&det[\overline{\psi}_{f}\frac{(1-\gamma_{5})}{2}\psi_{f^{\prime}}]\\}-\mathcal{U^{\prime}}(\Phi[A],\overline{\Phi}[A];T)\end{split}$
(2)

where, $D_{\mu}=\partial_{\mu}-iA_{\mu}$ and $A_{\mu}$ is the gauge field
which absorbs the strong interaction coupling and is given by
$A_{\mu}(x)=g_{S}A_{\alpha}^{\mu}\lambda^{\alpha}/2$. $f$ denotes the flavors
$u$, $d$ or $s$ respectively. The matrices $P_{L,R}=(1\pm\gamma_{5})/2$ are
respectively the left-handed and right-handed chiral projectors, and the other
terms have their usual meaning, described in details in Refs.
bhattacharyya2010investigation ; deb2009mesonic . The confinement/
deconfinement properties of the quarks are described by the effective
potential $\mathcal{U}^{\prime}(\Phi,\bar{\Phi};T)$ (with the Vandermonde term
$J[\Phi,\overline{\Phi}]=(27/24\pi^{2})[1-6\Phi\overline{\Phi}+4(\Phi^{3}+\overline{\Phi}^{3})-3(\Phi\overline{\Phi})^{2}]$)
expressed as:

$\displaystyle\frac{\mathcal{U}^{\prime}(\Phi[A],\overline{\Phi}[A];T)}{T^{4}}=\frac{\mathcal{U}(\Phi[A],\overline{\Phi}[A];T)}{T^{4}}-\kappa
ln[J(\Phi,\overline{\Phi})]$

where $\mathcal{U}(\phi)$ is a Landau-Ginzburg type potential given by,

$\frac{\mathcal{U}(\Phi,\overline{\Phi};T)}{T^{4}}=-\frac{b_{2}(T)}{2}\overline{\Phi}\Phi-\frac{b_{3}}{6}(\Phi^{3}+\overline{\Phi}^{3})+\frac{b_{4}}{4}(\overline{\Phi}\Phi)^{2}$
(3)

with
$b_{2}(T)=a_{0}+a_{1}\left(\frac{T_{0}}{T}\right)+a_{2}\left(\frac{T_{0}}{T}\right)^{2}+a_{3}\left(\frac{T_{0}}{T}\right)^{3}$.
For the effective potential, we can choose the following parameters :
$T_{0}=0.19$ GeV,
$a_{0}=6.75,a_{1}=-1.95,a_{2}=2.625,a_{3}=-7.44,b_{3}=0.75,b_{4}=7.5,\kappa=0.13$.

Now, we know that the energy of a quark can be written as,
$E_{i}=\sqrt{p_{i}^{2}+m_{i}^{2}}$. Using this and Eq. (2), we can obtain the
Grand potential density for the PNJL model in MFA (mean field approximation)
as,

$\begin{split}\Omega^{\prime}&=\mathcal{U}^{\prime}[\Phi,\overline{\Phi},T]+2g_{S}\sum_{f=u,d,s}\sigma_{f}^{2}-\frac{g_{D}}{2}\sigma_{u}\sigma_{d}\sigma_{s}-6\sum_{f}\int_{0}^{\Lambda}\frac{d^{3}p}{(2\pi)^{3}}E_{pf}\Theta(\Lambda-|\vec{p}|)\\\
&-2\sum_{f}T\int_{0}^{\infty}\frac{d^{3}p}{(2\pi)^{3}}ln[1+3(\Phi+\overline{\Phi}e^{-(E_{pf}-\mu)/T})e^{-(E_{pf}-\mu)/T})+e^{-3(E_{pf}-\mu)/T})]\\\
&-2\sum_{f}T\int_{0}^{\infty}\frac{d^{3}p}{(2\pi)^{3}}ln[1+3(\overline{\Phi}+{\Phi}e^{-(E_{pf}+\mu)/T})e^{-(E_{pf}+\mu)/T})+e^{-3(E_{pf}+\mu)/T})]\\\
&=\Omega-\kappa T^{4}ln[J(\Phi,\overline{\Phi})]\end{split}$ (4)

where, $\sigma_{f}^{2}=(\sigma_{u}^{2}+\sigma_{d}^{2}+\sigma_{s}^{2})$
,$\sigma_{f}^{4}=(\sigma_{u}^{4}+\sigma_{d}^{4}+\sigma_{s}^{4})$ and $\Lambda$
is the three momentum cut off. The pressure of the strongly interacting matter
can be written in terms of the potential density as,

$P(T,\mu)=-\Omega(T,\mu)$ (5)

and the number density of the quarks can be written as,

$n_{q}(T,\mu)=-\frac{\partial{\Omega(T,\mu)}}{\partial{\mu}}$ (6)

where T is the temperature and $\mu$ is the quark chemical potential. So far,
we have discussed about PNJL model for infinite volume. Now we will discuss
about its implementation in finite volume. To do this, the first step is to
choose the proper boundary condition (Periodic for bosons but anti-periodic
for fermions). This would give an infinite sum over the discrete momentum
values, given by, $p_{i}=n_{i}\pi/R$ (where i=x,y,z and $n_{i}$ are all
positive integers and R is the lateral size of a cubic volume). Then, the
effects of surface and curvatures should be incorporated. However, in this
calculations, we have taken up some simplifications, such as, (i) neglecting
surface and curvature effects, (ii) treating the infinite sum as an
integration over a continuous variation of momentum, (iii) not using any
modification to the mean field parameters due to finite size effects, as was
taken in deb1 .

##  III Results

In this present paper, we have discussed the location of critical point using
PNJL model calculation considering both infinite volume and finite volume
(with R=2 fm) system. Using PNJL model we have obtained data sets for the
number density of quark in a system as a function of the temperature of the
system at a fixed quark chemical potential ($\mu_{q}$) for different beam
energies (7.7 , 11.5, 14.5, 19.6, 27, 39, 62.4, 130, 200 GeV). The quark
number density ($n_{q}$) as a function of Temperature (T) and $n_{q}/T^{3}$ as
a function of Temperature (T) are plotted. Using the PNJL model parameters we
have also obtained similar datasets for strange and baryons number densities
which are shown in a similar manner. The list of the value of constant $\mu$’s
($\mu_{S}$ is strange quark number density and $\mu_{B}$ is the number density
of baryon) for different beam energies are shown in Table:1

Table 1: Value of $\mu_{q},\mu_{S}$, and $\mu_{B}$ for different beam energies Beam Energy $(\sqrt{s}_{NN})$ | 7.7 GeV | 11.5 GeV | 14.5 GeV | 19.6 GeV | 27.0 GeV | 39.0 GeV | 62.4 GeV | 130.0 GeV | 200.0 GeV
---|---|---|---|---|---|---|---|---|---
$\mu_{q}$ (MeV) | 140.33 | 105.33 | 88.00 | 68.67 | 52.00 | 37.33 | 24.33 | 12.00 | 8.00
$\mu_{S}$ (MeV) | 127 | 117 | 112 | 106 | 101 | 97 | 94 | 90.14 | 90.34
$\mu_{B}$ (MeV) | 421 | 316 | 264 | 206 | 156 | 112 | 73 | 36 | 24

In the graphs, we have highlighted the temperature points corresponding to the
beam energy values used in BES. To get these temperature (T) and $T/\mu_{B}$
from the value of beam energies ($\sqrt{s_{NN}}$), we have used the following
expressions, as was used by Cleymans cleymans2006comparison :

$T(\mu_{B})=a-b\mu_{B}^{2}-c\mu_{B}^{4}$ (7)
$\mu_{B}(\sqrt{s})=\frac{d}{1+e\sqrt{s}}$ (8)

where, $a=0.166\pm 0.002$GeV, $b=0.139\pm 0.016$ GeV-1 , $c=0.053\pm 0.021$
GeV-3 , and $d=1.308\pm 0.028$GeV, $e=0.273\pm 0.008$ GeV-1 . The values of
critical temperature, as obtained from calculations, is given in the following
Table: 2

Table 2: Value of $T_{c}$ for different beam energies Beam Energy $(\sqrt{s}_{NN})$ | 7.7 GeV | 14.5 GeV | 19.6 GeV | 27.0 GeV | 39.0 GeV | 62.4 GeV | 130.0 GeV | 200.0 GeV
---|---|---|---|---|---|---|---|---
Light Quark (MeV) | 228.70 | 233.50 | 234.65 | 235.50 | 235.95 | 236.20 | 236.25 | 236.25
Strange Quark (MeV) | 234.20 | 234.55 | 234.70 | 234.75 | 234.75 | 234.90 | 234.95 | 234.95
Baryon (MeV) | 153.15 | 211.75 | 219.40 | 226.50 | 231.25 | 234.15 | 235.65 | 235.90

###  III.1 Strange Quark

The strange quark number density ($n_{S}$) as a function of the Temperature
(T) is shown in Figure. 1. The coloured lines in each case, corresponding to
nine beam energy values (7.7 GeV, 11.5 GeV, 14.5 GeV, 19.6 GeV, 27.0 GeV, 39.0
GeV, 62.4 GeV, 130 GeV, 200 GeV). The $n_{S}$ vs.T graph have Poisson like
nature, with its maxima at around T $\sim$ 160 MeV in case of infinite volume
and around T $\sim$ 190 MeV for finite volume case.

Figure 1: $n_{S}$ vs. $T$ plot for different beam energies. Left: infinite
volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

Now, at the two sides of the maxima, the first order derivative or the slope
of these curves (positive or negative) has maxima. The region between these
two points is the critical region. So, we can expect that the critical end
point lies near the maxima of these graphs. The green points refer to the
values of critical temperature for each beam energy obtained from theoretical
calculations. For each line, the green points are at the base of the Poisson
curve. The red points corresponds to the beam energy points obtained using Eq.
(7) and Eq. (8). The black points refer to similar calculation, correspond to
the HRG model hrg_data . For the infinite volume case, both the red and black
points lie near the maxima of the lines. But for the finite volume case, the
points lie at far left from the maxima. PNJL and HRG model based calculations
for BES energies are close to each other and the calculated $T_{c}$ values are
close to 235 MeV for both the finite and infinite cases.

Figure 2: $n_{S}/T^{3}$ vs. $T$ for different beam energies. Left: infinite
volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

The $n_{S}/T^{3}$ vs. T graphs in Figure. 2 have the same Poisson like nature.
So, the nature of the slope of these curves, changes at the neighbourhood of
the maxima. Therefore, we can expect that the critical point is present near
the maxima of these curves too. In case of infinite volume, for lower beam
energy (7.7 GeV, 11.5 GeV), the red and black points lie near the maxima of
the lines. But, as the energy value increases, the points shift to the right
of the maxima more and more. For finite volume case, these points are at the
left of the maxima for each line (energy). The critical temperature points
(green) lie near the base of the lines, for each case.

Figure 3: $n_{S}$ vs. $T/\mu_{q}$ plot for different beam energies. Left:
infinite volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

The x-axis (temperature axis) of the $n_{S}$ vs. T graphs can be scaled by
diving it by the value of constant $\mu_{S}$ for each beam energies to get
$n_{S}$ vs. $T/\mu_{S}$ plots in Figure. 3. From the $n_{S}$ vs. $T/\mu_{S}$
plots, the peaks for different beam energies can be seen clearly. From these
peaks, we can estimate the temperature at which, the critical end point may
present. Red and black points lie near the maxima for lower beam energies in
case of infinite volume, but far away from the peaks for finite volume system
in all energies.

Figure 4: $n_{S}/T^{3}$ vs. $T/\mu_{q}$ plot for different beam energies.
Left: infinite volume, right: finite volume with R=2 fm. Points corresponds to
the temperature associated with BES energies from PNJL (red) and HRG (black)
model.

The $n_{S}/T^{3}$ vs. $T/\mu_{S}$ plots in Figure. 4 have the same nature that
of $n_{S}$ vs. $T/\mu_{S}$ plots. So, similar type of conclusions can be made
from these graphs too. The position of the three types of points are similar
to that of above plots.

###  III.2 Light Quark

The variations of $n_{q}/T^{3}$ as a function of $T$ for light quark for
finite and infinite volume are shown in Figure. 5.

Figure 5: $n_{q}/T^{3}$ vs. $T$ plot for different beam energies. Left:
infinite volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

This $n_{q}/T^{3}$ vs. T plots [5] do not have a Poisson like nature. These
plots have a rapid change in their slope, in the small T region and the slope
decreases slowly after the peak for the infinite volume case. An inflection
point is located near to left of the peak, and another one is little bit
further from the peak at right. So the critical point is assumed to be at the
peak of the curves. The red and black points are very far away (at left) from
the peak for every energy value. The $T_{C}$ points (green points) are located
very near to the peak of the curves. In case of finite volume system size, the
peaks are not present (may be due to unavailability of data). Both the
calculated $T_{c}$ values are similar for finite and infinite cases.

The $n_{q}/T^{3}$ vs. $T/\mu_{q}$ plots in Figure. 6 have the same kind of
nature as the $n_{q}/T^{3}$ vs. T graphs, as the only difference is that, the
T axis is scaled by a factor of $\mu_{q}$ . So, the position of the CEP can
also be estimated from these graphs. The calculated $T_{c}$ values (green
points) are very near to the peak of the curves. Other two type of points
(points corresponding to temperature associated with BES energies, calculated
in PNJL and HRG model) are far away from the peak.

Figure 6: $n_{q}/T^{3}$ vs. $T/\mu_{q}$ plot for different beam energies.
Left: infinite volume, right: finite volume with R=2 fm. Points corresponds to
the temperature associated with BES energies from PNJL (red) and HRG (black)
model.

The variations of $n_{q}$ as a function of $T$ for light quark for finite and
infinite volume are shown in Figure. 7. The $n_{q}$ vs. T graphs have no
maxima as such for both the finite and infinite volume. Finding CEP from this
type of plot is difficult. In the region between the black and green points,
there is an inflection point, where the slope of the curves have changed it’s
nature. Critical point lies near to that point.

Figure 7: $n_{q}$ vs. $T$ plot for different beam energies. Left: infinite
volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

Figure. 8 shows the variation of $n_{q}$ as a function of the temperature
scaled with $\mu_{q}$ for finite and infinite volume. The $n_{q}$ vs.
$T/\mu_{q}$ curves for different beam energies have a clear point of
inflection for infinite volume case. For finite volume case, this point is not
so clear. The green points lie at the right, and the black (and red) points
lie at the left of that point. The HRG and PNJL calculations are in good
agreement for all beam energies.

Figure 8: $n_{q}$ vs. $T/\mu_{q}$ plot for different beam energies. Left:
infinite volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

###  III.3 Baryon

For baryon number density the variations of $n_{B}/T^{3}$ as a function of $T$
for infinite and finite volume are shown in Figure. 9 for different beam
energies. The variations has similar dependency as in quark number density.
These curves have a rapid change in their slope, in the small T region and the
slope decreases slowly after the peak for the infinite volume case. A
inflection point is located near to left of the peak, and another one is
little bit further from the peak at right. So the critical point is assumed to
be at the peak of the curves. The red and black points are very far away (at
left) from the peak for every energy value. Calculation based on PNJL (red
point) and on HRG model (black points) are also shown on the corresponding
energy curves.

Figure 9: $n_{B}/T^{3}$ vs. $T$ plot for different beam energies. Left:
infinite volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

Figure. 10 shows $n_{B}/T^{3}$ as a function of $T/\mu_{B}$ for different RHIC
energies for infinite and finite volume. The position of the CEP can also be
estimated from these graphs. The points corresponding to temperature
associated with BES energies, calculated in PNJL and HRG model are also shown
and are far away from the peak.

Figure 10: $n_{B}/T^{3}$ vs. $T/\mu_{B}$ plot for different beam energies.
Left: infinite volume, right: finite volume with R=2 fm. Points corresponds to
the temperature associated with BES energies from PNJL and HRG model.

The variations of $n_{B}$ as a function of $T$ for different energies for
infinite and finite volume are shown in Figure. 11. Here also, the variations
has similar dependency as in quark number density. These curves have no maxima
as such for both the finite and infinite volume. Finding CEP from this type of
plot is difficult.

Figure 11: $n_{B}$ vs. $T$ plot for different beam energies. Left: infinite
volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

The variations of $n_{B}$ as a function of $T/\mu_{B}$ for different energies
for infinite and finite volume are shown in Figure. 12. The $n_{B}$ vs.
$T/\mu_{B}$ curves for different beam energies have a clear point of
inflection for infinite volume case. For finite volume case, this point is not
so clear. The black (and red) points lie at the left of that point. The HRG
and PNJL calculations are in good agreement for all beam energies. All of
these four kinds of plots for the baryons show the exactly same nature as the
quarks. So, the location of the critical points in the QCD phase diagram can
be estimated from these plots similarly.

Figure 12: $n_{B}$ vs. $T/\mu_{B}$ plot for different beam energies. Left:
infinite volume, right: finite volume with R=2 fm. Points corresponds to the
temperature associated with BES energies from PNJL (red) and HRG (black)
model.

##  IV Summary

In this work, we have tried to find out the critical end point of the QCD
phase diagram in PNJL model. This is a very important work at hand, for high
energy physicists, all over the world. This is so important because, by
finding the CEP of the QCD phase diagram, we can get an idea about the
critical point of the early universe, i.e, about the temperature and state of
matter just after the Big Bang. To find the CEP, many different approaches are
taken by different scientists at different times. With the advancement of
technology in the world, we are not just stuck with the theoretical
approaches. We have found paths to reach more close to the CEP in the
experimental world. One of these aspects in the experimental world is the
heavy ion collider. Inside these colliders, we can collide ions, protons, or
so on, to get the elementary constituents of a particle. But as these
collisions are instantaneous and the elementary constituents produced, are
evanescent, it is very difficult to observe them directly. So, we need
different models to estimate the values of different observables. The PNJL
model is such an effective model. In our present study, we have discussed the
formalism of the PNJL model briefly in section 2. In the following section, we
have demonstrated our work to find the location of the critical end point. We
have plotted the number density and number density/$T^{3}$ as a function of
temperature and $T/\mu_{x}$ , where x stands for light quarks, strange quarks,
and baryons. We have analyzed the nature of the curves and tried to find out
the point of inflection (The point of inflection is the point on a curve,
where a curve changes its nature.) and the maxima of the curves. The
inflection point of the graph gives an estimation for the critical region and
the maxima give an estimation for the CEP of the QCD phase diagram.

To conclude, we have studied the finite volume number density fluctuations at
RHIC energies with system size both as finite and infinite in PNJL model. The
variation of $n_{q}/T^{3}$ with T has the same nature in both these cases. The
peak value for $n_{q}/T^{3}$ is around 0.44 for infinite volume and is around
0.18 in case of finite volume for $\mu$ =105.33 MeV (for $\sqrt{s}_{NN}$ =11.5
GeV) which is lesser than 0.44. For other beam energies this difference
between peak values of $n_{q}/T^{3}$ can also be seen. The difference in peaks
can also be seen in other kinds of plots too. So, we can conclude that the
number density of the quarks and baryons are dependent on the size of the
system. The nature of our $n_{S}$ vs. T (or $T/\mu_{q}$ ) and $n_{S}/T^{3}$
vs. T (or $T/\mu_{q}$ ) plots for strange quarks are like Poisson
distributions. So, it has points of inflections at the two sides of the
maxima. We have highlighted the temperature points corresponding to the BES
energies on the curves. Some of these points, lie at the maxima, and some
points lie at the two sides. These points lie in between the inflexion point
and the maxima. Also, these points do not lie very near to the $T_{c}$ values.
Therefore, the critical end point may not occur at the measured RHIC beam
energy values exactly. So, at the energies, at which Beam Energy Scans are
taken, we may not observe all the behaviors which occur at the critical point.
Although, these measured points are located at the neighbourhood of the
desired points, the measurement will have influence of critical region and the
region near the critical end point can be explored by RHIC BES.

###### Acknowledgements.

P.Deb would like to thank Women Scientist Scheme A (WOS-A) of Department of
Science and Technology (DST) funding with grant no SR/WOS-A/PM-10/2019 (GEN).
The authors would also like to thank Prof. Raghava Varma (IITB) and Ratnamay
Koley (RKMRC) for the valuable discussion.

## References

  * (1) J. Adams et.al., Nucl. Phys. A 751 102–183 (2005).
  * (2) M. M. Aggarwal et. al., Phys. Rev. C 82 024905 (2010).
  * (3) D. Adamova et. al., Phys. Rev. Lett. 90, 022301 (2003).
  * (4) H.-T. Elze and W. Greiner, Phys. Lett. B 179, 385 (1986).
  * (5) M. Luscher, Commun. Math. Phys. 104, 177 (1986).
  * (6) J. Gasser and H. Leutwyler, Phys. Lett. B 188, 477 (1987).
  * (7) Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961); 124, 246 (1961).
  * (8) O. Kiriyama and A. Hosaka, Phys. Rev. D 67, 085010 (2003).
  * (9) J. Braun, B. Klein and P. Piasecki, Eur. Phys. Jr. C 71, 1576 (2011).
  * (10) J. Braun, B. Klein and B.-J. Schefer, Phys. Lett. B 713, 216 (2012).
  * (11) K. Fukushima, Phys. Lett. B 591 277 (2004).
  * (12) C. Ratti, M.A. Thaler, and W. Weise, Phys. Rev. D 73 014019 (2006).
  * (13) R.D. Pisarski, Phys. Rev. D 62 111501 (2000); A. Dumitru and R.D. Pisarski, Phys. Lett. B 504 282 (2001); 525 95 (2002); Phys. Rev. D 66 096003 (2002).
  * (14) K. Fukushima, Phys. Rev. D 77 114028 (2008).
  * (15) H. Hansen, W.M. Alberico, A. Beraudo, A. Molinari, M. Nardi, and C. Ratti, Phys. Rev. D 75 065004 (2007).
  * (16) M. Ciminale, R. Gatto, N.D. Ippolito, G. Nardulli, and M. Ruggieri, Phys. Rev. D 77 054023 (2008).
  * (17) S.K. Ghosh, T.K. Mukherjee, M.G. Mustafa, and R. Ray, Phys. Rev. D 73 114007 (2006); S.K. Ghosh, T.K. Mukherjee, M.G. Mustafa, and R. Ray, Phys. Rev. D 77 094024 (2008).
  * (18) A. Bhattacharyya, P. Deb, S. K. Ghosh, R. Ray, S. Sur, Phys.Rev. D 87 054009 (2013).
  * (19) A. Bhattacharyya, R. Ray and S. Sur, Phys. rev. D 91 051501 (2005).
  * (20) A. Bhattacharyya, P. Deb, S. K.Ghosh, R. Ray, Phys. Rev. D 82 014021 (2010).
  * (21) P. Deb, A. Bhattacharyya, S. Datta, and S.K. Ghosh,Phys. Rev. C 79 055208 (2009); A. Bhattacharyya, P. deb, A. Lahiri, R. Ray, Phys. Rev. D D 82 114028 (2010); A. Bhattacharyya, P. Deb, A. Lahiri, R. Ray, Phys. Rev. D 83 014011 (2011).
  * (22) A. Bazavov and B. A. Berg, Phys. Rev. D 76, 014502 (2007).
  * (23) G. Boyd et. al., Nucl. Phys. B 469 419 (1996).
  * (24) J. Engels et. al., Nucl. Phys. B 558 307 (1999).
  * (25) Z. Fodor and S.D. Katz, Phys. Lett. B 534 87 (2002); Z. Fodor, S.D. Katz, and K.K. Szabo, Phys. Lett. B 568 73 (2003).
  * (26) C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, Ch. Schmidt, and L. Scorzato, Phys. Rev. D 66 074507 (2002); C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F.Karsch, E. Laermann, and Ch. Schmidt, Phys. Rev. D 68 014507 (2003); C.R. Allton, M. Doring, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, and K. Redlich, Phys. Rev. D 71 054508 (2005).
  * (27) P. de Forcrand and O. Philipsen, Nucl. Phys. B642 290 (2002); B673 170 (2003).
  * (28) Y. Aoki, Z. Fodor, S.D. katz, and K.K. Szabo, Phys. Lett.B 643 46 (2006); Y. Aoki, G. Endrodi, Z. Fodor, S.D. Katz, and K.K.Szabo, Nature (London) 443 675 (2006).
  * (29) E. Megias, E. Ruiz Arriola, and L.L. Salcedo, Rom. Rep.Phys.58, 081 (2006); E. Megias, E.R. Arriola, and L.L. Salcedo, Proc. Sci., JHW2005 (2006) 025; E. Megias, E.R. Arriola, and L.L. Salcedo, Nucl. Phys. B, Proc. Suppl.186 256 (2009); E. Megias, E.R. Arriola, and L.L. Salcedo, Phys. Rev. D 81, 096009 (2010).
  * (30) D. Tlusty, arXiv:1810.04767.
  * (31) G. Odyniec and STAR collaboration, PoS CORFU2018 151 (2019).
  * (32) H-S. Zong and W-M. Sun, Phys. Rev. D 78 054001 (2008).
  * (33) J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, Phys. Rev.C 73 034905 (2006).
  * (34) F. Karsch and K. Redlich ,“Moments of charge fluctuations, pseudo-critical temperatures and freeze-out in heavy ion collisions ” Phys. Lett. B 695 (2011)


# On performance of thin-film meso-structured perovskite solar cell through
experimental analysis and device simulation

F. Bonnín-Ripoll1 Ya. B. Martynov2 R. G. Nazmitdinov 3,4 K. Tabah5 C. Pereyra5
M. Lira-Cantú5 G. Cardona6 R. Pujol-Nadal1 1Departament d’Enginyeria
Industrial i Construcció, Universitat de les Illes Balears, E-07122 Palma,
Spain.  2State Scientific-Production Enterprise “Istok” Fryazino, Russia.
3Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear
Research, 141980 Dubna, Russia.  4Dubna State University, 141982 Dubna, Moscow
region, Russia.  5Catalan Institute of Nanoscience and Nanotechnology (ICN2),
CSIC and the Barcelona Institute of Science and Technology (BIST), Building
ICN2, Campus UAB, E-08193 Bellaterra, Barcelona, Spain.  6Departament de
Matemàtiques, Universitat de les Illes Balears, E-07122 Palma, Spain.

###### Abstract

In the last few years there is an unprecedented progress in the increase of
the power conversion efficiency of perovskite solar cells. Evidently, further
advances of the efficiency of these devices will depend on the constraints
imposed by the optical and electronic properties of their constituents. Quite
apparently that during the manufacturing process of a solar cell, there is an
inevitable variation in the thicknesses of various functional layers, which
affects the optoelectronic characteristics of the final sample. In this work a
possible strategy of the analysis of the solar cell performance is suggested,
based on statistically averaging procedure of experimental data. We present a
case study, in which the optoelectronic properties of the meso-structured
perovskite solar cell (with a mesoporous TiO2 layer) are analysed within the
method providing a deeper understanding of the device operation. This method
enables an assessment of the overall quality of the device, pointing pathways
towards the maximum efficiency design of a perovskite solar cell by material
properties tuning.

Keywords: mesoporous, optoelectronic measurements and analysis, ray tracing,
transfer matrix method, transport equations, power conversion efficiency

##  I INTRODUCTION

Nowadays organic-inorganic halide perovskites, as solar cell light-absorbers,
attract substantial attention. Solar cells with such photoactive layers are
considered as one of the most promising competitors to silicon-based
photovoltaics, reaching efficiencies of $25.8\%$ only in two decades of
development [1]. Low-cost production of these materials with high power
conversion efficiency (PCE) pledge fast development of flexible solar cells
[2]. Moreover, there is challenging future in improvement of their performance
in order to reach their maximum efficiency (cf.[3]) by means: engineering of
perovskite composition [4] or/and perovskite solar cell (PSC) constituents [5,
6, 7, 8, 9]; study and improvement of perovskite stability [10]; altering
characteristics of optoelectronic behavior of PSC constituents [11] by
incorporation, for example, light trapping mechanisms [12]; to name just a
few.

Evidently, for further efficiency improvements of PSCs, the analysis of
optical and electric properties of different layers in the device architecture
is pivotal for optimizing light-harvesting in the absorber, whilst allowing
coherent effects and parasitic absorption to be accounted in the stack design.
In addition, this analysis should be based on statistically reliable estimates
of noticeable dispersion of various sample parameters during the manufacturing
process. Partially, some attempts in this direction have been done in [13, 14,
15], where various type of simulations of optoelectronic properties of PSCs
have been conducted recently. It is worth noting that for the analysis of
optical properties various research groups (see, for example, [16, 17, 18, 19,
20]) use the software SCAPS-1D [21] as a main tool. This code does not,
however, consider interference, scattering, or intermediate reflectors (see
details in [21]), thus preventing the comprehensive analysis of the
antireflection layer’s interference and the influence of thin films
thicknesses.

Also, a broad range of electron (hole) transport materials and back metal
contacts have been considered to examine the performance of lead-free PSCs in
[22]. It also became common to introduce a mesoporous TiO2 layer between TiO2
electron collecting contact and perovskite (the main absorber layer). This
process decreases defect density at the absorber boundary, while may introduce
additional resistance. Nevertheless, the effect of a mesoporous TiO2 layer and
its influence on the PSC performance remains in infancy. Here, we employ a ray
tracing-based optoelectronic computational model [23] to reproduce the
behaviour of a perovskite solar cell with a mesoporous TiO2 layer and
characterize experimental PSC samples. To this aim we further developed this
model, incorporating a few novel elements to treat properly thin-film meso-
structured perovskite solar cells.

Figure 1: (a) Structure of the layers present in each manufacturing stage. A
total of five stages in which the optical behavior has been studied. (b)
Typical cross-section image of the sample, obtained by means of the scanning
electron microscope.

The experimental device configuration (see Fig.1) is selected as the well
established experimentally perovskite solar cell [24]. The measurements of
optical transmittance have been performed in each stage of the manufacturing
process. In this n-i-p junction device structure, apart the complex glass
substrate, a thin and compact TiO2 (c-TiO2) is manufactured as an electron
selective contact to facilitate the collection of photogenerated electrons
from perovskite absorber. Next, to improve the efficiency of carrier
collection a mesoporous TiO2 (mp-TiO2) is introduced to provide a suitable
perovskite deposition [25]. In addition, the partial infiltration of halide
perovskite into the mp-TiO2 layer facilitates higher quantum yield for photo-
excited charge separation. Spiro-OMeTAD has been utilized in establishing a
good p-contact between the gold (Au) counter electrode with a perovskite
absorber in our stack design. As mentioned above, we will elaborate the model
[23] and discuss in detail the effective medium approach for analysis of the
optical behaviour of mp-TiO2. Our approach is experimentally verified with
measurements of the device reflection and transmission spectra. These
characteristics have been measured by means of spectrophotometry techniques
applied at each manufacturing stage. Finally, we make use of our approach to
evaluate the dependence of the PSC performance on the properties of the
mesoporous structure of the considered solar cell architecture, based on
statistically averaged physical and geometrical parameters of the device.

##  II Methodology

###  II.1 Device Characteristics

Since electronic and optical measurements cannot be performed on the same PSC
sample, we use two batches of 24 solar cells, manufactured at the same
conditions. One batch is used for optical measurements with interruptions on
each phase of the process, while the second one is used for J-V characteristic
measurements without any interruptions on each phase of the process.

Considering the deposition of a material as a new phase of the process, we
have performed the measurements of optical transmittance of each PSC layer in
each stage of the manufacturing process (see details in Appendix A). Stage 1
of the manufacturing process consists of the glass substrate, assuming the
material as a multilayer system itself, formed by sodalime glass, SnO2, SiO2
and fluorine doped tin oxide (FTO). In the stage 2 a layer of c-TiO2 of the
order of tens of nanometers is added by the spray pyrolysis technique. Next, a
mp-TiO2 layer is deposited in stage 3, followed by the subsequent deposition
of a perovskite (stage 4) and Spiro-OMeTAD (stage 5), the latter three using
the spin coating technique.

Next, using a solar simulator, the J-V characteristics of the PSCs in
operation a second batch of $24$ solar cells has been measured. The result of
this measurements provides the average value of the short-circuit current
$J_{sc}=20.07$ mA/cm2 with a standard deviation $0.99$ mA/cm2 at irradiance
$1000$ W/m2. In this study we will focus on reproducing the behaviour of the
reverse scan of the J-V characteristics. The scanning process conditions are
discussed in [24].

###  II.2 Optoelectronic Model

To analyse the amount of transmitted, absorbed and reflected light in each of
the materials present in the PSCs, we employ ray tracing simulations with the
aid of the ray trace OTSun python package [26]. It is a Monte Carlo ray
tracing where the optics is implemented through the Fresnel optics equations
in their most general form, complemented by the transfer-matrix method (TMM)
to account for the interference phenomenon [27, 28]. This software allows us
to perform the simulations by loading a geometry file designed with the aid of
FreeCAD [29], an open-source software for the parametric modeling of geometric
figures in 3D.

Input parameters in our calculations (References)
---
|  c-TiO2 |  mp-TiO2 |  Perovskite |  Spiro-OMeTAD
$\chi$ [eV] |  4 [30] |  4 [30] |  3.75 [31] |  2.12 [32]
Eg [eV] |  3.05 [33] |  3.05 [33] |  1.66 [34] |  3.1 [35]
$m_{e}^{\star}$ |  1 |  1 |  1 |  1
$m_{h}^{\star}$ |  1 |  1 |  1 |  1
$N_{d/a}$[cm-3] |  $1\times 10^{18}$ |  variable |  0 |  $1\times 10^{21}$
$\mu$ [cm2/Vs] |  0.006 [36] |  0.006 [36] |  variable |  0.0001 [37]
$\varepsilon$ |  60 [38] |  42.45** |  60 [39] |  3 [40]
$\ell_{\it D}$[nm] |  4 |  4 |  variable |  0.5
Table 1: Characteristics of the considered semiconductors: TiO2, FAPbI3
(Perovskite), Spiro-OMeTAD. We use the following notations: $\chi$ is an
electron affinity; Eg is an energy gap; $m_{e/h}^{\star}$ is an effective
electron/hole mass; $N_{d/a}$, is a donor/acceptor concentration; $\mu$ is a
mobility; $\varepsilon$ is a permittivity; $\ell_{\it D}$ is a diffusion
length.(**) The definition of the permittivity of the mp-TiO2 is discussed in
Sec.II; see Eq.10 and the accompanying discussion.

The simulations require knowledge of the layer thicknesses as well as the
complex refractive indexes of the corresponding materials as a function of
wavelength ($\lambda$). In particular, the complex refractive indexes of the
following materials have been taken from literature: FTO [41], TiO2 [42],
Spiro-OMeTAD [43]. In the case of SnO2 and SiO2 the data are taken from the
material manufacturer itself. According to the experimental setup [24] the
chemical structure of the perovskite layer has been associated with FaPbI3. In
this case the complex refractive index data are taken from [34], where the
real and imaginary parts of the refractive index have been measured for a
similar perovskite. The physical characteristics of all layers that form the
architecture of our device are displayed in Table 1.

On the other hand, if the layer contains a few constituents there is a need in
some effective approach. The refractive index of the glass has been calculated
from transmittance experimental measurements, while for the mp-TiO2 material
we have developed the effective approach to determine its complex refractive
index (see below).

### Glass complex refractive index

Since the glass substrate consists of a few thin components, to determine the
glass optical properties itself the oxides films have to be removed. To this
aim, on the glass multilayer system the etching process is done by using Zn
powder and HCl solution. As a result, the optical properties of the glass
substrate has been determined with the aid of the method [44], based on the
ideas discussed below. If the slab of material has plane-parallel faces, and a
material is absorbing partially ($k^{2}\ll n$), the real part $n$ and the
imaginary part $k$ of the refractive index of the slab can be calculated from
the intensity of reflectance $R(\lambda)\equiv R$ and transmittance
$T(\lambda)\equiv T$ spectra (for the sake of convenience we omit the symbol
$\lambda$). In particular, the following analytical formulas have been used:

$k=\frac{\lambda}{4\pi d}\ln\Bigr{[}\frac{R_{F}\cdot T}{R-R_{F}}\Bigr{]}\,,$
(1)

$n=\frac{1+R_{F}}{1-R_{F}}+\Bigr{[}\frac{4R_{F}}{(1-R_{F})^{2}}-k^{2}\Bigr{]}^{1/2}\,,$
(2)

where the slab geometrical thickness, $d$, is assumed as known, and

$R_{F}=(2+T^{2}-(1-R)^{2}-\\{[2+T^{2}-(1-R)^{2}]^{2}-4R(2-R)\\}^{1/2})(2(2-R))^{-1}\,.$
(3)

Experimental results for the transmittance $T(\lambda)$ and reflectance
$R(\lambda)$ spectra, and the complex refractive index, obtained by means of
Eqs.(1)–(3), are shown in Fig.2.

Figure 2: Left panel: the results of transmittance and reflectance
measurements of the glass substrate as a function of the wavelength are shown.
Right panel: the real ($n$) and imaginary ($k$) components of the calculated
complex refractive index of the glass substrate as a function of the
wavelength are displayed.

### Mesoporous layer as an effective medium: an optical approach

Mesoporous TiO2 is one of the commonly used electron transport materials
(ETMs) in PSCs. The mp-TiO2 film serves not only as a scaffold for the
perovskite layer but also as a route for electron transport. However, in many
simulations (see, for example, [45, 46, 47, 48]) the effect of the porosity of
the mesoporous TiO2 layer on optical and electronic properties meso-structured
perovskite solar cells is ignored.

The presence of mesoporous thin films requires a proper approach to determine
its optical characterization. Several effective medium approximations have
been suggested to describe the optical properties of heterogeneous thin films
by considering them as a homogeneous media with an effective refractive index,
along with effective absorption coefficient. Most of the models have been
developed for the effective refractive index omitting the absorption
coefficient (see for a review [49]). However, the absorption can affect the
optoelectronic properties of ETM, and, consequently, the properties of the
whole system.

To consider the complex refractive index we apply the Volume Average Theory
(VAT) model, discussed among others in [49]. It provides the effective complex
refractive index (subindex eff) of a material (see Eqs.(4)–(7)) composed of a
continuous phase material matrix (subindex c) and a dispersed phase material
matrix (subindex d), with a porosity $p$:

$n_{eff}^{2}=\frac{1}{2}(A+\sqrt{A^{2}+B^{2}})$ (4)
$k_{eff}^{2}=\frac{1}{2}(-A+\sqrt{A^{2}+B^{2}})$ (5)
$A=p(n_{d}^{2}-k_{d}^{2})+(1-p)(n_{c}^{2}-k_{c}^{2})$ (6)
$B=2pn_{d}k_{d}+2(1-p)n_{c}k_{c}$ (7)

The properties of the mp-TiO2 material can be characterized using an effective
medium that contains a mixture of the c-TiO2 material and air, as suggested in
[50, 51, 52]. In the latter paper, however, the effect of the extinction
coefficient $k$ is missing. Evidently, it may play important role, since the
absorption in the mp-TiO2 layer depends on the degree of its porosity. To find
thickness configurations that match best the experimental results we developed
the special procedure to quantify the minimal deviation of the theoretical
estimations from the experimental values for the transmittance and the short-
circuit current (see Figs.9, 10 and the accompanying text). Once the above
analysis has been completed, to gain much more accurate results on the optical
properties of the PSC we use the OTSun python package.

### Impact of porosity and Lambertian scattering

The knowledge of the thickness configuration enables to us to simulate the
photogeneration rate for its application into the electronic transport
equations. For the sake of convenience we recapitulate briefly the basic steps
of this procedure (see details in [23]).

According to the Beer-Lambert law for a given absorption coefficient
$\alpha(\lambda)=4\pi k(\lambda)/\lambda$ of a material, the loss of light
intensity is a function of the path length $\ell$ of the light beam through a
material:

$N(\ell,\lambda)=N_{0}(0,\lambda)e^{-\alpha(\lambda)\mathsf{\ell}}\,.$ (8)

Here, $N_{0}(0,\lambda)$ is the initial number of photons of the incident ray
on a material per unit area, per wavelength, and per unit time (for each time
when a ray impacts onto a material according to the path trajectories
determined by the ray tracing simulation). With the aid of OTSun and Eq.(8),
we calculate the generation rate produced by a single ray

$G(z)=\int_{0}^{\infty}\frac{\alpha(\lambda)N_{0}(0,\lambda)}{\cos\theta}e^{-\alpha(\lambda)z/\cos\theta}d{\lambda}\,,$
(9)

where $z=\ell\cos\theta$ is the depth from a semiconductor surface, and
$\theta$ is the ray refracted angle. The function $G(z)$ provides the number
of electrons per unit volume and per unit time, generated at each point in the
device due to the photon absorption. The total function $G(z)$ is the sum of
all rays (transmitted or reflected) passed at the point $z$. As a result,
obtained PSC thickness configuration and the rate function $G(z)$ allow to
calculate the J-V characteristics by means of the transport model [23]. In
this model there is a set of parameters that are relatively well established
either from literature or from measurements. Besides, there is another set
that consists of some uncertain parameters. All of them are listed in Table 1,
where uncertain parameters are named as “variable”. The latter ones are
treated as fitting parameters to explain the measured data.

There are a few comments in order, however. At this stage we have to elaborate
the model [23], taking into account the presence of the mp-TiO2 layer.
Consequently, to carry out the calculations within the transport model we have
to define the permittivity and the conductivity of the mp-TiO2 layer
$\sigma_{mp}=q\mu_{mp}N_{mp}$; here $\mu_{mp}$ is a majority carriers
mobility, and $N_{mp}$ is the doping concentration.

From the optical properties of different layers, we obtain the fact that the
mp-TiO2 layer includes air-filled holes (porosity). As the concentration of
these holes is small ($\sim 20\%$), we apply the electrodynamic approach
developed for weakly mixed systems (see the textbook [53]). In this case the
permittivity of the mp-TiO2 material is provided by the expression for
permittivity of the mixture $\varepsilon_{mix}$:

$\varepsilon_{mix}=\varepsilon_{TiO_{2}}+c\frac{3(\varepsilon_{a}-\varepsilon_{TiO_{2}})\varepsilon_{TiO_{2}}}{\varepsilon_{a}+2\varepsilon_{TiO_{2}}}\,,$
(10)

where $\varepsilon_{TiO_{2}}$ and $\varepsilon_{a}$ are the permittivities of
TiO2 and air, respectively; $c$ is the concentration of holes. As for the
conductivity, we make use of the fact obtained from the optical properties
fitting. Due to small concentration of air-filled holes in the mp-TiO2 layer
one can suppose that the conductivity in this layer should be reduced on 20%
as well. As a result, it seems reasonable to assume that for the donor
concentration holds the relation $N_{mp}=0.8\cdot N_{c}$, while the mobility
of the carriers remains the same $\mu_{mp}=\mu_{c}$, as in the c-TiO2 layer
(see Sec.III) To obtain the better agreement with the J-V characteristics of
our experimental device, we propose to vary only the unknown parameters
$N_{mp}$, the mobility of carriers $\mu$ and their relaxation time $\tau$ in
the perovskite layer, keeping the other parameters fixed.

It is important to stress that OTSun provides a more accurate form of the
photogeneration rate $G(z)$ compare to the TMM approach used in Appendix A,
since non-specular reflection and multiple reflection are considered. In
particular, from images obtained with the aid of the scanning electron
microscope (see Fig.8) it can be seen that the gold layer (the top layer) has
nano-sized roughness, of the order of magnitude of its thickness ($\sim 80$
nm). With a high degree of certainty, it can be argued that these textures
provide a Lambertian reflection to the surface. Such textures randomly reflect
the light that impacts on the gold layer, and, consequently, increase the
optical path in the perovskite, and the photocurrent (cf.[11]). As a result,
taking into account all facts above, we obtain by means of Eq.(9) that
$J_{sc}=20.15$ mA/cm2 (see Fig.10), i.e., it is slightly higher than the TMM
result and closer to the experimental value.

Figure 3: Experimental (solid line) and theoretical (dotted line)
transmittance curves for the different manufacturing stages of the PSCs under
study. For the stage identifications see also Fig.1

##  III RESULTS AND DISCUSSION

Begin with, we recall the main result on the analysis of the layer thicknesses
(see details in Appendix A). The layer configuration chosen as the theoretical
model of the PSCs studied is (in nm): $11$ (SiO2), $16$ (SnO2), $565$ (FTO),
$24$ (c-TiO2), $240$ (mp-TiO2, porosity = $20\%$), $500$ (Perovskite), $250$
(Spiro-OMeTAD). The result of this final disposition that produces the best
consistency between theoretical and experimental values for the short-circuit
current and the optical transmittance in each stage of the manufacturing
process is displayed in Fig.3.

There is a rather well agreement between experimental and theoretical results,
except two cases: stage 3 and stage 4. The visible difference between the
experimental and theoretical curves is located in the wavelength region
$380-500$ nm. This result implies some inaccuracy in estimation of the mp-TiO2
layer absorption coefficient. Theoretical result exhibits the larger
absorption for the perovskite FAPbI3 (stage 4) in comparison with that of the
experimental sample in the wavelength range $770-800$ nm. This disagreement is
due to the difference between the experimental value of the perovskite
absorption coefficient and that adopted from [34]. Since this difference is
located at a small interval $760<\lambda<800$ nm of the right boundary of the
considered wavelength region, it may be safely ignored in the description of
the J-V characteristics discussed below.

The values of theoretical layer thicknesses and material parameters (see Table
1) are used as the input parameters in the electronic transport model to fit
the simulated and measured J-V characteristics of the PSC. The fitting
procedure is described in [23] but here we mention some peculiarities:

  1. 1

The main source of difference between the simulated and the measured short-
circuit current has been related in [23] to the reflection properties of the
PSC layers. In present paper the reflection properties are already taken into
account by our optical model. It results in the remarkable accord between
theoretical and experimental values of the short-circuit current (see below).

Figure 4: RMSE as a function of carrier mobility in perovskite for different
doping concentration in the mp-TiO2 layer. The solid curves are associated
with the results of calculations without any additional resistances. The
dotted line is associated with the results obtained by means of the additional
resistances (see details in the text).

  2. 2

To fit the simulated and “measured averaged” open circuit voltage ($V_{oc}$)
one needs to vary the carrier lifetimes ($\tau$) in the perovskite [23]. Once
we vary the donor concentration $N_{mp}$ and the mobility $\mu$, it is
required to fit a $\tau$-value to reproduce the experimental value of
$V_{oc}$. To control the accuracy of the fitting procedure of J-V
characteristics, we introduce the root mean square error

$RMSE=\frac{1}{N}\sqrt{\sum_{i}[J_{sc}^{exp}(V_{i})-J_{sc}^{th}(V_{i})]^{2}}\,,$
(11)

where the number of steps $N$ is defined by the voltage range interval
$V=0-V_{oc}$ passed with a step size $\Delta V$. At each new value of the
voltage $V_{i}=i\cdot\Delta V$, we solve the system of nonlinear transport
equations (see Sec.2.4 in [23]) to calculate the current and compare its value
with the experimental data. The results of calculations manifest the
dependence of the RMSE on values of $N_{mp}$ and the mobility $\mu$ (see
Fig.4).

  3. 3 Figure 5: J-V characteristics. In panels (a) and (b) solid line connecting solid circles is associated with the “measured averaged” values. The panel (a): dashed (dot-dashed) line connecting solid triangles (squares) is associated with the results obtained without (with) additional resistances at $20\%$ porosity of mp-TiO2 layer. The panel (b): dashed line connecting solid triangles is associated with the results obtained without additional resistances (see for details the text). Light blue band around all these curves marks the range in which the measured J-V characteristics of different samples lie. The insert in the panel (a) displays the circuit layout of the generator.

We recall that the analysis of the transmittance resulted in a degree of
porosity of the order of 20% for the mp-TiO2 layer. Evidently, the presence of
the porosity in the mp-TiO2 layer should affect the donor concentration in
this material in comparison to that in the c-TiO2 layer, which is $1\cdot
10^{18}$ cm-3 (see Table 1). Consequently, it is natural to suppose that the
donor concentration in the mp-TiO2 layer will decrease to the value $8\cdot
10^{17}$ cm-3, proportionally to the degree of its porosity. Taking into
account that the porosity affects the permittivity (see Eq.10) and the donor
concentration of the mp-TiO2 layer, we solve the transport equations and reach
the RMSE minimum (0.18 mA/cm2) at $\mu=5$ cm2/(V$\cdot$s) (see Fig.4;
triangles connected by solid yellow line). There is a visible difference
between the simulated and the “measured averaged” J-V characteristics (see
Fig.5(a)).

A fairly common trick in the calculation of J-V characteristics is to use the
analogy of the PSC to a current generator in parallel with a diode (see, for
example, the analysis in [11]). Hereafter we follow this procedure, since such
an approach enables to us to consider the parasitic resistance of the contacts
and feeding conductors, R; the leakage resistance, related to the local short-
circuit of electron/hole transport material layers, Rsh (see the insert in
Fig.5(a)).

Once all necessary additional resistances are generated, we obtain a good
agreement of our results with the experimental J-V characteristics. This
agreement (see Fig.5(a)) is reached with the aid of the following values of
the variable parameters: $\mu$=15 cm2/(V$\cdot$s) (see Fig.4); R=4
Ohm$\cdot$cm2; R${}_{sh}\rightarrow\infty$; $\tau=2.47$ ns that yields the
diffusion length $\ell_{D}=\sqrt{\mu k_{B}T\tau/|e|}=0.3$ $\mu$m.

It turns out, however, that we can obtain even better agreement varying the
donor concentration in the mp-TiO2 layer. The variation of the RMSE between
simulated and “measured averaged” J-V characteristics yields the RMSE absolute
minimum (see Fig.4) at the “optimal” mobility $\mu$=50 cm2/(V$\cdot$s), and
the donor concentration $N_{mp}=1\cdot 10^{17}$ cm-3. The optimal carrier
lifetimes $\tau$=0.7 ns yields again the diffusion length $\ell_{D}=0.3$
$\mu$m. At these parameters we obtain a remarkable agreement between the
simulated and the “measured averaged” J-V characteristics, without any
additional resistances (see Fig.5(b)). It seems that in the mp-TiO2 layer pore
holes generate all necessary additional resistance in the measured device.
Consequently, it appears that in the real device the pore dimension occupying
20% of the mp-TiO2 layer is sufficient to affect optical properties of the
PSC. It turns out, however, that only small part of the remaining TiO2 take
part in the conductance due to large amount of defects on the pore boundaries.

It is noteworthy, the parameters values obtained by means of our fitting
procedure are located within the parameter range considered in literature. In
particular, several studies suggest that the doping values of the mp-TiO2
material are in the range of $10^{16}$ to $10^{19}$ (cm-3) [54]. Moreover, the
carrier mobility range measured for perovskite is in the range from 5 to 50
(cm2/(V$\cdot$s))[55]. The degree of the accuracy of our analysis is
summarized in Tables 2,3. It demonstrates that the short-circuit current, the
open circuit voltage, the fill factor and the power conversion efficiency for
both simulated and experimental PSC are in a remarkable agreement.

|  $J_{sc}$ (mA/cm2) |  $V_{oc}$ (V) | FF | PCE(%)
---|---|---|---|---
Experiment | $20.072$ | $0.985$ | $0.695$ | $13.735$
Model (a) | $20.120$ | $0.981$ | $0.772$ | $15.233$
Model (a) + R | $20.101$ | $0.981$ | $0.705$ | $13.903$
Model (b) | $20.076$ | $0.984$ | $0.703$ | $13.882$
Table 2: Main PSC parameters for the “measured averaged” device and simulated results obtained by means of our approach. Model (a) is associated with results displayed on the the panel Fig.5(a), without any resistance; the corresponding parameters are presented in the first column of Table 3. Model (a) + R is associated with results displayed on the the panel Fig.5(a), with resistance; the corresponding parameters are presented in the second column of Table 3. Model (b) is associated with results displayed on the the panel Fig.5(b); the corresponding parameters are presented in the third column of Table 3. Parameters | Model (a) | Model (a) + R | Model (b)
---|---|---|---
$N_{mp}$ (cm-3) | $8\cdot 10^{17}$ | $8\cdot 10^{17}$ | $1\cdot 10^{17}$
$\mu$ (cm2/(V$\cdot$s)) | $15$ | $15$ | $50$
$\tau$ (ns) | $2.47$ | $2.47$ | $0.7$
R (Ohm$\cdot$cm2) | $-$ | $4$ | $-$
Rsh (Ohm$\cdot$cm2) | $-$ | $\infty$ | $-$
Table 3: Summary of the parameters used to represent Fig.5 and Table 2.

##  IV CONCLUSIONS

In summary we have successfully utilized the optoelectronic model based on
experimentally obtained optical features for planar meso-structured perovskite
solar cells. Our optoelectronic model incorporates ray tracing simulations,
the transfer-matrix method and the carrier transport equations. The ray
tracing allows the exploration of any geometry, since the solar collector
device is designed by means of the FreeCAD software. The combination of the
ray tracing with the transfer-matrix method allows to perform the optical
simulations considering Fresnel optics, internal reflections of the multilayer
system, scattering produced by nanoroughness layers, and effective layers with
some peculiarities (such as porous materials).

Large number of transmittance-reflectance measurements on the several stages
of manufacturing process have been used to fix the uncertainties of the PSC
layers thicknesses and the mesoporous TiO2 layer porosity. The same number of
measurements have been done to define the average J-V characteristics of
experimental devices. The obtained results within the optoelectronic model
have been further used to analyse J-V characteristics of the manufactured PSCs
by means of the optimization procedure based on the multiple-use of the
solution transport equations [23]. The developed model demonstrates good
agreement with both optical and electrical measurements leading to extraction
of the layers thicknesses, the porosity degree and the permittivity of the
mesoporous TiO2 material, the perovskite carriers mobility and lifetimes. In
addition, we have found that even small degree of the porosity can affect the
performance of the experimental devices, reducing essentially the effective
donor concentration (i.e., the conductance) of the mp-TiO2 layer in comparison
to that in the c-TiO2 layer.

Our approach allows the study of devices with different architectures.
Reversing the order of the materials in the solar cell structure does not
affect the development of the optoelectronic simulations. Thus our work
provides a useful platform for device optimization towards increasing PCE as
well as appropriate tools for the analysis of various optoelectronic
properties of meso-structured perovskite solar cells.

## Acknowledgements

This work is part of the project TED2021-132758B-I00, funded by MCIN/ AEI
/10.13039/501100011033/ and the European Union “Next Generation EU”/PRTR. It
has also been co-funded by the predoctoral contracts call of the Vice-
presidency and Innovation, Research and Tourism Department of the Government
of the Balearic Islands and the European Social Fund (ESF) [Grant
FPI/2144/2018]. Part of this work is under Physics Ph.D. programme for F.
Bonnín-Ripoll of the Universitat de les Illes Balears (UIB, Spain). We thank
Dr José F. González Morey (Serveis Cientificotècnics, UIB) for his technical
assistance. This work is partially supported by the Russian Science Foundation
(grant no. 23-19-00884). We thank the Spanish State Research Agency for the
grant Self-Power (PID2019-104272RB-C54 / AEI / 10.13039/501100011033), the
ProperPhotoMile Project (PCI2020-112185) and the OrgEnergy Excelence Network
(CTQ2016-81911- REDT), and to the Agència de Gestió d’Ajuts Universitaris i de
Recerca (AGAUR) for the support to the consolidated Catalonia research group
2021 SGR 01617 and the Xarxa d’R+D+I Energy for Society (XRE4S). Part of this
work is under Materials Science Ph.D. Degree for K.T. and the Chemistry Ph.D.
programme for C.P. of the Universitat Autònoma de Barcelona (UAB, Spain). ICN2
is supported by the Severo Ochoa program from Spanish MINECO (grant no.
SEV-2017-0706) and is funded by the CERCA Programme/Generalitat de Catalunya.”

##  Appendix A Supplementary data

Figure 6: Image of the perovskite solar cell samples at the final stage of the
manufacturing process.

### Device Characteristics

The device fabrication and the experimental work, presented in this paper,
have been carried out at the facilities of the Institut Català de Nanociència
i Nanotecnologia (ICN2) in collaboration with the Nanomaterials for
Photovoltaics Energy research group. The perovskite solar cells have been
fabricated with the FTO/c-TiO2/mp-TiO2/Perovskite/Spiro-OMeTAD/Au
configuration using multi-cationic perovskite. Each of these materials has its
own synthesis and deposition process; all fabrication details are discussed in
[24]. Example of final samples is shown in Fig.6.

Figure 7: Mean values and standard deviations of the optical transmittance
measurements for each manufacturing stages of the perovskite solar cell.
Standard deviations have been calculated at each wavelength for six
experimental samples (see also the main text).

### Transmittance

Transmittance measurements have been done by means of the UV-Vis CARY 4000
spectrophotometer at normal incidence, using the solid sample holder
accessory. The incident radiation consists of $200-800$ nm wavelength range of
the optical spectrum. The results of the mean value of the optical
transmittance measurements for each stage as a function of the wavelength are
shown in Fig.7. Our measurements make manifest that it is enough to fabricate
six samples at each stage to obtain a reliable result for the mean value.
Additionaly, we indicate the standard deviation from the mean value for each
curve, corresponding to each stage. The largest deviations are obtained for
the c-TiO2 and mp-TiO2 depositions (stages 2 and 3, respectively). Noteworthy,
up to now there is no any available criterium in literature regarding
definition of the effective border line between c-TiO2 and mp-TiO2
depositions. In fact, one of our goals is to propose such the effective
methodology.

Figure 8: Cross-section image of the samples, obtained by means of the
scanning electron microscope. It is quite visible the degree of the gold layer
roughness. The extracted approximate thickness of each layer is presented in
Table 4.

### Layer Thicknesses

The approximate thickness of each material layers has been estimated using
scanning electron microscopy (SEM) techniques (see Fig.8). These images are
essential for modelling the PSC, since each layer thickness assesment is a key
parameter in the following model calculations. The heterogeneity observed in
the layer thicknesses is due to the experimental process itself, subject to
variations in the initial conditions, the experience and skill of the
technician, the correct stoichiometric balance of the compounds, etc. The less
accuracy for the estimation of the c-TiO2 layer thickness is due to its
extension diapason in the range of a few tens of nanometers.

Material | Thicknesses (nm)
---|---
SiO2 | 10-30
SnO2 | 10-30
FTO | 500-600
c-TiO2 | 10-30
mp-TiO2 | 180-300
Perovskite | 400-550
Spiro-OMeTAD | 200-250
Table 4: Thickness range of the material layers present in the PSC samples
measured by scanning electron microscopy.

Note that the substrate is formed by four different materials. The thickness
of each component in the compound is provided by the manufacturer itself:
$2.2$ mm of glass, $25$ nm of SnO2 and SiO2, and $540$ nm of FTO. Even so, in
the images (see Fig.8) it can be seen that the FTO is not completely
homogeneous. Consequently, we consider a certain range of thicknesses, in
which the calculations for the characterization are carried out. From the
images obtained by SEM, the thickness ranges considered for the materials in
the successive calculations, are shown in Table 4.

### Detailed elaboration of thickness configuration

In order to fit theoretical and experimental transmittances, a concrete
thickness configuration of the PSC structure is needed. Since scanning
electronic microscope images show that the thicknesses of the layers is not
completely homogeneous, we must consider a range between which each layer
thickness lies. Before to resolve this problem we have to analyse the optical
properties of each layer.

One of the most common theoretical tools used to analyse the optical
properties of complex solar systems with a high accuracy is the so called
Monte-Carlo ray tracing simulation. It consists of a set of techniques that
determine the path of light through matter in a three-dimensional environment
with computer simulations [56]. We implement in our approach the Monte-Carlo
ray tracing technique, adding the transfer-matrix method (TMM) to characterize
the optical response of the PSC with the aid of the ray trace OTSun python
package [57, 58, 59]. In order to carry out simulations in a simpler and
faster way, first, we use the TMM only to determine the thicknesses
configuration. Based on the results discussed hereafter, the OTSun will be
used (see the main text) as the subsequent step where more accurate
calculations are needed. Consequently, using the ranges shown in Table 4, each
possible configuration is considered to provide the transmittance with the aid
of the TMM calculation, to obtain those that are closer to the experimental
curves. To reach our goal we introduce the root mean square error (RMSE) of
each thickness configuration for a given material thickness

$RMSE=\frac{1}{N}\sqrt{\sum_{\lambda=300}^{\lambda=800}[T^{exp}(\lambda)-T^{th}(\lambda)]^{2}}\,.$
(12)

Here, $N$ is the total number of wavelength steps; the step is equal to $1$ nm
in all cases; $T^{exp}$($T^{th}$) is the experimental (theoretical) optical
transmittance. Hereafter, it is convenient to introduce the range limits of
the thickness variation for each layer as [a,b] (in nm).

First, we carry out the calculations for the stage 1. To this aim we consider
the following thicknesses: [10,30] for SnO2 and SiO2; [500, 600] for FTO. As a
result, we choose those configurations that lie within 1% of the lowest RMSE.

Next, we proceed to the stage 2 (repeating similar calculations), using the
range limits obtained from the stage 1: [10,20] for SnO2 and SiO2; [550,600]
for FTO; [10,30] for c-TiO2. These calculations reduce the starting range
limits to [10,16], [14,20], [550,600], [22,30] for SnO2, SiO2, FTO, c-TiO2,
respectively.

Figure 9: The short-circuit current $J_{sc}$ as a function of mp-TiO2 porosity
(a), mp-TiO2 thickness (b), perovskite thickness (c) and Spiro-OMeTAD
thickness (d).

Afterwards, we move to the stage 3, where both the thickness and the porosity
of the mp-TiO2 are added to the system. An initial sweep is made considering
the previous ranges plus the mp-TiO2 [180,300] layer and its porosity
considered between 1% and 50%. Considering the configurations that lie within
1% of the lowest RMSE, the material thicknesses are varied in the following
intervals: [10,12] for SnO2; [14,18] for SiO2; [555,580] for FTO; [23,25] for
c-TiO2; [200,300] for mp-TiO2, and its porosity reduced to the 10%–20% range.

From our previous experience (see details in [11, 23]) it follows that small
changes in layer thicknesses (in order of a few nanometers) at the stage 2 do
not affect the optical transmittance of the device. Therefore, we choose the
mean values 11, 16, 565, and 24 nm (associated with SnO2, SiO2, FTO, c-TiO2
layer thickness, respectively), which are defined within the ranges imposed by
the stage 3. By fixing these values, we also reduce significantly the
calculation time, enabling us to explore the mesoporous material and its
subsequent influence on the properties of the absorbing layer.

We recall that the stage 4 includes the perovskite layer considered in the
[400,550] interval, while keeping the variability of the porosity and
thickness of the mp-TiO2 layer as mentioned above. Since there is a complex
interplay of the transmittance through the mp-TiO2 layer and its effect on the
transmittance through the perovskite, we consider the ranges of three
parameters simultaneously: i) the mp-TiO2 layer thickness; ii) its porosity;
iii) the perovskite thickness within 5% of the lowest RMSE selection. As a
result, we obtain that: the mp-TiO2 porosity oscillates between 10% and 20%,
while its thickness between 260 and 300 nm; the perovskite layer thickness
oscillates between 450 and 500 nm. In the stage 5 the Spiro-OMeTAD layer is
added with a thickness defined in [200,250] nm interval. Our calculations
yield that the configurations within 5% of the lowest RMSE lie in the
following interval (in nm): [270,300] for the mp-TiO2 layer; the porosity is
defined between 15% and 20%; [450,470] for the perovskite layer; [200, 250]
for the Spiro-OMeTAD layer.

Notice that the configurations that fit better the optical behaviour of the
experimental samples are not necessarily those that fit better the electronic
properties of the experimental samples. Consequently, our analysis is
supplemented by one more parameter: the short-circuit current $J_{sc}$.
Therefore, considering the ranges of parameters for each material obtained
under the previous procedure and making them less restrictive, we proceed to
the calculations of the short-circuit current $J_{sc}$, which will involve
adding the gold metal layer to the system. Again, we have to vary the
discussed parameters within the above limits (see below). Taking this fact
into account, the maximum possible photocurrent is calculated in the
assumption that all absorbed photons generate electron and holes and
contribute to the current [40]. In an efficient PSC this current is closely
related to the short-circuit current, $J_{sc}$. As a result, to reach a
consistently between optical and electronic properties we calculate the short-
circuit current for every configuration, taking into account that the mean
value of the measured current $J_{sc}=20.07$ mA/cm2 (see details in the main
text).

We recall that all the effects regarding the radiation-matter interaction
mentioned above, are needed to elucidate the number of absorbed photons inside
active materials. To this aim the photogeneration rate $G$ of the creation of
pairs of electrons/holes should be considered a position-dependent function
within the material thickness. To this aim we employ the method proposed by
Ball et al. [60] and the TMM library [27] to calculate the short-circuit
current $J_{sc}$ as:

$J_{sc}=q\int G(\lambda,z)d\lambda dz\,,$ (13)

where $q$ is the electron charge, and $z$ is the depth from the semiconductor
surface.

Figure 10: The short-circuit $J_{sc}$ as a function of the RMSE, calculated by
fitting the experimental and theoretical transmittance of stage 5. The optimal
configurations (red solid points) are located in the range of $\pm 5\%$
deviation from the experimental value of the short-circuit current. Black
point corresponds to the selected case, used for comparison of experimental
and theoretical J-V characteristics (see main text).

Turning to the calculation details of the discussed procedure, first, we fix
the multilayer glass thicknesses at $11$, $16$, $565$ (nm) for SnO2, SiO2,
FTO, respectively; and $24$ nm for c-TiO2. Second, it is observed that the
photocurrent is very sensitive to the thickness and porosity of the mp-TiO2
layer. Indeed, it should be noted that at stages $2$ and $3$ that correspond
to the deposition of the c-TiO2 and the mp-TiO2 layers, respectively, a
standard deviation $10$ times higher relative to those for other stages (see
Fig.7). Including the $J_{sc}$ calculations to the RMSE result obtained from
the transmittance fitting, we weakend the boundaries found above, in
particularly, for the porosity and for the mp-TiO2 layer thickness.
Consequently, we vary the porosity and thickness of the mp-TiO2 layer (from 2%
to 20%, from 200 to 300 nm, respectively); the perovskite thickness (from 450
to 550 nm); the Spiro-OMeTAD thickness (from 200 to 250 nm). As a result, we
do calculations by varying the following parameters:

  1. 1

porosity takes 10 possible values;

  2. 2

mp-TiO2 thickness takes 11 possible values;

  3. 3

perovskite thickness takes 11 possible values;

  4. 4

Spiro-OMeTAD thickness takes 6 possible values,

i.e. we have a set of $10\times 11\times 11\times 6=7260$ combinations
determined by the parameters variations. Altogether the chosen parameters
provide the basis for the calculation of the set $J_{sc}$ values (see Fig.9).

Among all the possible configurations, we have selected those that hold the
following criteria: i) be $\leq|5|\%$ with respect to the experimental value
of the short-circuit current $J_{sc}$; ii) be $\leq 5\%$ of the lowest RMSE
(red region in Fig.10). From this selection, the configuration with the
highest $J_{sc}$ has been chosen denoted by the black point on Fig.10. The
results for the short-circuit current $J_{sc}$ demonstrate: i) its less
sensitivity to the variation of the perovskite and Spiro-OMeTAD thicknesses;
ii) while the close correspondence to the experimental value of the short-
circuit current $J_{sc}$ takes place at the thickness $240$ nm and the
porosity value = $20\%$ of the mp-TiO2 layer.

Summarising, according to our analysis the theoretical configuration, chosen
to model the experimental samples, consists of the following layer thicknesses
(in nm): $11$ (SiO2), $16$ (SnO2), $565$ (FTO), $24$ (c-TiO2), $240$ (mp-TiO2,
the porosity = $20\%$), $500$ (Perovskite), $250$ (Spiro-OMeTAD). This
configuration is obtained as a result of optimization of the optical
transmittance of the device configuration at the value of the short-circuit
current $J_{sc}=19.83$ mA/cm2.

## References

  * Min  _et al._ [2021] H. Min, D. Y. Lee, J. Kim, G. Kim, K. S. Lee, J. Kim, M. J. Paik, Y. K. Kim, K. S. Kim, M. G. Kim, T. J. Shin, and S. Il Seok, Nature 2021 598:7881 598, 444 (2021).
  * Babu  _et al._ [2020] V. Babu, R. Fuentes Pineda, T. Ahmad, A. O. Alvarez, L. A. Castriotta, A. Di Carlo, F. Fabregat-Santiago, and K. Wojciechowski, ACS Applied Energy Materials 3, 5126 (2020).
  * Martynov  _et al._ [2017] Y. B. Martynov, R. G. Nazmitdinov, A. Moià-Pol, P. P. Gladyshev, A. R. Tameev, A. V. Vannikov, and M. Pudlak, Phys. Chem. Chem. Phys. 19, 19916 (2017).
  * Urzúa-Leiva  _et al._ [2020] R. Urzúa-Leiva, A. Narymany Shandy, H. Xie, M. Lira-Cantú, and G. Cárdenas-Jirón, New journal of chemistry 44, 14642 (2020).
  * Frontera  _et al._ [2019] A. Frontera, Y. Martynov, R. G. Nazmitdinov, and A. Moià-Pol, in  _Perovskite Solar Cells: Properties, application and efficiency_ (Nova Science Publishers/ Inc., 2019) Chap. 3, pp. 117–174.
  * Agresti  _et al._ [2019] A. Agresti, A. Pazniak, S. Pescetelli, A. Di Vito, D. Rossi, A. Pecchia, M. Auf der Maur, A. Liedl, R. Larciprete, D. V. Kuznetsov, D. Saranin, and A. Di Carlo, Nat. Mater. 18, 1228 (2019).
  * Martynov  _et al._ [2021] Y. B. Martynov, R. G. Nazmitdinov, P. P. Gladyshev, and A. Moià-Pol, Mendeleev Communications 31, 459 (2021).
  * Li  _et al._ [2021] T. Li, S. Wang, J. Yang, X. Pu, B. Gao, Z. He, Q. Cao, J. Han, and X. Li, Nano Energy 82, 10.1016/J.NANOEN.2021.105742 (2021).
  * Zhang  _et al._ [2021] X. Zhang, J. Wu, S. Wang, W. Pan, M. Zhang, X. Wang, Z. Lan, and J. Lin, ACS Applied Energy Materials 4, 10.1021/ACSAEM.1C01142 (2021).
  * Pereyra  _et al._ [2021] C. Pereyra, H. Xie, and M. Lira-Cantu, Journal of Energy Chemistry 60, 599 (2021).
  * Bonnín-Ripoll  _et al._ [2019] F. Bonnín-Ripoll, Y. B. Martynov, G. Cardona, R. G. Nazmitdinov, and R. Pujol-Nadal, Solar Energy Materials and Solar Cells 200, 110050 (2019).
  * Deng  _et al._ [2019] K. Deng, L. Li, K. Deng, and L. Li, Small Methods 10.1002/smtd.201900150 (2019).
  * Haidari [2019] G. Haidari, AIP Advances 10.1063/1.5110495 (2019).
  * Mishra and Shukla [2020] A. K. Mishra and R. K. Shukla, Materials Today: Proceedings 49, 3181 (2020).
  * Husainat  _et al._ [2020] A. Husainat, W. Ali, P. Cofie, J. Attia, J. Fuller, and A. Darwish, American Journal of Optics and Photonics 8, 6 (2020).
  * Al-Hattab  _et al._ [2021] M. Al-Hattab, L. Moudou, M. Khenfouch, O. Bajjou, Y. Chrafih, and K. Rahmani, Solar Energy 227, 10.1016/J.SOLENER.2021.08.084 (2021).
  * Hussain  _et al._ [2021] S. S. Hussain, S. Riaz, G. A. Nowsherwan, K. Jahangir, A. Raza, M. J. Iqbal, I. Sadiq, S. M. Hussain, and S. Naseem, Journal of Renewable Energy 2021, 10.1155/2021/6668687 (2021).
  * Yadav  _et al._ [2022] R. K. Yadav, P. S. Pawar, R. Nandi, K. R. E. Neerugatti, Y. T. Kim, J. Y. Cho, and J. Heo, Solar Energy Materials and Solar Cells 244, 10.1016/J.SOLMAT.2022.111835 (2022).
  * Ahmad  _et al._ [2022] K. Ahmad, M. Q. Khan, R. A. Khan, and H. Kim, Optical Materials 128, 10.1016/J.OPTMAT.2022.112458 (2022).
  * Karthick  _et al._ [2022] S. Karthick, S. Velumani, and J. Bouclé, Optical Materials 126, 10.1016/J.OPTMAT.2022.112250 (2022).
  * Burgelman  _et al._ [2021] M. Burgelman, K. Decock, A. Niemegeers, J. Verschraegen, and S. Degrave,  _SCAPS-1D Manual_ (2021).
  * Deepthi Jayan and Sebastian [2021] K. Deepthi Jayan and V. Sebastian, Solar Energy 217, 40 (2021).
  * Bonnín-Ripoll  _et al._ [2021] F. Bonnín-Ripoll, Y. B. Martynov, R. G. Nazmitdinov, G. Cardona, and R. Pujol-Nadal, Physical Chemistry Chemical Physics 10.1039/D1CP03313A (2021).
  * Xie  _et al._ [2021] H. Xie, Z. Wang, Z. Chen, C. Pereyra, M. Pols, K. Gałkowski, M. Anaya, S. Fu, X. Jia, P. Tang, D. J. Kubicki, A. Agarwalla, H. S. Kim, D. Prochowicz, X. Borrisé, M. Bonn, C. Bao, X. Sun, S. M. Zakeeruddin, L. Emsley, J. Arbiol, F. Gao, F. Fu, H. I. Wang, K. J. Tielrooij, S. D. Stranks, S. Tao, M. Grätzel, A. Hagfeldt, and M. Lira-Cantu, Joule 5, 1246 (2021).
  * Liu  _et al._ [2016] H. Liu, Z. Huang, S. Wei, L. Zheng, L. Xiao, and Q. Gong, Nanoscale 8, 6209 (2016).
  * Cardona and Pujol-Nadal [2020] G. Cardona and R. Pujol-Nadal, PLOS ONE 15, e0240735 (2020).
  * Byrnes [2017] S. Byrnes, tmm 0.1.7 : Python Package Index (2017).
  * Byrnes [2020] S. J. Byrnes, Multilayer optical calculations (2020), arXiv:1603.02720 [physics.comp-ph] .
  * Fre [2001] FreeCAD: Your Own 3D Parametric Modeler (2001).
  * Liu  _et al._ [2014] F. Liu, J. Zhu, J. Wei, Y. Li, M. Lv, S. Yang, B. Zhang, J. Yao, and S. Dai, Appl. Phys. Lett. 104, 253508 (2014).
  * Hima [2019] A. Hima, Int. J. Energetica 4, 56 (2019).
  * Zhang  _et al._ [2017] F. Zhang, Z. Wang, H. Zhu, N. Pellet, J. Luo, C. Yi, X. Liu, H. Liu, S. Wang, X. Li, Y. Xiao, S. M. Zakeeruddin, D. Bi, and M. Grätzel, Nano Energy 41, 469 (2017).
  * Farhana  _et al._ [2017] A. Farhana, R. Mahbub, S. S. Satter, and S. M. Ullah, Inter. J. Photoenergy 2017, 9846310 (2017).
  * Chen  _et al._ [2019] H. W. Chen, D. P. Gulo, Y. C. Chao, and H. L. Liu, Scientific Reports 2019 9:1 9, 1 (2019).
  * Wang  _et al._ [2017] X. Wang, L.-L. Deng, L.-Y. Wang, S.-M. Dai, Z. Xing, X.-X. Zhan, X.-Z. Lu, S.-Y. Xie, R.-B. Huang, and L.-S. Zheng, J. Mater. Chem. A 5, 1706 (2017).
  * Wojciechowski  _et al._ [2014] K. Wojciechowski, M. Saliba, T. Leijtens, A. Abate, and H. J. Snaith, Energy Environ. Sci. 7, 1142 (2014).
  * Snaith and Grätzel [2006] H. J. Snaith and M. Grätzel, Appl. Phys. Lett. 89, 262114 (2006).
  * Dutta  _et al._ [2015] S. Dutta, Leeladhar, A. Pandey, O. P. Thakur, and R. Pal, J. Vac. Sci. and Tech. A 33, 021507 (2015).
  * Green  _et al._ [2014] M. Green, A. Ho-Baillie, and H. Snaith, Nat. Photonics 8, 506 (2014).
  * García-Cañadas  _et al._ [2006] J. García-Cañadas, F. Fabregat-Santiago, H. J. Bolink, E. Palomares, G. Garcia-Belmonte, and J. Bisquert, Synthetic Metals 156, 944 (2006).
  * Ching-Prado  _et al._ [2018] E. Ching-Prado, A. Watson, and H. Miranda, Journal of Materials Science: Materials in Electronics 29, 15299 (2018).
  * Siefke  _et al._ [2016] T. Siefke, S. Kroker, K. Pfeiffer, O. Puffky, K. Dietrich, D. Franta, I. Ohlídal, A. Szeghalmi, E.-B. Kley, and A. Tünnermann, Advanced Optical Materials 4, 1780 (2016).
  * Chen  _et al._ [2015] C. W. Chen, S. Y. Hsiao, C. Y. Chen, H. W. Kang, Z. Y. Huang, and H. W. Lin, Journal of Materials Chemistry A 3, 9152 (2015).
  * Nichelatti [2002] E. Nichelatti, J. Opt. A: Pure Appl. Opt 4, 400 (2002).
  * Kim  _et al._ [2021] K. P. Kim, W. H. Kim, S. M. Kwon, J. Y. Kim, Y. S. Do, and S. Woo, Nanomaterials 2021, Vol. 11, Page 1233 11, 10.3390/NANO11051233 (2021).
  * Widianto  _et al._ [2021] E. Widianto, Shobih, E. S. Rosa, K. Triyana, N. M. Nursam, and I. Santoso, Optical Materials 121, 10.1016/J.OPTMAT.2021.111584 (2021).
  * Mehrabian  _et al._ [2023] M. Mehrabian, O. Akhavan, N. Rabiee, E. N. Afshar, and E. N. Zare, Environmental Science and Pollution Research 30, 10.1007/S11356-023-26497-1/FIGURES/6 (2023).
  * Hossain  _et al._ [2023] M. K. Hossain, G. F. Toki, A. Kuddus, M. H. Rubel, M. M. Hossain, H. Bencherif, M. F. Rahman, M. R. Islam, and M. Mushtaq, Scientific Reports 2023 13:1 13, 10.1038/s41598-023-28506-2 (2023).
  * Hutchinson  _et al._ [2010] N. J. Hutchinson, T. Coquil, A. Navid, and L. Pilon, Thin Solid Films 518, 2141 (2010).
  * Raoult  _et al._ [2019] E. Raoult, R. Bodeux, S. Jutteau, S. Rives, A. Yaiche, D. Coutancier, J. Rousset, and S. Collin, in  _EU PVSEC Proceedings_, Vol. 36th EU PVSEC 2019 (2019) pp. 757–763.
  * Raoult  _et al._ [2022] E. Raoult, R. Bodeux, S. Jutteau, S. Rives, A. Yaiche, A. Blaizot, D. Coutancier, J. Rousset, and S. Collin, Optics Express 30, 9604 (2022).
  * Hernández-Granados  _et al._ [2019] A. Hernández-Granados, A. N. Corpus-Mendoza, P. M. Moreno-Romero, C. A. Rodríguez-Castañeda, J. E. Pascoe-Sussoni, O. A. Castelo-González, E. C. Menchaca-Campos, J. Escorcia-García, and H. Hu, Optical Materials 88, 695 (2019).
  * Landau  _et al._ [1995] L. Landau, E. Lifshitz, and L. Pitaevskii,  _Electrodynamics of Continuous Media: Volume 8_ , Course of theoretical physics (Elsevier Science, 1995) p. 460.
  * Sellers and Seebauer [2011] M. C. Sellers and E. G. Seebauer, Thin Solid Films 519, 2103 (2011).
  * Herz [2017] L. M. Herz, ACS Energy Letters 2, 1539 (2017).
  * Bos-Coenraad  _et al._ [2017] J. Bos-Coenraad, L. Bunthof, and J. Schermer, Solar Energy 155, 1188 (2017).
  * Cardona and Pujol-Nadal [2018] G. Cardona and R. Pujol-Nadal, OTSun Python Package (2018).
  * OTS [2018] OTSunWebApp (2018).
  * Bonnín-Ripoll  _et al._ [2018] F. Bonnín-Ripoll, G. Cardona, and R. Nadal-Pujol, Tutorial 2 OTSun Web App (2018).
  * M. Ball  _et al._ [2015] J. M. Ball, S. D. Stranks, M. T. Hörantner, Sven Hüttner, Wei Zhang, E. J. W. Crossland, Ivan Ramirez, Moritz Riede, M. B. Johnston, R. H. Friend, and H. J. Snaith, Energy & Environmental Science 8, 602 (2015).