Non-destructive determination of phase, size, and strain of individual grains in polycrystalline photovoltaic materials

We demonstrate a non-destructive approach to provide structural properties on the grain level for the absorber layer of kesterite solar cells. Kesterite solar cells are notoriously difficult to characterize structurally due to the co-existence of several phases with very similar lattice parameters. Specifically, we present a comprehensive study of 597 grains in the absorber layer of a 1.64% efficient Cu2ZnSnS4 (CZTS) thin-film solar cell, from which 15 grains correspond to the secondary phase ZnS. By means of three dimensional X-ray diffraction (3DXRD), we obtained statistics for the phase, size, orientation, and strain tensors of the grains, as well as their twin relations. We observe an average tensile stress in the plane of the film of ~ 70 MPa and a compressive stress along the normal to the film of ~ 145 MPa. At the grain level, we derive a 3D stress tensor that deviates from the biaxial model usually assumed for thin films. 41% of the grains are twins. We calculate the frequency of the six types of $\Sigma$3 boundaries, revealing that 180{\deg} rotations along axis<221>is the most frequent. This technique can be applied to polycrystalline thin film solar cells in general, where strain can influence the bandgap of the absorber layer material, and twin boundaries play a role in the charge transport mechanisms.


Introduction
Photovoltaic thin-film technology is increasingly targeting alternative materials to meet the triple challenge of sustainability, low energy payback time, and scalability. Current technologies include polycrystalline CdTe [1] and Cu(In, Ga)Se2 (CIGS) [2], both with power conversion efficiencies that surpass 20%. A relatively new but promising candidate is Cu2ZnSnS4 (CZTS), with an efficiency of 11% [3], and the selenized version, Cu2ZnSn(S, Se)4, where efficiency has reached 12.6% [4]. All of these materials still perform below the Shockley-Queisser limit [5].
The performance of these materials is strongly dependent on their complex microstructures. One limiting factor shared among these semiconductors is the deficient open-circuit voltage (VOC) attributed, among other reasons, to the structural defects in the absorber layer. For example, a small grain size is associated with an increased amount of grain boundaries, which, if poorly passivated, can contribute to carrier recombination [6]- [9]. Secondary phases can cause other deficiencies. For CZTS with a typically Cu-poor and Zn-rich composition, secondary phases with different bandgaps form, such as the high bandgap ZnS, increasing series resistance when situated in the back contact of the solar cell or acting as a barrier to the charge carriers at the p-n junction [10]- [12]. In CIGS absorbers, which usually have a Cu-poor composition, a Cu(In, Ga)3Se5 phase can occur at the surface with a high density of indium or gallium appearing as copper antisite (In, Ga)Cu defects and acting as recombination centers [13]. A lattice mismatch between CIGS and Cu(In, Ga)3Se5 can also create structural defects and an increased density of recombination centers [14]. Moreover, the lattice spacing changes when modifying the material composition while tuning the band gap e.g., with the variation of Ga/In and Se/S ratios.
Furthermore, the structure of the grains and the local strain change inevitably, as the multicomponent materials undergo different treatments from the deposition and the annealing of the absorber layer to the post-treatment methods for the coating of the subsequent layers of the device. The change in lattice parameters, as a result of fabrication stresses, can affect not only the mechanical properties of the film (adhesion to the substrate, elastic modulus, and deformation [15]) but the electronic properties as well. As an example, theoretical calculations demonstrate the reduction of bandgap due to a tensile biaxial in-plane strain. In contrast, a compressive strain increases the bandgap [16], [17].
CIGS and CdTe exhibit high efficiency, but indium and telluride scarcity is a concern for scaling up module production. Moreover, recycling of systems is complicated because of Cd toxicity. In comparison, CZTS has ideal optoelectronic properties and is made of earth-abundant, non-toxic, and low-cost constituent elements. However, to improve the device efficiency, the structural characterization, such as the identification and quantification of secondary phases, and depiction of grain structural properties, such as strain and twinning, need further work.
An additional complication arising with CZTS is that the crystallographic structures of some of the phases involved have nearly identical lattice parameters, which makes it challenging to identify and quantify the phases. For instance, ZnS with a face-centered cubic crystal structure (F-43m), and the kesterite CZTS with a tetragonal body-centered structure (I-42m), are closely related. Doubling the a, b, or c axis of the cubic ZnS structure yields a unit cell corresponding to kesterite with a small tetragonal deformation | (2 ) ⁄ − 1 | < 0.0026 (with lattice constants from [18]).
A progress beyond the "trial-and-error" approach is vital to visualize the microstructure and local strain within the thin films in 3D and preferably also record their evolution during processing under conditions that are deemed "realistic". Such information can guide theoretical understanding and the development of models for quantitative prediction, thereby accelerating the design efforts. Moreover, physical parameters may be determined by comparing 3D models and 3D experimental data (e.g. [19]). However, the techniques currently employed to characterize the structural properties and the local stress have significant limitations:  Electron microscopy (EM) can provide atomic-scale insight [13] but is confined to studies at the surface or films of a few hundred nanometer thickness. Three-dimensional resolved mapping may be accomplished by a combination of Electron Backscattered Diffraction (EBSD) and serial sectioning using either FIB [20] or laser ablation [21]. However, the destructive procedure prohibits studies of dynamics and direct coupling to functionality. Moreover, the angular resolution achieved by EM makes quantitative stress determination difficult and does not allow for a distinction between the phases mentioned above with nearly identical lattice parameters.
 X-ray powder diffraction (XRD) and grazing incidence X-ray diffraction (GIXRD) provide bulk information about phases, orientations, and strain, but only about average properties. Typically these techniques can identify secondary phases at the level of a volume percentage of 1, but the lack of well-separated peaks in the powder diffraction patterns imply that, e.g. quantification of ZnS is not possible [22].
 X-ray nanoprobe and forward scattering ptychography methods relying on the use of a synchrotron can reveal the local elemental composition in 3D [23] but does not provide structural information. Moreover, sample must be quite small (<10 µm), and dynamics studies representing bulk conditions are, therefore, excluded.
 Spectroscopic methods like X-ray Fluorescence Spectroscopy (XRF), X-ray Absorption Near-Edge Structure analysis (XANES) [24] can reveal the elemental composition but does not reveal anything about the microstructure of the film.
