Quasicubic model for metal halide perovskite nanocrystals

Colloidal metal halide perovskite (MHP) nanocrystals (NCs) developed over the last few years exhibit unprecedented, technologically interesting properties. Cesium lead halide (CsPbX₃) NCs have been demonstrated to show bright emission covering the full visible spectral range by virtue of compositional variation and quantum size effects, with 90%+ quantum efficiencies without surface passivation.¹ The exciton radiative emission in these NCs at a low temperature is very fast, with measured radiative decay rates in the sub-nanosecond range,² 2–3 orders of magnitude faster than in other common semiconductor NCs (e.g., CdSe³ or InP⁴). Moreover, the lifetime increases with increasing temperature;⁵ at room temperature, they emit with nanosecond-scale lifetimes,¹ more than an order of magnitude faster than other semiconductor NCs (e.g., CdSe and InP NCs) that have been optimized for over two decades.

The phenomenology described above does not seem to be consistent with the standard model of exciton fine structure in NCs that worked out over the last 25 years, where, owing to the exchange splitting of the exciton, it is expected that there is always a dark, or dipole inactive, exciton ground state.^6,7 In the MHPs, owing to the two-fold degeneracy of the conduction and valence bands, each of which can be described as corresponding to angular momentum 1/2 states in the cubic phase, the exchange interaction is expected to split the ground exciton into a dipole-active (bright) triplet state with total angular momentum 1 and a lower energy, dipole-inactive (dark) singlet with total angular momentum 0. Such a fine structure was considered first in CuCl.⁸ Since the exciton would be expected to cool into the lowest state at low temperature, the lifetime should increase as the temperature is reduced. This discrepancy has spurred a series of experimental investigations directed at the detection of the ground dark exciton conducted on perovskite NC ensembles^9–11 and in single NCs.^12–14 Indirect evidence for a dark ground-state exciton has been reported based on three-level modeling of the temperature dependence of the photoluminescence decays in conjunction with magnetic brightening either with external fields^9,10 or via magnetic impurities,¹¹ and direct observation of the dark ground exciton state via magnetic activation has been reported in single formamidinium lead bromide(FAPbBr₃) NCs.¹⁵ However, similar efforts to directly detect a dark ground exciton in single CsPbBr₃ NCs through magnetic brightening have not been successful.^12,13 The failure to directly detect the expected dark ground exciton state in CsPbX₃ (X = Cl, Br, I) NCs, in conjunction with the unusual temperature dependence of the radiative lifetime in these systems, led to the proposal that the fine structure in CsPbX₃ NCs is determined not by electron-hole exchange but by the Rashba effect, leading to the prediction that the ground exciton state in CsPbX₃ NCs is dipole active.¹⁶ 

Apart from the question of the order of levels of the bright vs dark exciton states, there remains the question of what mechanism—short-range (SR) vs long-range exchange, intrinsic crystal field and/or shape anisotropy, or the Rashba effect—determines the order of levels amongst the optically active triplet exciton states in NCs of CsPbX₃ and other MHPs, and the magnitude of their splitting. Low temperature investigations of single CsPbX₃ NCs show either a two- or a three-line structure with millielectronvolt-scale splittings (with no observation of lower lying dark exciton states).^12–14,16 Three-line spectra observed by several groups are associated with emission from three bright excitons, consistent with orthorhombic symmetry.^2,12,14,16 CsPbBr₃ NCs of apparent tetragonal symmetry consistently show fine structure comprising a low lying linearly polarized line and an upper doublet, which is split in a magnetic field.¹² The observed level order of a singlet below a doublet has been variously ascribed to being due to intrinsic crystal field splitting,^12,14 splitting due to shape anisotropy,^17–19 or the Rashba effect.¹⁹ In the case of models involving intrinsic crystal field splitting, the magnitude and sign of the crystal field have been inferred from the level structure observed in single NC measurements assuming no other effects such as shape anisotropy; the resulting values estimated from CsPbBr₃ NCs^12,14 have the opposite sign from the value calculated for tetragonal CsPbBr₃ using density functional theory (DFT).¹⁹ Another issue is that the magnitude of the measured bright exciton splitting, that is, the splitting within the bright triplet manifold, which is on the order of 1 meV in CsPbBr₃ NCs,^2,12,14 is larger than the full singlet-triplet splitting expected due to the short-range (SR) exchange interaction splitting calculated using density functional theory (DFT).^16,19 This has led to several different proposals to explain the magnitude of the bright triplet level spacing in MHP NCs: that it is due to the confinement enhancement of the short-range exchange interaction;¹⁴ that the contribution of long-range exchange exceeds that of the short-range contribution;^15,17–19 or that the Rashba effect^16,19 determines the bright triplet level spacing in MHP NCs. A quantitative understanding of these various contributions to the exciton fine structure is critical if experimental measurements are to be tied to a microscopic theory of the exciton fine structure. For example, a quantitative exchange model is required to determine whether measurements of the temperature dependence of ensemble photoluminescence, which have been interpreted in conjunction with phenomenological three-level system modeling,^10,11 can be explained solely within an exciton fine structure model.

