Theory for Light Scattering from Powders

The theory underlying the Python package PDielec is described. PDielec calculates the infrared absorption characteristics of a crystalline material supported in a non-absorbing medium. PDielec post processes solid-state quantum mechanical and molecular mechanical calculations of the phonons and dielectric response of the crystalline material. The molecular and solid-state quantum mechanical (QM) calculations of response properties, such as the frequencies and intensities of infrared (IR) and terahertz (THz) radiation absorption, have become generally available in many molecular and solid-state computer programs. A common approach is to assume the harmonic approximation and calculate the mass-weighted force constant matrix (for molecules) or the dynamical matrix at the gamma point (for periodic solids). Diagonalisation of the matrix gives the frequencies for absorption and the normal modes (molecules) or phonon displacements (periodic solids). The calculation of the absorption intensity for each mode requires the calculation of the change in dipole moment caused by the displacement of the atoms for that mode. For solids where there is a large separation of charge, there can be a large coupling between a phonon mode and the internal field within a particle resulting from its morphology.

PDielec post processes the output of solid-state quantum mechanical and molecular mechanics (MM) based codes such as; VASP [1], CASTEP [2], CRYSTAL [3], Abinit [4], Quantum Espresso [5], Phonopy [6], and GULP [7], to predict the infrared absorption of crystalline insulator materials whose crystal size is small compared with the wavelength of the absorbing radiation. The package is suited for the calculation of the complex, frequency-dependent permittivity and its associated absorption of infrared radiation for a finely ground crystalline material dispersed in a low-loss dielectric medium such as KBr or Polytetrafluoroethylene (PTFE). A particular feature of the program is its ability to take into account the constant permittivity of the supporting medium and the particle shape of the material of interest through an effective medium theory. The documentation outlines the theory used by the program and gives some examples of the application of the program for ionic and molecular materials.

Equations (3) and (4) describe Beer-Lambert’s law [18] where \(\alpha\) is the (decadic) absorption coefficient (usually given in cm^-1), \(I\) and \(I_0\) are the intensities after and before absorption respectively and \(d\) is the path length.

(3)\[\frac{I}{I_{0}} = 10^{- \alpha d}\]

(4)\[log\left( \frac{I}{I_{0}} \right) = -\alpha d\]

It is common, especially in the chemistry community, when reporting infrared spectra to use a decadic molar absorption coefficient (\(a\)), which has units of L·mol^-1·cm^-1. The relationship between the absorption coefficient and the molar absorption coefficient [18] is;

(5)\[\alpha = aC\]

where \(C\) is the concentration of the absorbing species.

Molecular approach to absorption intensity

For molecules, the transition intensity \(I_k\) of the \(k^{th}\) vibrational mode (calculated from the change in dipole moment along the mode displacement) can be converted to an integrated molar absorption coefficient, \(A_k\), which can then be more readily compared with experiment. The theory for this is described by Wilson, Decius and Cross [19] and results in expressions such as those given in equations (6) and (7). The first expression shows the relationship between the integrated molar absorption coefficient and the transition intensity and uses the number of molecules per unit volume (\(N\)), the velocity of light (\(c\)) and the degeneracy of the mode (\(g_k\)). The second expression shows the appropriate conversion factors if the units for the integrated molar absorption coefficient are L·mol^-1·cm^-2 (1 L·mol^-1·cm^-2 = 0.01 km·mol^-1) and the units for the transition intensity are D²Å^-2amu^-1, where D represents the Debye unit of dipole moment and amu is an atomic mass unit. The factor \(log_e10\) arises due to the choice of a decadic Beer’s law.

(6)\[ A_k = \frac{N\pi}{3c^{2}\log_e10}g_kI_k\]

(7)\[A_k = \frac{Na\pi}{3000c^{2}2.302585}g_kI_k = 4225.6I_k\]

The derivation of the above expressions assumes that the rotational levels are not quantised and that the vibrational levels are thermally occupied according to a Boltzmann distribution. To use the calculated molecular intensities to predict a spectrum it is usual to assume [19] that each transition is associated with a Lorentzian line shape with a full width at half maximum (FWHM) of \(\sigma_k\). It is common, when reporting comparisons between theoretical and experimental spectra, to assume that the line widths are the same for all modes [20, 21]. Recent work on terahertz absorption in crystalline pentaerythritol tetranitrate (PETN) using molecular dynamics calculations [22] in combination with the direct calculation of the cubic anharmonic couplings of the normal modes [22] has shown that the FWHM of the intense absorptions may vary between 10 and 25 cm^-1. Assuming a Lorentzian line shape, the molar absorption coefficient for the \(k^{th}\) mode at wavenumber, \({\bar{\nu}}_{k}\), can be written as a function of frequency or wavenumber (\(\bar{\nu}\));

(8)\[a_k(\bar{\nu}) = \frac{2A_k}{\pi}\frac{\sigma_k}{4\left( \bar{\nu} - {\bar{\nu}}_k \right)^{2} + \sigma_k^2}\]

(9)\[a_k^{max} = \frac{2A_k}{\pi\sigma_k}\]

The maximum height of the Lorentzian, \(a_{k}^{\max}\) depends upon the value of \(\sigma_k\). As can be seen in Equation (10), the choice of normalisation for the Lorentzian means that integration of the molar absorption coefficient over wavenumber returns the integrated molar absorption coefficient and a sum over all the bands provides the total molar absorption coefficient \(a^{mol}(\bar{\nu})\) as a function of wavenumber, calculated from the intensities of each band. Equation (12) shows the relationship between the absorption and the molar absorption coefficients. \(C\) is the concentration usually expressed in mol·L^-1.

