Dielectric function decomposition by dipole interaction distribution: application to triclinic K2Cr2O7

Here we present a general approach for the description for the frequency dependent dielectric tensor coefficients for optically anisotropic materials. Based on symmetry arguments we show that the components of the dielectric tensor are in general not independent of each other. For each excitation there exists an eigensystem, where its contribution to the dielectric tensor can be described by a diagonal susceptibility tensor. From the orientation of the eigensystem and the relative magnitude of the tensor elements, the dipole interaction distribution in real space can be deduced. In the limiting cases, the oriented dipole approach as well as the tensor of isotropic and uniaxial materials are obtained. The application of this model is demonstrated exemplarily on triclinic K2Cr2O7 and the orientation and directional distribution of the corresponding dipole moments in real space are determined.


Introduction
The optical response of a material is determined by its dielectric function (DF) and the knowledge of that is important for the design, fabrication and understanding of optoelectronic devices. Also the DF gives a deeper understanding of the physical properties of the material, like transition energies and phonon energies [1]. In order to obtain insight into the nature of the polarizability, the DF is usually decomposed into its excitations and the contribution of each excitation is described by using a model dielectric function [2]. For crystals with an orthorhombic symmetry or higher, this decomposition of the DF can be easily done for each tensor component independently from each other, since for these crystals the DF can be diagonalized for all energies simultaneously. The dielectric axes, i.e. the axes given in the framework of the diagonalized tensor, coincide with those given by crystal axes. As a consequence, all excitations are aligned along directions defined by the crystal axes.
In contrast to that, for monoclinic and triclinic materials the DF cannot be diagonalized in general. Only in the transparent spectral region, the DF can be diagonalized but not for all energies simultaneously and the orientation of the eigenvectors (normal modes) depends on the wavelength. Furthermore, despite both the real as well as the imaginary parts of the DF are each represented by symmetric tensors, they have different eigensystems, and the excitations are no longer aligned along directions defined by the crystallographic axes. We note, that the DF of a monoclinic form can be block-diagonalized, i.e. into a 2 × 2 matrix and a scalar along the crystallographic symmetry-axis. The excitations along and perpendicular to the symmetry axis can be described independently of each other. In the last years, materials with a monoclinic and triclinic crystal structure went into focus of research, like Ga 2 O 3 [3][4][5][6], CdWO 4 [7,8], Y 2 SiO 5 [9], α-PTCDA [10], pentacene [11], BiFeO 3 [12], K 2 Cr 2 O 7 [13], CuSO 4 · 5H 2 O [14]. Due to the above mentioned problem in the decomposition, most of these works are limited to the determination of the line shape of the DF tensor elements only. Furthermore, the tensor elements are often treated independently of each other which can lead, from a technical point of view, to a large correlation between the tensor elements and to unphysical results.
In order to explore the nature of the polarizability in such materials in the infrared spectral range, Born proposed a model which takes into account the orientation of the dipole moment of the phonons with respect to the laboratory system [15,16]. This oriented dipole approach assumes that each excitation which contributes to the DF has a dipole moment of a well-defined direction which does not necessarily coincide with the crystallographic axes. The magnitude and dispersion of the DF tensor elements can be then deduced by a rotation of the susceptibility of the framework of the excitation into the laboratory system [5] or by a projection of the dipole orientation [3]. Both approaches lead to the same result. The decomposition of the DF into its components and their orientation can for instance be used in order to explain the different optical properties of the polymorphs, e.g. the polymorphs of calcium carbonate, i.e. calcite and aragonite [17]. By means of the oriented dipole approach, the phonon properties of optically anisotropic materials were already investigated [3,5,8,9,13,14,18,19]. Recently, the generalization of this approach for electronic excitations, such as excitons, which can be used up to the UV spectral range was demonstrated on β-Ga 2 O 3 by using ellipsometry [5,6]. However, such model can only be applied for transitions which have a well-defined orientation of the dipole moment, as it is typically the case for phonons or excitonic contributions. This requirement is not fulfilled in general for all materials especially if excitations cannot be resolved separately due to the energetical overlap. Furthermore, as known from microscopic band theory, different k-points and different bands can contribute in general to a definite transition energy in real crystals which leads to an angular smearing of the dipole orientation of the corresponding excitations.
In this work we present a general approach for the description of the DF by using the dipole interaction distribution (DID) of each excitation. Based on symmetry arguments and without making any assumption on the nature of the polarizability we show that the components of the dielectric tensor cannot be independent of each other. For each excitation there exists an eigensystem, where its contribution to the DF can be described by a diagonal susceptibility tensor. The relative magnitude of the diagonal elements are a measure of the corresponding DID in real space. In the limiting cases, the oriented dipole approach as well as the DF for isotropic and uniaxial materials is obtained. The application of this model is exemplarily demonstrated for triclinic K 2 Cr 2 O 7 and the DID of the involved excitations in real space are deduced.