Complicating the analysis of the experimental data pertaining to exciton fine structure in MHP NCs is the fact that experimental values for key materials parameters needed to build quantitative models are often unavailable. Indeed, effective masses and dielectric constants, which determine the exciton radius and the NC-size dependence of the exciton transition energies in NCs, have commonly been taken from density functional theory (DFT) calculations in order to model size-dependent exciton structure.^1,16,17 Taking CsPbX₃ as an example, however, the reported effective mass values based on DFT^1,16 are significantly different than the recent experimentally measured values in bulk CsPbX₃²⁰. Additionally, due to the absence of direct experimental values, the exchange constants have so far been estimated from DFT^16,19 for the short-range exchange, and for the long-range exchange, the constant has been estimated from effective masses obtained either from DFT calculations^15,17 or experimental measurement.^18,19 In the case of CsPbBr₃, the resulting estimates of the LR exchange parameter ℏω_LT range from 3.5 meV¹⁹ to 6.3 meV,¹⁷ while a recent direct measurement gives 5.4 meV.²¹ 

These uncertainties outlined above underscore the need for a comprehensive analysis of the size-dependent quantum confinement energy and exciton fine structure in MHPs. To this end, in this paper, we develop a quasicubic (QC) model for excitons in MHP NCs. In a quasicubic model, we assume cubic symmetry to start with and consider departures from cubic to tetragonal or orthorhombic symmetry (either associated with a lattice distortion or a shape asymmetry) via perturbation theory. We begin with a coupled-band k ⋅ P model for the carrier confinement energies derived by Kang and Wise²² for spherical NCs, and from it, we derive an isotropic, nonparabolic effective mass model which can be applied to carriers and excitons in MHP NCs. We apply the model to CsPbBr₃ NCs whose size is comparable to or a few times larger than the exciton Bohr radius—that is, in an intermediate confinement regime where correlated electron-hole motion is important—as distinct from the usual strong confinement approximation where NC size is much less than the exciton radius. To establish the validity of the model, we test the coupled-band theory and the effective nonparabolic mass model against absorption data measured in CsPbBr₃ NCs reported by Brennan et al.²³ Using parameters derived from the measured bulk effective mass data of Ref. 20, we find excellent agreement with the Brennan NC data. Next, we develop a quasicubic model for the effect of tetragonal and orthorhombic lattice distortions and apply this model to calculate the exchange-related exciton fine structure in cube-shaped NCs of CsPbBr₃ of tetragonal and orthorhombic lattice symmetry. We use the model to predict the exciton bright triplet level order for several hybrid MHPs of interest and compare the results against DFT calculations. Finally, we introduce the effect of NC shape anisotropy and possible Rashba terms into the model. We conclude by calculating the fine structure of tetragonal CsPbBr₃ NCs, with and without Rashba terms, and discussing the results in the context of available experimental data.

We begin consideration of the quantum-confined carrier states in metal halide nanocrystals by noting that both the conduction and valence band-edge states in the cubic phase perovskites can be characterized by angular momentum J = 1/2.²⁴ As such, the conduction and valence bands are both simple bands without the orbital degeneracy that complicates analysis, for example, of the valence bands in the more familiar III-V and II-VI semiconductors. The simplest model for the quantum confined levels in MHP nanocrystals can therefore reasonably be based on a simple single-band effective mass approximation.

In this approximation, the wavefunction of the ground quantum size level in a cube-shaped NC with edge length L is given by

Assuming parabolic band dispersion, we compute the ground state confinement energy for this state as

where m^* is the effective mass of the carrier. We note that exactly the same ground state energy is obtained for a spherical NC with radius $aeff=L/3$⁠, which has wavefunction in spherical coordinates of the form

where $Y00(θ,ϕ)$ is the spherical harmonic for the state with angular momentum l = 0 and angular momentum projection m = 0 and j_l is the spherical Bessel function with index l = 0. We will therefore refer to the state described by Eq. (1) as s-like. We give the corresponding expressions for the excited p-like state with its nodal plane perpendicular to z,