(10)\[A_k = \int{a_k(\bar{\nu})d\bar{\nu}}\]

(11)\[a^{mol}(\bar{\nu}) = \sum_k{a_k(\bar{\nu})}\]

(12)\[\alpha^{\text{mol}}(\bar{\nu}) = Ca^{\text{mol}}(\bar{\nu})\]

A comment should be made about the various units which can be used for these quantities. A common unit for the transition intensity is D²Å^-2amu^-1, and another is km·mol^-1. However, it should be pointed out that strictly speaking the latter unit refers to the integrated molar absorption coefficient as defined above in Equation (10) and therefore relies on the assumptions made in its derivation. ( 1 D²Å^-2amu^-1 is equivalent to 42.256 km·mol^-1 ).

Solid state approach to absorption intensity

The optical properties of a solid are determined by its complex, frequency-dependent relative permittivity, \(\boldsymbol{\varepsilon}(\bar{\nu})\), and in particular the real and imaginary components, \(\boldsymbol{\kappa}\) and \(\mathbf{n}\), of the complex refractive index, \(\tensorbf{N}{(\bar{\nu})}\) where;

(13)\[\tensorbf{N}({\bar{\nu}})^{2} = {\tensorbs{\varepsilon}}({\bar{\nu}})\]

(14)\[\tensorbf{N}\left( \bar{\nu} \right) = {\tensorbf{n}}\left( \bar{\nu} \right) + i{\tensorbs{\kappa}}\left( \bar{\nu} \right)\]

The intensity of absorption is given by the imaginary component of the refractive index which for an isotropic material is [23];

(15)\[\begin{split}I = I_{0}e^{- \frac{4\pi\kappa(\bar{\nu})d}{\lambda}} \\ I = I_{0}e^{- 4\pi\bar{\nu}\kappa(\bar{\nu})d}\end{split}\]

(16)\[\begin{split}- ln\left( \frac{I}{I_{0}} \right) = 4\pi\bar{\nu}\kappa(\bar{\nu})d \\ - log\left( \frac{I}{I_{0}} \right) = 4\pi\bar{\nu}\kappa(\bar{\nu})d \cdot log(e)\end{split}\]

Comparison with the definition of the absorption coefficient from Beer-Lambert’s law, Equation (3), and using Equation (12) gives;

(17)\[\alpha^{sol}(\bar{\nu}) = 4\pi\bar{\nu}\kappa(\bar{\nu}) log(e)\]

(18)\[a^{sol}(\bar{\nu}) = \frac{\alpha^{\text{sol}}(\bar{\nu})}{C}\]

Since the refractive index is dimensionless, the absorption coefficient, \(\alpha^\text{sol}\) is specified in cm^-1. The superscripts ‘sol,’ for solid, and ‘mol,’ for molecular, are used here to distinguish between the two methods of calculating the absorption (\(\alpha\)) and molar absorption coefficients (\(a\)). In the calculation of the imaginary component of the refractive index it is necessary to choose the solution which gives a positive value. This is consistent with the Kramers-Kronig relationship between the real and imaginary components [24].

To calculate the relationship between absorption and molar absorption coefficients it is necessary to know the concentration. For solid-state calculations, the required unit is; moles of unit cells per liter. One of the drawbacks of this molar absorption coefficient unit is that the number of molecules in a unit cell can change depending on whether a supercell, primitive or non-primitive unit cell is being used. A more natural unit would be to use a mole of formula units, or a mole of molecules. To aid comparison between such calculations PDielec can calculate concentration in both moles of atoms and moles of molecules. However for the rest of this paper Equation (19) will be used, where \(V\) is the volume of the unit cell, and therefore the concentration \(C\) is moles of unit cell/liter.

(19)\[C = \frac{f \cdot 1000cm^{3}}{VN_{a}}\]

The volume fraction, \(f\), of the dielectric material in a supporting matrix of non-absorbing material is included in the expression for the concentration as it will be useful when the theory for mixtures is developed.

For a periodic system, the permittivity tensor can be calculated as a sum over Lorentz oscillators, incorporating an imaginary loss component through the damping factor \(\sigma_k\) [25]. The frequencies of the oscillators are the transverse optic (TO) phonon frequencies of the system.

(20)\[{\tensorbs{\varepsilon}}(v) = {{\tensorbs{\varepsilon}}}_{\infty} + \frac{4\pi}{V}\sum_{k}\frac{{{\tensorbf{S}}}_{k}}{\nu_{k}^{2} - \nu^{2} - i\sigma_{k}\nu}\]

(21)\[\tensorbf{S}_{k} = {\bar{\mathbf{Z}}}_{k}{\bar{\mathbf{Z}}}_{k}^\text{T}\]

(22)\[\bar{\mathbf{Z}}_{k} = \sum_{a}{\tensorbf{Z^{a}}}{\bar{\mathbf{U^{a}}}}_{k}\]

(23)\[\begin{split}\tensorbf{D} {\bar{\mathbf{U}}}_{k} = \Lambda_{k}{\bar{\mathbf{U}}}_{k}\\ \nu_{k}^{2} = \Lambda_{k}\end{split}\]

(24)\[I_{k} = tr\left( \tensorbf{S}_{k} \right)\]

\(V\) is the volume of the unit cell, \({\tensorbf{S}}_{k}\) is the dipole oscillator strength tensor for the \(k_{th}\) transition, with a TO frequency of \(\nu_k\) and \(\tensor{\varepsilon}_{\infty}\) the optical permittivity tensor, which represents the electronic contribution to the permittivity. The intensity of a transition, \(I_k\), is given by the trace of the oscillator strength tensor, Equation (24). The damping factor \(\sigma_{k}\) removes any discontinuities at the TO frequencies. Since the oscillator strengths and phonon frequencies can be calculated routinely in solid-state quantum mechanical packages, the calculation of the frequency-dependent complex permittivity using Equation (20) is straightforward.