In this paper, we propose three-dimensional X-ray diffraction (3DXRD) as a tool for studying the microstructure and local stress in the photovoltaic polycrystalline films. This non-destructive technique combines highly penetrating hard X-rays from a synchrotron source and the application of 'tomographic' approach to the acquisition of diffraction data [25]- [32]. For grains with a known phase and sizes larger than a few micrometers, 3DXRD can generate 3D maps of several thousands of grains, revealing their shape, orientation, and type II stress (as averaged over each grain) and their variation with time [19]. For grains with a size in the 0.1-1 µm rangeas is typical for photovoltaic polycrystalsshape information is not available, but one can still determine the position, size, orientation, and strains of each grain as a function of time, and thereby generate statistics on the dynamics at the grain scale [27], [32]- [35].
However, the application of standard 3DXRD software to thin-film solar cells is hampered by the complication of phase identification. In principle, a standard single crystal crystallographic analysis can be applied to each grain, a method known as multigrain crystallography [29], [36].
Here we present an approach where a priori information about the photovoltaic materials is used to facilitate the generation of comprehensive statistics of phase, grain size, strain, and twinning relations by standard 3DXRD software. We discuss the importance of such data for R&D in photovoltaics and outline how this work can be generalized to the generation of 3D in situ movies of the microstructure.
The method will be presented with reference to and demonstrated on a specific example: a CZTS (kesterite) solar cell device. We examine the crystallographic properties of this semiconductor on the grain level and the mechanical deformation in the film that the experimental data reveal. Moreover, we present approaches to get around the crystallographic challenges that this absorber layer imposes in order to identify and quantify secondary phases, stress values, and twin boundaries in the material. In our view, other chalcogenide thin-film systems such as CIGS and CdTe could also benefit from this type of 3DXRD analysis gathering statistical information about the absorber film microstructure buried in the multilayer device structure.
In the first part of this study, we present the crystallography related to the CZTS absorber layer and the challenges of identifying the secondary phases. Next, we introduce the principles of the 3DXRD technique within the context of the absorber layer microstructure and present an appropriate data analysis pipeline. Subsequently, we present an experimental 3DXRD study of CZTS including the sample details, and the results. In the final part, we discuss the connection of the results to photovoltaics properties and how recent developments of 3DXRD can advance the characterization of thin films even further.