The top line in the equation above gives the wavefunction and energy for a cube-shaped NC, while the lower line shows the expressions for the corresponding state in a spherical NC. In the last line, the term α_1,1 ≈ 4.493 is the first root of the spherical Bessel function j₁(x). Using the effective radius $aeff=L/3$ in the expression for the spherical NC gives a confinement energy within 2.2% of the result calculated for the cube-shaped NC.

While the parabolic, simple band approach above is certainly valid for sufficiently large NCs, growth and characterization of CsPbBr₃ NCs with sizes down to edge lengths of ∼3 nm, comparable to the typical exciton radius in these materials, have recently been reported.²³ For such small NCs, the confinement energies are appreciable and we expect band nonparabolicity and band coupling effects to become important. A model for treating this problem for spherical NCs with isotropic bands, which utilizes the fact that quantum size levels in a spherical NC are eigenstates of total angular momentum, was developed by Kang and Wise²². Of course a complication in applying this model to MHP nanocrystals is that the fact that MHP nanocrystals have a nanocube morphology with cubic or nearly cubic shape.^1,25 The difference between a cubic shape and a spherical shape can be described in terms of a perturbation; writing the perturbation in terms of a spherical harmonic expansion, the lowest order terms in the perturbation potential have angular momentum L = 4. These terms can have no effect in first order on the lowest energy quantum size levels and only a weak effect in second order since the coupling is via highly excited states.²⁶ It is therefore a good approximation to model cube-shaped MHP nanocrystals as spherical with an effective radius $aeff=L/3$⁠. As noted above, this correspondence is exact for the ground state and within ∼2% for the excited p states within the simple effective mass approximation.

The energy of quantum confined levels is determined by the boundary conditions. Here, we assume impenetrable boundary conditions. For the lowest quantum confined levels, which have total angular momentum F = 1/2, the energies, E, are determined by the relations (see the supplementary material)

where for the s-like ground quantum size level l = 0 and for the first excited p-like states l = 1. Note that these expressions represent the fact that the s-like state in the conduction band is coupled to a p-like state in the valence band and vice versa. The quasimomenta k in Eq. (5) are determined by the bulk band dispersion relation,

In Eqs. (5) and (6), the Kane energy E_p is given by

where the Kane momentum matrix element $P=−i⟨s|P^|Z⟩$ and m₀ is the free electron mass. To connect this picture to the simple effective mass theory, we note that in the 2-band model the band edge effective masses for the conduction and valence bands are given by¹⁹ 

In Fig. 1, we show the confinement energy for the ground, s-like, and first excited, p-like, quantum size levels in the conduction band of CsPbBr₃ nanocrystals, using the experimentally derived effective mass parameters from Ref. 20 which are summarized in Table I. The comparison of the confinement energies calculated using the parabolic band model, Eq. (2), with the 2-band result, shown in Fig. 1, reveals that the two-band model has significantly lower confinement energy than that calculated in the parabolic model. In part, this is due to the nonparabolic dispersion of the bulk bands in the two-band model. This can be accounted for in principle by employing an energy-dependent effective mass. We illustrate this for the conduction band states. Starting with the 2-band dispersion relation, Eq. (6), we rearrange to obtain the expression

where the energy-dependent effective mass m₀/m_c(E) ≡ E_p/3(E − $Ev$⁠) reflects the hyperbolic bulk band dispersion. The energy-dependent effective mass defined in the last equation reflects the hyperbolic bulk band dispersion. When implemented in single band effective mass theory, it leads to the following hyperbolic expression for the confinement energy, $Es=E−Ec$⁠, of the ground state developed by Wang et al. in Ref. 27,

The corresponding expression for the excited p-state is, with $Ep=E−Ec$⁠,

Using these hyperbolic expressions results in confinement energies which approach the full 2-band result more closely than the simple parabolic band approximation, as shown in Fig. 1, but still systematically overestimate the confinement energy calculated in the full two-band model.

The reason for the remaining discrepancy is that in the full two-band model, the ground conduction band state, which is predominantly s-like, is coupled to the p-like valence band state and vice versa for the ground valence band state [see Eq. (5)]. As a result, the dispersion relation which governs the quantum confined levels differs from the bulk dispersion relation, Eq. (9). In the case of the conduction band s-state, we approximate the difference by replacing the term $Ev$ in Eq. (9) by the energy of the lowest p-like valence band quantum size level, namely, $Ev−Ep$⁠, where $Ep$ is taken from Eq. (11). For the excited p-like conduction band state, which couples to the s-like valence band state, we use the replacement $Ev−Es$⁠, where $Es$ is taken from Eq. (10). This leads to the modified nonparabolic expressions for the confinement energies of the ground and first excited states, respectively, denoted as $Es′$ and $Ep′$⁠,

