Monte Carlo Modelling of Single-Crystal Diffuse Scattering from Intermetallics

Abstract

Single-crystal diffuse scattering (SCDS) reveals detailed structural insights into materials. In particular, it is sensitive to two-body correlations, whereas traditional Bragg peak-based methods are sensitive to single-body correlations. This means that diffuse scattering is sensitive to ordering that persists for just a few unit cells: nanoscale order, sometimes referred to as “local structure”, which is often crucial for understanding a material and its function. Metals and alloys were early candidates for SCDS studies because of the availability of large single crystals. While great progress has been made in areas like ab initio modelling and molecular dynamics, a place remains for Monte Carlo modelling of model crystals because of its ability to model very large systems; important when correlations are relatively long (though still finite) in range. This paper briefly outlines, and gives examples of, some Monte Carlo methods appropriate for the modelling of SCDS from metallic compounds, and considers data collection as well as analysis. Even if the interest in the material is driven primarily by magnetism or transport behaviour, an understanding of the local structure can underpin such studies and give an indication of nanoscale inhomogeneity.

1. Introduction

Short-range order (SRO) is present in almost all families of crystalline compounds, from metals to proteins [1,2,3,4,5,6,7,8]. SRO can influence electrical, magnetic and most other physical properties, including ferroelectricity, superconductivity and multiferroic behaviour.

SRO manifests in the diffuse scattering, the coherent scattered intensity which is not localised on the reciprocal lattice; in other words, it is found throughout reciprocal space, not just on the Bragg reflections at integer $hkl$. Thus, to best investigate the diffuse scattering it is necessary to survey a large region (area or volume) of reciprocal space with low noise and high dynamic range. This is not a trivial exercise, and much effort has gone into data collection and reduction [6,9,10,11].

Data are typically presented as reciprocal space cuts or sections, which essentially plot diffracted intensity as a function of position in reciprocal space.

Metals were an early test-bed for ways of modelling SRO, in particular chemical SRO as modelled by, for example, Cowley SRO parameters [12,13,14,15]. Cowley realised that Fourier transforming the diffuse intensity could give atomic pair correlations when the scattering admitted a direct interpretation, for example when looking at a diffuse peak that would sharpen to a Bragg spot on going through a phase transition. Warren and co-workers showed how the atomic size effect (the dependence of interatomic spacing on species, most simply conceptualised as thinking about atoms as being of different radii) caused asymmetries in the scattering [16]. When the system is relatively simple, sometimes an analytical form can be found to yield the distribution of scattering.

If the underlying crystallography is simple, it may be possible to use an essentially analytic analysis, as for example can be obtained by expanding the diffraction equations [17] and using conditional probabilities to express the various terms. These probabilities can then be adjusted and the expected scattering calculated.

However, in more complex cases, in particular systems containing many atomic species and/or in which the atoms form into clusters with their own structure factors that then conflate with the scattering from the defects and the local ordering, it is often difficult or impossible to interpret the scattering directly or to meaningfully invert it to get the real space structures. These, and cases where we must allow for displacive relaxation around defects, require a more model-based approach. When contrast between scatterers is weak (atoms nearby on the periodic table will have very similar X-ray scattering factors), it may be necessary to use neutron and/or X-ray single-crystal diffuse scattering (SCDS) data. Neutron diffraction requires larger crystals, which may be difficult to obtain, so it may be that X-ray SCDS is coupled with neutron pair distribution function analysis (PDF; [18,19,20]), obtained from polycrystalline specimens.

A wide range of local structures have been observed in metallic compounds, from classic examples like chemical substitution and resulting clustering or anti-clustering in alloys, through to subtle phenomena related to the atomic size effect and even the rotation of large motifs, such as the cages of atoms seen in complex intermetallics [21]. For relatively simple systems, recent advances allow almost direct interpretation of the diffuse scattering, while developments in detailed calculation methods, like density functional theory and molecular dynamics, allow direct calculation of low energy short-ranged order configurations when not too many atoms are required [22,23,24,25].