In some cases, using Equations (21) and (22), PDielec calculates the oscillator strengths from the Born charge matrix for atom \(a\) and its contributions to the \(k_{th}\) phonon mode [25]. As shown in Equation (23), at the \(\Gamma\) point the \(k_{th}\) phonon mode is described by the eigenvector, \(\bar{\mathbf{U}}_{k},\) and eigenvalue, \(\Lambda_{k}\), of the mass-weighted, dynamical matrix, \(\tensorbf{D}\), which is a 3Nx3N matrix, where N is the number of atoms in the unit cell. The eigenvalues are the squared frequencies of the phonon modes, Equation (23). The displacement of each atom in the \(k_{th}\) is proportional to \(m_{a}^{-1/2}\), where \(m_a\) is the mass of atom \(a\). The dynamical matrix has 3N eigenvectors and eigenvalues, of which three should be zero due to translational invariance. If there are any negative eigenvalues the system is unstable to some displacement and therefore not at an energy minimum.

For ionic systems, it is common practice in solid-state QM and MM programs to include a long wave-length, non-analytic correction to the mass-weighted dynamical matrix at the \(\Gamma\) point, which describes the coupling of the longitudinal optic (LO) modes to the induced field resulting from the vibration. This may be written for atoms \(s\) and \(t\) and their Cartesian components \(\alpha\) and \(\beta\) as [25];

(25)\[\left( \tensorbf{D}^\text{LO}_{\mathbf{q \rightarrow 0}} \right)_{s,\alpha;t,\beta} = \left( \tensorbf{D} \right)_{s,\alpha;t,\beta}\mathbf{+}\frac{4\pi}{V\sqrt{M_{s}M_{t}}}\mathbf{\ }\frac{\left( {\bar{\mathbf{q}}}^{\text{T}}{{\mathbf{\ }\tensorbf{Z}}}_{s} \right)_{\alpha}\left( {\bar{\mathbf{q}}}^{\text{T}}{\tensorbf{Z}}_t \right)_\beta}{{\bar{\mathbf{q}}}^{\text{T}}\mathbf{\cdot}{{\tensorbs{\varepsilon}}}_{\infty} \cdot \bar{\mathbf{q}}}\]

The mass weighting has been incorporated through the mass of the atoms, \(M_s\) and \(M_t\). The correction depends upon the direction, \(\bar{\mathbf{q}}\), that the long wave-length limit is approached. Diagonalisation of the corrected matrix gives the squared frequencies of N-1 LO modes and 2N-2 TO modes, Equation (23). In some of the examples given below the LO frequencies will be given for comparison with the TO frequencies.

Effect of particle shape on infrared absorption

It has long been recognised that, especially for ionic materials, the local field within a crystal and its coupling with the transverse optical phonons has an important effect on the position and intensity of the absorption. Fröhlich [26] was one of the first to point out that the frequency of absorption of a small ionic sphere embedded in a low dielectric medium is shifted to lie between the transverse and longitudinal optical frequencies of the material making up the sphere.

In the development of the theory used in PDielec an important assumption is that the particle size of the crystallites in the sample is small compared with the wavelength of light. Using this approach Genzel and Martin [27] were able to explain the observed infrared absorption of small spheres of MgO crystallites and the effect of the permittivity of the supporting medium on the spectrum. Studies of the infrared absorption by small particles of α-Fe₂O₃ using an effective medium theory and an absorption/scattering theory [28, 29] showed that it was necessary to adjust not only the damping factors in Equation (20) but also the permittivity of the matrix and the volume fraction of the dielectric medium to fit the experimental data. The volume fraction was used to account for aggregation effects as the volume fraction increased. It was also shown that effective medium theories were only applicable to particles smaller than the wavelength of light. For larger particles, the scattering from the particles becomes increasingly important.

More recently Balan and others in a series of papers [30, 31, 32] used density functional calculations together with an effective medium theory to calculate the infrared absorption of several minerals incorporating information about the crystallite shape. In an experimental and theoretical study of irradiated kaolinite [32] it was shown that exposure to radiation resulted in shifts in the infrared spectrum which could be accounted for by increasing the polarisability of the particles through an increase in the optical permittivity tensor.

The underlying theory adopted by PDielec is based on similar premises to the work described above, namely that the dielectric response of small spherical, ellipsoidal, slab-like or needle-like crystallites randomly distributed in a non-absorbing medium such as PTFE, KBr or Nujol, is the same as that of an effective medium material whose frequency-dependent dielectric response can be calculated from the frequency-dependent permittivity tensor of the crystal (as calculated by solid-state QM or MM calculations), the shape of the crystallites and the permittivity of the non-absorbing medium (taken to be a constant over the frequency range of interest).

The development of the theory reported here closely follows the work of Sihvola [33]. It will be assumed that the inclusion particles, which may be non-isotropic, ellipsoidal (including spherical, needle-like and plate-like), are randomly orientated in an embedding, non-absorbing medium such as PTFE, KBr or Nujol. It should be emphasized that, whilst PDielec can take account of particle shape, particle and matrix permittivity many additional aspects of infrared absorption that need to be considered when comparing calculated and experimental results. Most notable of these are; the coupling between phonons and mobile electrons or holes (so-called phonon-polariton coupling) [34], the scattering which starts to dominate as the particles get larger [32] and the agglomeration of particles as the volume fraction increases.

The polarisability of an isolated particle