where the energy-dependent effective mass, m′, is given in this “s/p” nonparabolic approximation as

In the last expression, we use Eqs. (10) and (11) for $Es$ and $Ep$ to obtain

Note that, unlike in the hyperbolic approximation, Eqs. (10) and (11), the expressions for the energies of the s- and p-like states in Eqs. (12) and (13) now both involve the same energy-dependent effective mass, m′, defined in Eq. (15). Remarkably, using this expression for the nonparabolic effective mass reproduces the full two-band result very closely in contrast to either the simple parabolic band approximation or the hyperbolic approximation, as shown in Fig. 1.

Having developed expressions for the quantum size levels of electrons and holes in MHP nanocrystals in Sec. II, we now consider the problem of excitons. We begin with a parabolic band approximation as before.

As discussed previously, cube-shaped MHP nanocrystals are now being fabricated with edge lengths L ranging from a few nanometers to ∼15 nm, that is, the grown sizes range from the exciton Bohr radius to several multiples of it. Such NCs are in an intermediate confinement regime,^16,19 where confinement effects on the carrier energies are significant but the motion of the electron and hole are correlated. To capture the correlated motion of the electron and hole, we write the envelope function of the confined excitons using a one-parameter ansatz function,^16,19

Here, ψ_gr is the confined wavefunction found previously for carriers in the strong confinement limit, Eq. (1). This function satisfies the boundary condition that the envelope function vanishes at the NC surface, while the term involving β, the variational parameter, builds correlated electron-hole motion into the ansatz wavefunction. In the expression, N(β) is a normalization factor. In the variational approach, we calculate the expectation value of the two-particle Hamiltonian,

where ϵ_in is the dielectric constant inside the NC. The procedure involves computing the kinetic energy [the expectation value of the first and second terms on the right-hand side of Eq. (17)] and the potential energy [the expectation value of the third term on the right-hand side of Eq. (17)]. Expressing these as dimensionless integrals (given in Ref. 19), we can write the exciton energy E_x relative to the bandgap, E_g, as

For a given NC size L, we minimize the energy with respect to the dimensionless variational parameter b = $β$ L to obtain the exciton energy E_x(L, b_opt), as described in Refs. 16 and 19. We find that the ground state exciton energy smoothly varies between the strong confinement result when L < a_x,

and the weak confinement result when L ≫ a_x,

In these expressions, the exciton binding energy, B_x, is given in terms of the hydrogen Rydberg, Ry, as

μ = m_em_h/(m_e + m_h) is the reduced effective mass given in terms of the electron and hole masses, m_e and m_h, and M_x = m_e + m_h is the total exciton effective mass.

The exciton energy calculated in these limits for parabolic band dispersion, with energy independent electron and hole masses m_e and m_h, is shown in Fig. 2. In the figure, energy is given in units of the exciton binding energy B_x relative to the bandgap E_g and is plotted against the dimensionless quantity $(ax/L)2$⁠. The energy in the intermediate confinement limit is shown in the figure as the solid black line. For small NC sizes, this converges as expected with the result calculated in the strong confinement limit, Eq. (19), plotted as the blue dashed line. For large sizes, it converges to the weak confinement limit, Eq. (20), shown as the red dashed line, approaching the exciton binding energy B_x for large sizes.

From the plot, it is clear that for sizes as large as L < ∼3a_x, the intermediate confinement energy is very close to that calculated in the strong confinement limit. Indeed, a reasonable approximation to the energy in the intermediate confinement limit is to use the kinetic energy term associated with the strong confinement limit in Eq. (19) but to replace the Coulomb energy term [the third term on the right hand side in Eq. (19)] by the square root of the sum of the squares of the exciton binding energy B_x and the Coulomb energy calculated in first order perturbation theory,

The resulting approximate expression for the exciton energy in intermediate confinement is given by

The exciton energy calculated in Eq. (23) is shown in Fig. 2 as the black dashed line and is seen to represent an excellent approximation to the full variational calculation shown in the figure.