Crystallographic aspects of kesterite
First, we must distinguish between kesterite and disordered kesterite, the latter the most frequently observed structures for CZTS. X-ray and Neutron studies have demonstrated that the quaternary compound CZTS, crystallizes in the kesterite structure (I-4) [37], [38]. The "disordered" kesterite structure associated with space group (I-42m ) was first observed by Schorr [39]. In this phase Cu and Zn cations intermix in the Cu-Zn layers (z=1/4, 3/4) of stoichiometric CZTS [39], see Figure  1a. The critical temperature for the phase transition from the ordered to disordered kesterite is reported to be in the range Tc = 480-560 K [40], [18], [41]. These temperatures are all below the annealing treatments under which CZTS is usually grown (720-830 K). Therefore, disordered kesterite will form initially, and ordering among Cu and Zn can only be controlled during the cooling process.
Moreover, CZTS films are grown in Cu-poor, Zn-rich conditions to obtain high-efficiency devices [42], [43]. The off-stoichiometric CZTS maintains the kesterite-type structure with variations in the lattice parameters due to the altering composition and the cation disorder [44]. The pure-phase kesterite phase only exists in a narrow area of the ZnS-CuS2-SnS2 phase diagram [45], [46]. Thus, secondary phases tend to form in the off-stoichiometric films, for instance, ZnS with a facecentered cubic crystal structure (F-43m), Figure 1 b). The two structures are closely related: doubling the a, b, or c axis of the cubic structure of the ZnS yields a unit cell corresponding to kesterite and with nearly identical lattice parameters, see Table 1.   3DXRD is a well-established tool for non-destructive characterization of grains in 3D. Based on the use of a monochromatic beam from a synchrotron source, the experimental geometry is sketched in Figure 2. Similar to the rotation method, diffraction images are acquired while rotating the sample around an axis (ω) perpendicular to the incoming beam. It images the intensity of the diffraction spots originating from the individual grains. Figure 2 displays a stack of recorded diffraction images for a small rotation range, showing the evolution of the intensity within a region of interest comprising the diffraction spot from one reflection. Typically, a focused line beam is used, thereby characterizing one slice in the sample. To provide 3D information, the sample is then and translated along z, and the data acquisition is repeated. In this way, one characterizes multiple slices that correspond to consecutive z-positions in the sample. In the far-field version of 3DXRD, which is of interest here, the sample-detector distance is relatively large (tens of centimeters to meters), and the size of the detector pixels (a few hundred µm) similar to the size of the sample. This geometry is optimized for high angular resolution. For the relation between experimental observables (position of diffraction peaks on the detector and corresponding rotation angle ) and reciprocal space, we shall follow the FABLE conventions [50]. Figure 3 shows a plot of the position of all the harvested reflections from all the slices and their azimuthal angle η, according to the FABLE protocols. Let Gl be the reciprocal lattice vector corresponding to lattice planes (h,k,l) in a particular grain of interest, as defined in the laboratory system (see Figure 2). Let Gs be the same diffraction vector defined in the sample systemfixed with respect to the sample. The diffraction condition is fulfilled when Here  is the rotation matrix corresponding to rotation around the -axis, U is a matrix representing the orientation of a grain of interest, () are polar coordinates characterizing the direction of the diffracted beam, see Figure 2, and | ⃗ | = 2sin ( ) = 1 is given from Bragg's law, with d representing the spacing between crystallographic planes. B is a matrix that comprises information about the unit cell as expressed by reciprocal lattice constants (a*,b*,c*,       ): As usual in crystallography, the matrix A = (B -1 ) T , where (…) T symbolizes transposing, comprises the corresponding information about the direct space unit cell of the grain of interest, expressed by the direct space lattice constants (a, b, c, ). Notably, the matrix inverse of (UB) -1 gives the real space unit cell vectors ( ⃗, ⃗ ⃗ , ⃗) of the grain in the sample frame.
The grain elastic strain can be expressed in terms of the unit cell of a reference (unstrained) crystal A0 and a strained crystal A. We determine the deformation gradient tensor of the grain as: For the small strain levels of relevance to this study and in the absence of rotation, the infinitesimal strain tensor is applicable and is, by definition, given by the symmetric tensor where I is the identity matrix. Ultimately the strain is transformed in the sample coordinate system by applying the orientation of the grain: 4.2 Conventional 3DXRD data analysis and its relation to polycrystalline photovoltaic materials 3DXRD methods are usually applied to studies of polycrystalline materials. It requires knowing the space group and unit cell lattice parameters for the unstrained material (that is, with a known matrix B0). The position of diffraction spots on the detector are given by the grain orientation with only small perturbations due to strain. In this case, one may initially assume that all grains are associated with the undistorted matrix B0. The process of identifying grains, multigrain indexing, then becomes a question of identifying orientations, U, that complies with Eq.(1), for a set of known (h,k,l) indices. As a result, the reflections determined are sorted into groups, where each group represents one grain. The main limitation is the overlap of diffraction spots. For inorganic materials exhibiting weak textures, up to around 3000 grains can be indexed [28] from a single rotation scan. Our approach is to utilize a line beam that limits the number of simultaneously illuminated grains to avoid spot overlaps.
Following this indexing step, all the tools of single-crystal crystallography can be applied to each grain. The relative grain volume can be estimated from the integrated intensities of the assigned reflections. A least-square fit can be performed to determine all nine U and B components, by minimizing the angular distance between the predicted reflections, cf. Eq. (1), and the experimentally determined ones. Next the strain tensor can be determined from B by Eqs. (3)-(4), where the grain unit cell, determined during indexing, is compared to an unstrained reference B0 [26]. When required, crystal structure refinement can also be used to optimize the position of the atoms within the unit cellwith a quality in the results that can match that of refinement based on single-crystal diffraction [51], [28].
In principle, the 3DXRD formalism, as expressed by Eq. (1), allows for indexing without any prior information by operating in the 9D space, spanned by U and B. In this way, 3DXRD could handle any number of arbitrary unknown phases, strained or unstrained. However, brute force procedures are too slow to be operational. A general-purpose method involving searching only in 3D has been suggested [36], but this algorithm still lacks a sufficiently robust software implementation. In this work, phase identification from a database search could provide sufficiently accurate for the unit cell parameters of the phases in the sample.
The polycrystalline photovoltaic materials, and particularly the kesterite solar cell, pose a special challenge as several complications are present simultaneously:  Several phases are present, and some might not be known a priori.  Some phases may exhibit a doubling of the unit cell in one direction, and their lattice parameters give rise to -angles that are nearly indistinguishable.  Twins may appear, leading to a large fraction of reflections being shared by more than one grain.  The specimens are subject to mechanical stress, originating in the thin film preparation.  The grains are sub-micron in size leading to signal to noise ratios (S/N) in the diffraction data.
These CZTS data comprise additional information, as the doubling of a unit cell leads to superstructure peaks. These may, however, be weak, and spurious peaks from other phases can cause unexpected overlaps.
In the following, we present an approach that overcomes these challenges and generates a list of grains. Each grain is associated with an orientation, a size, and a phase related to a unit cell. The unit cell parameters represent a strained state, caused by stoichiometry changes and an externally imposed mechanical strain. We demonstrate how to calculate an overall strain for the film and subtract it from the oriented unit cells of the grains.
In section 4.3, we present our approach to indexing the grains and identifying their unstrained unit cells. Section 4.4 describes how to exploit the results for the statistical description of phases, grains size distributions, stresses, and texture with a special focus on potential twin relations.
In an initial exploratory phase, we discovered that a data analysis based on the existence of two phases, a cubic and a tetragonal, was consistent with the data. Figure 3b displays the lines associated with the cubic and tetragonal phases. Hence, we shall assume two phases in the following. The data analysis pipeline used is sketched in Figure 4. It is structured in three parts. Its implementation is based on existing 3DXRD software, throughout, primarily the ImageD11 software [52].

