Chemical ordering beyond the superstructure in long-range ordered systems

To describe chemical ordering in solid solutions systems Warren-Cowley short-range parameters are ordinarily used. However, they are not directly suited for application to long-range ordered systems, as they do not converge to zero for large separations. It is the aim of this paper to generalize the theory to long-range ordered systems and quantitatively discuss chemical short-range order beyond the superstructure arrangements. This is demonstrated on the example of a non-stoichiometric B2-ordered intermetallic alloy. Parameters of interatomic potentials are taken from an embedded atom method (EAM) calculations and the degree of order is simulated by the Monte Carlo method. Both on-lattice and off-lattice methods, where the latter allows individual atoms to deviate from their regular lattice sites, were used, and the resulting effects are discussed.


Introduction
In elemental systems consisting of one type of atom, the structural configuration on the atomic scale in principle is defined by the crystal lattice. In contrast, in multi-component systems aspects of order give rise to additional degrees of freedom. On the one hand, this pertains to long-range order, where the crystal structure's sublattices have different elemental make-ups. On the other hand, short-range order corresponds to energetically preferred local arrangements within the stochastic occupations of the lattice sites by chemically distinct atoms.
While the state of long-range order is determined by the sublattice compositions, the number of degrees of freedom of short-range order is in principle unbounded. In the simplest case of a binary system on a Bravais lattice, the Warren-Cawley short-range order parameters [Cowl 50,Warr 51] quantify the probabilities of pairs of sites being occupied by a given combination of elements. These pair-level correlations are accessible directly in the intensity variations of diffuse scattering, and both elastic neutron scattering (for a review see Ref. [Schw 98]) and X-ray scattering (see Ref. [Scho 99b]) in the diffuse regime have been used for determining shortrange order in a range of systems.
Determination of long-range order via the intensities of Bragg and super-structure peaks is a routine experimental procedure [Xiao 95,Inde 90]. In contrast, diffuse scattering experiments are much more labourious. This pertains even more for studies of short-range order in long-range ordered compounds, where after the pioneering work by Georgopoulos and Cohen [Geor 81] on β -NiAl not much information has been published. The richness of aspects of short-range order that can generally be expected in ordered compounds, including both correlations of the occupations of sites within a given sublattice and on distinct sublattices, might be a reason for this. Note that even in the above-mentioned example of β -NiAl the occupational disorder is restricted to a single sublattice, which allowed it to be described within the conventional Bravais lattice formalism.
The outstanding properties of ordered intermetallics, including high yield strengths and high Young's moduli at low mass densities, corrosion resistance and high-temperature creep resistance [Saut 95, Pfei 07], to name just a few, are due to the ordered arrangement of the atoms. However, lacking a detailed picture of local atomic arrangements, this has hitherto been considered only in terms of long-range order parameters. It is the aim of this work to take the first step towards a full model by providing the necessary formalism for describing short-range order correlations within longrange ordered systems, and to elucidate the accompanying phenomena by way of simulating temperature-dependent long-and short-range order. Specifically, we will consider the evolution of short-range order over an order-disorder transition in a non-stoichiometric B2 (CsCl) intermetallic alloy. We will assume a realistic potential in the vein of the embedded-atom method (EAM) [Daw 84]. Such potentials are non-discrete, which allows us to compare simulations that explore the full classical parameter space including atomic displacements to simulations where the atoms are restricted to positions on an ideal lattice.
The paper will be organized in the following way: In Section 2, we will review the classical short-range order parameters and give a generalization for defining ordering on sublattices. We will detail the applied simulation technique in Section 3, which will be followed by a discussion of the results in Section 4 and a brief summary in Section 5.