Figure Fig. 12 shows a schematic of the field and polarisation inside an inclusion with non-isotropic permittivity \(\tensor{\varepsilon}_{i}\) embedded in a supporting medium with permittivity, \(\varepsilon_e\). The internal field within the inclusion is indicated by \(\fieldbf{E}_{i}\) the external, applied field is indicated by \(\fieldbf{E}_{e}\) and the induced polarisation in the inclusion is shown by \(\fieldbf{P}\).

Fig. 12 Schematic showing the various terms used to define the effective medium of an inclusion

The figure shows the field and polarisation inside an inclusion with non-isotropic permittivity \({\tensorbf{\varepsilon}}_{i}\) embedded in a supporting medium with permittivity \(\varepsilon_e\). The internal field within the inclusion is indicated by \(\fieldbf{E}_i\), the external, applied field is indicated by \(\fieldbf{E}_e\) and the induced polarisation in the inclusion is shown by \(\fieldbf{P}\)

The electric field internal to the inclusion gives rise to a polarisation density which is no longer necessarily aligned with the field because the material is non-isotropic. The polarisation density in the inclusion can be expressed as the tensor product of the permittivity contrast between the inclusion and the supporting medium and the (as yet unknown) internal field.

(26)\[\fieldbf{P} = \left( \tensorbs{\varepsilon}_{i} - \varepsilon_{e}{\tensorbf{1}} \right){\fieldbf{E}}_{i}\]

For any ellipsoidal shape (including sphere, slab and needle) with volume V, the polarisation density throughout the particle is uniform and integrating over all space gives the field-induced dipole moment of the inclusion, \(\fieldbf{p}\).

(27)\[\fieldbf{p} = V\fieldbf{P} = V\left( \tensorbs{\varepsilon}_{i} - \varepsilon_{e}\tensorbf{1} \right)\fieldbf{E}_{i}\]

The dipole and the external field, \({\fieldbf{E}}_{e}\), are related by the polarisability tensor, \(\tensorbs{\alpha}\).

(28)\[\fieldbf{p} = \tensorbf{\alpha}\fieldbf{E}_{e}\]

Equations (27) and (28) allow the determination of the polarisability, once the field internal to the inclusion has been expressed in terms of the shape of the inclusion and its permittivity. The polarisation within the inclusion gives rise to a depolarisation field, \(\fieldbf{E}_{d}\), which depends on the shape of the inclusion through the symmetric and unit-trace depolarisation tensor, \(\tensorbf{L}\).

(29)\[\fieldbf{E}_d = - \frac{1}{\varepsilon_e}\tensorbf{L} \cdot \fieldbf{P}\]

The internal field is the sum of the external field and the depolarisation field.

(30)\[\fieldbf{E}_{i} = \fieldbf{E}_{e} + \fieldbf{E}_{d}\]

The depolarisation matrix is determined from the shape of the crystallite. For an ellipse, slab or needle, a unique axis [hkl] using the crystal axes of the unit cell is specified. PDielec transforms these to a cartesian coordinate system using the unit cell lattice vectors. In the case of a slab morphology, the unique direction is normal to the surface specified by its Miller indices (hkl). The definitions of the various depolarisation tensors are indicated in Table below.

Table 5 Definitions used for the depolarisation tensor

Shape	Depolarisation Tensor
Sphere	\(\tensorbf{L} =\frac{1}{3}\left( \fieldbf{V}_1 \fieldbf{V}_1^T + \fieldbf{V}_2 \fieldbf{V}_2^T+\fieldbf{V}_3 \fieldbf{V}_3^T \right)\)
Slab	\(\tensorbf{L}=\fieldbf{V}_1 \fieldbf{V}_1^T\)
Needle	\(\tensorbf{L}=\frac{1}{2}\left( \fieldbf{V}_2\fieldbf{V}_2^T + \fieldbf{V}_3 \fieldbf{V}_3^T \right)\)
Ellipsoid	\(\tensorbf{L}=a\fieldbf{V}_1 \fieldbf{V}_1^T + b\fieldbf{V}_2 \fieldbf{V}_2^T + b\fieldbf{V}_3 \fieldbf{V}_3^T\)

The three directions defined by \(\fieldbf{V}_1\), \(\fieldbf{V}_2\) and \(\fieldbf{V}_3\) are mutually orthogonal cartesian vectors calculated from [hkl] for an ellipse or needle and (hkl) for a slab. In the case of a slab, needle or ellipsoid, \({\bar{V}}_{1}\) defines the unique direction and the other vectors are orthogonal to it. For the case of an ellipsoid, the parameters a and b in the Table depend on the ratio, \(z\), of the unique axis length to the length of an axis perpendicular to it [33].

For z > 1 the ellipsoid is prolate;

(31)\[e = \sqrt{1 - z^{- 2}},\ a = \frac{\left( 1 - e^{2} \right)}{2e^{3}}\left( \log\frac{1 + e}{1 - e} - 2e \right),\ b = \frac{1}{2}\left( 1 - a \right)\]

For z < 1 the ellipsoid is oblate

(32)\[ e = \sqrt{z^{- 2} - 1},\ a = \frac{\left( 1 + e^{2} \right)}{e^{3}}\left( e - \operatorname{}e \right),\ b = \frac{1}{2}\left( 1 - a \right)\]

The depolarisation tensor acts as a projection or screening operator describing the effect of the geometry of the inclusion on the depolarisation field which results from its polarisation. For instance, in the case of a needle, only polarisation perpendicular to the needle axis contributes to the depolarising field, whilst for a slab only polarisation perpendicular to the slab face may contribute. Similarly, for a sphere, all directions contribute and so the depolarisation matrix is diagonal with a value of 1/3 for each diagonal element, as the trace of the depolarisation tensor must be 1. It follows from Equations (26), (29) and (30) that;