Identifying grains and their crystal structure.
In the first part, the experimental data from different slices of the sample (different z-positions) are analyzed independently. For each slice, initially, a background is subtracted from the raw images, and the diffraction spots are identified ("peak search"). Based on the statistics of these reflections, several global parameters related to the experimental setup are refined, including wavelength, sample-to-detector distance, and tilts of the detector. We assume all grains belong to the same phase with a cubic symmetry and an "average" lattice constant, a0 * , corresponding to a common B0 = a0 * I. Excluding superstructure peaks and using only diffraction spots positioned at anglescorresponding to the cubic phase (the three green lines marked in Fig 3b), grains are found by the classical monophase 3DXRD indexing program ImageD11. The result is a list of grains, each with an associated (UB) -1 matrix and a list of reflections.
Next, we assume that the mechanical stress gives rise to comparable strains in the grains if asserted in the sample coordinate system. The grain strain can then be expressed in terms of an "average" contribution and a residual that is specific to the grain. The average strain tensor of the film may be determined in several ways. The approach used here is to focus on the diffraction spots belonging to a specific (hkl) family. For each diffraction spot, i, we can determine the shift in  position, ∆2 ℎ . The corresponding normalized scattering vector in the sample system is ⃗⃗ = ⃗⃗⃗⃗⃗⃗ | | ⁄ . From the differentiation of Bragg's law, the  shift corresponds to an axial strain (a strain in the direction ⃗⃗ ) of Here 2θ 0 hkl is the average two-theta angle of the (hkl) Debye-Scherrer ring. By definition, ) ( ).
From this follows that the average strain tensor elements for the hkl ring, ij, can be determined by a linear least-squares fit of experimental data to Eqs. (6) and (7). Let the resulting matrix be  mat . Subsequently, for each grain,  mat is subtracted to correct the unit cell parameters and obtain the "strain-free" lattice parameters.
Next, the grain strain tensor is calculated applying Eq. (4) and the obtained "strain-free" unit cell as the unstrained reference. Then, we calculate the stress tensor using Hook's law Subsequently, the stress in the film is obtained by transforming the grain into the sample reference system.
In the second part, the superstructure peaks are taken into consideration. These appear at positions in the reciprocal space that corresponds to a doubling of the direct space unit cell. To study this systematically, for each grain, we form the supercell (2 ⃗, 2 ⃗ ⃗ , 2 ⃗) in direct space, a doubling in all directions of the unit cell ( ⃗, ⃗ ⃗ , ⃗). The reciprocal space unit cell is correspondingly halved in all directions. A search is now performed within the full set of reflections from the original peak search of reflections positioned at the nodal points of the supercell. The reflections appearing with an odd number in any of the three directions are "superstructure peaks"; they do not belong to the original cell. Searching for grains with a double unit cell in the ⃗ direction, the number of reflections appearing at (h, 2k, 2l) with h odd is compared to the number of reflections assigned to the supercell. If the ratio is above a certain threshold, defined by S/N and spurious background, the grain is defined to have a double cell with the preferred axis along ⃗. Likewise, searches are performed on unit cells that are doubled along ⃗ ⃗ or ⃗. It is generally observed that the shortest unit cell lattice and the axis with odd reflections coincide. No occurrence is found of cells being doubled in more than one direction. Figure 5 shows a plot of the ratio of the number of superstructure diffraction spots to the total number of diffraction spots for a given grain as a function of the derived grain volume. For large grains, most the weak diffraction spots at low scattering angles from the superstructure are recorded. Hence, it is straightforward to classify a grain as tetragonal or cubic based on the existence of superstructure peaks. As expected, with a decreasing grain volume, more of the superstructure peaks become too faint to be recorded. For relatively small grains with few or no superstructure reflections, the distinction between tetragonal or cubic was therefore based on the unit cell geometry. This classification scheme turns out to be very robust.
For grains that are classified as tetragonal, the shortest axis is doubled. Moreover, the unit cell axes were permuted, such that this short axis becomes the ⃗-axis.
The (illuminated) volume of the grain, assuming proportionality of the volume with the reflection intensities, is given by Eq. (11), where ∑ I ̅ g is the sum of the average intensities of all indexed grains, and Vsample is the illuminated sample volume.
= ̅ ∑ ̅ × (11) Figure 5. The ratio of the superstructure peaks to the total number of peaks for a given grain is compared to the grain volume.
In the third part of the flow diagram, the slices are combined. Grains in neighboring slices with a similar (UB) -1 within given tolerances are considered identical. For these, the weighted average of (UB) -1 by the integrated intensities is determined. Moreover, the (sub) volumes determined in the different slices are added. Following this, for each grain, a final optimization step is performed. Here all parameters in (UB) -1 and the 3D center-of-mass position of the grain (which will cause minor translations of diffraction peaks not described by Eq. 1) are refined.
The result of the data analysis is, therefore, a division of the grains into two "phases": one with no cell doubling and with a unit cell that is close to cubic, and another with a unit cell that is close to a tetragonal unit cell with unit cell parameters (a, a, 2a, 90 o , 90 o , 90 o ). We shall term these the "cubic" phase associated with ZnS and the "tetragonal phase" associated with CZTS. Each grain is associated with an orientation, U, a B matrix, a center-of-mass position, and a list of reflections. Moreover, we have determined the average strain tensor within the (illuminated part of the) sample.