Theory
Warren and Cowley [Cowl 50, Warr 51, Cowl 60] introduced parameters to quantify short-range order in binary systems. They can be defined as where c A (c B ) is the overall concentration of atoms of type A (B) and p AB lmn is the probability that a given pair of sites separated by a lattice vector r lmn is occupied by distinct atoms. These parameters were originally introduced to describe SRO in systems where no LRO is present. Without symmetry breaking, short-range interactions lead to correlations that decay with distance, so that for large separations the occupation of the sites is independent and therefore p AB lmn = c A c B and α lmn = 0. Generally, the scattered intensity due to the sites of a Bravais lattice being occupied by different elements can be described by the expression [Schw 13] I( q ) = lmn α lmn cos( q · r lmn ). (2) Here q denotes the position in reciprocal space, r lmn are the lattice vectors, and the intensity is understood in terms of the Laue unit where f X is the q-dependent scattering factor. If the temperature compared to the ordering forces in the system is sufficiently low, a long-range superstructure lattice is formed. To give an example one may consider a B2 ordered system (ClCs structure), where the underlying bodycentred cubic Bravais lattice dissociates into two simple cubic sublattices denoted as α and β. The long-range order parameter is given by the difference of the occupations of the sublattices where c X χ is the probability for a given site on sublattice χ to be occupied by an atom of element X.
With long-range order, the definition given in Eq.
(1) will give non-zero values even for large r lmn . For a stoichiometric system with perfect order and no defects the short-range order parameter α lmn is either -1 or 1, depending on whether r lmn connects sites on the same or on different sublattices. In this case Eq. (2) gives δ-like superstructure peaks. The same is true for a long-range ordered system where the atomic concentration is non-stoichiometric, only that the excess atoms lead to the absolute value of the short-range order parameters being smaller than 1.
Specifically, assuming no short-range order we havē for atoms within different sublattices and for atoms within the same sublattice. With the definition of η and Eq. (1) we get where the positive sign is for same sublattices (ξ = χ) and the negative sign for different sublattices.
Under the presence of short-range order, the Warren-Cowley parameters for small r lmn will reflect the actual short-range order, but for large separations they will converge to the appropriate limiting value of the above two. The short-range ordering can now be quantified by the difference of the actual value to the uncorrelated expression of Eq. (5) which indeed will converge towards zero after a small number of shells for non-critical temperatures. This allows the classical short-range order parameter to be written as a combination of a part that is due to long range order and a part that is strictly due to short-range ordering. Note that for the B2 case considered here, with two inequivalent sublattices, for given lmn there are either two intra-sublattice short-range order coefficients or one intersublattice coefficient, corresponding to the cases respectively.

Simulation
Configurations representative for the equilibrium configuration of the system were generated by the Metropolis-Hastings Monte Carlo method [Land 09]. The internal energies of the specific configurations were computed according to the embedded atom method (EAM)-like potentials with long-range atomic interactions proposed by Ouyang et al. [Ouya 12] for Fe-Al intermetallic alloys. For the two superstructures on the body-centred cubic lattice, the B2 phase with FeAl stoichiometry and the D03 phase with Fe3Al stoichiometry, the resulting lattice constants, elastic constants, heat of fusion, point defect formation enthalpies, and the phonon dispersions were reported to agree very well with available experimental data. [Ouya 12,Ouya] In the following we will identify Fe with A and Al with B.
We implemented a grand-canonical Monte Carlo simulation, where a lattice site is chosen stochastically and its occupation is flipped with a probability that depends on the difference in the internal energies of the original and the flipped configuration and on a chemical potential. In order to achieve high efficiencies in the face of highly differing sublattice occupations, we allowed for a bias in attempt frequencies: specifically, if for a given sublattice the trial frequency for A→B swaps is ωA→B and analogously for the reverse transition, then performing the exchanges with probabilities fulfills detailed balance.[Hast 70] Note that here ∆H incorporates both the differences in internal energy as well as the effect of the chemical potential. We implemented this bias by keeping lists of the sites on each sublattice occupied by either element. Choosing ωA→B/ωB→A approximately as c B ξ /c A ξ (but fixed for the duration of the simulation) maximizes the acceptance rate and therefore the efficiency. We ensured the correct composition of the system via a temperaturedependent chemical potential. We simulate a physically meaningful evolution of time by incrementing at each trial with an exponentially-distributed random variable with a parameter according to the instantaneous trial frequency, which depends on the sublattice occupations. We considered two statistical ensembles. The on-lattice approach constrains the atoms to their periodical lattice sites, so that the occupations of the sites are the only degrees of freedom. In the off-lattice approach, the atoms have the additional degree of freedom of deviating from their ideal positions. This was implemented by shifting random single atoms by random vectors and accepting the move with the Metropolis probabilities. As specific occupations of the lattice sites prefer specific off-lattice deviations (for instance, antistructure atoms of the larger element will lead to outward relaxations), allowing this degree of freedom slows down the occupational dynamics. We observed a slowingdown of about a factor of ten. Further, in the off-lattice simulations the lattice parameter was considered as an additional degree of freedom subject to Metropolis dynamics under zero external pressure, allowing the simulation cell to thermally fluctuate in size and to expand with increasing temperature.
We performed the averaging necessary to compute the simulated statistical parameters (such as the internal energy, lattice constant, sublattice occupations and pair distribution functions) by an exact temporal integral over the configurations as opposed to the conventional sampling approach. This was implemented by keeping records of the point in time a given quantity was changed the last time, and adding its value to the accumulating average with a weight given by the elapsed timespan when it is changed the next time. For the example of the pair distribution function, the neighbours of a flipped site are gone over, and the appropriate entry in the pair distributions is incremented proportional to the age of the pair, which is the minimum of the age of the two participating sites. Thanks to this approach and to the biasing discussed above, we can detect also correlations in defect species that become very small at low temperature, which would not be possible by a sampling approach.
The size of the system with periodic boundary conditions for the simulation cell was chosen to be 32 3 × 2 = 65536 atoms. We have verified that this size was large enough to avoid serious finite-size effects, which is also corroborated by the sharpness of the ordering phase transition we will report below. A composition of c A = 54% was used for the calculations. In the on-lattice case, we used a lattice constant of a = 2.978Å at all temperatures, which is the low-temperature zero-pressure equilibrium value for this composition. No vacant sites were allowed.
Ideally, in the off-lattice case atoms are kept on their position in the lattice only by the interactions with their neighbours. At elevated temperatures this would give a non-zero possibility for the system to transition towards a neighbouring local minimum of the high-dimensional potential landscape even under the absence of point defects, corresponding to diffusional jumps via some mechanism of direct exchange. This would break neighbourhood relations and prohibit the computation of order parameters. Therefore, during atom displacements we enforced the constraint that the position of each atom in relation to its eight direct neighbours has to conform to the order relations on the lattice, specifically, along each dimension the coordinate of a given central atom has to be larger than those of the four neighbours on one side and smaller that those of the neighbour on the other side.
Of course, this constraint suppresses melting and leads to a different ensemble to be simulated. However, it is plausible that significant probabilities of such large displacements would lead to a melting of the unconstrained crystal, and indeed they appear only at temperatures where the lattice parameter has expanded by 3% and root-mean-square displacements with respect to the nearest neighbours have reached 18% of the nearest-neighbour distance. Following the Lindemann criterion[Grim 74], we assume that the range of existence of the crystal has been covered and stop our simulation at this point. In contrast, in the on-lattice ensemble the atoms are fixed to their lattice sites so that crystals of arbitrarily large temperature can be considered and simulated.