(33)\[\fieldbf{E}_i = \fieldbf{E}_e - \frac{1}{\varepsilon_e} \tensorbf{L} \left( \tensorbs{\varepsilon}_i - \varepsilon_e{\tensorbf{1}} \right){\fieldbf{E}}_i\]

Rearrangement allows the internal field of the inclusion to be expressed in terms of the known permittivities, the shape of the inclusion and the external field.

(34)\[\fieldbf{E}_i \left( \varepsilon_e{\tensorbf{1}} + \tensorbf{L} \left( \tensorbs{\varepsilon}_i - \varepsilon_e{\tensorbf{1}} \right) \right) = \varepsilon_e\fieldbf{E}_e\]

(35)\[\fieldbf{E}_i = \varepsilon_e{\fieldbf{E}}_e\left( \varepsilon_e\tensorbf{1} + \tensorbf{L} \left( \tensorbs{\varepsilon}_i - \varepsilon_e \tensorbf{1} \right) \right)^{- 1}\]

Substituting the internal field expression, (33), into Equation (27) for the dipole moment and requiring the dipole moments calculated using the polarisation density to equal those calculated from the polarisability allows the polarisability to be written;

(36)\[\tensorbs{\alpha} = V\varepsilon_e\left( \tensorbs{\varepsilon}_i - \varepsilon_e\tensor{1} \right)\left( \varepsilon_e\tensorbf{1} + \tensorbf{L} \left( \tensorbs{\varepsilon}_i - \varepsilon_e\tensorbf{1} \right) \right)^{- 1}\]

Although it has not been specified explicitly the permittivity of the inclusion, and therefore the polarisability tensor, are frequency dependent through the oscillator strengths of each phonon mode contributing to the permittivity according to Equation (20). The calculation of the complex, frequency-dependent polarisability tensor of the composite material is the key step in the calculation of its effective permittivity.

The effective permittivity of a mixture

To extend this approach to include the effect of a number of inclusions we need to introduce the concept of an effective permittivity, \(\tensorbf{\varepsilon}_{eff}\), which describes the behaviour of an average field, \(\left\langle \fieldbf{E} \right\rangle\), where the angle brackets indicate an average over a volume of the composite material. It is required that the average electric flux density \(\left\langle \fieldbf{D} \right\rangle\) is the same in the effective medium as in the composite medium;

(37)\[\left\langle \fieldbf{D} \right\rangle = \tensorbs{\varepsilon}_{eff}\left\langle \fieldbf{E} \right\rangle = \varepsilon_e\left\langle \fieldbf{E} \right\rangle + \left\langle \fieldbf{P} \right\rangle\]

The averaging is necessary because the polarisation within a given inclusion has an effect on the field in other inclusions. The local field in the cavity left by a single inclusion embedded in the average polarisation field is given by;

(38)\[\fieldbf{E}_L = \left\langle \fieldbf{E} \right\rangle + \frac{1}{\varepsilon_e}\tensorbf{L}\left\langle \fieldbf{P} \right\rangle\]

The local field ‘excites’ the inclusion resulting in a dipole moment \(\mathbf{p}_{mix}\) that is related to the polarisation through the number density of inclusions, \(n\), and through the polarisability of the inclusion, which is already known from Equation (36).

(39)\[\left\langle \fieldbf{P} \right\rangle = n{\fieldbf{p}}_{mix} = n\left\langle \tensorbs{\alpha}\fieldbf{E}_{L} \right\rangle\]

The angle brackets around the product of the polarisability and the local field indicate that it is necessary to average the polarisation according to the distribution of alignments of inclusions. In this work, it will be assumed that the inclusions are randomly aligned. Substituting the expression for the local field, Equation (38), gives;

(40)\[\left\langle \fieldbf{P} \right\rangle = \left( \tensorbf{1} - \frac{n\left\langle \tensorbf{\alpha}\tensorbf{L} \right\rangle} {\varepsilon_e} \right)^{- 1} n\left\langle \tensorbs{\alpha} \right\rangle\left\langle \tensorbf{E} \right\rangle\]

Mixing rules

There are many mixing rules which have been proposed to describe the homogenization of composite materials and a lot of work has been done in comparing their accuracy. Here two methods will be used. The first and the most commonly used method is the Maxwell-Garnett mixing rule [33]. Indeed this has been implied by the use of Equation (37) to define the effective permittivity. The other commonly used method is the Bruggeman mixing rule [33], which differs substantially in the way the two components of the system are treated. It is usually stated that the Maxwell-Garnet mixing rule is good for low-volume fractions of the inclusion and the Bruggeman approach should be better for higher-volume fractions [35]. In addition to these mixing rules one other approach will be described, namely the Averaged Permittivity (AP) mixing rule, which calculates the absorption spectrum ignoring the effects of the internal field on the absorption and can therefore be used as an indicator of the shifts in frequency and intensity which have occurred owing to the effect of particle shape.

Maxwell-Garnett mixing rule

The Maxwell-Garnett approach for treating the properties of heterogeneous mixtures assumes that the average field and the average flux density result from volume fraction weighted sums. Substituting Equation (40) into Equation (37) gives the Maxwell-Garnett effective permittivity;

(41)\[\tensorbs{\varepsilon}_{mg} = \tensorbf{1} + \left( \tensorbf{1} - \frac{n\left\langle {\tensorbs{\alpha}}{\tensorbf{L}} \right\rangle}{\varepsilon_e} \right)^{- 1}n\left\langle \tensorbs{\alpha} \right\rangle\]