Statistics on the structural and mechanical properties of the grain ensemble
As already stated, once the grains have been identified, their properties can be subject to statistical analysis for understanding both average properties and the heterogeneity within the sample. The local texture (within the volume studied) can be derived from the grain orientations. The grain sizes and their distribution can be determined from the integrated intensities. Strain components and their distributions can be derived via Eqs. (3)-(4), changes in stoichiometry may be inferred directly from changes in unit cell parameters.
The spatial resolution of far-field 3DXRD did not allow us to identify neighboring grains, and therefore their misorientation angle is not accessible. But twins can be identified from angular relationships. In practice, reflections can be shared among grains because of coincidental overlap or due to the presence of twin relationships between grains, so thresholds must be introduced in the data analysis. Notably, twin relationships may also exist between the tetragonal and the cubic phase.
We calculate the number of shared reflections between grains by computing the scattering vectors of each grain to look for overlapping of Bragg peaks. Using equation (1), the error of the computed hkl should be below ~0.02 to be considered as part of the grain. We confirm the twin relation among pairs with 30% of reflections overlapping, if the pair has a certain misorientation angle associated to a symmetry operation.
When comparing a pair of grains, we compute a natural lattice mismatch via the deformation gradient tensor F. Based on Eq. (3), here we utilize the reciprocal lattice of a grain ( ) , considered as the reference lattice that has been deformed by the inverse transpose of F when compared to grain m, whose reciprocal lattice is ( ) .

= ( ) − • ( )
Following conventions in the field of continuum mechanics, we perform a polar decomposition of F to produce a pure rotation R and a pure stretch tensor Us, the right stretch tensor (not to be confused with U-rotation matrix in 3DXRD) Next we calculate the Biot strain tensor Eq. (14), which is equivalent to the infinitesimal strain tensor in the absence of rotations. The overall strain magnitude is given by the Euclidean norm of the Biot strain tensor ‖ ‖.
This rotation matrix R, is used to determine the angle and axis of rotation between the two grains being compared. The equivalent symmetry operations are applied to the reference grain, calculating the strain magnitude. The symmetry transformation with the lowest strain value describes the misorientation angle between two grains, where the lattices are also well matched. For the tetragonal structure, eight transformations are possible, whereas the cubic structure allows 24.