Results
As the most obvious aspect of chemical ordering in the system, long-range order was calculated as a function of temperatures in keeping with Eq. (3). Results are shown in Fig. 1. They display the prototypical behaviour of Ising-like systems with only finite-energy excitations, so that the order parameter converges to its zero-temperature value faster than any power law. Specifically, below about 0.5 Tc the system is practically fully long-range ordered, that is, the α sublattice is occupied exclusively by A atoms, while the surplus of A atoms leads to constitutional antisites on the β sublattice. Further, it is evident that the choice of the statistical ensemble has a significant effect on the bcc ↔ B2 phase transition temperature, with kBTc equal to 0.282 eV in the on-lattice case and about 0.183 eV in the off-lattice case. This is not surprising for two reasons: First, it is widely assumed [Fult 10] that oscillatory normal modes in a given system get softer with disorder, thereby increasing vibrational entropy and lowering the free energy of the disordered phase. Of course, the system has this freedom only in the off-lattice approach, which would imply a lowering of the critical temperature compared to the on-lattice case. Second, the relaxations around antistructure atoms mentioned above decrease their formation energy compared to the onlattice case. This results in disorder being also less costly in energy in case of the off-lattice approach.
The temperature-dependent short-range order parameters for the first four neighbouring shells were calculated according to Eq. (1). The results are given in Fig. 2. It can be seen that at low temperatures, where the system is well ordered, the classical short-range order parameters do not deviate appreciably from the appropriate non-correlated values according to Eq. (5). This holds true even more for the inter-sublattice parameters, as the vanishing antisite concentration on the α sublattice c B α does not allow for any correlations to develop. For the intra-sublattice parameters very small deviations are discernible at low temperatures, which have to be due to short-range ordering of the consti- tutional antisites on the β sublattice. In contrast, on the high-temperature side sizeable short-range correlations are stable to the highest simulated temperatures, which, as was to be expected, are largest for closest neighbours. While the classical short-range order parameters are sufficient to calculate the total occupational scattering according to Eq. (2), in the long-range ordered phase a large fraction of this intensity goes into the superstructure peaks. For discussing additional short-range order (or equivalently diffuse scattering), we also computed the sublattice-resolved short-range order coefficients according to Eq. (6). Results are shown in Fig. 3.
From the on-lattice results it is evident that above Tc the two sublattices are equivalent, and therefore the pertaining short-range order parameters are equal, and, as can be deduced from Eqs. (7), they are also equal to the classical parameters in Fig. 2. Obviously also above Tc, the shortrange order already hints at the eventual B2 symmetry breaking, as the inter-sublattice parameters are consistently negative, corresponding to a preference for pairs of unlike atoms if the sites are on different sublattices. Long-range order as illustrated in Fig. 1 obviously is not sensitive to these issues. Mutatis mutandi, the same holds for the intrasublattice parameters. Here it is worth pointing out that, perhaps contrary to intuition, over a wide temperature range correlations over 220 are stronger than those over 200 , even though the latter correspond to smaller separations.
The general behaviour of negative inter-and positive intra-sublattice short-range order parameters persists also some way below Tc, while the antisite concentration on the α-sublattice is becoming increasingly smaller and correlations between these antisites vanish. The constitutional antisites on the β-sublattice, on the other hand, develop a characteristic signature corresponding to an increasing preference for nearest-neighbour pairs of unlike atoms within the sublattice. This suggests that the system would undergo D03 ordering at even lower temperatures. Further, it is remarkable that up to a non-linear rescaling of temperature, the off-lattice short-range order parameters behave quantitatively very similar to the on-lattice parameters.
The short-range order diffuse scattering intensity as given by Eq. (2) is depicted in Fig. 4. For the two chosen temperatures both above and below Tc, the diffuse scattering looks qualitatively very similar. However, while the intensity stays roughly constant at the Γ point, it decreases pronouncedly over the rest of the Brillouin zone during the ordering of the system. Note that the integrated scattering that is due to the sites of a lattice being occupied by idealized point-like scatterers (as assumed in the definition of short-range order intensity) does not depend on the arrangement of the scatterers. When long-range ordering is present, the δ-like superstructure peaks compensate the reduction in short-range order diffuse scattering intensity. This is a consequence of a version of Parseval's theorem, which relates the squared l2-norm of a function on a lattice (here the scattering lengths) to the squared L2-norm of its Fourier transform (here the intensity). The Laue unit is chosen just so that the mean intensity is unity. If there is no long-range order, this directly holds also for the diffuse short-range order intensity. Under non-zero long-range order, on the other hand, the super-structure peaks account for a part of this scattered intensity, which in turn decreases the mean diffuse intensity.
Apart from the decrease in intensity, the figures also demonstrate that the short-range order diffuse scattering in general has only the periodicity of the reciprocal lattice of the Bravais lattice. In the present case, only at very low temperatures, where because of the vanishing antisite concentration on the α-sublattice all inter-sublattice parameters go to zero, the short-range order intensity would increasingly display the simple-cubic reciprocal-space symmetry of the B2 structure.