The fact that the polarisability tensor has a volume term in it (Equation (36)) means that the terms in Equation (41) containing \(n\tensorbs{\alpha}\) depend on the volume fraction \(f\). Although written as a tensor, because the assumption has been made that the inclusions are randomly orientated, the effective permittivity has to be diagonal with equal complex values. Since the polarisability is complex and frequency dependent, the effective permittivity is also and its calculation using Equations (41) and (36) need to be calculated over the frequency range of interest.

Bruggeman mixing rule

In the Maxwell-Garnett mixing formalism there is a distinction between the inclusion and the supporting medium which results in an asymmetry in the treatment of the two species in the mixture. Instead, the Bruggeman mixing rule assumes that each species is polarised against the background of the effective medium and therefore the polarisation in the two components cancels each other out;

(42)\[\left\langle \fieldbf{P}_1 \right\rangle + \left\langle \fieldbf{P}_2 \right\rangle = 0\]

where the components are now labelled 1 and 2 rather than external and internal. The polarisation for species 1 and 2 with a number density of species represented by \(n_1\) and \(n_2\) can be obtained from the polarisability of the species (Equation (42));

(43)\[\left\langle \fieldbf{P}_1 \right\rangle = n_1 \left\langle \tensorbs{\alpha}_1 \right\rangle\fieldbf{E}\]

Substituting Equation (43) into Equation (42) leads to the requirement that;

(44)\[n_1\left\langle \tensorbs{\alpha}_1 \right\rangle + n_2\left\langle \tensorbs{\alpha}_2 \right\rangle = 0\]

Taking Equation 19 and generalising it for species \(i\), (where \(i\) is 1 or 2) embedded in an effective permittivity given by \(\tensorbs{\varepsilon}_{br}\);

(45)\[\tensorbs{\alpha}_i = V_i \tensorbs{\varepsilon}_{br} \left( \tensorbs{\varepsilon}_i - \tensorbs{\varepsilon}_{br} \right) \left( \tensorbs{\varepsilon}_{br} + \tensorbf{L} \left( \tensorbs{\varepsilon}_i - \tensorbs{\varepsilon}_{br} \right) \right)^{- 1}\]

Equation (44) has to be solved for \(\tensorbs{\varepsilon}_{br}\). Since the systems considered here are isotropic with random inclusions, a solution has to be found for a complex value of the Bruggeman permittivity at each frequency considered. An issue in the use of Equation (45) is that the same depolarisation matrix is being used for both species, which is not always appropriate. The solution of Equation (44) can be achieved either by iteration or by casting the equation as a minimisation problem. The iterative approach implemented in PDielec involves repeated application of (46) until convergence [36]. The starting point for the iterations is taken as the Maxwell-Garnett solution for the first frequency and then the solution at the previous frequency is used to start the next iteration.

(46)\[\tensorbs{\varepsilon}_{br} = \frac{f_1{ {\tensorbs{\varepsilon}}}_1\left\lbrack {\tensorbf{1}} + {\tensorbf{L}}\left( {\tensorbs{\varepsilon}}_1 - \tensorbs{\varepsilon}_{br} \right) \right\rbrack^{- 1} + f_2\tensorbs{\varepsilon}_2\left\lbrack \tensorbf{1} + \tensorbf{L}\left( \tensorbs{\varepsilon}_{2} - \tensorbs{\varepsilon}_{br} \right) \right\rbrack^{- 1}} {f_{1}\left\lbrack \tensorbf{1} +{\tensorbf{L}}\left( \tensorbf{\varepsilon}_1 - \tensorbs{\varepsilon}_{br} \right) \right\rbrack^{- 1} + f_{2}\left\lbrack \tensorbf{1} + \tensorbf{L}\left( \tensorbs{\varepsilon}_2 - \tensorbs{\varepsilon}_{br} \right) \right\rbrack^{- 1}}\]

Although the Bruggeman permittivity is written here as a tensor, the polarisabilities in Equation (44) have to be averaged over the random orientation of the inclusions and therefore the homogenized material is isotropic with a single complex value for the diagonal tensor. Also, as with the Maxwell-Garnett mixing rule, since the polarisability is complex and frequency dependent, the effective permittivity is also, and its calculation using Equation (41) needs to be performed over the frequency range of interest.

The choice between using the Bruggeman or Maxwell-Garnett model is often governed by the assumption that the Maxwell-Garnett model works well at low concentrations and the Bruggeman model works better at higher concentrations. Work by Karkkainen et al., using a finite difference method for random mixtures of non-absorbing materials, indicated that the Bruggeman approximation works best when there is some clustering of the inclusions and the Maxwell Garnett model works best when there is no clustering [37].

The Bruggeman solution has been shown to be unphysical in certain circumstances [38]. In particular when the real components of the permittivity have different signs and when the absolute values of the real components are much larger than those of the imaginary components. Unfortunately, it may be that these conditions will apply to modelling infrared absorption. As a result, only a few of the examples discussed below will include results using the Bruggeman mixing rule; the majority will use the Maxwell-Garnett mixing rule.

Averaged-Permittivity mixing rule

It is useful to be able to compare the effective medium theories with the absorption predicted using no shape information, that is using only the TO frequencies.

(47)\[\tensorbs{\varepsilon}_{TO} = f\left\langle \tensorbs{\varepsilon}_i \right\rangle + \left( 1 - f \right)\varepsilon_e\]

Equation (47) defines an isotropic permittivity which can be used to calculate such an absorption coefficient. The angle brackets indicate an average of orientation. This mixing rule provides a useful comparison between the absorption calculated without any shape effects and that calculated including shape effects using the other mixing rules presented above. At low concentrations, the peak positions of the AP mixing rule will be at the TO frequencies.

Particle size effects