Experimental 5.1 Sample description and preparation
The investigated kesterite solar cell architecture consists of a stack of layers deposited on a molybdenum coated soda-lime glass substrate (Mo-SLG). The CZTS absorber layer is fabricated by pulsed laser deposition (PLD) and annealed in a high-temperature sulfurized atmosphere to form the polycrystalline kesterite film (~400 nm thick). The subsequent coatings are a CdS buffer layer (60 nm), an intrinsic ZnO window layer (50 nm), an Indium Tin Oxide (ITO) contact layer (200 nm), and an MgF2 anti-reflection coating (100 nm). The fabrication details can be found in [57].
We cut a 40 (W) x 300 (H) µm 2 piece of the solar cell, as shown schematically in Figure 6a. To maximize the signal to noise ratio of the diffracted intensity originating from the 1 µm sized grains, we reduced the 1 mm thickness of the Mo-SLG by mechanical polishing and milling by a focused ion beam (FIB) down to 4 µm thickness (see Figure 6b).

3DXRD experiment
The experiment was carried out at the Advanced Photon Source synchrotron at the 1-ID beamline.
A monochromatic X-ray beam (52 keV

Results
Following the procedure outlined above, a total of 597 grains were identified, 582 tetragonal and 15 cubic. As a figure-of-merit we note that 33% of the diffraction spots identified by the peak search were assigned to these grains. It is possible to identify more grains, but they will be associated with larger errors.

Averaged strain and stress in the sample and the grains.
Shown in Fig 7a is the variation of the 2θ angle with the azimuthal angle, , and with the rotation angle,  for all reflections in the (103) lattice plane. We define the reference angle 2 0 103 as the average angle. There is a systematic displacement of the 2θ angle with both  and , which we attribute to an external stress field. Hence, we can determine a strain tensor, corresponding to this external stress, using the fitting procedure formulated in Eqs. (6)-(7).
The average strain tensor elements for each of the 7 slices obtained by this least-square approximation are listed in Table 2. By subtracting the strain in the tetragonal and cubic grains, we can observe the improvement of the lattice parameters in both phases, as the distributions of the corrected lattice parameters become narrower (see Figure 7b, c).  In both b) and c) the term the "measured" denotes the original data while the term "corrected" are the results after subtracting the effect of an external stress acting on the film.
Subsequently, for each grain we can derive its strain tensor elements using the corrected unit cells as the reference unstrained state and Eqs. (3)-(4). Figure 8a defines the elements of the strain tensor in the tetragonal crystal structure. The normal strains directions of ε1, ε2, ε3 are perpendicular to the (100), (010), (001) planes, whereas the shear strains ε4, ε5, and ε6 are coplanar to the planes where the normal strain is applied. The shear strains are in equilibrium, implying that ε4: ε23= ε32, ε5: ε13=ε31, ε6: ε12= ε21, and therefore they appear in two planes. Figure 8b shows the resulting correlation between the normal strain components ε1 and ε3. We observe the typical behavior of a deformed object where an increasing strain along the a-lattice will result in a decreasing strain along the c-lattice. The slope is determined to be -0.83 by a linear fit to the data. However, we also note a substantial scatter in these data, caused by grain-grain interactions.
Next, we use Eq. (9) and the corresponding elastic constants to calculate the stress components for each grain. Histograms of these components are shown in Fig 8c for the tetragonal phase. The normal stress along the a-direction σ1 has a slightly right-skewed distribution suggesting predominant compressive stresses, whereas the normal stresses along b (σ2) and c (σ3) directions have left-skewed distributions that correspond to tensile stresses. The shear stress σ4 has an almost normal distribution, whereas σ5 has a bimodal distribution, and σ6 a left-skewed one. We can now calculate the macroscopic stresses in the sample by averaging over the grains while taking into consideration their orientation. The results are listed for each slice in Table 3 and displayed in Figure 9a. We visualize the normal stresses σ1, σ2, and σ3 along the x-, y-, and z-axes correspondingly. The shear stresses are coplanar to the planes where the normal stresses are applied. Similar to the strain depiction, the shear stresses appear in two planes as they are in equilibrium, meaning that σ4: σ23=σ32, σ5: σ13=σ31, and σ6: σ12= σ21.
The film shows a compressive strain in the normal direction of the film plane that corresponds to an averaged compressive stress of -144.7 MPa. The corresponding tensile strain within the film plane results in an average stress of 53.6 MPa along the y-axis and 81.9 MPa along the z-axis (Figure 9). The shear strains ε4 and ε5 are zero within experimental error. On the other hand, the non-zero component ε6, coplanar to the (xz)-and the (yz)-planes, is associated with a compressive shear stress σ6 = -38.1 MPa. This can be explained by a misalignment of the sample in the diffractometer.  (Eq. (7)).