The consideration of excitonic effects above was conducted in the parabolic band approximation. To account for band nonparabolicity, we can take the reduced effective mass as energy-dependent, or equivalently as dependent on the NC size, L, using μ(L)⁻¹ = 2/m′(L), where Eq. (15) defines the effective nonparabolic mass in the Kane model. This approach allows us to use Eq. (23) for nonparabolic bands with μ(L)⁻¹ = m_e(L)⁻¹ + m_h(L)⁻¹ and introducing the size dependence of exciton binding energy $Bx(L)=μ(L)/(m0ϵin2)⋅Ry$⁠.

We now compare the models developed above to experimental data. In Fig. 3(a), we show the room temperature absorption peak energy of cube-shaped CsPbBr₃ NCs measured as a function of nanocrystal edge length, L, reported by Brennan et al. in Ref. 23. Panel (b) of the figure shows the same data plotted vs 1/L² to emphasize the departure of the measured data from parabolic dispersion, which would show a linear dependence on 1/L². In the figure, we show the theoretical exciton energy calculated within the intermediate confinement approximation, corresponding to Eq. (18), for several approximations. First, the solid black line shows the intermediate confinement energy, Eq. (18), calculated within the nonparabolic model, using the size-dependent effective mass given in Eq. (15). The T = 300 K parameters used in this calculation are based on the measurements of Yang et al.,²⁰ summarized in Table I. The only parameter adjusted from their measured low temperature bulk values is the room temperature bandgap E_g(300 K) which was adjusted to achieve best fit to the measured room temperature CsPbBr₃ NC data of Ref. 23. To underscore the failure of the parabolic band approximation, we also show as a blue dashed line in the figure the exciton energy calculated for parabolic band dispersion, using Eq. (18) with the band-edge reduced exciton effective mass. This calculation only matches the data for the largest NC sizes; owing to significant band nonparabolicity, a parabolic dispersion model does not match the experimentally measured size dependence. Finally, for reference, we show the result of a full 2-band model, using Eqs. (5) and (6), to determine the kinetic energy and the approximate expression for the potential energy in intermediate confinement, given by Eq. (22). This calculation matches both the data of Brennan et al. and the nonparabolic effective mass model very well with only a slight departure in energy at the smallest NC sizes.

Having developed a model for the quantum confinement of excitons in MHP NCs in Sec. III, we now turn to the description of the exciton fine structure, that is, the millielectronvolt scale structure created by spin-dependent interactions between the electron and the hole. We start the discussion with the short-range (SR) exchange.

The SR exchange can be written in the form of a spin-dependent contact interaction,^19,28

In this expression, $I$ is the 4 × 4 unit matrix, σ_e and σ_h are Pauli operators representing the electron and hole spin (true spin, not their total angular momenta), Ω is the volume of the crystal unit cell, while C is the exchange constant. To obtain the size-dependent SR exchange energies, we simply average the exchange interaction over the total exciton wavefunction, including the underlying periodic part as well as the slowly varying envelope. We write the general form of the total exciton wavefunction, first, for a NC of cubic symmetry,

Here, f(r_e, r_h) is the envelope function for the exciton which describes the electron and hole motion within the NC. Within the intermediate confinement regime, f(r_e, r_h) is given by Eq. (16). The Bloch functions $uJec(re)$ and $uJhv(rh)$ represent the band-edge periodic basis functions for the exciton which contain the spin functions. For a cubic perovskite, as noted previously, the conduction and valence band states each have angular momentum J = 1/2. For the valence band edge, which has s-symmetry, these are the even parity states of angular momentum J = 1/2, which we write $u1/2,±1/2v$⁠, given by¹⁶ 

where the spinor functions ↑ and ↓ are the eigenfunctions of the electron spin projection operator s_z = ±1/2. The spin-orbit interaction splits the conduction band, which has p-symmetry, into lower, band-edge states with angular momentum J = 1/2 and upper states with J = 3/2. The conduction band edge states $u1/2,±1/2c$ are given by¹⁶ 

where X, Y, and Z are the three orbital Bloch functions which transform as the functions having p symmetry. Since there are two conduction band states and two valence band states, it is clear that the exchange interaction, described by Eq. (24), will take the form of a 4 × 4 matrix. Diagonalization of the resulting matrix shows that exchange splits the exciton into an upper, optically active triplet and a ground optically passive dark singlet. The size dependence of the exciton fine structure splitting can be obtained by averaging Eq. (24) over the wavefunction of the confined exciton. Since the envelope function in Eq. (25) is not itself spin-dependent, the result of this averaging takes the form

Here, C is the exchange constant, $M̃SR$ is a dimensionless matrix determined by the Bloch functions, Eqs. (26) and (27), while Θ(L) is the overlap factor representing the probability of the overlap of the electron and hole in the same unit cell,

