Xray and neutron scattering on disordered nanosize clusters: a case study of leadzirconatetitanate solid solutions
Abstract
Defects and frequently used defect models of solids are reviewed. Signatures for identifying the disorder from xray and neutron scattering data are given. To give illustrative examples how technologically important defects contribute to xray and neutron scattering numerical method able to treat nonperiodical solids possessing several simultaneous defect types is given for simulating scattering in nanosize disordered clusters. The approach takes particle size, shape, and defects into account and isolates element specific signals. As a case study a statistical approximation model for leadzirconate titanate [Pb(ZrTi)O, PZT] is introduced. PZT is a material possessing several defect types, including substitutional, displacement and surface defects. Spatial composition variation is taken into account by introducing a model in which the edge lengths of each cell depend on the distribution of Zr and Ti ions in the cluster. Spatially varying edge lengths and angles is referred to as microstrain. The model is applied to compute the scattering from ellipsoid shaped PZT clusters and to simulate the structural changes as a function of average composition. Twophase coexistence range, the so called morphotropic phase boundary composition is given correctly. The composition at which the rhombohedral and tetragonal cells are equally abundant was . Selected xray and neutron Bragg reflection intensities and line shapes were simulated. Examples of the effect of size and shape of the scattering clusters on diffraction patterns are given and the particle dimensions, computed through Scherrer equation, are compared with the exact cluster dimensions. Scattering from two types of 180 domains in spherical particles, one type assigned to Tirich PZT and the second to the MPB and Zrrich PZT, is computed. We show how the method can be used for modelling polarization reversal.
Contents
1 Introduction
The paper is organized as follows. Section 1.1 summarizes the basic concepts of xray and neutron scattering, sections 1.2 and 1.3 review the models of disordered materials and summarize the pair distribution function method, respectively, and finally a brief review of the crystal structures of leadzirconatetitanate [Pb(ZrTi)O, PZT] is given in section 1.4. PZT was chosen as an example material due to its wide use in applications, which are often based on the controlled use of defects. Section 2 describes the numerical method applied in the present study. Section 2.1 describes the method after which selected case studies on PZT clusters are given in sections 2.2, 2.3 and 2.4.
1.1 Xray and neutron diffraction.
Diffraction techniques, notably xray diffraction, are the most commonly applied materials characterization methods. The theory of xray scattering and diffraction is well established and below key principles (see, e.g., ref. [1]) are given. The farfield amplitude is given by Eq. (1)
(1) 
where is the number of atoms in the scattering volume, and are the scattering amplitude and the position vector of the atom n, respectively, and is the scattering vector, depicted in Fig. 1. By denoting the unit propagation vectors of the incoming and scattered radiation of wavelength by and respectively, is . In crystals the intensity maxima correspond to the reciprocal lattice points , where , and are the reciprocal lattice vectors and , and are integers. Eq. (1) can be generalized by replacing the discrete atomic densities by a continuous electron density :
(2) 
so that is the Fourier transform of . Inversely, electron density is given by Eq. (3)
(3) 
The scattered intensity is given by Eq. (4)
(4) 
and, in terms of the convolution obtained by substituting Eq. (1) to Eq. (4), which reads that the Fouriertransform of the autocorrelation function (the Patterson function) equals to the intensity. In the case of an infinite crystal the electron density can be expressed as a convolution between the function representing the electron density inside a unit cell and a series of Dirac functions representing the crystal lattice. Fouriertransform of the convolution results in the structure factor and further gives the wellknown result according to which the intensity of the reflection is . Section 2 focuses on the nanosize clusters which possess shortrange order but lack translational symmetry. The nanosize itself does not remove the translational symmetry: we consider cases in which each ’unit cell’ has own ’lattice parameters’ (i.e., the assumption of translational symmetry is abandoned).
If the absorption, extinction, thermal displacements, angle dependent polarization corrections and instrument related factors are not considered, the elastic scattering intensity can readily be computed from Eq. (5) once the atomic scale structure is known:
(5) 
Neutron diffraction intensity is modelled similarly, the scattering amplitudes are replaced by nuclei specific neutron scattering lengths . Since xrays interact with electrons and neutrons interact with nuclei the techniques are in many ways complementary. In contrast to xrays, neutrons frequently scatter strongly from light nuclei, such as oxygen and hydrogen, which is essentially why neutrons suit for the determination of light elements positions. The neutron scattering lengths of different isotopes are often quite different, even possessing different signs, which has made isotope substitution a technique for pinpointing structural details. Since neutrons possess a magnetic moment they interact with the magnetic moments of electrons, the cross section being the same order of magnitude as the neutronnuclear interaction. Correspondingly, neutron scattering is the most common method for determining magnetic ordering. Magnetic form factors used in the calculations of the cross sections for magnetic scattering of neutrons are given in terms of sums of exponential functions whose coefficients are tabulated in ref. [2]. The form factors decay with increasing so that often the minimum spacing included in the magnetic scattering model is Å. The form factors depend on the valence state of the ions and thus it is frequently necessary to complement neutron scattering data by alternative measurements to clarify the valence state(s) of the ions. Since the scattering power of atoms for neutrons is not dependent (), in contrast to xrays where the atomic scattering factors fall away rapidly at high, strong diffuse scattering can be observed well beyond the range where Bragg peaks occur [3].
1.2 Common Models for Disordered Systems
Computation of diffraction intensities is straightforward once the positions of atoms are given (see, for instance refs. [1] and [4]). The most challenging problem is to find the atomic scale structure (structural model) corresponding to the measured intensity. Numerous recipes to solve the problem have been developed, for instance see refs. [4] and [5]. Space group symmetry determination from the Laue symmetry and the reflection conditions, as obtained from the diffraction patterns, is given in ref. [6]. At the final stage of the structure solving task one introduces a model which is refined by adjusting model parameters so that the difference between the computed intensity and the measured intensity is minimized. A wellknown refinement technique for powders is the Rietveld method [7, 8]. In the case of known average symmetry, equation (5) is not applied directly, but one computes the squared value of the absolute value of the structure factor, of the Bragg reflection . , however, tells nothing about the linewidths or the shape of the profile function of the reflection. Numerous line shapes were derived for different instruments and sometimes the choice of the profile function for describing the sample and instrument originating line broadening is challenging. The broadening parameters are refined as a part of the structural model. More recently, instrument related line broadening has also been computed from the known diffractometer properties as summarized in ref. [9]. In principle, direct application of Eq. (5) gives the sample contribution to the lineshape.
Direct application of equation (5) is a heavy computational task even for relatively small atom clusters. Generally, the problem is challenging once the material lacks periodicity at least in one dimension as even a straightforward computation of the scattered intensity from a collection of atoms with known positions becomes computationally formidable task (see, however, Model #5 in section 1.2). Thus, to make structure solving or even refinement possible approximations are required. The most complex structures require tailored solutions as they do not possess the high symmetry required by the commonly available programs.
Defects in solids.
Defects in solids are classified as zerodimensional (point defects), onedimensional (line defects) and twodimensional (surfaces), see Fig.2.
Point defects introduce a change in the scattering power by interfering the coherence between scattered waves by introducing a disturbance to the scattering amplitudes (or lengths) and scatterer positions. In diffraction experiments this effect is hardly seen if the defect concentration is small. However, once the concentration increases not only is the scattering power changing via change in average (substitutional disorder) but also bond lengths and angles are spatially varying (displacement disorder). Mathematical treatment of displacement disorder correlated with substitutional disorder in solidsolutions is given in ref. [1]. Phenomenologically the displacement disorder following the substitutional disorder can be described using the concept of microstrain introduced in ref. [10] and it is commonly used to model the dependent line broadening in Rietveld refinement: for instance GSAS [11] Rietveld refinement software has an option to use profile functions with microstrain broadening formulated in ref. [10]. Table 1 summarizes common disorder types observed in solids and the characteristic signatures in xray scattering and diffraction.
The models given in Table 1 may appear rather specific. However, the tabulated signatures are also found in more general cases. For instance, the signal related to the displacement disorder increases with increasing reflection index and vanishes for small values of . The signatures can be used for distinguishing different structural models. For example, both sinusoidal composition modulation and sinusoidal displacement modulation produce satellite reflections at the same spacing, if the propagation vector is the same. The two cases can be identified by their different intensity dependencies. Also the case of correlated substitutional and displacement disorder can be distinguished from the plain substitutional and plain displacement disorder by the characteristic asymmetric intensity of the satellite peaks of each pair. We note that the pure sizerelated line broadening gives the same particle size estimate, no matter which reflection is chosen. However, if the broadening is not solely due to the particle size it is important to correctly model the disorder. Section 2 considers several cases in which nonperiodical displacement and substitutional disorder are correlated and the results are found to be consistent with the observations seen in periodically modulated cases.
Model #1 is a common and straightforward way to describe scattering from a wellmixed solid solution formed by nottoodifferent size constituent elements. The success of the Model #1 is essentially due to the fact that diffraction tends to emphasize average structure and to suppress deviations from it. However, if the composition and atomic scattering factors vary periodically, for instance as a sine wave form with a propagation vector , then the Bragg peaks (principal nodes) are surrounded by satellites at distances of each Bragg peak.
In the case of a pure sinusoidal composition variation (in contrast to displacement disorder Models #3 and #4 which also exhibit characteristic satellite peaks) the ratio of the intensity of the satellite reflections to that of the corresponding Bragg peak intensity is a constant for all the nodes. Following ref. [1], if atomic scattering factor varies as a sine wave, , then the diffracted wave amplitude is . Thus, diffraction pattern exhibits two satellites at distances of each node, with intensities which are times that of the principal node.
Model #2 is a wellknown description of the atomic displacement originated diffuse scattering implemented in every standard Rietveld refinement software. Model #5 is based on the generation of a scattering object from two types of layers and with proportions and , respectively. Layers are added to the lattice one at a time and the probability of the new layer being an or type is dependent only on the preceding layer [3]. The scattering object is constructed using conditional probabilities , , and , where (0) denotes that the site is occupied by an ()type layer. In the case of an infinite number of layers rather simple expression is obtained for the diffuse intensity, see Model #5 in Table 1.
Model #5 is an example of a model suitable for describing a layerbylayer crystal growth. Contemporary growth chambers are often equipped with the in situ monitoring facilities, such as Reflection High Energy Electron Diffraction facility allowing to ensure the growth of correct type of layers [12]. Optimization of growth parameters includes the adjustment of substrate temperature, gas atmosphere (e.g., oxygen gas pressure) and deposition rate related parameters (such as laser beam fluence in pulsed laser ablation deposition or sputtering power). During in situ monitoring one observes the diffraction pattern and, if required, conducts parameter adjustment until a correct phase is formed. In the case of a single layer it is rather straightforward to observe when a correct diffraction pattern emerges. Model #5 and their extensions are amenable for in situ modelling of multilayer thin film structures consisted of stacked thin layers possessing different types of layers.
Model  Disorder type  Mathematical model  Signature  Examples 

#1  Substitutional  Diffraction pattern is similar to the  Solidsolutions of similar size  
one possessing a translational  atoms. The most common  
symmetry: instead of elemental  way to model solidsolutions  
scattering factors certain sites  by the Rietveld refinement.  
possess an average scatterer .  
#2  Displacement,  Part of the diffraction line intensity  Isotropic thermal motion of  
DebyeWaller  is shifted to the background. FWHM  an atom.  
value of the peak is unaffected.  
#3  Displacement  A pair of satellite peaks of equal  Nonhomogeneous CuNiFe  
intensity are surrounding  alloys.  
each reflection .  
#4  Correlated  ,  Ratio of the intensities of the satellites  Solid solutions of different 
Substitutional and  to the normal node is asymmetric:  size atoms.  
Displacement  ,  
, .  
#5  Substitutional  Probability of the new layer  Layerbylayer crystal growth  
an or depends only on the  .  when only short range  
immediately preceding layer.  forces are important.  
#6  Dislocation  At high dislocation densities peak  
broadening . 
Dislocations come in many forms as do their structural models. Model #6 summarizes the angle dependent features for high dislocation densities. Dislocations possess technologically interesting problems as they are often strongly responding to external stimulus and thus are timedependent. The atomic scale structure is not only different in dislocation core but the atomic positions close to a dislocation are different from the surrounding matrix [13, 14], making them visible through scattering techniques. In thin films edge dislocations have a crucial role as they are favoured over the stressed thin film state once the mismatch between the substrate and thin film exceeds a threshold value [15, 16], see Fig. 3. The plasticity of metals is explained in terms of dislocations. Dislocation cores also serve as a diffusion path in many materials.
Fig. 4 shows schematically why the models applied to data collected on commonly used BraggBrentano geometry can give misleading information in the case of multilayers. First, the net intensity is not a simple superposition of the intensities scattered from individual layers and substrate but involves interfacial layers. Second, especially in the case of epitaxial thin films the scattering should be considered to take place in a large entity formed by different layers: in terms of the Eq.(5) the summation involves all the atomic pairs in the multilayer structure.
Insitu scattering measurements on crystal growth is becoming a routine experimental techniques and thus models for layerbylayer growth are required. Modelling techniques based on the Markov chain and Ising models are described in ref. [3].
Similarly the structure in the vicinity of the domain wall separating the different domains also has an atomic scale structure different from either of the domains. In 180 domain boundary the polarization direction is smoothly reversed [18] which is accompanied by a spatial variation of atomic positions. The modelling task of multilayer thin films and domain structures is similar and is given in section 2.4.
1.3 Pair distribution function method
Even though isotropic disordered materials, such as a glass or an amorphous material, do not have a longrange order, they nearly always possess shortrange order. Neutron and xray scattering techniques are commonly used for determining the distribution function giving the statistics of atom pairs. Following the treatment given in ref. [1], the principle of the pair distribution function (PDF) method can be formulated by expressing the interference function in terms of the scattered intensity and structure factor , forming the average over all orientations of the vector connecting the atoms of the pair with respect to the scattering vector (the sample is assumed to be isotropic), and expressing the interference function in terms of the pair distribution functions, Eqs. (6):
(6) 
where is the number of atoms of all types per unit volume.
The reverse Fourier transform of the interference function gives , which thus would require that is experimentally determined in all of reciprocal space. Instead, a common way is to construct a structural model and compute the PDF function, compare it with the experimental data and optimize the model parameters to minimize the difference between the computed and measured intensity. Rather recently, a new method for the calculation of xray and neutron powder diffraction patterns from the Debye scattering equation was given in ref. [19]. PDF functions were computed as an intermediate stage for computing the diffraction patterns. The method is based on the splitting of pairwise atomic interactions into two contributions, the first from latticepair vectors and the second from cellpair vectors. Illustrative application examples of the PDF method can be found from ref. [20], with an emphasis on the neutron scattering studies of silica.
In order to contrast the method described below and the PDF method we note that the latter suits well for extracting information from homogeneously disordered materials, whereas the present method is aimed to provide information about spatially confined defects, such as domain boundary. In the present work we also pinpoint the spatial location of the defect.
1.4 Solidsolution with correlated substitutional and displacement disorder: leadzirconatetitanate
Despite its longhistory [21] PZT is still an intensively studied ferroelectric oxide [22] exhibiting exceptionally high piezoelectric properties when the amount of titanium and zirconium is roughly equal [21, 23]. Numerous space group assignments for nominally similar samples can be found in the literature (for reviews, see refs. [24, 25]), which is essentially due to the localscale disorder modelled by different lowsymmetry structures.
Average symmetries.
PZT has a perovskite O crystal structure in which the cation site is statistically occupied by Ti and Zr ions. Metrically, the structure is close to a cube in which Pb cations are approximately at the cube corner, Zr and Ti cations are at the cube centre and oxygen anions are close the cube face centres. A wellknown summary of the composition and temperature dependent phases of PZT is given by the phase diagram of ref. [21]. Titaniumrich structure is traditionally modelled by assuming an ideal space group symmetry ( for ), which is achieved by placing a compositionally averaged ’pseudoatom’ at the site. The anions and cations are displaced from the centrosymmetric positions along positive and negative axis direction, respectively. For Zrrich concentrations the average structure is rhombohedral, the anion and cation displacements being along the cube diagonal. Notably challenging is the composition at which the two phases coexist [26]. The twophase coexistence is typically modelled by refining the total intensity by the superposition of the intensities from the pseudotetragonal and pseudorhombohedral phases.
The average crystal structures of PZT as a function of can be summarized thus: at room temperature the crystal structure remains tetragonal up to at which composition (termed the morphotropic phase boundary, MPB) a rather complex set of phases emerges, including a monoclinic phase [27, 28, 29] and coexisting low and hightemperature rhombohedral phases [28, 29] . As pointed out in ref. [30] there is no real boundary and the phase coexistence region is extended close to the PbZrO end of the solidsolution system. It is evident that the phase exists as judged by several highresolution powder diffraction studies, but the phase transition mechanism resulting in an average phase, and even the stability of the phase, is still under investigation as it is often linked to the extraordinary piezoelectric properties of PZT at the MPB composition. A thermodynamical study showed that an 8thorder Devonshire theory is required to explain the monoclinic phase [31], suggesting that the transition is quite unusual. This was also discussed in ref. [32], where a twoorderparameter thermodynamic model was developed for PZT to account for its peculiar features. An early report of the phase in a ferroelectric perovskite oxide is about PbNbFeO (PNF) compound [33, 34]. In the case of PZT the phase has often been claimed to favour polarization rotation and further to be responsible for the good piezoelectric properties. We note that also in the case of PNF the Fe and Nb cations are disordered. The fact that the monoclinic symmetry only tells that the polarization vector is within a mirror plane does not mean that the polarization can rotate within . Nevertheless, it has been proposed that a continuous polarization rotation between the and directions along the plane the two direction vectors span would be energetically preferable [35]. However, we do not adopt that view as detailed in refs. [25] and [36]. Instead, to understand the polarization reversal we focus on the domain wall and wall motion, which involve fairly complex timedependent structural changes [23].
At higher values of the dominant phase can, to a good approximation, be described as rhombohedral, there being two variants, with (, at low temperatures) and without (, at high temperatures) octahedral tilts. More precisely, a recent PDF and Rietveld refinement study has shown that Zrrich Pb(ZrTi)O powders possess mixed phases, described by model for and model for , where and refer to two polarization direction variants of the phase [37]. Neutron diffraction [38] and highresolution xray diffraction [39] studies verified the rhombohedral symmetry and also revealed the presence of the monoclinic phase (assigned to symmetry in ref. [38]), providing further support to the idea that the monoclinic phase is not due to the presence of adaptive phases. At the highest values of , an antiferroelectric orthorhombic phase is formed.
Deviations from average symmetries: Signatures of disorder in PZT.
Though the crystal structure models are sufficient for explaining an average structure, they are insufficient when local scale structure is considered. For instance, Raman measurements reveal that the spectra are not consistent with the average crystal structure[40, 41, 42]: the number of Raman active modes is about twice the number corresponding to an ideal structure. In diffraction experiments disorder affects the intensity distribution: The fraction of inelastically scattered intensity is increased (as discussed below in the context of atomic displacement parameter (ADP)) and the tail regions of the Bragg reflection gain intensity. The displacement disorder can be divided into Pb displacements (offsite cations) and different positions of the cations, Zr and Ti [30, 43, 44, 45, 46, 47]. Large offcentre Pb displacements are common in perovskite oxides [33, 48, 49, 50, 51]. Randomly distributed Zr and Ti atoms cause locally varying bond lengths and angles, frequently approximated by introducing so called microstrain. In contrast to the ADP induced diffuse scattering the reflection widths are dependent in microstrained samples. A random distribution of Zr and Ti results in corresponding distribution of Pb displacements. If the point defect concentration is sufficiently small (of the order of or less) so that the longrange order characteristic to the crystal symmetry is not changed, defects mainly contribute to the intensity of the tail areas around the Bragg reflections: diffuse scattering is increased, whereas elastic scattering is decreased. The effect is most clearly seen at small spacing area (see, e.g., ref. [17]). Correspondingly, point defects are frequently seen as abnormal ADP’s, which are either anomalously large or unphysically small: even diagonal components can be negative. This issue is addressed in ref. [5]. The negative values can be due to the fact that ADP’s try to model two separate types of disorder, dynamic (thermal vibrations) and static (substitutional and/or positional disorder). Small substitute atomic concentrations are seldom capable of introducing large changes in xray scattering intensities (unless there is a large difference in the numbers of electrons of the atoms) though they may result in detectable changes in neutron scattering. Point defect concentrations of a fraction of at.% may, however, result in observable changes in bond lengths [52].
At large point defect concentration it is better to abandon the concept of wellordered host crystal with point defects and to model the system from the beginning. Fig. 5 schematically illustrates a binary solid solution and the commonly used approximation in which one does not make a distinction between the two types of atoms, but only considers a pseudoatom taken to possess a scattering amplitude formed as a composition weighted average of the scattering amplitudes of the two atoms. Threedimensional periodicity is commonly introduced in a similar manner. Thus, the disorder is averaged away as it is a heavy task to compute the scattering intensity corresponding to a huge unit cell. Below we introduce an approximation which keeps the essential features of the solidsolution, namely the inhomogeneous distribution of two types of atoms and the corresponding variation in bond lengths, corresponding to the case illustrated in Fig. 5(c).
2 Scattering from PZT clusters
A method for modelling the local structure under the constraint that the average structure remains intact was recently developed [53] and applied for addressing the cation displacements as a function of hydrostatic pressure [54]. Central part of the modelling work is the parameterization of the disorder. As in the case of crystals with welldefined space group symmetries there is a requirement to fill the space exactly once. This sets limits to the cells and their mutual connectivity. A straightforward way to do this is given below. Focus is on the construction of a model which takes the disorder into account yet is computationally sufficiently simple. We first describe how a single scattering cluster is constructed. The material to be modelled can entirely be consisted of a single scattering cluster, or as a special case the cluster can be a unit cell in a crystallographical sense.
Next section introduces the parameters used in this paper after which application examples are given.
2.1 Generation of the scattering cluster
Below a construction of a model of a combined substitutional and displacement disorder is given. The influence of the crystal size and shape on the peak profiles is taken into account. Scattering intensity is computed for single clusters. The program was written using the Clanguage and a message passing interface (MPI). Computational platform was provided by the CSC (Finnish IT Center for Science Ltd., administered by the Ministry of Education, Science and Culture).
We first generate a scattering cluster in which each ion is placed on a specific site (described below). Ti and Zr are statistically distributed after which the atomic positions are adjusted. Specifically, each cell in the cluster is either tetragonal or rhombohedral and the ion positions are relaxed accordingly. The relaxed structural parameters are tabulated in Table 2. Atomic positions are adjusted by bondvalencesum method [55].
Initial parameters.
Experimental roomtemperature structural values of PbTiO and PbZrO are used for generating the initial values for the scattering cluster. In this stage, also alternative methods could be used. The lattice parameters for PbTiO are: Å and Å and for PbZrO the lattice parameters are Å. Though PbZrO is orthorhombic, we average the structure to be cubic (space group ). For the computation of the scattering power the positions of the ions are referred to points which give the origins of the cells. In the special case of a scattering volume possessing translational symmetry the points form a crystallographical lattice. Here the focus is on the cases lacking translational symmetry.
The structure of the cells depends on the Zr concentration: at high titanium concentrations the cells are tetragonal, at two phases coexist and at high Zrconcentrations rhombohedral cells are dominant. This feature is embedded in the model as described in flow chart 2.1.
Table 2 lists the sites in tetragonal and rhombohedral cells and Table 3 lists the neutron scattering lengths and coefficients used for the computation of xray scattering factors.
Tetragonal structure  

Atom  
Pb  (Pb)  (Pb)(Pb)  (Pb)(Pb) 
Zr  (Zr)  
Ti  (Ti)  
I  (O)  
II  (O)  
III  (O)  
Rhombohedral structure  
Atom  
Pb  (Pb)  (Pb)(Pb)  (Pb) (Pb) 
Zr  (Zr)  (Zr)  (Zr) 
Ti  (Ti)  (Ti)  (Ti) 
I  (O)  
II  (O)  
III  (O) 
Parameter  Pb  Zr  Ti  O 

9.405  7.16  3.370  5.803  
13.4118  2.06929  1.28070  0.250800  
31.0617  17.8765  9.75950  3.04850  
0.690200  1.27618  7.85080  13.2771  
13.0637  10.9480  7.35580  2.28680  
2.35760  11.9160  0.500000  5.70110  
18.4420  5.41732  1.69910  1.54630  
8.61800  0.117622  35.6338  0.323900  
5.96960  3.65721  1.90210  0.867000  
47.2579  87.6627  116.105  32.9089  
4.075  0.186  0.219  0.049  
8.506  2.245  1.807  0.032 
Flow chart.
Flowchart describing the steps involved in the scattering cluster construction.

The size of the scattering cluster is given in terms of parallel epipeds, or cells. In the present example a rectangular parallel epiped cluster, consisted of cells, is constructed. Each cell is referred to by indices , and .

Atoms are placed into the initial sites according to Table 2. Either Ti or Zr is inserted into each cell with probabilities and , respectively. A random number generator is used for this purpose.

The disorder to be modelled is correlated substitutional and displacement type caused by the random distribution of Ti and Zr ions. For modelling purpose we use three sum functions , and . gives the number of Zr atoms in planes , being . Functions and posses similar meaning.

Cell edges are readjusted so that the compatibility condition is fulfilled. Fig. 6 illustrates the compatibility in twodimensional case.
The cell edges , , and are relaxed either to be rhombohedral or tetragonal:
if and and
else
.
Thus, both rhombohedral and tetragonal cells can coexists within a same cluster. To further simplify the model, the largest rhombohedral cell is determined after which all rhombohedral cells are constrained to have this size. This also ensures that all cells are continuously connected to adjacent cells. A parameter is introduced. This is related to the fact that average valence exceeds 4 in Zrrich PZT . To obtain nominal valence the cell volume is gradually expanded by increasing the value of till atoms can be positioned so that they possess nominal valences.For simplicity, linear relationship between lattice constants and sum functions is assumed. We also assume that the same composition variation causes the same structural variation in the and directions. This also implies that for a single phase cluster the average values of and , respectively denoted as and , are equal. This is because , where is the number of Zr atoms in the cluster. The relationships are no longer true for clusters possessing coexisting tetragonal and rhombohedral cells. By giving up these conditions one could approach more complex cases. We also assume that , and directions are orthogonal.

Cation positions are adjusted to fulfil the nominal bondvalence sums, described below.

Parallel epiped shaped particles are hardly seen in real materials, so an ellipsoidal cut of the cluster is formed. The semiaxes of ellipsoids, including spheres, are chosen to be parallel to the cell edge directions. Rounded cluster has a feature that the subsidiary maxima are diminished, in contrast to the parallel epiped clusters.

The scattering power of the cluster is computed by Eq. (5) for selected directions. Directions are given in terms of , where , and , where , , and are unit vectors parallel to the positive , and axis directions, respectively. The scattering power is multiplied by Lorentz factor (Eq.(8)) in the case of neutrons, and by the factor (Eq.(9)) in the case of xrays, see Appendix I.
Bondvalencesum based atomic position adjustment.
Bondvalencesum (BVS) method [55] is applied to calculate the cation positions with respect to the oxygen polyhedra. Alternative computational techniques, such as using empirical potentials for structure optimization, or experimental techniques, could also be used. In the model Zr and Ti cations are displaced along the axis direction if the cell is tetragonal, otherwise the displacement is along the cell diagonal (the rhombohedral cells). Lead cations, consistently with the experimental observations, are always displaced along the cell diagonal.
As the cells are too tight for Zr (there is no position at which the Zr valence would be +4), the Zr positions are first adjusted so that Zr valences are minimized, after which the Ti positions are adjusted so that the average valence is 4. If the average valence is above 4, the cell size is increased by increasing the parameter , see item 4 in flow chart 2.1.
2.2 Single domain spherical clusters
To address the local structural changes and microstrain spherical clusters with a radius of 18 cell lengths () are generated.
Correlated substitutional and displacement disorder.
PZT is evidently a material in which displacement disorder follows substitutional disorder: volumes with high Zrconcentration have larger cells. The distribution of cations itself is assumed to be random as described in flow chart 2.1. Figures 7 and 8 show the  and cation displacements as a function of .
Model indicates rather small displacement of Zr and Ti cations from the oxygen octahedra centre at and in the vicinity of the MPB region. The model predicts that the ferroelectric polarization is essentially due to the Pb and Ti cation displacements. Above it is essentially the Pbdisplacements which are responsible for polarization in PZT. The increasing width of the cation displacement distribution with increasing indicates larger deviation from the average symmetry: In the case of perfect translational symmetry all  and all cation displacements would be peaked at a single value. Also the number of displaced cations is strongly decreasing with increasing : The cation displacements are centred close to the oxygen octahedron centre when . As Fig. 7 shows, the width of the Pbcation displacements is largest in the twophase region (green data). It is also evident that the Pbion displacements are increasing with increasing . As was discussed in ref. [54], the cation displacements are probably underestimated: If one gives up the constraint that the average valence of cations should be and that the Pbcation should have a valence of and replaces it with a less severe constraint that the sum of the two types of cations should be then the cations can be displaced by a larger amount also at Zrrich areas. This would be compensated by a smaller Pbdisplacements. However, such a computation would require an energy minimization. The BVS values do not necessarily correspond to the energy minimum, though it is reasonable to assume that they are not too far off as the BVS parameters are based on a fit to a vast number of experimental data. It is reasonable to assume that most reported structures correspond to the energy minimum.
Cluster and cell dimensions.
Fig. 9(a) shows the number of rhombohedral cells as a function of . The model is seen to be consistent with the known MPB composition (green shadows in Fig. 9(a)(d)), as the strong increase onset at shows. This implies that relatively straightforward statistical approach (given in flow chart 2.1) is capable of explaining the twophase coexistence. The twophase coexistence within a same cluster affects the lineshapes in the twophase regions, resulting in less accurate lattice parameter values in the MPB region, see Fig. 9(d). Fig. 9(c) plots the cluster dimensions and compares the dimensions to the values obtained through the Scherrer equation:
(7) 
where and are the Bragg peak centre position and wavelength values and is the fullwidthathalfmaximum of the Bragg peak. All angles are given in radians. The constant is characteristic to the spherical shaped particles. The maximum dimension of the cluster is slightly larger than the value obtained from the Scherrer equation. This is partially due to the fact that the cluster is not exactly a sphere but is consisted of pseudocubic cells so that the maximum dimension, given by blue markers in Fig. 9(c) are slightly larger than the diameter of the spherical surface fit to go as close to the cluster exterior as possible. We note that the line broadening due to the particle size effect alone yields an apparent size which is independent of the order of the reflection, while it depends on the reflection order in correlated substitutional and displacement disorder. This is also illustrated in Fig. 10 which plots the particle size, estimated through Eq. (7) for PbTiO and and PZT compositions from 14 reflections. The effect of combined substitutional and displacement disorder is most evident in the sample, seen as a large number of satellite peaks (though there is no simple modulation of a periodic structure, we still use the term satellite to distinguish the subsidiary minima peaks due to the limited cluster size and disorder generated peaks).
Correlated displacement and substitutional disorder is most clearly seen in the and type reflections, whereas the reflections possess nearly symmetric distribution of satellite peaks, see Fig. 12. The remarks found in the case of periodical composition and displacement disorder (see Table 1) are seen in Figs. 11 and 12: The intensities of the satellite peaks belonging to the same pair centred at each Bragg reflection are not identical. The signatures are generic to PZT nanoclusters, though the details depend on the cluster statistics.
2.3 Ellipsoid shaped single domain clusters
Figure 13 shows the xray scattering profiles along three directions for and clusters. Table 4 gives the average lattice parameter and cluster size estimates, obtained through Bragg and Scherrer equations, respectively. The widths correspond to the values estimated from the Scherrer equation, though they are systematically smaller than the cluster dimensions, being about 96 % from the maximum ellipsoid dimension. The error is very similar to the error obtained in the case of spherical clusters for which the Scherrer equation (7) holds. As in the case of spherical clusters the Bragg equation values are rather close to the known average values. Partially the error is due to the reasons discussed in the context of spherical clusters. The impact of the width of the scattering intensity on the average lattice parameter estimation is rather small.
0  3.900000, 3.900300  3.900000, 3.900015  4.150000, 4.150105 

0.514485  4.079813, 4.078121  4.074603, 4.074910  4.090201, 4.092251 
0  93.600000, 89.503502  421.200000, 405.758631  74.700000, 71.396015 
0.514485  97.700121, 93.289568  439.059234, 422.897821  73.485345, 70.458609 
2.4 Domains in spherical clusters
As an example of a complex defect system domain walls and domains in spherical clusters are modelled. Fig. 14 illustrates the structural model. The cases treated below possess complex combination of displacement disorder (when the boundary halves the particle, lefthand side domain can be obtained from the righthand side domain by atomic displacements) and substitutional disorder.
The cluster for simulations is constructed as explained in section 2.1, except that the cation displacements directions are reversed at different domains. In each domain the cation positions are relaxed so that the bondvalence sum criteria is fulfilled. In contrast to many standard models, the two domains are not considered separately; instead the scattering intensity is computed for the entire cluster. The purpose of the simulation is to show that (i) domain structure has a significant impact on the scattering intensity, (ii) the method suits for insitu modelling (e.g., polarization reversal studies), (iii) the impact of different elements can be isolated and (iv) certain reflections can no more be described by a single asymmetric peak, even if they would originally correspond to a single symmetric peak. Two domains types, referred to as T and Rdomains, are considered. In the Tdomain the domain wall is perpendicular to the direction and in the case of the Rdomain the wall is perpendicular to the direction. The presence of T and Rdomains is most evidently revealed as a split of the type reflections. The split itself depends on the domain wall position in the cluster, as is seen from Figs. 15 and 16. In the case of Rdomain and reflection also the scattering intensity in the tail regions of the Bragg peak is considerable, consistently with the known strong diffuse scattering observed in many Pbperovskites. As reflection indicates, see Fig. 16, the strong diffuse scattering is due to the PbPb scattering and is easily seen when the particle is divided into domains.
The most significant contribution to the split is from Pbions, which are displaced from the corners to a cell diagonal direction. The parallel epipeds are not identical and thus there is a variation in the Pbdisplacement directions within a domain. In Tdomains the type reflections are left intact by the domain wall, whereas the intensity depends strongly on the domain wall position in the case of the and reflections. In Rdomains all reflections intensities and line shapes strongly depend on the domain wall position.
Applications.
By adjusting domain wall position and orientation and cation displacement direction different cases (e.g., 71 domain walls) can be constructed in a straightforward manner. Domain wall motion takes place in polarization () switching. Also structural changes occurring in the domain wall region can be implemented into the model by constructing a domain wall with a finite thickness. Tdomain boundary addressed above fulﬁls the electrical boundary condition, , so that the boundary is not charged. The mechanical compatibility condition is satisfied as there is no abrupt change in the average lattice parameters across the boundary (i.e., the boundary is stressfree). Also more complex cases, such as the formation of impurity phase or cases not fulfilling the electrical and mechanical boundary conditions (headtohead Rdomain is an example) can be constructed. A plausible application would be a modelling of insitu measurements of domain walls in ferroelectric materials. Due to the high brightness xray synchrotron radiation based experimental techniques are evident methods for addressing timedependent phenomena. As an example, the nonlinear effects in the coupling of polarization with elastic strain and the initial stage of polarization switching were addressed in refs. [56, 57]. In these studies capacitors containing 35 nm thick epitaxial Pb(ZrTi)O ferroelectric thin films were studied by timeresolved xray microdiffraction technique in which highelectric field up to several hundred MV/m pulses were synchronized with the synchrotron xray pulses. Laboratory scale measurements can also be used to address the in situ structural changes in thin films. An example is given in ref. [58] which reports the changes in lattice parameter (chemical expansivity) and its further use for quantifying oxygen reduction reaction processes and vacancy concentration changes in LaSrCoO thin films under chemical and voltage stimuli.
An example of the use of an laboratory xray diffractometer to address timedependent ferroelectric domain reversal is given in ref. [59], where the changes in the volume fractions of the 90 domains parallel to the electric field direction were calculated from the intensities of the diffraction peaks.
Also magnetic scattering can be treated in a manner analogous to the nuclei scattering. In magnetic materials the domain boundary region often have spatially large extent (e.g., Néel and Bloch walls) also in a direction perpendicular to the domain wall. In multiferroic materials the ferroic properties are not necessarily well coupled, and thus it is crucial to be able to isolate different contributions to the scattering intensity. For instance, magnetization reversal may not be accompanied by apparent changes in nuclei arrangements.
All above disorder cases can be combined to model the often complex structures of multilayer thin films. Especially in structures formed from very thin layers (say, of the order of few nanometers) the scattering intensity should not be treated as originating from independent layers. A better way is to model the entire structure and to impose the required boundary conditions for interfaces.
Conclusions
Applications of xray and neutron scattering techniques for analysing defects were reviewed. Focus was on the common approaches applied for modelling defects. A method for analysing scattering data collected on nanoparticles was described with necessary compatibility conditions. The method is capable of isolating different contributions to the scattering intensity, such as element specific scattering, microstrain in solidsolutions, particle size and shape effects and domains: structural disorder is not averaged away. Scattering measurements provide information without destructive sample preparation. A case study on leadzirconatetitanate nanoparticles was given. Potential applications include nanoparticles, disordered nanoparticles and insitu studies of structural changes in ferroic materials.
Acknowledgments
CSC  IT Center for Science Ltd. is acknowledged for providing a computing environment.
References
 [1] A. Guinier. Xray Diffraction in Crystals, Imperfect Crystals, and Amorphous Bodies. Dover Publications, Inc., New York (1994).
 [2] International Tables for Crystallography, Volume C: Mathematical, Physical and Chemical Tables. E. Prince, (Ed.), Kluwer Academic Publishers, Dordrecht, The Netherlands (2002).
 [3] T. R. Welberry. Diffuse XRay Scattering and Models of Disorder. Oxford University Press Inc., New York (2010).
 [4] M. M. Woolfson. An introduction to Xray crystallography. Cambridge University Press, Cambridge (1997).
 [5] W. Massa. Crystal Structure Determination. SpringerVerlag, Berlin (2000).
 [6] International Tables for Crystallography A: Spacegroup Symmetry. Th. Hahn, (Ed.), Kluwer Academic Publishers, Dordrecht (2005).
 [7] H. M. Rietveld. J. Appl. Crystallogr. 65, 2 (1969).
 [8] The Rietveld Method. R. A. Young, (Ed.), Oxford University Press Inc., New York (1996).
 [9] M. Birkholz. Thin Film Analysis by XRay Scattering. WileyVCH Verlag GmbH & Co. KGaA, Weinheim (2006).
 [10] P. W. Stephens. J. Appl. Crystallogr. 32 281 (1999).
 [11] A. C. Larson and R. B. Von Dreele, General Structure Analysis System (LANSCE MSH805, Los Alamos National Laboratory, Los Alamos, NM, 2000).
 [12] M. D. Biegalski, D. D. Fong, J. A. Eastman, P. H. Fuoss, S. K. Streiffer, T. Heeg, J. Schubert, W. Tian, C. T. Nelson, X. Q. Pan, M. E. Hawley, M. Bernhagen, P. Reiche, R. Uecker, S. TrolierMcKinstry, and D. G. Schlom. J. Appl. Phys. 104, 114109 (2008).
 [13] R. Phillips. Crystals, Defects and Microstructures: Modeling Across Scales. Cambridge University Press, Cambridge (2001).
 [14] A. Kelly and K. M. Knowles. Crystallography and Crystal Defects. Second Edition. A John Wiley & Sons, Ltd., Publication, Malaysia (2012)
 [15] K. W. Kolasinski. Surface Science: Foundations of Catalysis and Nanoscience. Wiley, Chippenham (2008)
 [16] H. Lüth. Solid Surfaces, Interfaces and Thin Films. SpringerVerlag, Berlin (2001).
 [17] L. S. Zevin and G. Kimmel. Quantitative Xray Diffractometry. SpringerVerlag, New York (1995).
 [18] B. A. Strukov and A. P. Levanyuk. Ferroelectric Phenomena in Crystals: Physical Foundations. SpringerVerlag, Berlin (1998).
 [19] N. W. Thomas. Acta Cryst. A66, 64 (2010).
 [20] M. T. Dove. Structure and Dynamics: An atomic view of materials. Oxford University Press, Inc., New York (2003).
 [21] Jaffe, B.; Cook, W. R. and Jaffe, H. (1971). Piezoelectric Ceramics, Academic Press, New York.
 [22] A basic search by a keyword ”lead zirconate titanate” in a topic field and timespan from 2012 to 2014 gave 1114 records. When the time span was expanded to cover years 19602014, the number of records was 6874, indicating that a large number of the research records are recent.
 [23] M. E. Lines and A. M. Glass, Principles and Applications of Ferroelectric and Related Materials (Clarendon Press, Oxford, 2001).
 [24] B. Kocsis, J. M. PerezMato, E. S. Tasci, G. de la Flor and M. I. Aroyo. J. Appl. Cryst. 47, 1165 (2014).
 [25] J. Frantti. J. Phys. Chem. B, 112, 6521 (2008).
 [26] W. Cao and L. E. Cross. Phys. Rev. B 47, 4825 (1993).
 [27] B. Noheda, D. E. Cox, G. Shirane, J. A. Gonzalo, L. E. Cross, and S.E. Park, Appl. Phys. Lett. 74, 2059 (1999).
 [28] H. Yokota, N. Zhang, A. E. Taylor, P. A. Thomas, and A. M. Glazer, Phys. Rev. B 80, 104109 (2009).
 [29] J. Frantti, J. Lappalainen, S. Eriksson, V. Lantto, S. Nishio, M. Kakihana, S. Ivanov, and H. Rundöf, Jpn. J. Appl. Phys., Part 1 39, 5697 (2000); J. Frantti, S. Eriksson, S. Hull, S. Ivanov, V. Lantto, J. Lappalainen, and M. Kakihana, J. Eur. Ceram. Soc. 24, 1141 (2004); J. Frantti, S. Eriksson, S. Hull, V. Lantto, H. Rundöf, and M. Kakihana, J. Phys.: Condens. Matter 15, 6031 (2003); J. Frantti, Y. Fujioka, J. Zhang, S. Wang, S. C. Vogel, R.M. Nieminen, A. M. Asiri, Y. Zhao, A. Y. Obaid, and I. A. Mkhalid, J. Appl. Phys. 112, 014104 (2012).
 [30] A. M. Glazer, P. A. Thomas, K. Z. BabaKishi, G. K. H. Pang, and C. W. Tai, Phys. Rev. B 70, 184123 (2004).
 [31] D. Vanderbilt and M. H. Cohen, Phys. Rev. B 63, 094108 (2001).
 [32] A. J. Bell and E. Furman, Jpn. J. Appl. Phys., Part 1 42, 7418 (2003).
 [33] N. Lampis, P. Sciau and A. Geddo Lehmann. J. Phys.: Condens. Matter 11, 3489 (1999).
 [34] V. Bonny and M. Bonin and P. Sciau and K. J. Schenk and G. Chapuis.Solid State Communications. 102, 347 (1997).
 [35] H. Fu and R. E. Cohen. Nature. 403, 281 (2000).
 [36] J. Frantti, Y. Fujioka and R. M. Nieminen. J. Phys.: Condens. Matter 20, 472203 (2008).
 [37] N. Zhang, H. Yokota, A. M. Glazer, Z. Ren, D. A. Keen, D. S. Keeble, P. A. Thomas and Z.G Ye. Nature Communications, DOI:10.1038/ncomms6231
 [38] D. Phelan, X. Long, Y. Xie, Z.G. Ye, A. M. Glazer, H. Yokota, P. A. Thomas, and P. M. Gehring, Phys. Rev. Lett. 105, 207601 (2010).
 [39] S. Gorfman, D. S. Keeble, A. M. Glazer, X. Long, Y. Xie, Z.G. Ye, S. Collins, and P. A. Thomas, Phys. Rev. B 84, 020102 (2011).
 [40] J. Frantti, V. Lantto, S. Nishio and M. Kakihana. Phys. Rev. 59, 12 (1999)
 [41] J. Frantti, Y. Fujioka, A. Puretzky, Y. Xie, Z.G. Ye and A. M. Glazer. Journal of Applied Physics, 113, 174104 (2013).
 [42] J. Frantti, Y. Fujioka, A. Puretzky, Y. Xie, Z.G. Ye., C. Parish and A. M. Glazer. J. Phys.: Condens. Matter, 27, 025901 (2015).
 [43] D. L. Corker, A. M. Glazer, R. W. Whatmore, A. Stallard, and F. Fauth, J. Phys.: Condens. Matter 10, 6251 (1998).
 [44] J. Ricote, D. L. Corker, R. W. Whatmore, S. A. Impey, A. M. Glazer, J. Dec, and K. Roleder, J. Phys.: Condens. Matter 10, 1767 (1998).
 [45] C. Muller, J.L. Baudour, V. Madigou, F. Bouree, J.M. Kiat, C. Favotto, and M. Roubin, Acta Crystallogr., Sect. B: Struct. Sci. 55, 8 (1999).
 [46] J. Frantti, S. Ivanov, S. Eriksson H. Rundöf, V. Lantto, J. Lappalainen, and M. Kakihana, Phys. Rev. B 66, 064108 (2002).
 [47] D. J. Goossens, Acc. Chem. Res. 46, 2597 (2013).
 [48] N. de Mathant, E. Hussont, S. G. Calvarint, J. R. Gavarris, A. W. Hewat and A. Morel. J. Phys.:Condens. Matter 3, 8159 (1991).
 [49] C. Malibert, B. Dkhil, J. M. Kiat, D. Durand, J. F. Bérar, and A. Spasojevicde Biré. J. Phys.: Condens. Matter 9, 7485 (1997) .
 [50] B. Dkhil, J. M. Kiat, G. Calvarin, G. Baldinozzi, S. B. Vakhrushev, and E. Suard. Phys. Rev. B 65, 024104 (2001).
 [51] M. Paściak, T.R. Welberry, A.P. Heerdegen, V. Laguta, T. Ostapchuk, S. Leonid and J. Hlinka. Phase Transitions, DOI: 10.1080/01411594.2014.981266.
 [52] J. Frantti and V. Lantto. Phys. Rev. B. 221 56 (1997).
 [53] J. Frantti and Y. Fujioka, “Method and system for analysing data obtained using scattering measurements from disordered material,” WO2011/018554  EP2464962 B1 (2014).
 [54] J. Frantti, Y. Fujioka, J. Zhang, J. Zhu, S. C. Vogel, and Y. Zhao. Rev. Sci. Instrum. 85, 083901 (2014).
 [55] I. D. Brown. The Chemical Bond in Inorganic Chemistry: The Bond Valence Model. Oxford University Press Inc., New York (2009).
 [56] A. Grigoriev, R. Sichel, H. N. Lee, E. C. Landahl, B. Adams, E. M. Dufresne, and P. G. Evans. Phys. Rev. Lett., 100, 027604 (2008).
 [57] A. Grigoriev, R. Sichel, J. Y. Jo, S. Choudhury, LQ. Chen, H. N. Lee, E. C. Landahl, B. Adams, E. M. Dufresne, and P. G. Evans. Phys. Rev. B, 80, 014110 (2009).
 [58] M. D. Biegalski, E. Crumlin, A. Belianinov, E. Mutoro, Y. ShaoHorn, and S. V. Kalinin. Appl. Phys. Lett. 104, 161910 (2014).
 [59] A. Pramanick and J.L. Jones. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control. 56, 1546, (2009).
3 Instrument specific corrections
Scattering geometry and the nature of the incident radiation (for instance, xrays or polarized or unpolarized neutron beams) affects the observed scattering intensity. The simplest correction is the diffraction geometry dependent Lorentz factor , a trigonometric factor due to the fact that different reflections stay different times at the reflection position. In the case of a constant angular speed and single crystal sample is (see, e.g. ref. [5]). For powders the prevailing measurement geometry is the BraggBrentano, or geometry, in which the sample is rotated at constant speed. If one assumes that the incident beam probes large set of randomly oriented crystals a further geometrical factor is required as only a fraction of the crystals are in a reflection position with respect to the detector. The Lorentz factor is the same for xrays and neutrons. Thus, the neutron scattering intensities computed through Eq.(5) were multiplied by a factor
(8) 
In the case of xrays one needs a further correction to take the polarization of xrays into account. Typically xrays emerging from a xray tube are unpolarized, and after reflecting from the diffracting plane the component of the electric field parallel to the plane is not attenuated, whereas the electric field component perpendicular to the diffracting plane is attenuated by a factor of , so that the total intensity is reduced by a polarization factor . Monochromators affect the polarization factor and modified expressions for are required: when a monochromator for an incident beam is used [5]. is a constant typically close to unity.
The effect of Lorentz and polarization corrections is typically put together and expressed as the factor, which multiplies the scattering power or, in the case of crystals, . In this study the xray scattering intensities were multiplied by a factor
(9) 
4 Line shapes in diffraction experiments
Instrument also affects the observed lineshape. When the scattering power is computed by equation (5) as a function of for a chosen direction of , the cluster size and shape dictates the linewidths and, if the cluster is small, is almost entirely responsible for the broadening. In contrast, when the intensity for a specific reflection is computed as , a specific peak profile must be assumed. Typically, one needs to consider the sample size and shape, possible defects causing symmetric and asymmetric broadening, and the instrumental contribution as is done in conventional Rietveld refinement. The functions describing different contributions are convoluted, which is computationally heavy task. Correspondingly, a large number of different types of profile functions have been introduced to model different cases.