Photovoltaic effect in multi-domain ferroelectric perovskite oxides
Ying Shi Teh and Kaushik Bhattacharya
Division of Engineering and Applied Science
California Institute of Technology
Pasadena, CA 91125 USA
Abstract
We propose a device model that elucidates the role of domain walls in the photovoltaic effect in multi-
domain ferroelectric perovskites. The model accounts for the intricate interplay between ferroelectric
polarization, space charges, photo-generation and electronic transport. When applied to bismuth ferrite,
results show a significant electric potential step across both 71
◦
and 109
◦
domain walls, which in turn
contributes to the photovoltaic (PV) effect. We also find a strong correlation between polarization and
oxygen octahedra tilts, which indicates the nontrivial role of the latter in the PV effect. The domain
wall-based PV effect is further shown to be additive in nature, allowing for the possibility of generating
above-bandgap voltage.
1 Introduction
In conventional photovoltaics, electron-hole pairs are created by the absorption of photons that are then
separated by an internal field in the form of heterogeneous junctions such as p-n junctions. Less than a
decade ago, Yang
et al.
[1] reported large photovoltages generated in thin films of multi-domain bismuth
ferrite (BFO), and suggested a new mechanism where electrostatic potential steps across the ferroelectric
domain walls drives the photocurrent. This discovery has since revitalized the field of photoferroics. Among
many ferroelectric oxides, BFO has particularly attracted considerable interest due to its high ferroelectric
polarization and relatively small bandgap.
Many novel experiments were subsequently devised to investigate the role of domain walls in the observed
photovoltaic (PV) effect in ferroelectric perovskites. Alexe and Hesse [2] performed measurements of the local
photoelectric effect using atomic force microscopy (AFM). They found that the photocurrent is essentially
constant across the entire scanned area, hence indicating the absence of the domain wall (DW) effect. The
nanoscale mapping of generation and recombination lifetime using a method combining photoinduced tran-
sient spectroscopy (PITS) with scanning probe microscopy (SPM) points to a similar conclusion [3]. This
led to the hypothesis that the bulk photovoltaic (BPV) effect, which arises from the noncentrosymmetry of
perovskites, is the key mechanism instead [4, 5]. However further recent studies focused on characterizing
both the BPV and DW effects show that the latter effect is much more dominant [6, 7]. The lack of clear
concensus among the scientific community on the key mechanism in the PV effect in perovskites as well as
on the role of domain walls could be understood from the inherent difficulties in the experimental techniques.
The nanoscale order of ferroelectric domain walls makes it difficult to probe into and separate the effects
from the bulk domains and the domain walls. Other issues such as defect formation and grain boundaries in
perovskite crystals further complicate the analysis.
First-principles calculations have provided a detailed understanding of the structure of domain walls
[8, 9, 10], and have established the drop in electrostatic potential across it. However, they are limited to a
few nanometers, and cannot examine the interaction of domain walls with other features. On the other hand,
models at the device scale provide understanding at the scale [11], but assume
a priori
the polarization and
other aspects of the domain wall. Finally, various phase field models provide understanding of the domain
pattern [12, 13], but in the absence of space charge and photocurrent. Thus, there is a gap in our modeling
of the intricate interplay between space charge, ferroelectric polarization and electronic transport.
This paper seeks to fill this gap by building on prior work of Xiao
et al.
[14, 15] and Suryanarayana and
Bhattacharya [16] who developed a continuum theory of semiconducting ferroelectrics including electron and
1
arXiv:1811.07948v1 [cond-mat.mtrl-sci] 19 Nov 2018
Figure 1: Crystal structure of bulk BFO. The two O
6
octahedra rotate out-of-phase about the polarization
axis marked by the dotted line.
hole transport. We extend their work to include photogeneration due to illumination and study photovoltaic
effect in ferroelectric perovskite oxides. We investigate the photovoltaic response of BFO films with different
domain wall configurations by solving the model at the device scale. At a 71
◦
or 109
◦
domain wall, we
observe a change in the component of polarization perpendicular to the domain wall. This in turn results in
a relatively large electrostatic potential step across the wall which allows for separation of photogenerated
electron-hole pairs. Thus this model supports the hypothesis of domain walls contributing to the photovoltaic
effect.
We emphasize that this model does not
a priori
assume the domain wall structure or the electrostatic
potential step across it. Instead, this is a prediction of the model that is based on well-established Devonshire-
Landau models of ferroelectrics and lumped band models of semiconductors.
The rest of the paper is organized as follows: In Section 2, we briefly review the structure of BFO and
discuss the classifications of ferroelectric domain walls in BFO. We develop the theoretical framework in
Section 3. We apply the theory to examine PV effect in ferroelectrics with domain walls in Section 4. We
conclude with a brief discussion in Section 5.
2 Bismuth Ferrite
In this work, we focus on bismuth ferrite (BFO), though we note that the same framework can be applied to
other ferroelectrics and the results are expected to be similar qualitatively.
At room temperature, BFO has a rhombohedral phase with space group
R
3
c
(Figure 1). The displace-
ments of the atoms from the ideal cubic structure in this phase lead to a spontaneous polarization pointing
in the [111] pseudocubic direction. Another distintive feature is the network of O
6
octahedra surrounding
the Fe ions that rotate or tilt out-of-phase about the polarization axis. This is also commonly known as the
antiferrodistortive (AFD) mode and has been found to play an important role in the ferroelectric phase of
the material [10].
Electric polarization in rhombohedral BFO can take one of the eight variants of the [111] pseudocubic
direction which gives possible domain wall orientations of 71
◦
, 109
◦
and 180
◦
. On each domain, there can be
two possible orientations of oxygen octahedra. We follow Lubk
et al.
’s method of classifying oxygen octahedra
tilts (OTs) across domain walls as either
continuous
or
discontinuous
[9]. In the continuous case, the direction
of oxygen octahedra tilt remains the same along the polarization vector field. In the discontinuous case, the
direction reverses across the domain wall.
3 Theory
We consider a metal-perovskite-metal (MPM) structure that is connected to an external voltage source to
form a closed electrical circuit (see Figure 1). The multi-domain ferroelectric perovskite film occupying the
2
Figure 2: Schematic of a device model in a metal-perovskite-metal configuration.
space Ω is subjected to light illumination. The two metal-perovskite interfaces are denoted by
∂
Ω
1
,
∂
Ω
2
∈
∂
Ω. All the processes are assumed to occur at constant temperature
T
. We present the equations and their
physical meanings here. Readers may refer to Appendix A for the thermodynamically consistent derivation.
3.1 Charge and electrostatic potential
The total charge density (
x
∈
Ω) is given by
ρ
=
q
(
p
v
−
n
c
+
z
d
N
+
d
−
z
a
N
−
a
)
,
(1)
where q is the electronic charge,
n
c
is the density of electrons in the conduction band,
p
v
is the density of
holes in the valence band,
N
+
d
is the density of ionized donors,
N
−
a
is the density of ionized acceptors,
z
d
is
the valency number of donors, and
z
a
is the valency number of acceptors. The polarization and space charge
in the ferroelectrics together generate an electrostatic potential. This is determined by Gauss’ equation
∇·
(
−
ε
0
∇
φ
+
p
) =
ρ,
(2)
where
ε
0
is the permittivity of free space, subject to appropriate boundary conditions.
3.2 Transport equations
In the presence of light illumination, an incident photon may be absorbed in the semiconductor to promote an
electron from the valence band to the conduction band, thus generating an electron-hole pair in the process
of photogeneration. The reverse may also occur such that an electron and a hole recombine. Electrons and
holes may also move from one point to another point, as represented by the electron and hole density flux
terms,
J
n
and
J
p
. With conservation of electrons and holes, we can relate the time derivatives of densities
of electrons and holes to the aforementioned processes via the following transport equations,
̇
n
c
=
−∇·
J
n
+
G
−
R,
(3)
̇
p
v
=
−∇·
J
p
+
G
−
R,
(4)
where
G
is the rate of photogeneration which can be taken to be proportional to the intensity of light
illumination, while
R
is the recombination rate. Here we assume that the only form of recombination present
is radiative recombination, which involves the transition of an electron from the conduction band to the
valence band along with the emission of a photon.
R
takes the form of
B
(
n
c
p
v
−
N
2
i
) with the intrinsic
carrier density being given by
N
i
=
√
N
c
N
v
exp(
−
E
c
−
E
v
2
k
B
T
) [17]. The radiative recombination coefficient
B
is
a material property, independent of the carrier density.
The electron and hole fluxes
J
n
,
J
p
are taken to be proportional to the gradient in its electro-chemical
potential [17] via
J
n
=
−
1
q
ν
n
n
c
∇
μ
n
,
(5)
J
p
=
−
1
q
ν
p
p
v
∇
μ
p
,
(6)
where
ν
n
and
ν
p
are the electron and hole mobilities respectively.
In this work, the diffusion of donors and acceptors are neglected.
3
3.3 Free energy
The free energy of the ferroelectric is postulated to be of the form
W
=
W
DL
(
p
,θ
) +
W
G
(
∇
p
,
∇
θ
) +
W
n
c
(
n
c
) +
W
p
v
(
p
v
) +
W
N
d
(
N
+
d
) +
W
N
a
(
N
−
a
)
.
(7)
The various terms are currently explained.
W
DL
refers to the Devonshire-Landau free energy of bulk ferroelectrics. In addition to the typical primary
order parameter of electric polarization
p
, we include a second order parameter — oxygen octahedral tilts
θ
.
We adopt the following energy form for BFO [12]. The corresponding coefficients can be found in Table 1.
W
DL
=
a
1
(
p
2
1
+
p
2
2
+
p
2
3
) +
a
11
(
p
4
1
+
p
4
2
+
p
4
3
) +
a
12
(
p
2
1
p
2
2
+
p
2
2
p
2
3
+
p
2
1
p
2
3
)
+
b
1
(
θ
2
1
+
θ
2
2
+
θ
2
3
) +
b
11
(
θ
4
1
+
θ
4
2
+
θ
4
3
) +
b
12
(
θ
2
1
θ
2
2
+
θ
2
2
θ
2
3
+
θ
2
1
θ
2
3
)
+
c
11
(
p
2
1
θ
2
1
+
p
2
2
θ
2
2
+
p
2
3
θ
2
3
) +
c
12
[
p
2
1
(
θ
2
2
+
θ
2
3
) +
p
2
2
(
θ
2
1
+
θ
2
3
) +
p
2
3
(
θ
2
1
+
θ
2
2
)]
+
c
44
(
p
1
p
2
θ
1
θ
2
+
p
1
p
3
θ
1
θ
3
+
p
2
p
3
θ
2
θ
3
)
.
(8)
The energy stored in the ferroelectric domain walls is accounted for through the gradient or Ginzburg
energy term
W
G
which includes the energy cost associated with rapid change in polarization and octahedral
tilts.
W
G
(
∇
p
,
∇
θ
) =
1
2
a
0
|∇
p
|
2
+
1
2
b
0
|∇
θ
|
2
.
(9)
Here we assume that the gradient terms are isotropic for simplicity but can easily be modified.
W
n
c
,W
p
v
,W
N
d
,W
N
a
in equation (7) are the free energies of electrons in the conduction band, holes in the
valence band, donors and acceptors respectively. The explicit expressions of these energies can be determined
by considering each system as a canonical ensemble in the framework of statistical mechanics[16]
W
n
c
(
n
c
) =
n
c
E
c
+
k
B
T
[
−
N
c
log
N
c
+
n
c
log
n
c
+ (
N
c
−
n
c
) log(
N
c
−
n
c
)]
,
(10)
W
p
v
(
p
v
) = (
N
v
−
p
v
)
E
v
+
k
B
T
[
−
N
v
log
N
v
+
p
v
log
p
v
+ (
N
v
−
p
v
) log(
N
v
−
p
v
)]
,
(11)
W
N
d
(
N
+
d
) =(
N
d
−
N
+
d
)
E
d
−
(
N
d
−
N
+
d
)
k
B
T
log(2
z
d
)
+
k
B
T
[
−
N
d
log
N
d
+
N
+
d
log
N
+
d
+ (
N
d
−
N
+
d
) log(
N
d
−
N
+
d
)]
,
(12)
W
N
a
(
N
−
a
) =
N
−
a
E
a
−
(
N
a
−
N
−
a
)
k
B
T
log(2
z
a
)
+
k
B
T
[
−
N
a
log
N
a
+
N
−
a
log
N
−
a
+ (
N
a
−
N
−
a
) log(
N
a
−
N
−
a
)]
.
(13)
3.4 Polarization and tilt equations
The polarization and tilt evolve according to the (time-dependent) Landau-Ginzburg equations
1
ν
p
̇
p
=
∇·
∂W
G
∂
∇
p
−
∂W
DL
∂
p
−∇
φ,
(14)
1
ν
θ
̇
θ
=
∇·
∂W
G
∂
∇
θ
−
W
DL
∂θ
.
(15)
where
ν
p
,ν
θ
are the respective mobilities. They are subject to natural boundary conditions
ˆn
·
∂W
g
∂
∇
p
= 0
,
(16)
ˆn
·
∂W
g
∂
∇
θ
= 0
.
(17)
4
3.5 Electrochemical potentials
The electrochemical potentials are obtained from the energy to be
μ
n
=
∂W
n
c
∂n
c
−
qφ,
(18)
μ
p
=
∂W
p
v
∂p
v
+
qφ,
(19)
μ
N
+
d
=
∂W
N
d
∂N
+
d
+
qz
d
φ,
(20)
μ
N
−
a
=
∂W
N
a
∂N
−
a
−
qz
a
φ.
(21)
At thermal equilibrium,
μ
n
=
−
μ
p
=
−
μ
N
+
d
=
μ
N
−
a
=
E
F
m
where the magnitude of
E
F
m
is the workfunction
of the metal electrode. Further, using equations (10) to (13) we can invert the relations to obtain
n
c
=
N
c
1 + exp(
E
c
−
E
F
m
−
qφ
k
B
T
)
,
p
v
=
N
v
1 + exp(
E
F
m
−
E
v
+
qφ
k
B
T
)
,
N
+
d
=
N
d
[
1
−
1
1 +
1
2
z
d
exp(
−
E
F
m
+
E
d
−
qφz
d
k
B
T
)
]
,
N
−
a
=
N
a
[
1 + 2
z
a
exp(
−
E
F
m
+
E
a
−
qφz
a
k
B
T
)
]
−
1
,
(22)
consistent with Fermi-dirac distribution [18]. Finally, assuming
N
c
>> n
c
and
N
v
>> p
v
, equations (5) and
(6) become
J
n
=
−
ν
n
k
B
T
q
∇
n
c
+
ν
n
n
c
∇
φ,
(23)
J
p
=
−
ν
p
k
B
T
q
∇
p
v
−
ν
p
p
v
∇
φ.
(24)
Equations (23) and (24) show that each of
J
n
and
J
p
can be resolved into two contributions: (1) a diffusion
current, driven by concentration gradient of carriers, and (2) a drift current, driven by an electric field. By
applying the Einstein relation which relates diffusion constant
D
to mobility
ν
via
D
=
νk
B
T/q
, we recover
the equations that are typically written to describe the flow of electrons and holes in solar cells.
3.6 Ohmic boundary conditions
We prescribe ohmic boundary conditions at the contacts with metal electrodes following [19] for convenience.
We have
n
c
=
N
c
e
−
(
E
c
−
E
f
m
)
/k
B
T
p
v
=
N
v
e
−
(
E
f
m
−
E
v
)
/k
B
T
}
on
∂
Ω
1
∪
∂
Ω
2
.
This is equivalent to assuming that the Fermi level in the semiconductor aligns with that of the metal, hence
giving rise to electron and hole densities that are independent of the applied voltage.
5
3.7 Steady-state model
At steady state, all the fields of interest do not vary with respect to time. Further, we are interested in domain
walls, and therefore can assume that things are invariant parallel to the domain wall. This means that we
have one independent space variable which we denote
r
. We denote the components of polarization and tilt
parallel (respectively perpendicular) to the domain wall to be
p
s
,θ
s
(respectively
p
r
,θ
r
). With
z
d
=
z
a
= 1,
we have a coupled system of differential equations for region
x
∈
(0
,L
), where
L
is the length of the film.
a
0
d
2
p
r
dr
2
−
∂W
DL
∂p
r
−
dφ
dr
= 0
,
(25)
a
0
d
2
p
s
dr
2
−
∂W
DL
∂p
s
= 0
,
(26)
b
0
d
2
θ
r
dr
2
−
∂W
DL
∂θ
r
= 0
,
(27)
b
0
d
2
θ
s
dr
2
−
∂W
DL
∂θ
s
= 0
,
(28)
−
ε
0
d
2
φ
dr
2
+
dp
r
dr
=
q
(
p
v
−
n
c
+
N
+
d
−
N
−
a
)
,
(29)
−
dJ
n
dr
+
G
−
B
(
n
c
p
v
−
N
2
i
) = 0
,
(30)
−
dJ
p
dr
+
G
−
B
(
n
c
p
v
−
N
2
i
) = 0
,
(31)
J
n
=
−
ν
n
k
B
T
q
dn
c
dr
+
ν
n
n
c
dφ
dr
,
(32)
J
p
=
−
ν
p
k
B
T
q
dp
v
dr
−
ν
p
p
v
dφ
dr
,
(33)
where
N
+
d
=
N
d
[
1
−
1
1 +
1
2
exp(
−
E
F
m
+
E
d
−
qφ
k
B
T
)
]
,
N
−
a
=
N
a
[
1 + 2 exp(
−
E
F
m
+
E
a
−
qφ
k
B
T
)
]
−
1
,
N
i
=
√
N
c
N
v
exp(
−
E
c
−
E
v
2
k
B
T
)
,
with boundary conditions
dp
r
dr
(
r
= 0) =
dp
r
dr
(
r
=
L
) = 0
,
dp
s
dr
(
r
= 0) =
dp
s
dr
(
r
=
L
) = 0
,
dθ
r
dr
(
r
= 0) =
dθ
r
dr
(
r
=
L
) = 0
,
dθ
s
dr
(
r
= 0) =
dθ
s
dr
(
r
=
L
) = 0
,
φ
(
r
= 0) = 0
, φ
(
r
=
L
) = 0
,
n
c
(
r
= 0) =
n
c
(
r
=
L
) =
N
c
e
−
(
E
c
−
E
f
m
)
/k
B
T
,
p
v
(
r
= 0) =
p
v
(
r
=
L
) =
N
v
e
−
(
E
f
m
−
E
v
)
/k
B
T
.
6
Symbols
Values
Units
a
1
−
1
.
19
×
10
9
V m C
−
1
a
11
9
.
93
×
10
8
V m
5
C
−
3
a
12
3
.
93
×
10
8
V m
5
C
−
3
b
1
−
1
.
79
×
10
10
V m
−
3
C
b
11
1
.
14
×
10
11
V m
−
3
C
b
12
2
.
25
×
10
11
V m
−
3
C
c
11
1
.
50
×
10
10
V m C
−
1
c
12
7
.
50
×
10
9
V m C
−
1
c
44
−
1
.
60
×
10
1
−
V m C
−
1
Table 1: Coefficients of Laudau-Devonshire energy for BFO
3.8 Numerical issues
The model derived above comprises of differential equations that are nonlinear and coupled, which can prove
troublesome numerically. So we non-dimensionalize the problem as in Appendix B. Further, we notice that
the coupling between the first five governing equations and the rest of the model is weak. This is especially
so when the length of the simulated device is much smaller than the Debye length, or when the dimensionless
quantity
δ
is small, which is generally the case in the simulations in this paper. Therefore we treat them
as two subproblems, that are then solved self-consistently until convergence occurs. Each subproblem is
constructed within the finite difference framework, and the resulting system of nonlinear equations is solved
using the trust-region dogleg method.
4 Application to Bismuth Ferrite
4.1 Material constants
The coefficients of the Devonshire-Landau energy for BFO in equation (8) are presented in Table 1. They
are derived to match the values of spontaneous polarization, tilt angles, and dielectric constant [20, 10, 21].
Other material parameters including band structure information [22] and carrier mobility values [23] are
listed in Table 2. The values of
a
0
and
b
0
are chosen to match a ferroelectric domain wall width of 2 nm.
Typically BFO exists as a n-type semiconductor due to oxygen vacancies. It can also become p-type with Bi
deficiency. Here we restrict our simulations to n-type semiconductors.
4.2 Two-domain ferroelectrics
4.2.1 71
◦
and 109
◦
domain walls
We consider a device comprising of a BFO film with two ferroelectric domains separated by either a 71
◦
or
109
◦
domain wall, with continuous or discontinuous oxygen octahedra rotations across the DW. This gives a
total of four different cases, as illustrated in Table 3.
Figures 3 and 4 show the variation of various field quantities when the perovskite film is exposed to light
illumination and shorted. Notice that in all cases, the perpendicular component of the polarization
p
r
is not
constant in the vicinity of the domain wall. In other words, the polarization is not divergence free, and we
see a voltage drop across the domain wall. The polarization profile (i.e.
p
r
) of 71
◦
and 109
◦
domain walls
with continuous OT are qualitatively similar to those obtained from first-principles calculations [9]. This
voltage drop across the domain wall leads to charge separation of photogenerated electron-hole pairs, and a
non-zero photocurrent. This is evident in the current-voltage plots shown in Figure 5 and is consistent with
the mechanism proposed by Yang
et al.
[1].
Figure 5 shows that the magnitude and direction of photocurrent due to the domain wall effect hinge
greatly upon the changes in the crystallographic structure across the domain wall. The case with 109
◦
DW
and continuous OT gives a positive short-circuit current, which is in the same direction as net polarization
7