Meier and Wokaun [39] outlined an approach to treating large (metal) spherical particles, where particle size is incorporated in terms up to 3^rd order in the wave vector k. Using Equations (26) and (30) we can write;

(48)\[\fieldbf{P} = \left( \tensorbs{\varepsilon}_i - \varepsilon_e \tensorbf{1} \right)\left( \fieldbf{E}_e + \fieldbf{E}_d \right)\]

(49)\[\begin{split}\begin{aligned} \fieldbf{E}_d &= - \frac{\left( 1 - x^2 - i\frac{2}{3}x^3 \right)}{\varepsilon_e}{\tensorbf{L}} \fieldbf{P} \\ x &= ak = \frac{2\pi a}{\lambda} \end{aligned}\end{split}\]

Here \(a\) is the radius of the spherical particle and \(x\) is the dimensionless ‘size’ of the particle with respect to the wavelength of the incident light. The first term relating the depolarisation field to the polarisability is the same as that used above in Equation (29). The second term is a dynamic depolarisation term and the third term is a radiation damping correction to the electrostatic solution [39].

A slightly different, but related, approach is presented by Shivola [33]. Starting from Equation (33), a size-dependent term \(G(x)\) is introduced as indicated by the work of Peltoniemi [40];

(50)\[\begin{split}\begin{aligned} G\left( x \right) &= G_1\left( x \right) + G_2\left( x \right) \\ G_{1}\left( x \right) &= \ \frac{2}{3}\left\lbrack \left( 1 + ix\right)e^{- ix} - 1 \right\rbrack \\ G_{2}\left( x \right) &= \ \left ( 1 + ix - \frac{7}{15}x^{2} - i\frac{2}{15}x^{3} \right)e^{- ix}-1 \end{aligned}\end{split}\]

The modified equation for the relationship between the internal field and the external field based on equation (33) becomes;

(51)\[\fieldbf{E}_i = \fieldbf{E}_e - \frac{\left( 1 - G\left( x \right) \right)}{\varepsilon_e}\tensorbf{L} \left( \tensorbs{\varepsilon}_i - \varepsilon_e\tensorbf{ 1} \right)\fieldbf{E}_i\]

This leads to a modified equation for the polarisability of spherical particles;

(52)\[\tensorbs{\alpha}\left( x \right) = V\varepsilon_e\left( \tensorbs{\varepsilon}_i - \varepsilon_e\tensorbf{1} \right)\left( \varepsilon_e\mathbf{1} + \left( 1 - G\left( x\right) \right)\tensorbf{L} \left( \tensorbs{\varepsilon}_i - \varepsilon_{e}\mathbf{1} \right) \right)^{- 1}\]

Using the modified, sized dependent polarisability all the Bruggeman and Maxwell mixing schemes can be implemented in a way that incorporates size effects. Generally, the onset of changes in the calculated absorption is an indication that size effects are important and should be treated properly.

Mie Theory of Light Scattering

For spherical particles that are comparable in size to the wavelength of light, the theory developed by Mie and described fully by van de Hulst [23] can be used. Unfortunately, this theory is only applicable to spherical, isotropic, particles where the separation between the particles is large compared with the wavelength of light.

PDielec applies Mie theory using the Python library PyMieScatt [8]. To treat anisotropic systems, PDielec first calculates the Maxwell-Garnett permittivity for the anisotropic sphere embedded with random orientation in the supporting matrix. The program then calculates the required permittivity of an isotropic sphere which would give the same Maxwell-Garnett permittivity as the anisotropic system. The Mie scattering is then performed on this isotropic sphere.

An alternative method that diagonalises the real permittivity and then transforms the complex permittivity into the new basis was found to cause anomalous absorptions at the TO frequencies and so the method is no longer used. However, the theory is given here in case it ever has to be revived.

(53)\[\begin{split}\begin{aligned} & {\tensorbf{U}^{\text{T}}} \tensorbs{\varepsilon}^{real} \tensorbf{U} = \tensorbs{\varepsilon}^{diagonal} \\ & {\tensorbs{\varepsilon}}^{new} = \tensorbf{U}^\text{T} \tensorbs{\varepsilon}^{full}\tensorbf{U} \end{aligned}\end{split}\]