Conclusion
In the theoretical part of this work, we have shown how for long-range ordered systems the classical short-range order coefficients can be split into a term that depends only on the degree of long-range order and a term that is due to actual short-range order. The first term depends only on the sublattices the respective sites are on and does therefore not decay with distance. It is responsible for the sharp superstructure peaks. The second term represents deviations in the correlations in pair occupations from the long-range order term. For vectors in the structure's Bravais lattice, i.e. vectors that connect sites within a given sublattice, this short-range order term can further be written as a sum of parameters of correlations within the distinct sublattices.
Going beyond standard Metropolis simulations, we have presented an implementation of a Monte Carlo simulation in Hastings' framework, with a site-age based computation of the desired statistical quantities. These two complications allowed for an efficient study also of the very rare defects on the majority sublattice of an off-stoichiometric system. We simulated according to an EAM model of the Fe-Al system and demonstrated an increase in disordering temperature Finally, as an exemplary use of the formalism introduced in the first part of this work, we computed the evolution of the distinct short-range order coefficients with temperature, showing an increase of short-range order correlations with decreasing temperature in the disordered phase, with a maximum at the ordering temperature, and a successive decrease as the correlations become long-ranged. We also illustrated the corresponding diffuse scattering.
Additionally we could show that even tough on-lattice and off-lattice simulation approaches show qualitatively similar results for ordering parameters, the differences can not be neglected.
Apart from exemplifying the concepts as done here, it is also conceivable to use simulations of diffuse scattering as substitutes for experimental data when short-range order information is necessary. Specifically, this pertains to the problem of data modeling in atomic-scale X-ray Photon Correlation Spectroscopy (aXPCS) [Leit 12,Stan 14], where in the simplest approximation[Leit 11] the measured q-dependent correlation times are a product of the shortrange order diffuse scattering and a factor depending on the jump geometry. Only in rare cases[Leit 09] experimental short-range order data is available for the same composition and temperature range[Scho 99a], or simple qualitative models such as nearest-neighbour site exclusions[Stan 13] can describe the data. In the general case, simulations of the diffuse intensity as presented here can serve as a starting point for interpretation [Stan 16]. Further, studies as presented here can also be used to test potentials against each other or against experimental data, as the diffuse shortrange order intensity is sensitive to minute details of the atomic interactions.