The integral above is taken over the volume V = L³ of the cube-shaped NC with edge length L. Diagonalization of $M̃SR$ for cubic perovskites yields a singlet level with eigenvalue 0 which corresponds to the dark exciton and 3 degenerate upper levels with eigenvalues equal to 2/3CΘ which are optically allowed.¹⁹ The size dependence of the splitting is given by the overlap factor Θ(L). Using the ansatz function, Eq. (16), in the intermediate confinement regime, we have

where N(b_opt) is the normalization calculated for each L for the variational optimum parameter b = b_opt(L). We show this function in Fig. 4, normalized to the value of the overlap calculated in the bulk, $Θbulk=Ω/πax3$⁠, and plotted vs the ratio L/a_x of the edge length, L, of a cube-shaped NC to the exciton radius a_x. The function N(b_opt) → 1 for L ≪ a_x recovering the well-known 1/L³ dependence of the overlap Θ in the strong confinement limit. However, the exchange splitting calculated this way is seriously underestimated for L > a_x because the correlated motion of the electron and hole is neglected in this limit since the overlap factor approaches the bulk value in the weak confinement limit L ≫ a_x: $Θ→Ω/πax3$⁠. Analytical expressions for the overlap in the strong, weak, and intermediate approximations are derived in Ref. 19.

For NCs with tetragonal and orthorhombic crystal structure, we expect that the 3 upper triplet levels will be split. The splitting is described by matrix $M̃SR$⁠, which will be different for tetragonal and orthorhombic lattice symmetry due to a modification of the Bloch functions of the conduction band edge, relative to those shown in Eq. (27). In Sec. IV A 1, we review the quasicubic model, which shows how the modified Bloch functions depend on the perovskite crystal symmetry for the tetragonal and orthorhombic phases, which serves as a starting point for our understanding of the fine structure.

In the quasicubic model, we consider the effect of the noncubic lattice structure as a simple distortion of the unit cell from the cubic phase.²⁹ The approach is motivated by the observation that the lattice constants of the tetragonal and orthorhombic phases are typically within ≈1 − 2% of those of the cubic α phase, viewed in a nonprimitive supercell aligned to the tetragonal or orthorhombic unit cell, as appropriate. Alternately, the unit cell parameters of the noncubic phases can be expressed in “normalized” form in terms of an equivalent distortion of the lattice parameters of the cubic lattice.³³ Either way, the differences can be quantified by defining strains along the principle axes relative to the cubic lattice. Having defined the strains, we then parameterize the effect on the conduction band edge states using a deformation potential model.³⁰ The strain deformation Hamiltonian is constructed using the theory of invariants as

where U_d is a deformation potential; e_ii are the components of the strain tensor with i running over x, y, z, taken parallel to the crystal principal axes; and L_x,y,z are the matrices representing the x, y, and z projections of angular momentum l = 1. Together with the 3 × 3 unit matrix, $I$⁠, the matrices L_x,y,z serve as base matrices for the deformation potential Hamiltonian for the p-like conduction band with no spin-orbit coupling. Since we are only interested in the splitting caused within the conduction band manifold, we have subtracted off the volume dilatation in Eq. (31). In the basis |l,m⟩ with $l=1$⁠, taken in the order |1,1⟩, |1,0⟩, |1,−1⟩, the deformation potential Hamiltonian then has the form

Equation (32) can be rewritten via symmetry-adapted strain combinations and we define the tetragonal, δ, and orthorhombic, ζ, crystal field parameters accordingly,

For the valence bands, which are known to have s symmetry,²⁴ the only effect of the strain is a volume dilatation term—since the valence bands are orbitally quenched, there is no splitting as in the case of the conduction band; as a result, the valence band Bloch functions are unmodified in the quasicubic model.

To find the modified Bloch functions, we need to add spin-orbit coupling, $ĤSO=23ΔSOL⋅S$⁠, to the quasicubic strain deformation model, where L is the orbital angular momentum and S is the spin. Then, it is straightforward to find the modified Bloch functions for the conduction and valence bands within a six-band model for the conduction electrons. For perovskites of tetragonal symmetry, these functions were found in Refs. 31 and 32. This theory was extended in first order perturbation theory to the case of NCs with orthorhombic symmetry by Fu et al.¹² In the supplementary material, we follow Fu et al. and show that the Bloch functions for the conduction band edge can be written approximately for all three perovskite crystal symmetries (cubic, tetragonal, and orthorhombic) as