Theory
The DF is in general complex-valued and the real and imaginary part of each tensor component is connected with each other by the Kramers-Kronig relation [20]. This can be easily seen by the fact that the macroscopic polarisation (P) induced by the electromagnetic wave is given by with E, and 0 being the electric field, the DF, which is in general a tensor, and the vacuum permittivity, respectively. Equation (1) has to be fulfilled for all directions of the electric field. If we choose the electric field in such a way, that it is aligned along an axis of the coordinate system, then each component of the polarisation vector depends only on one tensor component of the DF, i.e.
for k = j and δ ij being the Kronecker delta. From causality it follows that the induced polarisation cannot appear before the stimulus, i.e. P(t) = 0 and thus for each component ij − δ ij = 0 holds for t < τ. If we further transform the polarization from the time domain into frequency domain, we obtain and The index 1 and 2 indicates the real and imaginary part of the tensor component of the DF, respectively. From the application of Cauchy's theorem, which connects real and imaginary part of the Fourier transformation of analytical functions, and equation (2) we obtain the validity of the Kramers-Kronig transformation for each component individually.
In the following we will deduce the relation between the components of the DF. In order to explore the nature of the polarisability and to extract the properties of the involved transitions, the susceptibility and thus the contribution of each transition to the DF is described by a model function χ i (see e.g. references [2,21,22]), which represents the susceptibility of the excitation. The exact line shape of χ i depends on the nature of the involved transitions. Without making an assumption of the relation between the individual tensor components of the DF, we define as the experimentally determined tensor of the DF referred to the coordinate system defined by x, y and z and assume that each component can be described by a sum of model functions, i.e.
where δ ij describes the vacuum contribution to the DF and N ij the number of contributions to ij . Note, at this point we do not explicitly imply that each model function has to appear in each component. Summing up transitions which have a similar line-shape and differ only in their magnitude of the individual components of the susceptibility tensor, equation (5) can be written in matrix form as withχ k being a normalized complex-valued scalar function and A k a matrix, the so-called amplitude matrix, describing the line-shape and the magnitude of the susceptibility tensor of the nth transition, respectively. The symbol represents the unit matrix and N the number of contributions to the susceptibility which differ in their line-shape. Note, the amplitude matrix A k does not depend on the energy. Since in monoclinic and triclinic crystals, due to the non-orthogonal angle between two or all crystallographic axes, there does not exist a well defined Cartesian coordinate system, so a coordinate system defined by x, y and z is chosen arbitrarily. In case of monoclinic crystals, only one axis of the Cartesian coordinate system can be aligned along the symmetry axis. Therefore, the determined DF could be also referred to a coordinate system given by x , y and z . The DF in this system is denoted by . Similar to , in the (x , y , z )-coordinate system the components of can also be described by model functions, resulting in a similar expression as given in equation (5) for . Since the dielectric response must be independent of the choice of coordinate system and thus fulfills transformation properties, the tensors for the two coordinate system are connected by each other by where R is the rotation matrix which transforms a vector from one into the other coordinate system. If we now compare the tensor components of and with each other and take into account that they are related with each other by equation (7) we can deduce: (a) The components of or rather cannot be described independently of each other. This can be seen by the fact, that if we assume that the components of can be described independently of each other, then the components of consist of a linear combination of all components of , i.e. ij = nkm a ijkm χ n,km and the coefficients of a ijkm are then determined by the rotation of the two coordinate systems against each other. Thus, the components of ij cannot be independent of each other. As we assumed that the components of are independent of each other, this would imply that there is a well-defined Cartesian coordinate system defined by (x, y, z) where all tensor components of the DF are independent of each other. However, this is in contradiction to the fact that in triclinic and monoclinic crystals, due to the non-orthogonal angles between the crystallographic axes, no well-defined Cartesian coordinate system exists and thus the components of the DF are not independent of each other in any Cartesian coordinate system. (b) The amplitude matrix A n is symmetric, since the DF is a symmetric tensor (in the absence of magnetic fields) [23], and thus the contributions of each excitation to the DF as well. (c) In general, each transition which is contributing in a given framework to an off-diagonal element of the DF contributes also to at least two diagonal elements in the same framework, i.e. at least two diagonal elements of A n are non zero. The amplitudes of the contribution can differ for each component. (d) As mentioned above, the amplitude matrix A n is a real-valued symmetric matrix. Since a real-valued symmetric matrix can be diagonalized, for each transition there exists a coordinate system in which the amplitude matrix A n is a diagonal matrix. Note, only one A n,ii has to be non-zero. Thus, for each transition exist only 3 independent amplitude parameters describing the magnitude of the contribution to the DF. This is also reflected by the fact that a 3 × 3 matrix has three invariants with respect to a similarity transformation, namely its trace, determinant and the sum of the minors. As mentioned in equation (6), the contribution of each excitation to the susceptibility tensor is given by χ n = A nχn with χ n being a scalar function, the susceptibility tensor χ n is diagonalizable for all energies as well. If we choose the xx-element as the non-zero component, the susceptibility tensor of the nth excitation in its eigensystem can be written as with A n being a real-valued scalar, describing the magnitude of the contribution of the nth excitation to the χ D n,xx in its corresponding eigensystem. The abbreviations f n,y and f n,z are real-valued and describe the different weights of the components for the nth transition. Since the contribution to the susceptibility tensor of each transition is described by a Kramers-Kronig consistent function, the ratio between the imaginary and real part for all components of the susceptibility tensor for a given transition is the same and thus f n,i is a real valued number. If we further choose the coordinate system in such way that the xx-component has the highest amplitude and take into account that the imaginary part of χ D ii has to be positive, since otherwise gain would be present, which is impossible without the presence of additional excitation, not considered here. Thus it holds 0 f i 1.
Taking these considerations into account, the entire DF in the x-y-z-system is given by with χ D i representing the susceptibility of the nth excitation in the eigensystem of the excitation. The rotation matrix R n transforms the susceptibility tensor of the nth excitation from its eigensystem into the laboratory system and depends on the Euler angles φ n , θ n and ψ n . Note, although the susceptibility tensor for each excitation is diagonalizable for all energies simultaneously, the entire DF is not diagonalizable in real and imaginary part simultaneously.
In the following we show that the coefficients f i defined by equation (8) can be interpreted as a measure for the DID in real space of the corresponding excitation. We consider, that each transition contributing to the DF consists of an ensemble of microscopic excitations with similar transition energy, amplitude and broadening but different orientation. If these microscopic excitations are discretely distributed in real space, e.g. if they are aligned along crystallographic axis, the coefficients f i are given by with r i,j being the jth component of the direction vector of the ith microscopic excitation.
In the case that the microscopic orientations are continuously distributed in real space, e.g. the excitations cannot be spectrally resolved due to an spectral overlap the microscopic excitations or different k-points and bands contribute to the excitation, then the sum in equation (10) has to be replaced by an integral. In the excitation framework the orientation of the microscopic excitations can be described by where α and β represents the polar and azimuthal angle. If g = g(α, β) represents the distribution of the orientation of the microscopic excitations, the weighting factors are given by with α m,β being the maximum polar angle for a given azimuthal angle β. Since the exact distribution of the microscopic excitation is not known and to reduce the number of required parameters as well as the correlation between these parameters, we assume for simplicity a uniform distribution around the x -axis, i.e. g(α, β) = 1 and that the microscopic dipole orientations are within a cone which describes in the y -z -plane a superellipse with its major axis parallel to the y -and z -axis. The relation between α m,θ and β is then The angles α y and α z represent the opening angles of the cone described by the microscopic excitations in the x -y and x -z -plane, respectively (figure 1). For k = 2 the cone would describe an ellipse in the y -z -plane. The weighting factors f i of the susceptibility tensor in the eigensystem are then given by (2 + cos(α m,β ))sin 4 α m,β 2 cos 2 (β)dβ (14a) If α y = α z = π/2, i.e. the microscopic transitions are equally distributed in all directions, we get f y = f z = 1. In this case the susceptibility tensor reduces to a scalar as it is the case for isotropic materials. For uniaxial materials the transitions are equally distributed in one plane only, e.g. the x-y-plane, but have no contribution to the other direction, i.e. f y = 1 and f z = 0. For an excitation which is well defined along the x -direction α y = α z = 0 holds and thus f y = f z = 0 which is then equivalent to the dipole ansatz approach [3,5,[14][15][16]. Hence, the quantities f y and f z are a measure for the DID of the transition in the x-y-and x-z plane, respectively, and P = (f y + f z )/2 is a measure for the distribution of the dipole interaction in the entire space. For P = 1 the dipole moments are equally distributed in the entire space as expected for isotropic materials whereas P = 0 represents the case where the dipole moment is oriented along a single well-defined direction.
By using R(φ n , θ n , ψ n ) = R z (φ n )R x (θ n )R z (ψ n ) [24,25], where the index denotes the mobile rotation axis, the symmetry axis of the dipole moment distribution of the nth transition in the laboratory system is given by cos(φ n ) cos(ψ n ) − cos(θ n ) sin(φ n ) sin(ψ n ) cos(ψ n ) cos(θ n ) sin(φ n ) + cos(φ n ) sin(ψ n ) sin(φ n ) sin(θ n ) For the case f y = f z = 0 this vector r represents the orientation of the dipole in the oriented dipole approach. It is interesting to note, that although the orientation of a dipole pointing in an arbitrary direction can be represented by two angles, for tensors describing the rotation of the respective tensorial property in space three angles are needed. Here, the third angle ψ is required for the additional revolution around the dipole direction. Thus, by using the orientated dipole approach, three Euler angles have to be considered for triclinic materials. Note, for monoclinic crystals due to symmetry considerations two Euler angles are sufficient. In this case the DF consists of a preferred direction and the excitations are oriented either parallel or perpendicular to this direction.

Experimental
We exemplarily demonstrate the application of the model presented in section 2, on K 2 Cr 2 O 7 bulk single crystals, which were grown at room temperature by evaporation of a saturated K 2 Cr 2 O 7 solution. The infrared properties of this crystal and details of the crystal preparation can be found in [13]. K 2 Cr 2 O 7 has a triclinic crystal structure and belongs to the C 1 i space group. The lattice constants and the angles between the crystallographic axes are [26] [2,28]. The MM can be represented by a matrix of 2 × 2 block matrices. The diagonal block matrices are mainly connected to the reflection coefficients for the s-and p-polarized light, whereas the off-diagonal ones are a measure for the above mentioned mode conversion. For the analysis, the MM is normalized to M 11 since this element is related to the total reflected intensity.
By using a dual rotating compensator ellipsometer all 15 elements of the normalized MMs in the spectral range from 0.747 eV up to 6.30 eV were determined for the above mentioned growth faces. For each face we used angles of incidence of 50 • , 60 • and 70 • and rotated the sample around the surface normal of each face from 0 • to 315 • with an increment of 45 • . Additionally, we took measurements for in-plane rotation angles of 30 • and 60 • . The orientation of the sample with respect to the laboratory system was described by using Euler angles (φ, θ, ψ) in the zxz-notation [24]. In the following we choose the sample coordinate system in such way, that the x-y plane coincides with the (100) plane of K 2 Cr 2 O 7 and the x-axis is within the (001) plane. A schematic of the sample and the chosen coordinate system is given in figure 2(b).
In order to extract the DF from the recorded MM, all sample orientations were analyzed simultaneously by using a layer stack model consisting of a semi infinite substrate, representing the bulk properties of K 2 Cr 2 O 7 , and a surface layer. The DF of the surface layer was described by using a Bruggeman effective medium approach [29] where the DF of the substrate layer, i.e. K 2 Cr 2 O 7 , is mixed with that of void with the ratio of 1 : 1. For simplicity, only the diagonal elements of the DF of K 2 Cr 2 O 7 were used in order to calculate the optical properties of the surface layer. The thickness of the surface layer was estimated to be about 50nm.
For the determination of the DF we analyzed at first the transparent spectral region and described each component of the DF independently from the other one by means of a Cauchy approximation, i.e. √ 1,ij = A + B/λ 2 + C/λ 4 and 2,ij = 0, where the Cauchy parameter A describes the low energy limit and B and C the curvature of the refractive index. In doing so, a first guess of the DF, the thickness of the surface layer and the Euler angles were estimated. In a second step, in the entire spectral range for each energy the DF was determined, while keeping all other wavelength independent parameters constant, like the thickness of the surface layer. In the final step, all parameters were analyzed simultaneously and the DF of K 2 Cr 2 O 7 was obtained by using a wavelength-by-wavelength analysis, i.e. the components of the DF were obtained for each wavelength independently from each other. The only global parameters were the constant Euler angles as well as the surface thicknesses for each face. In doing so, no assumption regarding the line shape and the involved electronic transitions had to be made. Note, a Kramers-Kronig relation (KKR) between the real and imaginary part of the tensor elements has not been enforced in this analysis. This is not required since intensity and phase information from the ellipsometric measurements allow to retrieve the real and the imaginary parts independently. For the minimization of the mean square error between the experimental and calculated Müller matrix, we use the differential evolution algorithm implemented in BlackBoxOptim package for Julia with a subsequent minimization of the best population by using a BFGS algorithm [30].
In order to get a deeper understanding of the line shape the DF deduced by wavelength-by-wavelength analysis was parameterized by model dielectric functions. Thereby, the onset of the absorption was described by using a critical point function for parabolic band-band transitions (CPM0) developed by Adachi [31] with and where A, E 0 and Γ being the amplitude, central energy and broadening of the transition. It was shown, that the optical properties of K 2 Cr 2 O 7 are mainly determined by the presence of the dichromate ions (Cr 2 O 7 ) which cause a charge transfer from the ligands to the metal ions [32,33]. The line shape of these transitions were described by Gaussian oscillators, i.e. the susceptibility was described by The integral of the real part of the susceptibility of χ Gauss represents the Kramers-Kronig relation and results in a Dawson function [34].

Results
The experimentally determined Müller matrix is exemplarily shown for selected orientations and for an angle of incidence of 60 • in figure 3. In the whole investigated spectral range of the MM spectra broad peaks are observable which indicate the presence of absorption in the investigated spectral range. Furthermore, for all orientations we observe non-vanishing block-off diagonal matrix elements. This demonstrates the presence of mode conversion and thus signifies the optical anisotropic character of the material. The calculated Müller matrix by using the wavelength-by-wavelength approach as described in section 3 is shown by red solid lines which yield a good agreement with the experimental spectra. From the best match between the calculated and experimental MM, the DF was deduced and is shown in figure 4(a) in the chosen laboratory coordinate system as indicated in figure 2(b). It should be mentioned, that although not directly included in the wavelength-by-wavelength analysis, a numerical analysis by using Kramers-Kronig consistent b-spline functions [35] yields that all obtained tensor elements fulfill the KKR ( figure 4(a)). On the one hand this corroborates the validity of the expected KKR for each tensor element as stated above and on the other hand demonstrates the quality of the fit. Note, the analysis of the deduced components of the DF by means of the b-spline functions was only used in order to verify the validity of the KKR. The b-spline functions were not considered in the further analysis.  As expected for a triclinic material, all six tensor components of the DF are non-zero and there exists no set of Euler angles where the complex non-diagonal tensor components vanish for all energies. This is in line to the results obtained in the infrared spectral range [13]. This is in contrast to monoclinic and orthorhombic structures where at least one set of Euler angles can be found so that two or all three off-diagonal elements vanish, respectively. Only for an energy of E ≈ 2 eV there exists a coordinate system where two off-diagonal elements vanish. Thus, at this energy the shape of the DF is similar to that of a monoclinic system. Interestingly, the diagonal elements xx and yy are quite similar. This indicates, that in the given laboratory system the anisotropy in the x-y-plane is only weakly pronounced. This is also Table 1. The parameters of the excitations deduced from the line shape analysis. The band-band transition and the Gaussian oscillators are denoted by CPM 0 and G i , respectively. The opening angles α i where calculated by equation (14) assuming an unweighted orientation distribution. reflected by the weak pronounced off-diagonal tensor element xy , which is almost one order of magnitude smaller than xz and yz ( figure 4(a)).
In order to obtain a deeper insight into the distribution of the electronic transitions, the DF deduced by the wavelength-by-wavelength analysis was described by means of model dielectric functions as described in section 2. A good description can be obtained by using a band-band transition and Gaussian oscillators. The entire model dielectric function is then given by equation (6). For the determination of the corresponding parameters we analyzed all tensor elements simultaneously and used for the minimization of the objective function a differential evolution algorithm. By using this approach a good agreement between the tensor components of the DF obtained by the wavelength-by-wavelength analysis and model dielectric functions is obtained ( figure 4(a)). This line-shape analysis yields that the DF is mainly determined by transitions described by Gaussian oscillators which we attribute to the transitions caused by the charge transfer from the ligands to the metal ions of the dichromate ions. Due to the spectral overlap of the band-band transitions with the higher excitations as well as of the broadening of band-band transition it was not possible to determine the band gap energy for a well-defined polarization. Thus an averaged band-gap energy of about E ≈ 2.45 eV was determined. This is also reflected by the amplitude matrix of the susceptibility tensor in the eigensystem of the excitation which holds similar values for all diagonal elements (f y ≈ 0.93 and f z ≈ 0.8) indicating that the polarization dependence of the band gap energy is quite small.
As mentioned above, the line shape of the diagonal elements of the DF are similar and instead of a well-defined dipole direction for the higher excitations only a direction of an ensemble of excitations can be determined. This is also reflected by the diagonal elements of the corresponding weighting matrix of the susceptibility in the eigensystem which are similar and indicate a broad directional distribution of the involved transitions. We attribute this finding to the relatively large broadening of these transitions and thus the large spectral overlap with the energetically neighbored excitations. Only for the oscillators G 3 , G 6 , G 9 and G 10 a significant inhomogeneous directional distribution was obtained, i.e. the directional distribution in the x-z-plane differs strongly from that one in the x-y-plane. Especially the oscillator G 3 and G 10 are distributed in a plane only. The determined properties of the susceptibility tensors in the eigensystem and its orientations with respect to the laboratory system for each excitation are summarized in table 1.
A similar line shape analysis can be obtained by using the dipole approach as it was demonstrated for monoclinic materials [3,5,6]. However, this would enhance the numbers of required oscillators and thus of the involved parameters. It also would lead to a strong correlation between the parameters and making the interpretation difficult. Note, a description of each component of the DF independently from the other one, although as explained above not physically correct, will not lead to a reduction of the number of free parameters, as e.g. for the description of xx and xy four and six excitations are required, respectively. By using equation (14) and the obtained weights f i of the susceptibility tensor, the DID for each oscillator in the chosen laboratory coordinate system was calculated. The deduced orientation of the corresponding dipole moment and its directional distribution are depicted in figure 4(b). This demonstrates quite nicely the aforementioned distribution of the DID of the excitations and shows that, as expected for triclinic materials, the excitations are not necessarily aligned to the crystallographic axes.
For completeness, the parameters of the high energy transitions, with a transition energy above the measured spectral range are given in table 2. A determination of the DID for these transitions is due to the limited spectral range not possible and thus the contributions for each tensor component was analyzed independently from each other.

Summary
We demonstrated that the dielectric function can be decomposed into its excitations by taking into account their dipole orientation distribution. From symmetry arguments it follows that all components of the dielectric function in general are not independent of each other. The contribution of each excitation to the dielectric function is given by a complex-valued diagonal susceptibility tensor and three Euler angles. The latter describes the orientation of the eigensystem of the excitation with respect to the laboratory coordinate system. The relative magnitude of the corresponding susceptibility tensor elements represent the dipole orientation distribution in real space. In the special case of a well-defined orientation of the dipole or if the dipole is equally distributed in one plane or in all directions, the oriented dipole approach and the tensor for uniaxial and isotropic materials are obtained, respectively. The application of our generalized model for the dielectric response was successfully demonstrated on triclinic K 2 Cr 2 O 7 and the experimentally determined dielectric function was very well reproduced. Based on the relative amplitudes of the susceptibility tensor of each excitation, its orientation and its dipole orientation distribution in real space was obtained.