Grain size distribution and orientations
The grain volumes were derived using Eq. (11). The histogram for the tetragonal phase shown in Fig 10 a) and b) are consistent with a log-normal size distribution with a cut-off at lower radii due to thresholding of the intensities. The grain size is 0.32 ±0.26 µm 3 , and the corresponding radius 0.47 ±0.18 µm (see Figure 10a, b). The radius is obtained assuming a cylindrical volume for the CZTS grain with a height of 0.45 µm, the film thickness. For the cubic grains, the grain statistics is scarce. The volume is 0.25 ±0.16 µm 3 , Figure 10c), and the corresponding equivalent-sphere radius is 0.43 ±0.13 µm, (see Figure 10d). The orientation distribution functions were computed using the MTEX MATLAB toolbox and a 5° resolution [58]. The pole figures of the tetragonal grains for the planes {100}, {110}, {001} and {112} are shown in Figure 11, with the yz-plane being the film plane on the diffractometer. The pole figures for the cubic grains are not shown because of poor grain statistics (15 grains).

Figure 11
Pole figures of the tetragonal grains in multiples of random distribution (mrd).

Twin boundaries
Twin boundaries are often described with the quantity Σ, which is defined by the ratio between the area enclosed by a unit cell of the coincident lattice sites and the standard unit cell.
To identify twins among the tetragonal grains a search was performed for pairs of grains sharing 30% of the reflections or more. We find that the misorientation angles of these in all cases correspond to one of four values, corresponding to the symmetry operations with the minimum strain between the compared grains.
120 twin pairs were identified. Their resulting misorientation angles are shown in Figure 12ad.The most frequent type of grain boundary is identified as Σ3, characterized by a rotation of 70.53° and the corresponding symmetrically equivalent misorientation angles 109.47°, 131.81°, and 180°. Σ3 boundaries are also detected by EBSD in chalcopyrite thin films [59]. Our measurements deviate from the mentioned angles due to the tetragonal distortion c/2a< 1 in the kesterite structure. We report the average rotation angles, the corresponding rotation axis and the transformed lattice plane in Table 4.  Among the cubic grains, only one pair of twinned grains was found with a rotation of 180°-<211>cubic transforming the plane (211)cubic. This transformation corresponds to a Σ3 twin boundary. The equivalent symmetric transformation is the Σ3 60°-<111>cubic, typically found in the diamond-type structure [59].
Special orientation relationships between grains of the different phases were not considered.