Here, the terms X, Y, and Z denote the orbital Bloch functions which transform as the functions having p symmetry. The angles θ and ϕ are defined via the spin-orbit splitting, Δ_SO, of the conduction band and the tetragonal^31,32 and orthorhombic^12,14 crystal fields, δ and ζ, respectively, as follows:

It is worth noting that in the case of cubic symmetry, δ = 0, ζ = 0 so that $tan⁡2⁡θ=22$ or θ ≈ 35.26° and ϕ = 0. In this case, the conduction band Bloch functions reduce to those in Eq. (27). In the tetragonal case, ζ = 0 so that ϕ = 0, and therefore, the conduction band Bloch functions can be represented as eigenstates of the projection of total angular momentum along the c-axis, taken to be along z.^31,32

The spin-orbit coupling parameter Δ_SO can be obtained from experiments (by measuring the spacing between the conduction band edge states and the higher lying states of angular momentum J = 3/2) or can be calculated using DFT (which is the approach taken later on in this paper). The crystal field parameters δ and ζ can also be determined by calculation using DFT. Absorption measurements can be used to determine the splitting as suggested in Ref. 14; however, if this approach is taken, polarization resolved measurements are required in order to distinguish which of the J = 3/2 upper conduction band states moves upwards and which moves downwards under the action of the crystal field, which is required to determine the sign of the crystal field splitting. For example, in the tetragonal case, a positive crystal field causes the J_z = ±3/2 states, which are inactive for light polarized along the ẑ axis, to move upwards in energy while the J = 3/2, J_z = ±1/2 states, which are active for ẑ polarized light, to move downwards.¹⁹ This trends become transparent in the no-spin-orbit coupling case, where the Z state is shifted downwards relative to the X and Y orbital states for a positive crystal field; see Eqs. (32) and (33).

Averaging Eq. (24) over the four conduction and valence band Bloch function products, |c₁⟩|$v↑$⟩, |c₁⟩|$v↓$⟩, |c₂⟩|$v↑$⟩, and |c₂⟩|$v↓$⟩, we obtain the short-range exchange operator in the same from as Eq. (28), but now the matrix $M̃SR=M̃SR(θ,ϕ)$ is a function of the Bloch function phase angles which encode the symmetry. This matrix is given by

Within the exchange model, the matrix $M̃SR(θ,ϕ)$ again determines the level order and the oscillator transition strengths of each level, independently of the envelope wavefunction of the confined exciton. Diagonalizing the matrix $M̃SR$ in Eq. (37) for angles θ and ϕ defined in Eqs. (35) and (36) results in energies for each exciton eigenstate X_i, where X_i = {X, Y, Z, D}, expressed in the general form

Here, the energy scale is set by CΘ, while $ẼXi(θ,ϕ)$ is a dimensionless eigenvalue of the matrix $M̃SR$⁠. Note that in the equation above, we have made the replacement, ℏω_st = 2/3CΘ_bulk, where ℏω_st is the bulk exciton singlet-triplet splitting and $Θbulk=Ω/(πax3)$ is the overlap factor for the bulk. The functions $ẼXi$ are easy to evaluate and are given within the k · P model in Table II for the cubic, tetragonal, and orthorhombic phases. The average of $ẼXi(θ,ϕ)$ taken over the 3 triplet states X, Y, Z is 2/3 in all cases, as shown in the table. As a result, the average singlet-triplet splitting always equals 2/3CΘ. In the cubic and tetragonal phases, ϕ = 0, while ϕ ≠ 0 in the orthorhombic phase.

The transition dipole matrix elements for each exciton sublevel are found as the matrix element of the momentum operator, $p^$⁠, between the exciton eigenfunction, $ψXi$⁠, and the crystal ground state, G. The result can be expressed for each exciton eigenstate, X_i, in the general form¹⁹ 

where P is the Kane momentum matrix element defined above, $p̃Xi(θ,ϕ)$ is a dimensionless function describing the dependence of the matrix element on the crystal structure, and $O$ is an overlap integral between the electron and hole wavefunctions, given by

where the integral is taken over the entire NC volume V. Analytical expressions for $p̃Xi(θ,ϕ)$ for NCs with different crystal symmetries—cubic, tetragonal, and orthorhombic—are shown in Table II. Using the dipole matrix elements defined in Eq. (39), we can calculate the oscillator strength for each exciton state X_i,