However, when many atoms are involved and the correlation lengths encompass many unit cells, the number of atoms involved is beyond the scope of such methods. Then, the ability to model a crystal of >${10}^5$ atoms becomes useful. Methods like 3D-ΔPDF [26] offer what are almost “direct methods” for such systems and are currently a fascinating field of development. The reverse Monte Carlo (RMC) approach [19,27,28] offers a means to directly fit the diffuse scattering data, but can be limited in the size of simulation that can be implemented because of the way in which a single atomic move must have a significant effect on the goodness of fit of the model.

Thus, at this time, the most flexible approach remains the forward Monte Carlo (MC), though it has its own weaknesses, in particular one must posit the nature of the disorder and then find a means of introducing that disordered structure into the model, before calculating the Fourier transform of the model and testing the theory. The process can be slow; models are difficult to optimise; and knowing what to include in the model (what forms of disorder and how to induce them) requires considerable insight. Further, since disorder can take on so many forms, it is often necessary to write bespoke computer code to tackle a given problem, something which is time consuming and not conducive to broad acceptance of the technique.

This paper aims to very briefly look at Monte Carlo analysis of diffuse scattering, particularly as it pertains to metallic materials, alloys and the like. The fascinating field of quasicrystals, many of which are metallic, will not be covered. This field has been surveyed in a range of detailed and high quality presentations, which need not be repeated here [29,30,31].

2. Data Collection

The experiments considered here use large slices of reciprocal space, rather than collecting intensity at a few key scattering vectors. This allows elucidation of SRO that is anisotropic or only affects small regions of reciprocal space. Similarly, the use of pair distribution function and powder diffraction is not discussed, though both are very important techniques [3,18,19,32].

The quality and quantity of data required depends, of course, on the experiment being undertaken. Ideally, the different scatterers will have well-differentiated cross-sections for the radiation being used. If the disorder is anisotropic, then data that extend in three dimensions are desirable. If local ordering is only significant in, say, the $ab$ plane, then collection of the $hk0$ section of reciprocal space may be sufficient. If quantitative comparison of the calculated SCDS with the observed is desired [33], the observed data must show low noise, few artefacts, and a background that can be removed either by subtraction of “blank” runs or some other method, like fitting a function to it. For qualitative comparison with calculations, showing whether features are present or not, for example, noisier data may be acceptable, and the less quantitative results of electron diffraction are also useful. Analysis of SCDS is often limited by the data that can be obtained, but as long as features in the scattering can be identified as “real”, then some insight can be gained.

2.1. X-ray

Assuming that the X-ray source is a constant wavelength, monochromated source, volumes of diffuse scattering are collected by rotating a sample in front of an area detector. Earlier work often made use of a line counter [34], but the modern prevalence of area detectors has rendered this approach largely redundant.

The main variation is in the choice of detectors. In particular, while much important data collection has made use of image plates [31,35,36,37,38,39,40,41,42], the use of electronic counters that can provide a high dynamic range has become possible [43,44,45]. These have a much improved duty-cycle. Experiments with image plates at synchrotrons, where beams are very intense, can follow an exposure of a few seconds, rarely more than 30 s, with a readout time of a minute or more, which is not good use of the intense and expensive beam.

Figure 1 presents a generic schematic diagram of a constant wavelength experiment. The main parameters include the sample to detector distance, the wavelength and whether the beam path is enclosed in a vacuum or He-filled vessel, which reduces noise, or is through air, which tends to result in intense forward scattering that requires careful correction and collections of “blank” runs, which can then be subtracted from the data. Other corrections may be required depending on the nature of the detector and the stability of the beam and the nature of the beam. If a laboratory source is being used, the compromise between intensity and quality of monochromation can result in the beam possessing a white component, which is much weaker than the characteristic radiation, but nevertheless results in a radial streak through the Bragg peaks, because of the long exposures required to reveal the diffuse scattering. Other artefacts that would not be apparent in an experiment using shorter exposure times may also be revealed. These include X-rays that pass through the image plate and scatter off components of the detector and re-enter the image plate from behind (this was discovered when the shadows of the image plate mounting screws were projected onto the detector(!)), as well as resolution streaks, discussed in Figure 15 of [46].