The kesterite solar cell
Identification of secondary phases. The 3DXRD analysis of the CZTS absorber layer revealed the presence of the ZnS phase representing 2.5% of the total number of grains in the film. This small amount of ZnS is not detrimental to CZTS devices. However, when present in large amounts, it can block the charge transport or increase the series resistance in the solar cell [60].
A similar value has been measured in a sputtered CZTS film with 3.1% of ZnS by X-ray absorption near-edge spectroscopy (XANES) at the sulfur K-edge [61]. A co-sputtered CZTS film, also measured by XANES, yielded 10.5% of ZnS [61]. These measurements are within the ZnS limit detection in XANES of 3% [24].
Film averaged stress. The stress and strain components vary slightly between slices with an average standard error of 2.23 MPa. These almost constant values are a testament to the robustness of the 3DXRD method.
Comparing to literature, Johnson et al. consider the formation of a biaxial tensile thermal stress at the Mo/CZTS interface during annealing as a result of the Mo deposition stress and the thermal expansion mismatch stresses [62]. In their study, wafer curvature measurements show compressive deposition stresses of about -400 MPa to -38MPa for optimized Mo sputtering deposition, and a compressive deposition stress of about -100 MPa for CZTS by co-sputtering. The normal tensile stresses over the PLD-CZTS film plane agree with the biaxial tensile stresses of the co-sputtered CZTS, although the compressive normal element and the shear components are missing in their model. Moreover, our calculated stresses are in the same order of magnitude as their deposition stresses.
In addition to the thermal stresses inflicted during annealing, we can also consider the entangled combination of the chemical distortion caused by the off-stoichiometric composition of the film and the inflicted mechanical deformation through the cutting of the solar cell and the removal of the glass through mechanical polishing. Likewise, FIB etching could have introduced artifacts such as ion implantation and structural damage [63], [64]. Unfortunately, sample preparation is unavoidable as the 5 mm thick amorphous glass substrate causes a background that buries the signal of the grains.
Grain averaged stress and strain. At the grain level, the slope derived from linear fitting of the normal strain components ε1 and ε3 shown in Fig 8b, -0.83, is larger than the biaxial relaxation coefficient of -1.23 reported by Li et al. [16]. The forces along the c-axis are set to zero in their model, whereas our stress measurements have shown non-zero elements. The same study by Li et al. calculated an increase in the bandgap with the increase of compressive biaxial strain for ε1 above -1.5% and a decrease for ε1 below -1.5%. Our ε1 strain values oscillate between compressive and tensile strain in a range of [-25; 20]×10 -4 , which implies that the bandgap is not homogeneous among the grains. The overall bandgap of the film is the result of contributions from all of these grain bandgaps. We also note that the strain is on the order of 10 -4 , which agrees with the strains measured in CuInSe2 thin-film solar cells [65].
Grain properties. The standard error of the grain variations in angle, length and volume among the measured tetragonal unit cell parameters is 0.008°, 8.8×10 -4 Å, and 0.012 Å 3 respectively, demonstrating a high accuracy in the 3DXRD measurements.
The estimation of the grain size based on the intensities of the reflections agrees with the scale shown in the SEM image (Figure 6c). The grain volumes of CZTS grains are larger than those of the ZnS grains.
The main conclusion from the texture analysis is that the <112> poles are preferably aligned to the normal direction of the film, whereas a faint discontinuous ring in the same pole figure could indicate a weak fiber texture. Moreover, the poles <001> and <001> are almost aligned parallel to the surface of the film. This out-of-plane (112) fiber texture has also been observed in a cosputtered CZTS film [66]. The link between this texture and efficiency is not clearly established, but many studies report the (112) preferred orientation in CZTS films deposited by PLD [67]. In CIGS films, a (200)/(204) preferred orientations yields higher efficiencies than CIGS films with (112) orientations [68].
Twin boundaries. The six variants of Σ3-type twin boundaries have also been observed with electron microscopy in CIGS, CGS, and CIS solar cells [59]. Σ3 boundaries have lower defect density compared to random grain boundaries according to electron-beam scattering diffraction (EBSD) and cathodoluminescence (CL) measurements [69].
The formation energies of Σ3 are low, and hence they are common in CZTS films. In our results, 41.2% of the total number of grains are Σ3 twins. The most frequent twin operation is the Σ3{112} that corresponds to the 180° around <221>. Characterization of grain Σ3 {112} twin boundaries in CIGS has been done extensively, revealing a rather benign electronic behavior [69]. Additionally, the formation of Cu vacancies and InGa antisite defects in Σ3 {112} have been experimentally confirmed [70]. Similarly, one could expect the development of defects in CZTS Σ3{112} twin boundaries. In a theoretical study by Wong et al. [71], Σ3{112} grain boundaries are constructed with ZnCu and SnZn defects based on anion-anion terminations. Such grain boundaries are detrimental and can acts as seeds for secondary phases. According to this model, we could speculate that ZnS grains might be lying close to the Σ3{112} twin boundaries. Moreover, first principle calculations have predicted Cu-poor anion terminated (-1-1-2) surfaces to situate VCu defects, which are benign for the solar cell performance [72].

3DXRD limitations and new horizons
This paper demonstrates that far-field 3DXRD is suitable for providing comprehensive statistical information about the ensemble of grains in the absorber layer. However, the position of the grain is not resolved. In outlook, one can make mapping of grains in 3D using the 3DXRD scanning modality that employs a smaller X-ray beam (200 nm) and records diffraction patterns at each yzposition of the sample [73]- [75]. Thus, grain positions and strain maps with a higher resolution can be achieved [76]. A drawback with scanning techniques is that the acquisition time increases. However, the next generation of synchrotron sources, such as the Extremely Brilliant Source, EBS, in Grenoble, which was successfully put into operation in Summer 2020, promises an increase in the data acquisition speed of all types of 3DXRD modalities by a factor of 10-50. Moreover, preliminary results on a new full field modality known as High-resolution 3DXRD suggests that 300 nm spatial resolution is within reach [32].
In outlook, combining scanning 3DXRD with X-ray Beam Induced Current (XBIC) and X-ray fluorescence (XRF) could reveal the relation between microstructure and photovoltaics properties of the device and localized elemental composition. By performing in-operando studies implementing the mentioned techniques, we could analyse the effects of grain boundaries on PV properties and identify elemental clusters that tend to populate the grain boundaries for passivation. These studies could improve current models for thin-film optimization [71], [77], [78].

Conclusion
We have characterized the microstructure of a PLD-deposited CZTS absorber layer buried within the stack of layers that constitute a full solar cell device. We demonstrate that 3DXRD can distinguish between phases with nearly identical unit cell parameters. As a result, we found 597 grains; 582 were identified as tetragonal and 15 as cubic. We extracted the strain and stress components both at the sample and at the grain level. We provided extensive statistics of the tetragonal and cubic grains, including the number of grains, sizes, orientations, and twin boundaries of each phase and discussed the relevance of this information for CZTS design.
More generally, the most common photovoltaic thin-film materials are chalcogenides with cubic and tetragonal structures. Structural characterization of these with traditional methods is hampered by the same issues as CZTS. Hence, we propose that the 3DXRD methodology may be applied to index grains of other absorber materials such as CdTe (F-43m) and CIGS (I-42d).