where f₀ = E_p/ℏω is the oscillator transition strength of a free carrier transition at the band edge, approximately equal to f₀ ≈ E_p/Eg, $O$ is the overlap integral [see Eq. (40)], and finally, $f̃Xi=|p̃Xi(θ,ϕ)|2$ is a reduced oscillator strength function that determines the relative oscillator strength of each exciton sublevel. Notably, within the six-band k · P/quasicubic model for the conduction bands, $f̃Xi$ is the same function of θ and ϕ as the function $ẼXi$ that appears in the expressions for the energy, Eq. (38).¹⁹ Later, in this paper, we will check the validity of this result for several materials (CsPbBr₃, FAPbBr₃, and methylammonium lead iodide, MAPbI₃) using DFT. We find that while the relationship is approximately true for the tetragonal phases of these materials, the relationship typically fails for the more complicated orthorhombic phases.

Having developed the quasicubic model to describe the exciton fine structure created by lattice distortion, we will estimate the SR exchange constant and the crystal field parameters for several materials of interest. DFT will be used for these calculations as well as for the testing of the quasicubic model. We start however with a qualitative estimation of the effect of the lattice distortion within the quasicubic model.

Using the expressions in Table II for the energy of each exciton state $ẼXi$⁠, and the definitions of the phase angles θ, ϕ in Eqs. (35) and (36) in terms of the crystal field parameters δ, ζ in Eq. (33), we can write expressions for the fine structure energies simply in terms of the lattice constants of each material. To do this, we note that the expressions for the exchange energies in Table II can be linearized for small crystal field parameters δ ≪ Δ_SO, ζ ≪ Δ_SO. The lattice distortion causes the dimensionless exciton energies to deviate from the average value 2/3 by an amount ΔẼ_Z,X,Y = Ẽ_Z,X,Y − 2/3 which can be written as

In Eq. (42), we have expressed the crystal field parameters δ, ζ in terms of the deformation potential and the lattice strains relative to the cubic phase. It bears emphasis that the assumption underlying the model is that the crystal field can be described in terms of a single deformation potential. To the extent that this assumption is valid, the result above shows that we can predict the order of levels solely on the basis on the quasicubic strain, which is determined by the lattice constants of a given material.

Given a material with tetragonal or orthorhombic lattice structure, we must first determine the strain associated with each principle axis. To consider the electronic structure of the tetragonal or orthorhombic phases in the quasicubic approximation, we need to transfer from the primitive unit cell of the cubic phase to a nonprimitive cell aligned to the tetragonal or orthorhombic primitive cell. For example, for the orthorhombic lattice structure, it is convenient to transfer to a nonprimitive cell of the cubic lattice spanned by vectors a, b, c with dimensions $2a0:2a0:2a0$⁠, where a_o is the cubic phase lattice constant, and the vectors a and b are rotated by 45° to the cubic phase x, y axes, while the vector c is aligned to the cubic phase z axis and twice the length of the primitive vector. Equivalently, we can define scaled quasicubic lattice parameters for the tetragonal and orthorhombic phases which represent a distorted primitive cubic cell.³³ These scaled cell parameters, a_i, are equal to b_i/γ_i, where b_i are the lattice parameters of the given noncubic phase (tetragonal or orthorhombic) and the factor γ_i is the dimension of the corresponding nonprimitive cubic cell in the absence of any lattice distortion from cubic symmetry. In the orthorhombic case, for example, $γx=γy=2$⁠, while γ_z = 2. Finally, to eliminate volume dilatation, rather than computing the strain relative to the cubic phase, we can simply compute the quasicubic strain relative to the parameter $āN=(axayaz)1/3$⁠, which is the lattice parameter of a primitive cubic cell with the same volume as the scaled noncubic cell volume,

In Table III, we show this analysis for tetragonal CsPbBr₃ and FAPbBr₃ both corresponding to space group P4/mbm, No. 127. Also for comparison, we analyze tetragonal MAPbI₃ (space group I4/mcm, No. 140). Assuming a positive deformation potential in Eq. (42) for all three materials, we expect that CsPbBr₃ and MAPbI₃ will have a fine structure with a singlet Z-polarized exciton higher in energy than the degenerate X/Y polarized doublet, while FAPbBr₃ is expected to have the opposite level order, namely, the Z-polarized exciton should lie below the X/Y polarized doublet. Anticipating results to be developed later, we show in the table the level orders calculated using DFT for tetragonal CsPbBr₃ as well as mimic structures corresponding to the experimental FAPbBr₃ and MAPbI₃ lattice structures in which the organic cations have been replaced by Cs. In each case, the level order agrees with expectations based on the quasicubic model, validating the qualitative applicability of the quasicubic model for the tetragonal structures as well as the assumption of a positive deformation potential constant.