The high intensities at a synchrotron can cause problems when the area detector intercepts a Bragg reflection; depending on the design of the detector, a wire or a pixel can become saturated. In CCD devices, charge can spill over and contaminate surrounding pixels (deep depletion devices overcome this somewhat); in a wire detector, a bright spot anywhere on the wire may force the removal from the data set of all “pixels” measured by that wire [11].

Other issues include ghosting, when a pixel value on a measurement is partly influenced by the previous measurement. This can happen in image plates, where a very highly exposed pixel may not be fully “reset” by the readout, and thus, its value on the next exposure is not correct.

Traditionally, flat reciprocal space cuts have been reconstructed from the curved sections collected in an experiment such as that in Figure 1. Flat sections generally admit to easier visual interpretation, as the normal is everywhere the same and corresponds to a particular reciprocal space direction. However, from a computational point of view there is little difference between calculating the scattering in a flat or curved section. Further, at high X-ray energies the radius of the Ewald sphere is so large that each exposure is almost a flat section in reciprocal space anyway. In such cases, it is sensible to align the crystal carefully, such that useful data can be obtained with relatively few exposures. This leads to the ability to do parametric studies of diffuse scattering, which is an area under-exploited at this time.

Figure 1. A schematic diagram of a diffuse scattering collection using a 2D detector. The sample angle is ω; incoming X-rays are of known wavelength, λ; and the scattering angle is as usual $2\theta$; but because we wish to transform the detector coordinates into $hkl$’s, we work with x and y coordinates on the detector. During a single exposure, the sample is typically rocked through an angle $d\omega \sim {0.25}^{\circ }$, then ω is incremented by $d\omega$ and the measurement repeated. After $180/d\omega$ such exposures, enough data points have been collected to reconstruct most of reciprocal space [10] out of the the maximum value, which is given by the radius of the detector, the sample-detector distance, and λ.

Figure 1. A schematic diagram of a diffuse scattering collection using a 2D detector. The sample angle is ω; incoming X-rays are of known wavelength, λ; and the scattering angle is as usual $2\theta$; but because we wish to transform the detector coordinates into $hkl$’s, we work with x and y coordinates on the detector. During a single exposure, the sample is typically rocked through an angle $d\omega \sim {0.25}^{\circ }$, then ω is incremented by $d\omega$ and the measurement repeated. After $180/d\omega$ such exposures, enough data points have been collected to reconstruct most of reciprocal space [10] out of the the maximum value, which is given by the radius of the detector, the sample-detector distance, and λ.

It may be noted that static and dynamic displacements cannot be distinguished with an X-ray experiment because, compared to the high energies of the X-rays, all atomic motions are of very low energy (seem very “slow”) and are seen as “static”; this is an area where neutrons may be preferable.

2.2. Neutron

Neutron diffraction comes in essentially two varieties: constant wavelength and time of flight. The latter is most commonly found at a spallation neutron source, while the former is found at reactor sources or a steady-state spallation source, like SINQ, the Swiss Spallation Neutron Source.

A constant wavelength experiment essentially uses the same configuration as for the X-ray case (Figure 1). A typical example is the Wombat instrument at the Bragg Institute at the Australian Nuclear Science and Technology Organisation (ANSTO) [47,48]. This instrument uses a two-dimensional detector to collect a sort of “cake slice” of diffraction space, such that data collected at multiple sample angles can be combined to give a volume from which sections can be extracted. Such an instrument does not select for neutron energy, so scattering from dynamic effects like phonons overlaps with that from static structures like chemical short-range order. This is much as for X-rays, except that the neutron energy is much lower, and inelastic effects may change the neutron wavelength substantially, which has the effect of “moving” the scattered beam around on the detector and, thus, shifting the inelastic scattered intensity to different positions in the reciprocal space map. Such effects can in some cases be interpreted usefully [49]. They do lead to a reduction of the symmetry of the pattern and may limit the ability to quantitatively model the scattering. If diffuse scattering is measured using an instrument that can select for neutron energy, for example a chopper spectrometer, then static can be separated from dynamic, although that depends on the energy resolution of the instrument; quasi-elastic scattering may be binned in with the “strictly elastic” scattering.