The core of Mie theory is the calculation of efficiency factors S1(θ) and S2(θ) (page 125 in reference [23]. These factors depend upon the size of the particle relative to the wavelength of the light (\(x\)) and the refractive index of infrared active material (\(m_1\)) relative to that of the supporting medium (\(m_2\)). Both refractive indices are scalars, \(m_1\) is complex as light is absorbed by the material and \(m_2\) is assumed to be real.

For each of the three diagonal elements of \(\tensorbs{\varepsilon}_{new}\) the relative refractive index, \(m\) is calculated;

(54)\[\begin{split}\begin{aligned} m_1 &= \tensorbs{\varepsilon}_{w,w}^{new} \\ m &= \frac{m_{1}}{m_{2}} \\ x &= \frac{2\pi a}{\lambda} \\ \lambda &= \frac{\lambda_{\text{ vac}}}{m_{2}} \end{aligned}\end{split}\]

The size parameter depends upon the refractive index of the supporting medium because the wavelength of light is smaller than its vacuum value in a medium with a refractive index greater than one. The subscript \(w\) indicates one of the three directions x, y or z in the basis which diagonalises the real permittivity tensor.

Using the routines available in PyMieScatt [8] the efficiency factor \(S(0)\), where \(S(0) = S1(0) = S2(0)\), is calculated from \(m\) and \(x\). The extinction (this includes absorption and scattering) along with retardation are reflected in the overall complex refractive index, page 129 in reference [23];

(55)\[\begin{split}\begin{aligned} \widetilde{m}_w &= m_2\left( 1 + iS\left( 0 \right)2\pi Nk^{- 3} \right) \\ N &= \frac{f}{V_{sphere}} \\ V_{sphere} &= \frac{4}{3} \pi a^{3} \\ k &= \frac{2\pi}{\lambda} \end{aligned}\end{split}\]

Here \(N\) is the number density of particles and \(k\) is the wave-vector of the light in the surrounding medium. The subscript \(w\) indicates that there are three values for the effective refractive index, depending on which diagonal value of the permittivity tensor is taken. There is a change in sign from the expressions given by Hulst, this is to make the imaginary component compatible with the effective permittivities. Once the effective refractive index is known in each direction the effective permittivity can be calculated using the one-third rule [41];

(56)\[\begin{split}\begin{aligned} \widetilde{m}_{eff} &= \frac{1}{3}\sum_w{\widetilde{m}}_w \\ \varepsilon_{eff} &= \widetilde{m}_{eff}^2 \end{aligned}\end{split}\]

This approach to taking anisotropy into account when the embedded particles are anisotropic but randomly oriented is approximate but has been shown to have a reasonably wide range of application [41].

Calculation of ATR Extinction Spectra

Attenuated total reflectance spectroscopy (ATR) uses a medium of high refractive index contrast to modify the reflectance and transmission of the infrared radiation. Balan has reported the calculation of the ATR spectrum of apatite [42] and chrysotile [43]. A similar approach is adopted by the PDielec package. Denoting the angle of incidence by \(\theta\) and the refractive index of the high index, non-absorbing material by \(n_0\) and that of the absorbing material by \(n_s\), the Fresnel amplitudes of reflection of the parallel polarised component (\(r_p\)) and perpendicular polarised component (\(r_s\)) is given as follows; [44]

(57)\[\begin{split}\begin{aligned} r_s &= \frac{ -n_0 \cos{\theta} + \sqrt{n_s^2 - n_0^2 \sin^2{\theta}}}{n_0 \cos{\theta} + \sqrt{n_s^2 - n_0^2 \sin^2{\theta}}} \\ r_p &= \frac{ \frac{n_s^2}{n_0} \cos{\theta} + \sqrt{n_s^2 - n_0^2 \sin^2{\theta}}}{\frac{n_s^2}{n_0} \cos{\theta} + \sqrt{n_s^2 - n_0^2 \sin^2{\theta}}} \end{aligned}\end{split}\]

Using \(f_p\) and \(f_s\) as the fractions of P and S polarised incident light, the extinction is expressed in absorbance units as;

(58)\[\begin{split}\begin{aligned} e_s &= -log_{10}( r_s^*r_s ) \\ e_p &= -log_{10}( r_p^*r_p ) \\ e_{total} &= -log_{10}( f_p r_p^*r_p + f_s r_s^*r_s ) \\ \end{aligned}\end{split}\]

The approach adopted by Balan and in the PDielec package is to use the calculated effective permittivity to determine the (complex) refractive index of the material being studied. ATR measurements are taken with compressed powders of the material, so the appropriate effective medium can be determined from a high volume fraction of the crystal in air.

Scattering Effects due to Air Inclusions

A common method for taking a terahertz spectrum is to dilute the active dielectric material in a transparent material such as polyethylene or PTFE spheres and then compress the composite material. This leads to the presence of large air inclusions, or bubbles in the sample material. The size of these inclusions can vary in size from 10 to 100 microns and is therefore capable of scattering terahertz radiation. To account for this the effective dielectric constant of the supporting matrix can include these scattering effects using Waterman-Truell scattering theory [45].

(59)\[K = \sqrt{ k^2 [ 1 + \frac{ 2 \pi N f(0)}{k^2} ]^2 - [ \frac { 2 \pi N f(\pi)}{k^2} ]^2 }\]

Where \(N\) is the number density of voids, \(K\) is the effective wavenumber of the wave in the material, \(k\) is the wavenumber of the wave in the support without air inclusions and \(f(0)\) and \(f(\pi)\) are the forward and backward amplitude scattering factors.

(60)\[f(\theta) = \frac{i}{k} S(\theta)\]

The scattering factors are determined using Mie theory, where it is assumed that there are spherical voids in the matrix, which are large enough to scatter light.

(61)\[\begin{split}\begin{aligned} m_2 &= \frac{K \lambda} {2 \pi} \\ \epsilon_e &= m_2^2 \\ \end{aligned}\end{split}\]

The effective permittivity of the matrix with voids is then used in the calculation of the effective medium theory for the dielectric material. The calculation proceeds as it did in the case with no scattering, only now a complex permittivity is used to describe the supporting matrix.

Removing Baseline Artefacts from Experiment

In the process of reproducing experimental spectra using calculated information, it is often necessary to remove any baseline artefacts from the experimental data. A relatively straightforward approach is adopted by PDielec using a Hodrick-Prescott filter. The approach is also called asymmetric least squares smoothing.

If there is a data sequence \(y\), sampled at regular intervals we want to determine the baseline sequence \(y\) such that x and y are as close as possible to each other with the requirement that x is as smooth as required.

It can be shown [46] that this can be achieved by an iterative solution for \(x\) of ;

(62)\[( \tensorbf{W} + \lambda \tensorbf{D}^T \tensorbf{D}) x = \tensorbf{W} y\]

where \(\tensorbf{W}\) is a diagonal matrix of weights modified during the iterations to make values of \(x\) lie below those of \(y\) and \(\tensorbf{D}\) is a 3rd order differencing matrix. Since the differencing matrix is very sparse, sparse matrix routines are used to solve the equations.

In PDGui the value of \(\lambda\) is supplied as \(\log \lambda\).