At a spallation neutron source, the time structure of the pulse collapses an entire diffraction pattern into a single pixel on a detector, meaning that such instruments, for example SXD (single crystal diffractometer) at ISIS [50,51] and TOPAZ at the Spallation Neutron Source [52], collect very large volumes of reciprocal space with a single sample setting. Rotating the sample leads to rapidly scanning a large volume, generally much larger than that accessible at a constant wavelength source. On the other hand, instrument resolution can vary dramatically from forward- to back-scattering detectors, and since the experiment is essentially imaging reciprocal space, this can affect the interpretability of some patterns. Further, such instruments are often “open” in geometry, without collimation between sample and detector. Thus, they effectively image the sample onto the detector, meaning that anisotropic sample shape can lead to odd-shaped features. This is not an issue when the feature is to be integrated up to get an intensity for conventional Bragg analysis, but when reciprocal space maps are being looked at, it can have an effect.

It is possible to use energy discrimination on spallation instruments [49,53], and again, this yields the possibility of separating dynamic from static effects.

Whether constant wavelength or spallation, polarisation analysis can be used to separate magnetic from structural diffuse scattering [54,55,56,57].

3. Basic Principles of Monte Carlo Modelling of SRO

This topic is dealt with in great detail elsewhere [17,58,59,60], so a simple outline will suffice; Figure 2 summarises the process.

Figure 2. The overall MC modelling procedure. The flow chart illustrated in Figure 3 is an expansion of the box labelled “Do a Monte Carlo simulation to equilibrate the structure”. This diagram assumes a least squares procedure based on calculating a ${\chi }^2$ statistic for the model (or perhaps a kind of R-factor [61,62]); but often, the comparison will be done heuristically by the investigator, and the results will be more qualitative. The initial model is based on the average structure from Bragg data.

Figure 2. The overall MC modelling procedure. The flow chart illustrated in Figure 3 is an expansion of the box labelled “Do a Monte Carlo simulation to equilibrate the structure”. This diagram assumes a least squares procedure based on calculating a ${\chi }^2$ statistic for the model (or perhaps a kind of R-factor [61,62]); but often, the comparison will be done heuristically by the investigator, and the results will be more qualitative. The initial model is based on the average structure from Bragg data.

At its simplest, the type of MC modelling considered here has just a few steps.

• Decide on a starting configuration for the model. This usually means creating (in a computer) a $M\times N\times P$ array of unit cells, typically 32 on a side, and populating it with atoms based on the average structure determined by conventional studies.
• Choose some interactions between atoms. To set up chemical SRO when there are two species, a typical interaction is a Ising-like potential for the energy associated with the occupancy of site i, $E_{\mathrm{occ}}^i$:

  $$E_{\mathrm{occ}}^i=-J_{\mathrm{NN}}\sum _{\mathrm{NN}}{S}_i{S}_j$$

  (1)

  where j indexes nearest neighbours and the sign of J determines whether a positive or negative nearest neighbour occupancy correlation, $C_{\mathrm{NN}}$, is energetically favourable. Further, such terms may be present for more distant neighbours. $S_j=\pm 1$.
  If it is displacements that are of interest, the simplest choice is to connect atoms with Hooke’s law springs The program ZMC [63] is designed to induce correlations amongst atomic and/or molecular displacements by causing the atoms to interact with surrounding atoms via Hooke’s law springs of the form:

  $$E_{\mathrm{inter}}=\sum _{\mathrm{cv}}{F}_i{(d_i-d_{0i}(1+{ϵ}_i))}^2$$

  (2)

  where $d_i$ is the length of vector i connecting atoms, $d_{0i}$ is its equilibrium length and $F_i$ is its force constant. The sum is over all contact vectors (cv). ${ϵ}_i$ is the “size-effect” term, which allows that the equilibrium length required for the calculation may not be the average length as determined from Bragg scattering; this is particularly likely to be the case in occupationally-disordered materials, where the Bragg-refined intermolecular distance is in fact an average over several different distances resulting from differing atomic or molecular species (or vacancies).
• The actual MC part happens as follows (summarised in Figure 3). An atom is chosen at random, and its energy is calculated. Its configuration is changed, and the energy calculation repeated. The new configuration is kept or rejected based on a simple criterion: if new energy is lower, it is kept, and it may be kept if new energy is higher, with some probability based on simulation “temperature”.
• Note that the configuration may be changed by adding small random variations to an atom’s variables (e.g., moving it slightly) or by swapping the variables of one site with those of another. Swapping is particularly useful as a means of maintaining an initial population of displacements or chemical species, while inducing correlations within that population.
• Once every site has been visited, on average, some large number of times, which could be ten, hundreds or thousands, depending on the needs of the simulation, the simulation is complete, and the atomic coordinates are read out.
• A Fourier transform program DIFFUSE [64] then calculates the diffuse scattering for comparison with the experiment.
• It is possible to embed this process within a procedure that automatically modifies the interaction parameters to try to improve the fit between calculated and observed diffuse scattering, although often useful results can be obtained by qualitative comparison, which can be used to reveal key aspects of the local order without comprehensive fitting.

Figure 3. A simple representation of a single forward MC step; a molecule may be a single atom or a more complex motif.

Figure 3. A simple representation of a single forward MC step; a molecule may be a single atom or a more complex motif.

The advantage of this approach is that the “energy” can be anything as long as it is quick to calculate. It may be relatively realistic or quite abstract, whatever suits the problem. However, knowing what disorder is present and then how to parameterise the interactions to induce it is not simple.

4. A Model System

In this section a model system, CePdSb, is considered from the point of view of inducing a range of local orderings and their resulting diffraction effects. No comparison with the observations is made, as we are looking simply to show how the disorder is modelled and some of the forms it can take.

CePdSb and related compounds form a family demonstrating a wide range of unusual magnetic phenomena, including the Kondo effect, heavy fermion behaviour and half-metal behaviour [65,66,67,68,69,70,71].

CePdSb itself shows a crystal structure in which the Pd and Sb lie on ordered sub-lattices at coordinates $(1/3,2/3,0.4684)$ and $(2/3,1/3,0.516)$ [66], with space-group $P{6}_{3}mc$ and lattice parameters approximately $a=4.935$Å and $c=7.890$Å. This is different from an earlier structure in which the z coordinates of both Pd and Sb were taken as 0.5 [67] and the Pd and Sb were considered to be randomly mixed across the Pd/Sb sites.

The Ce atoms lie at $(0,0,1/3)$ 2a positions, forming chains along the c axis. The structure is represented in Figure 4. For the purposes of demonstrating various diffraction effects, we will explore what happens when the Ce2 site ($z=2/3$) is occupied by approximately 67% Ce atoms and 33% vacancies.

Figure 4. A schematic diagram of the structure of CePdSb, showing the Ce layers and the Pb/Sb layers, the latter of which are not flat, but “puckered” [72].

Figure 4. A schematic diagram of the structure of CePdSb, showing the Ce layers and the Pb/Sb layers, the latter of which are not flat, but “puckered” [72].

If we take the average structure of CePdSb [66] and calculate the diffuse scattering, we of course see nothing of interest, as there are no short-range correlations. However, we may, for example, connect atoms with Hooke’s law springs (Equation (2)) and run a simulation. Figure 5 shows three sections through the diffuse scattering from CePdSb. The first row of images comes from a model in which there are no Ce vacancies and the atoms are connected by Hooke’s law springs. The interactions induce streaks, most apparent in the $hk5.5$ layer. The second row shows the same cuts, but for a model in which there are 33% vacancies on the Ce2 layers, and they are forced to cluster. In the third row, the vacancies anti-cluster, and we can see in Figure 5i that this induces sharp spots in the half-layer, $hk5.5$, where previously, there were only streaks (the streaks are in fact sections through planes of scattering that can also be seen in the $hk5$ layer, though being less obvious due to the bright spots). We can also see that the clustering has little effect on the $hk5$ layer, while in $hk0$ it causes the spots that are present in hexagonal motifs around each Bragg peak (one hexagon is noted by white lines in Figure 5a) to extend closer to the origin. These spots actually come from the fact that the Pd and Sb atoms are not on idealised positions, such as $(1/3,2/3,1/2)$. When the vacancies cluster, we have large regions of the crystal where the scattering from the Ce2 layer is absent (effectively, these are like crystallites of composition Ce${}_{0.5}$PdSb), giving different cancellation and allowing the spots to persist. When the vacancies anti-cluster (in the third row), the average scattering from Ce2 is preserved on the local scale, as well, and the cancellation is more like that seen in Figure 5a, though not identical.

Figure 5. Slices of calculated diffuse scattering from different models of CePdSb. Row 1: no vacancies. Row 2: 33% vacancies on Ce2 site, clustering. Row 3: 33% vacancies on Ce2 site, anti-clustering. $hk5.5$ layers are normalised more brightly to bring out the details. For details, see the text. h and k axes noted on (a) to indicate directions.

Figure 5. Slices of calculated diffuse scattering from different models of CePdSb. Row 1: no vacancies. Row 2: 33% vacancies on Ce2 site, clustering. Row 3: 33% vacancies on Ce2 site, anti-clustering. $hk5.5$ layers are normalised more brightly to bring out the details. For details, see the text. h and k axes noted on (a) to indicate directions.

In Figure 6, in rows 1 and 2, the displacive and occupancy effects are combined: the average distance atom-vacancy has been made 20% bigger than the average, while atom-atom is 10% smaller and vacancy-vacancy is 40% bigger. This is to mimic the effect sometimes seen where atoms move away from vacancies due to the lack of a bond, rather than moving into the gap. Row 1 is the model where the vacancies cluster; row 2 is where they avoid each other.

However, the third row of images in Figure 6 is the same as the second, but the size-effect signs have been reversed. Examining the two rectangles in Figure 6d,g shows how the brightness of consecutive spots is reversed by the change in size effect (white rectangles). Note however that other spots are relatively independent of this effect; this is one way in which this kind of modelling is useful, as it allows for the combined effects of the different structure factors (in a sense, each correlation has its own “structure factor”) and form factors and how they interact.

Figure 6. Slices of calculated diffuse scattering from different models of CePdSb, this time incorporating the atomic size effect. Row 1: Same as row 2 of Figure 5, but atoms move away from vacancies and vacancies away from each other. Row 2: Same as row 3 of Figure 5, but atoms move away from vacancies and vacancies away from each other. Row 3: Same as row 3 of Figure 5, but atoms move toward vacancies and vacancies toward each other. For details, see the text. h and k axes noted on one figure to indicate directions.

Figure 6. Slices of calculated diffuse scattering from different models of CePdSb, this time incorporating the atomic size effect. Row 1: Same as row 2 of Figure 5, but atoms move away from vacancies and vacancies away from each other. Row 2: Same as row 3 of Figure 5, but atoms move away from vacancies and vacancies away from each other. Row 3: Same as row 3 of Figure 5, but atoms move toward vacancies and vacancies toward each other. For details, see the text. h and k axes noted on one figure to indicate directions.

Note how the size effect is very different when applied to the clustering model (row 1) and the anti-clustering models (rows 2 and 3). Rows 2 and 3 of Figure 5 are different, but relatively subtly. Compare then rows 1 and 2 of Figure 6, which are the same two rows, now with the same kinds of size effects applied. Because the fraction of atom-vacancy bonds and atom-atom bonds is very different in the two models, the scattering is very different. This shows how strongly these effects can interact, something that can be difficult to disentangle without this kind of modelling to lean on.

Hence, even these relatively simple effects can have interesting and complex influences on the diffraction patterns of metallic systems. The MC model allows insight to be gained when the system is too complex to use direct inversion of the diffuse scattering to determine the correlations. In particular, exploring a range of representative models that look at various possible forms of SRO and their resulting diffraction is a useful guide to finding out what kinds of SRO are present in the real system.

5. Conclusions

Complex metallic systems, such as intermetallics, alloys, quasicrystals and Hume-Rothery phases, can all show detailed local ordering, which gives rise to highly structured and often very anisotropic single-crystal diffuse scattering. This paper reviews some of the issues associated with collecting and analysing such scattering and uses hypothetical calculation on the intermetallic CePdSb to illustrate some of the effects that may be observed in real systems.

Local order is important in determining many materials’ properties and should not be ignored when trying to relate structure to function, especially when phenomena on the nanoscale are to be considered.

By qualitatively inducing various orderings in an MC model, the signatures of these orderings can be determined and compared to the observed data, providing guidance as to what structures are present in the real material.