Abstract

Lattice distortions are intrinsic features of all solid solution alloys associated with varying atomic radii; this phenomenon facilitates the formation of single-phase solid solutions. Using high-entropy alloys (HEAs), as an example, we investigate the influence of variations in inter-atomic separations for stabilizing and controlling their structural, mechanical, and thermodynamic properties. This is done through a combination of statistical mechanics analysis and molecular dynamics simulations on simplified 2D systems, as well as a 3D crystals with harmonic and anharmonic inter-atomic bonds with varying natural inter-atomic separations. We demonstrate that the impact of this inter-atomic length disorder (representing static lattice distortion) and temperature fluctuations (representing dynamic lattice distortion) on fundamental and universal thermodynamic, structural, and elastic characteristics are similar and can be unified through effective temperature; i.e. a scaling law for HEAs that establishes a relationship between these factors. This scaling law reveals that different HEAs (i.e. varying degrees of local lattice distortions) collapse onto a single curve when plotted against the effective temperature. We demonstrate that lattice distortion significantly enhances the stability of solid solution alloys (relative to phase separation or ordering by effectively increasing the temperature of the system; this stabilization effect is particularly pronounced in HEAs).

Significance Statement

Compositionally induced atomic spacing fluctuations enhance the thermodynamic stability of multiprincipal component alloys. This is demonstrated through rigorous statistical mechanics of low-dimensional models and validated in 3D materials via classical molecular dynamics. We identify a robust scaling law that formally connects static bond length and thermal vibration-induced (dynamic) dispersion.

Introduction

High entropy alloys (HEAs) have garnered significant attention in recent years due to their fascinating range of physical and mechanical properties (e.g. see Refs (1–4)). These exceptional properties have been achieved by expanding the composition palette far beyond that of most conventional alloys. HEAs are often described as materials with nearly equimolar compositions of five or more principal elemental components (5). The term “high entropy” is associated with the observation that the entropy of mixing in an ideal, random solid solution increases with the number of components (n) according to the equation ΔSm=kBicilnci, where ci represents the atomic fraction of component i, the sum is over all n components, kB is the Boltzmann constant. This expression reduces to ΔSm=kBlnn for an equimolar alloy. The increase in the entropy of mixing with n contributes to the enhanced stability of the solid solution phase, reducing the tendency for formation of multiphase or ordered structures.

The entropy of an HEA encompasses contributions that extend beyond the straightforward entropy of mixing, including vibrational, electronic, and magnetic components (6). We note that due to the differing atomic sizes of individual alloy components, the atoms in an HEA do not sit on ideal lattice sites (defined by the Bravais lattice and basis vectors). Hence, the variation in atomic size and deviation of atom locations from the ideal lattice in a random solid solution provides an additional contribution to the entropy, which we denote noncomposition entropy ΔSnc. (This noncomposition entropy is distinct from the entropy of mixing ΔSm.) In a single-component glassy system, ΔSnc is referred to as “configurational” or “information” entropy or “complexity” (7). Random local lattice distortions (misfit) are responsible, in part, for the high strength of high entropy alloys (e.g. see Refs (8–11)). Here, we focus on the thermodynamic implications of such random distortions.

The effects of atomic-size-induced distortions on structure and properties are difficult to assess in real alloy systems since these are necessarily convoluted with variations in bonding between different types of atoms (12). One approach for simplifying the relationships between atom types and their bonding effects is to use average potentials, such as the average-atom embedded-atomic method (EAM) potential (13). However, this mean-field approach has limitations; it does not account for the influences of lattice distortion or atomic displacement and thus ΔSnc (14).

In this article, our objective is to disentangle these effects and specifically isolate the influence of atomic size-induced distortions in HEAs. To accomplish this, we concentrate on a theoretical model system, where all bond strengths are uniform, while variations in bond lengths persist. The general expression for the potential energy Φ of such a bond stretched to a length L from its natural length L0 is given by, Φ=12!K2(LL0)2+13!K3(LL0)3+  14!K4(LL0)4+. We normalize this expression by dividing both sides with K2L02 (dimensions of energy), yielding:

(1)

where ϕ=Φ/(K2L02), =L/L0, and kn=KnL0n2/K2 are the normalized quantities (see Table S1 for relevant normalized quantities).

Furthermore, to minimize the impact of bonding intricacies, we simplify our approach by idealizing atomic bonding as purely harmonic, i.e. set kn to zero for n>2. In fact, spring (harmonic-bonded) models have been utilized in various studies to investigate structural and mechanical properties for many years (15–20). Barron and Stacey (15, 17, 18) demonstrated that harmonic bond lengths are unaltered by heating, whereas harmonic-bonded crystals contract. Toda-Caraballo (19) employed both density functional theory (DFT) and a spring model to predict lattice displacements in HEAs, using synthetic pair potentials. Liu (20) introduced a method to manipulate disorder in an unstressed spring network through bond-stiffness disorder; randomly distorting an initially close-packed lattice structure showed that the network becomes auxetic. Of course, real materials are never truly harmonic and the anharmonic effects can be introduced by preserving the higher order stiffness terms kn,n>2 in Eq. 1.

We note that thermal vibrations and random atomic size effects share an important feature (6, 21–23). Much like “snapshots” of a crystal at finite temperatures show that atom positions fluctuate from their ideal lattice positions, a crystal that is randomly populated by atoms of varying sizes also exhibits atoms displaced from ideal lattice sites at 0 K. This structural analogy has been widely discussed in the literature (24–27). These discussions commonly focus on how both dynamic and static components contribute to the final atomic displacements from their ideal lattice positions. The observation that vibrational contributions to the free energy can be described by examining a time average of “snapshots” implies a fundamental connection between vibrational and noncompositional entropy. This connection, in turn, suggests an intimate relationship between lattice distortions in HEAs and their thermal behavior. Here, we examine the effect of this correlation between static lattice distortions and dynamic thermal vibrations on thermodynamic behavior.

In our atomic bonding model, we specifically narrow our focus to examine the influence of atomic bond length disorder, rather than the dispersion of atomic sizes. Although an n-component HEA possesses n(n+1)/2 distinct equilibrium bond lengths, in our model, we treat atomic size dispersion as a continuous and smooth bond length distribution, as firstly outlined in Refs (28, 29). We hypothesize that this relationship could potentially be formally articulated using the concept of dual lattices, which is often employed in lattice spin models. However, it is important to note that we have not pursued such a derivation. Consequently, in the subsequent sections of this article, we frequently refer to bond length disorder, which is a broader term than atomic size variation and provides a more comprehensive representation of alloy systems. Leveraging the analogy between vibrational effects and bond length disorder, we aim to predict fundamental properties of HEAs, including thermal expansion and the temperature-dependent behavior of elastic constants.

Results

Thermal behavior of a monatomic, harmonic-bond material

Before delving into the complexities of HEAs, we revisit the thermal behavior of simpler, single-component systems. As previously mentioned, understanding the thermal behavior of a single-component system sheds light on atomic size disorder in HEAs. We start by analyzing our investigation to a 2D system. These low-dimensional models are chosen for their amenability to rigorous theoretical analysis. Eventually, we apply our insights to a 3D crystal structure. In essence, we substitute a lattice of many atoms with varying bond lengths for an ensemble of smaller 2D diatomic models. This allows us to focus on ensemble averages rather than volumetric averages. While our approach is grounded in established principles of statistical mechanics, our primary aim here is to empirically demonstrate its validity. Detailed derivations and key findings for these simplified systems are provided in the Supplementary material.

2D: a harmonic-bond Λ model

The “Λ” model is designed to capture bond angle effects to represent an essential feature of crystal lattices; see Fig. 1(a). The Λ configuration is composed of a single “free” atom, bonded to two fixed atoms via harmonic springs. Here, we fix both the nondimensional spring stiffnesses and nondimensional unstretched spring length to 1.

Schematic of a) a symmetric Λ model and b) a centro-symmetric Λ model. The vertex atom is located at the origin. All geometric parameters are labeled in dimensionless forms. c) The vertex height ⟨y⟩ and d) vertical stiffness kyy vs. temperature T for different bond-length dispersions Σ. e) Vertex height ⟨y⟩ and f) vertical stiffness kyy vs. effective temperature Teff=T+αΣ2. The results for the Λ model and its centro-symmetric version are equivalent.
Fig. 1.

Schematic of a) a symmetric Λ model and b) a centro-symmetric Λ model. The vertex atom is located at the origin. All geometric parameters are labeled in dimensionless forms. c) The vertex height y and d) vertical stiffness kyy vs. temperature T for different bond-length dispersions Σ. e) Vertex height y and f) vertical stiffness kyy vs. effective temperature Teff=T+αΣ2. The results for the Λ model and its centro-symmetric version are equivalent.

The potential energy for the free atom at r=(x,y) is

(2)

where h, b, and γ0 are as in Fig. 1(a). The associated partition function is Z=(2π/β)eβUdrπ/(βbh) (see derivation of Eq. S9).

Expanding the potential energy about r=0 to the third order and the displacement in the partition function about r=0 to the fifth order (see derivation of S(S12)), we find the averaged position of the free (vertex) atom:

(3)

This result shows that the average position of the free atom shifts towards the base (analogous to a reduction of lattice parameter) with increasing temperature and the magnitude of the shift is directly proportional to temperature (at low temperature). The details of the calculations and the energy landscape are presented in Fig. S1. This finite-temperature contraction (y00) is consistent with the energy landscape.

It is important to note that the thermal contraction observed in the Λ model is not primarily associated with its lack of centro-symmetry. Consider the centro-symmetric double-Λ model Fig. 1(b) (where the base atoms are free to move in the vertical direction). Examination of this figure shows that horizontal (thermal) vibrations of the free vertex atom will induce a downward motion of the two upper base nodes/atoms; i.e. thermal contraction. Hence, thermal contraction is not associated with the asymmetry of the basic Λ model—rather it is a generic feature of such a 2D model (in fact for all dimensions, d2). This is clearly shown below for the centro-symmetric harmonic FCC lattice.

We now focus on the temperature dependence of the stiffness in the Λ model in the y-direction (i.e. x=0). The inverse stiffness for the Λ model is ky1=(y0/f)f=0, where f is a force on the atom in the y-direction. To obtain an explicit form of the temperature dependence of the stiffness (in the y-direction), we expand the potential for small displacements (see derivative of Eq. S27) and obtain

(4)

where b^b/h=tan(γ0/2) and γ0 is the bond angle (Fig. 1(a)). The temperature dependence of the stiffness is shown in Fig. 2 for several bond angles. For all γ0>0, the stiffness decreases with increasing temperature. This is associated with the displacement of the free atom mean position toward negative y with increasing temperature and the curvature of the potential energy landscape decreases with decreasing y (see Fig. S1). The stiffness reaches the minimum at γ0=2tan1(2/5)=83.62.

Temperature dependence of the stiffness of the Λ model for several different bond angles γ0 based upon Eq. 4. The inset shows the γ0-dependence of the slope at T=0, [d(ky(T)/ky(0))/dγ0]T=0.
Fig. 2.

Temperature dependence of the stiffness of the Λ model for several different bond angles γ0 based upon Eq. 4. The inset shows the γ0-dependence of the slope at T=0, [d(ky(T)/ky(0))/dγ0]T=0.

To validate the approximate theoretical predictions for the Λ model, we conduct a comparison between the predictions of the temperature dependence of the lattice parameter (Eq. 3) and the stiffness (Eq. 4) with the results obtained from MD simulations, as illustrated in Fig. S3. Remarkably, we observe excellent agreement at lower temperatures, although deviations became more pronounced at higher temperatures due to the inherent limitations of the small-displacement expansion employed in our analysis.

These results signify that, even in the scenario where bonding is perfectly harmonic but the bond angle γ0 deviates from 180, the Λ model experiences contraction and a reduction in the elastic constant as the temperature increases. Notably, a significant portion of the elastic constant’s softening as temperature rises can be attributed to thermal contraction. An examination of 4 reveals that, to the leading order, the elastic constant scales with the elevation of the free atom above its base as h2, with dh/dT<0.

3D: a harmonic-bond FCC crystal

We next turn our attention to the temperature-dependent properties of a 3D FCC crystal with harmonic bonds. Due to the increased geometric complexity inherent in a periodic 3D lattice—as opposed to the simpler diatomic and Λ models—we employ MD simulations (see Methods) to ascertain how temperature influences the lattice parameter a0 and the elastic constants Cij0.

In Fig. 3(a), it is evident that the time-averaged cubic lattice parameter a0 exhibits a decreasing trend as the temperature rises, as documented in Ref. (30). Figure 3(b) showcases a corresponding behavior, where all three cubic elastic constants (C110, C120, C440) also decrease with increasing temperature. These observations of contraction and softening in the elastic constants as temperature increases align seamlessly with the preceding analysis of the temperature dependence within the 2D Λ model. Note that existence of thermal contraction of the harmonic 3D FCC model is insensitive to model details (including second neighbors or not), as shown in the Fig. S7. It should be noted that the thermal expansion coefficient is nonzero near 0 K (in violation of the third law of thermodynamics), like in all classical (not quantum) MD simulations (31–33).

a) Lattice parameter a0 vs. temperature. b) Elastic constants C110, C120, and C440 vs. temperature.
Fig. 3.

a) Lattice parameter a0 vs. temperature. b) Elastic constants C110, C120, and C440 vs. temperature.

In the 2D Λ model, the bond angle undergoes temperature-dependent variations, which are associated with thermal contraction. In contrast, within the 3D FCC lattice, there is no alteration in the mean bond angle due to the presence of inversion symmetry within the FCC crystal structure. Nevertheless, the FCC lattice still experiences contraction as it is heated. This may be deduced from the simple, classical arguments of Stacey and Barron (15, 17) who conclude that transverse thermal vibrations lead to negative thermal expansion. This contraction scales proportionally to the square of the amplitude of the thermal vibrations, like say a Σ2. The amplitude, in turn, is directly proportional to the square root of the temperature, implying that, to the leading order, the contraction is linearly related to temperature—a behavior evident in Fig. 3(a).

Contrary to most materials, whose elastic constants decrease due to the anharmonicity of the crystal as temperature increases, our harmonic FCC crystal exhibits an interesting behavior: it contracts upon heating, but its elastic constants also decrease (34, 35). This seemingly paradoxical behavior can be understood by examining the variation in potential energy with applied strain, as detailed in Fig. S4. Elastic constants are directly linked to the second derivative of the potential energy with respect to strain. To leading order, this second derivative shows a linear decrease with compression. Given that the 3D lattice contracts upon heating, it is logical that the elastic constants would also decrease. This trend is consistent with our MD results (Fig. 3). The asymmetric nature of the potential energy curve (Fig. S4) suggests that thermal vibrations cause the thermal contraction; bonds are more easily compressed than stretched and spend more time in a compressed state. This asymmetric behavior, particularly the increased curvature with rising εxx, accounts for the decreasing trend in C11 as the system heats up and contracts (i.e. negative εxx).

A “multi”component harmonic-bond material

We introduce a simple, generic model for a high entropy alloy in which bonds between nearest neighbor atoms share the same stiffness, are harmonic, and have unstretched bond lengths chosen at random from a Gaussian distribution:

(5)

where the bonds have equilibrium lengths L0 distributed about a mean unstretched bond length L0, 0L0/L0 and the standard deviation of bond length distribution is Σ. This is a continuous version of an HEA with ideal bonds. While greatly simplified, this model allows us to focus solely on the effect of randomness in atomic size in a model amenable to analysis.

In the majority of solid solution HEAs, the atom sizes tend to exhibit close clustering. For instance, in extensively studied alloys like Cantor CrMnFeCoNi (6, 36, 37), the standard deviation in atomic sizes or bond lengths typically remains below 5%. It is noteworthy that at a finite temperature, even a single-component (harmonic) system also features a Gaussian distribution of instantaneous bond lengths. However, in this case, the standard deviation scales with temperature (following a T1/2 relationship) rather than being influenced by local composition variations.

2D: ensemble of random harmonic-bond Λ models

We first study a 2D model that captures bond angle effects present in crystal lattices; see Fig. 1(a). The Λ configuration is composed of a single free vertex atom, bonded to two fixed atoms via harmonic springs. Here, we fix the nondimensional stiffnesses of both bonds and choose each of the two unstretched bond lengths (1 and 2) at random from the Gaussian distribution, Eq. 5. The vertex atom located at r corresponds to the particular unstretched bond lengths 1 and 2. The probability that the vertex atom is at r is

(6)

where β is the inverse temperature, and we introduced an effective potential Ueff(r) and partition function Zeff. This effective potential energy is

(7)

where b and h describe the geometry of the Λ model (see Fig. 1(a)). By comparison with the case without bond length disorder (U(r)=12{[(x+b)2+(y+h)21]2+[(xb)2+(y+h)21]2}), we find that the finite-temperature, uniform bond length Λ model has the same energy as the ensemble averaged (random bond length) model if we set Σ2=β1=T. Since Ueff(r) differs from U(r) by a term independent of r, the equilibrium positions/standard deviations from Ueff(r) and U(r) are the same. So, the disordered system can be described as identical to the uniform-bond-length Λ model at effective temperature: Teff=T+Σ2 (see derivation of Eq. S34).

We conduct MD simulations to ascertain the ensemble-averaged thermal contraction and stiffness properties of the Λ model, examining their dependencies on both the standard deviation in bond lengths and temperature. These results are illustrated in Fig. 1. Similar to the behavior observed in the uniform-bond-length model, we observe that the vertex height y and the stiffness in the y-direction, denoted as kyy, exhibit a decrease as temperature increases. Analyzing the temperature dependencies of both the vertex atom position y and kyy, it is evident that altering the standard deviation in unstretched bond lengths primarily shifts the curves to higher temperatures. However, this effect is notably weaker for kyy and is not observed at the highest temperatures, as illustrated in Fig. 1(d).

Comparison of the MD data for y and kyy with the proposed

(8)

(where α is a constant) show two main features. Firstly, this scaling approach effectively collapses all lattice parameters and stiffness properties as functions of both T and Σ onto single curves, as vividly demonstrated in Fig. 1(d) and (e). The collapse is exceptionally close to perfect for lattice parameters and outstanding for stiffness, particularly at low temperatures. Secondly, it is noteworthy that the constant α is not equal to 1, as evident in Fig. S5. Moreover, α takes on distinct values for lattice parameters (α=0.69) and stiffness (α=0.10), emphasizing the differentiation between these properties.

As derived earlier, it is apparent that α=1 results in a poor agreement for both the vertex height and stiffness, as depicted in Fig. S5. The data strongly indicate that both y and kyy exhibit scaling behavior with respect to both T and Σ2. However, it is important to note that the scaling factors for these two properties are not identical. In essence, the effective temperature scales with both T and Σ2, with α serving as the scaling factor. The figures in Fig. 1(d) and (e) clearly demonstrate the accuracy of this scaling relationship for y and kyy concerning T and Σ2. Notably, α values that yield the best fit to y and kyy are found to be 0.69 and 0.10, respectively.

Concurrently, it is important to note that the derivation of thermal contraction is not the primary focus of our study. Instead, our emphasis lies in utilizing a simplified Λ model to derive the fundamental scaling law for lattice distortion, which constitutes our most significant finding.

3D: random harmonic-bond FCC crystal

We conduct atomistic simulations, encompassing both static relaxation and MD, to delve into a harmonically bonded FCC crystal where unstretched bond lengths are randomly drawn in accordance with Eq. 5. The resultant relaxed structure at 0 K is presented in Fig. 4(a), while Fig. 4(b) showcases the equilibrium bond-length distribution, offering insights into the structural distortions induced in both figures. The 0 K cubic lattice parameter a and the elastic constants (C11, C12, C44) are depicted in Fig. 4(c) and (d), respectively. These are plotted against the variance of the bond-length distribution, Σ2. Mirroring observations from the Λ models, the FCC lattice undergoes contraction with escalating bond length disorder. The effect of bond length dispersion at 0 K on the lattice parameter a and elastic constants mirrors that of the temperature in a scenario devoid of bond length dispersion. This observation is substantiated by comparing Fig. 4(c) and (d) and Figs. 3(a) and (b). Notably, elastic constants C11 and C12 diminish as bond disorder rises—a trend consistent with the Λ-model outcomes where the elastic stiffness kyy (akin to C11) weakens with mounting bond length dispersion. Conversely, the shear elastic constant C44 wanes with bond length disorder, albeit at a reduced rate than its counterparts. To put this in context, the value of kxx in the Λ model actually escalates with bond length disorder, as illustrated in Fig. S5(c). The decrease in elastic constants as lattice distortion increases aligns with Nöhring’s observation and statement (14) that the elastic constants of true random alloys are consistently lower than those of averaged-atom alloys.

a) Distorted harmonically bonded FCC crystal (a part of big sample obtained by MD) for Σ=0.071. b) Bond-length distribution after relaxation; the legend shows the values of its standard deviation Σ. c) The mean lattice parameter and b) elastic constants vs. Σ2 at T=0 and (c)–(f) vs. Teff=T+αΣ2 for several values of Σ. The best fit values of α are 3.6, 4.5, 7.0, and 3.3 for a, C11, C12, and C44, respectively.
Fig. 4.

a) Distorted harmonically bonded FCC crystal (a part of big sample obtained by MD) for Σ=0.071. b) Bond-length distribution after relaxation; the legend shows the values of its standard deviation Σ. c) The mean lattice parameter and b) elastic constants vs. Σ2 at T=0 and (c)–(f) vs. Teff=T+αΣ2 for several values of Σ. The best fit values of α are 3.6, 4.5, 7.0, and 3.3 for a, C11, C12, and C44, respectively.

The synergistic impact of temperature and bond length disorder on the lattice parameter and elastic constants is vividly demonstrated in Fig. 4(e)–(h); for the unscaled data, readers can refer to Fig. S6. Utilizing Eq. 8 to define an effective temperature Teff yields an almost flawless scaling or data collapse, reinforcing the validity of the proposed scaling relation. This lends credence to the idea that thermal and disorder effects in crystalline materials behave in a similar manner to those in the Λ models, a hypothesis initially supported by partition function analysis (applicable at least within the regime of low Teff). Both increasing disorder and rising temperature result in lattice contraction and a softening of the elastic constants.

Similar to the Λ model, the convergence of random bond lengths and finite-temperature FCC crystal data onto a unified curve underscores the previously observed “equivalence” between temperature and bond length disorder. The value of the parameter α in Teff is determined by fitting it to simulation data at low Teff for the parameters a, C11, C12, and C44; specifically, α is found to be 3.6±0.2, 4.5±0.2, 7.0±0.5, and 3.3±0.3, respectively. The distinct values of α associated with each quantity stem from the different geometries of the corresponding deformations: dilatation (a), uniaxial strain (C11), lateral strain (C12), and shear (C44).

The concept of an effective temperature, denoted as Teff, encapsulates the notion that structural disorder arising from thermal vibrations and local static lattice distortions is statistically indistinguishable. While we illustrate how T and Σ synergistically combine within the framework of an effective temperature, we have not yet provided an analytical prediction for the coefficients α. The value of α exhibits a complex dependence on multiple factors, including bond connectivity, lattice geometry, and the direction in which the properties of interest are measured, particularly in the case of tensor properties like Cij. As a result, α is not a universal constant; it does not possess the same value for all properties within a specific HEA, nor does it remain constant across all HEAs for a given property.

Additionally, we conduct MD simulations on a 3D lattice system characterized by bond strength dispersion while maintaining identical bond lengths. The outcomes of these simulations, as depicted in Fig. S22, reveal that bond strength distribution exerts only a minimal impact on the temperature-dependent behavior of lattice parameters or elastic constants. This outcome is not unexpected, as bond strength dispersion, under conditions of a constant equilibrium bond length, does not induce local lattice distortions at 0 K.

Discussion and conclusions

The primary objective of this study is to investigate the impact of bond length and atomic size disorder on structural and mechanical properties. The resulting lattice distortion plays a pivotal role in various aspects of compositionally complex and high-entropy alloys, as evidenced by prior research (38–41). To isolate the effects of atomic size variations and the resulting lattice distortions on thermodynamic behavior, we intentionally neglect the influence of bond strength or stiffness variations among alloy components. Our approach is rooted in a statistical mechanics treatment of both distortions arising from thermal fluctuations in bond length and those associated with bond length disorder. We explore these phenomena within a simple 2D model and subsequently validate the emergent features within a 3D crystal lattice using MD simulations.

The primary finding of this study is that the partition function for harmonic bonds, with a single unstretched bond length at a finite temperature (T), is equivalent to that of a system where unstretched bond lengths are randomly selected from a Gaussian distribution at T=0 in the 2D (Λ) model. This result can be understood by observing that the thermal distribution of bond lengths in a harmonic model follows a Gaussian distribution with a standard deviation proportional to T1/2. When the distribution of unstretched bond lengths is also Gaussian with a standard deviation of Σ, the two systems exhibit indistinguishable instantaneous configurations. The thermal average bond length and stiffness in a system where unstretched bond lengths are randomly selected are identical to those in a system with a single unstretched bond length at temperature T, provided we replace T with an effective temperature, denoted as Teff=T+αΣ2, where α is a constant. These findings have been validated through MD simulations in the 2D system (see Fig. 1).

While a direct analytical evaluation of the partition function for 3D FCC crystals is not feasible, we conduct MD simulations at various temperatures for both scenarios: the single unstretched bond length case and the case where unstretched bond lengths are randomly selected from a Gaussian distribution. Our simulations reveal that the equilibrium cubic lattice parameter (a) and the three cubic elastic constants (C11, C12, and C44) for both systems converge onto single curves at an effective temperature (Teff), as illustrated in Figs. 4 and S6. This outcome demonstrates the applicability of the 2D theoretical results to the 3D FCC lattice system. Additionally, similar effective temperature behavior is observed in the FCC model in the presence of the anharmonic term (k3). The effect of the strength of anharmonicity is studied for three values of k3 and the results are shown in Figs. S11–S13 (representatives for negative thermal expansion, Invar, and positive thermal expansion). Figures S14–S16 show how each of these families of curves collapse onto a single curve each with a different α corresponding to each physical quantity a, C11, C12, and C44 and for different values of k3. We further observe a similar convergence of data for specific properties with the introduction of an effective temperature (Teff) in other lattice structures, as demonstrated in the case of the 3D HCP lattice in Fig. S21.

The finding that temperature and bond length disorder are interchangeable holds intriguing thermodynamic implications for compositionally complex and high-entropy alloys. The contribution of compositional complexity to the system’s entropy is typically calculated using a mean-field or Bragg–Williams approximation based on random site substitution. In an equi-atomic, random solid solution comprising n components, the entropy of mixing is represented as ΔSm=kBlnn. This is often regarded as the primary contribution to the change in free energy when transitioning from a homogeneous single-component alloy to a high-entropy alloy, neglecting internal energy effects. Given that the key findings derived from the 2D Λ model align well with our MD simulations of the FCC lattice, we revisit the 2D Λ model to assess the influence of atomic size and unstretched bond length dispersion on the thermodynamic behavior of our harmonic-bond 3D FCC system.

The partition function of the Λ model (see Eq. S9) is Z=2π/βπ/(βeffbh)=(2π2T/bh)(T+αΣ2), where βeff=1/Teff. From this, we derive the Helmholtz free energy (F=TlnZ), as well as the entropy and internal energy. We can also explicitly add a term corresponding to the noncomposition entropy ΔSnc (for an n component system) in the free energy

(9)

where C is a numerical constant accounting for the detailed geometry of the model. The changes in free energy, entropy and internal energy on going from a mono-dispersed bond length system to one with bond length dispersion are

(10)
(11)
(12)

where we have set α=1 (for simplicity) and the final expressions are the leading order terms in Σ2/T. As a result, bond length distribution introduces a temperature-dependent entropy component, supplementing the traditional entropy term associated with the random distribution of atoms of different types on lattice sites (kBlnn). The contribution of noncompositional entropy (Snc) is also quantified by Nöhring and Curtin (14) using the average-atom EAM potential. Notably, this contribution is relatively small in scale. Their findings regarding the free energy difference (ΔF) align with Eq. 10, particularly in the low-temperature limit.

To contextualize the thermodynamic impact of bond length disorder, we have estimated its contribution to several HEAs. Our analysis involves MD simulations for a CuCrNiFeCo alloy system using EAM potentials (42). Figure 5(a) and (b) depict a typical relaxed structure and the associated bond length distribution. These results, as shown in Fig. 5(c)–(f), clearly illustrate that the scaling law involving the effective temperature (Teff) established for simple harmonic systems is also applicable to a real HEA.

An FCC (CuCr)x(NiFeCo)1−x HEA: a) Relaxed structure (a section of a much larger sample) for x=0.4 and Σ=0.0184 Å. b) Bond-length distribution after relaxation; the legend shows the concentration of CuCr and the corresponding standard deviation Σ. c) lattice parameter and (b)–(d) elastic constants (C11, C12, C44) vs. effective temperature for various Σ values.
Fig. 5.

An FCC (CuCr)x(NiFeCo)1−x HEA: a) Relaxed structure (a section of a much larger sample) for x=0.4 and Σ=0.0184 Å. b) Bond-length distribution after relaxation; the legend shows the concentration of CuCr and the corresponding standard deviation Σ. c) lattice parameter and (b)–(d) elastic constants (C11, C12, C44) vs. effective temperature for various Σ values.

We have also calculated the standard deviation of the nearest-neighbor distances for several experimentally realized HEAs. Specifically, for the FCC Cantor alloy (CoCrMnFeNi), we find Σ=0.0303, and for FCC CoCrCuFeNi, Σ=0.0294 (refer to Table S2 for more details). As previously discussed, bond length disorder effectively elevates the temperature above its actual value. At room temperature, this amounts to an increase of 42% for CoCrMnFeNi and 41% for CoCrCuFeNi (details are available in Table S2). At 1,000 K, the influence of bond length disorder is more modest, effectively raising the temperature by 12% for both alloys. It is crucial to highlight that while bond length disorder significantly impacts thermodynamic properties at lower temperatures, its effect on elastic constants and melting temperatures is much less pronounced. Moreover, these impacts are generally subdued in FCC HEAs owing to their lower levels of local lattice disorder (Σ5%) (29, 43) (see Figs S18 and S23). When comparing the magnitudes of atomic displacements caused by vibration and bond length effects, it is important to note that in real systems, the intrinsic lattice distortion scale is typically smaller than the vibrational scale. This observation applies to the variations in the atomic radii of the constituent elements, regardless of the presence of fewer or more atom types (cf. Fig. 5(a) and Fig. 6 in Ref. (39) and the Discussion in Ref. (25)).

It is intriguing to observe that one of the classic Hume–Rothery rules posits that stable substitutional solutions exist only when the atomic sizes of the solute and solvent atoms differ by no more than 15% (44, 45). This implies that there could be an upper bound on the extent of atomic size dispersion allowable in HEAs. The relevance of this rule to HEAs is an open question deserving further investigation; for example, see Ref. (46). We should note that in the HEAs we have referenced in Table S2, the maximum value of Σ/l¯ falls below this 15% threshold.

We demonstrate that temperature and lattice distortion jointly contribute to numerous thermodynamic, structural, and elastic properties via a unified effective temperature, encompassing both the actual temperature and the extent of lattice distortion. Lattice distortion notably enhances the stability of a single-phase solid solution, as commonly found in binary and multicomponent alloys, by elevating the effective temperature. Increasing temperature further reinforces the stability of a solid solution compared to phase separation or ordering. The established scaling law for lattice distortion is validated through simulations in a 2D Λ model, a 3D harmonic and anharmonic crystal model, and by employing MD simulations on a real HEA.

Methods

Molecular dynamics (MD) simulations of the 2D three-atom system depicted in Fig. 1(a) is performed using our in-house MD code, where the masses of the atoms are 1 gm/mole and the velocity Verlet integration time step is taken to be 0.05 fs. Temperature is maintained through a Langevin thermostat with a friction coefficient of 0.2 time units.

3D simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (47). In these simulations, all nearest atoms are connected by elastic harmonic bonds. The equilibrium lattice parameters at finite temperatures are obtained using the Nosé–Hoover thermostat and barostat (NPT) with a time step of 1 fs. The elastic constant tensor is obtained by applying strains to the box within NVT ensemble. The elastic constants at all temperatures are determined at fixed strains are calculated from Cij=Δσ/ε, where Δσ represents the stress change and ε is the applied strain. The reported elastic constants are obtained by extrapolating the finite strain results to zero strain at any temperature. Simulations cover a finite temperature range from 300 K to 1,500 K, with intervals of 200 K. The MD results are obtained through time and symmetry averaging in a system comprising 500,000 atoms, following the equilibration process.

The selection of the bond length distribution function, along with the extent of bond length disorder represented by Σ, is informed by computational analyses grounded in DFT as well as calculations based on EAM potentials. These analyses cover a wide array of HEAs and are corroborated by multiple studies (29, 43, 48). Empirical observations consistently reveal that the bond length distributions of both the first-nearest neighbor (1NN) and the second-nearest neighbor (2NN) conform to Gaussian patterns. Additionally, the standard deviations of these Gaussian distributions, denoted as Σ, have been observed to lie within the interval [0, 0.11]. For a more in-depth comparison between our findings and the DFT/EAM-calculated bond length distributions at 0 K, refer to Figs S18 and S23. Given these empirical and computational insights, our study adopts a Gaussian distribution for bond lengths, setting the standard deviation Σ within a slightly extended range of [0, 0.15]. For consistency with real HEA potential cases, Gaussian distributions with the standard deviation between 0 and 0.108 after relaxation at 0 K are used in manuscript. In simulations, we find when standard deviation exceeds 0.2, the structure would be unstable and show continuous shrinkage even at lower finite temperatures.

The core of this study revolves around the face-centered cubic (FCC) lattice, characterized by harmonic springs connecting nearest neighbor atoms. This lattice is simply chosen as an example in this article; the results are applicable to other lattices as well (we discuss nearest-neighbor hexagonal close-packed HCP and second nearest-neighbor bonded body centered cubic BCC, harmonically bonded lattices and see Figs. S19–S21, Table S3).

To provide a broader context, MD simulations are executed on a 3D HEA model employing extant EAM potentials. The objective of these simulations is to elucidate the convergences and divergences between the outcomes generated by our simplified model and those derived from more complex EAM-based calculations.

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

Z.W. and D.J.S. gratefully acknowledge the support of the Hong Kong RGC through the Collaborative Research Fund project number 8730054. S.P. acknowledges the support of the Hong Kong RGC General Research Fund (HKU 11211019). J.H. acknowledges support of the Early Career Scheme (ECS) grant from the Hong Kong Research Grants Council CityU21213921 and Donation for Research Projects 9229061.

Author Contributions

Z.W. conceptualization, software, formal analysis, investigation, visualization, methodology, writing original draft; A.S.L.S.P. conceptualization, software, investigation, methodology, writing review and editing; J.H. conceptualization, supervision, investigation, writing review and editing; D.S. conceptualization, resources, formal analysis, supervision, funding acquisition, writing original draft, project administration, writing review and editing.

Data Availability

All data required for main findings of this manuscript are included in the article and supplementary material.

References

1

Ikeda
 
Y
,
Grabowski
 
B
,
Körmann
 
F
.
2019
.
Ab initio phase stabilities and mechanical properties of multicomponent alloys: a comprehensive review for high entropy alloys and compositionally complex alloys
.
Mater Charact
.
147
:
464
511
.

2

George
 
EP
,
Raabe
 
D
,
Ritchie
 
RO
.
2019
.
High-entropy alloys
.
Nat Rev Mater
.
4
(
8
):
515
534
.

3

Daw
 
MS
,
Chandross
 
M
.
2021
.
Sluggish diffusion in random equimolar FCC alloys
.
Phys Rev Mater
.
5
(
4
):
043603
.

4

Roy
 
A
,
Singh
 
P
,
Balasubramanian
 
G
,
Johnson
 
DD
.
2022
.
Vacancy formation energies and migration barriers in multi-principal element alloys
.
Acta Mater
.
226
:
117611
.

5

Murty
 
BS
,
Yeh
 
J-W
,
Ranganathan
 
S
,
Bhattacharjee
 
PP
.
2019
.
High-entropy alloys
. Amsterdam, The Netherlands:
Elsevier
.

6

Ma
 
D
,
Grabowski
 
B
,
Körmann
 
F
,
Neugebauer
 
J
,
Raabe
 
D
.
2015
.
Ab initio thermodynamics of the CoCrFeMnNi high entropy alloy: importance of entropy contributions beyond the configurational one
.
Acta Mater
.
100
:
90
97
.

7

Nieuwenhuizen
 
TM
.
1998
.
Thermodynamics of the glassy state: effective temperature as an additional system parameter
.
Phys Rev Lett
.
80
(
25
):
5580
.

8

Okamoto
 
NL
,
Yuge
 
K
,
Tanaka
 
K
,
Inui
 
H
,
George
 
EP
.
2016
.
Atomic displacement in the CrMnFeCoNi high-entropy alloy–a scaling factor to predict solid solution strengthening
.
AIP Adv
.
6
(
12
):
125008
.

9

Song
 
H
, et al.
2017
.
Local lattice distortion in high-entropy alloys
.
Phys Rev Mater
.
1
(
2
):
023404
.

10

George
 
EP
,
Curtin
 
WA
,
Tasan
 
CC
.
2020
.
High entropy alloys: a focused review of mechanical properties and deformation mechanisms
.
Acta Mater
.
188
:
435
474
.

11

He
 
QF
, et al.
2022
.
A highly distorted ultraelastic chemically complex Elinvar alloy
.
Nature
.
602
(
7896
):
251
257
.

12

Liu
 
W
,
Lu
 
X-G
,
Hu
 
Q-M
.
2023
.
Comprehensive understanding of local lattice distortion in dilute and equiatomic FCC alloys
.
Mater Chem Phys
.
293
:
126928
.

13

Varvenne
 
C
,
Luque
 
A
,
Nöhring
 
WG
,
Curtin
 
WA
.
2016
.
Average-atom interatomic potential for random alloys
.
Phys Rev B
.
93
(
10
):
104201
.

14

Nöhring
 
WG
,
Curtin
 
WA
.
2016
.
Thermodynamic properties of average-atom interatomic potentials for alloys
.
Model Simul Mat Sci Eng
.
24
(
4
):
045017
.

15

Barron
 
THK
.
1957
.
Grüneisen parameters for the equation of state of solids
.
Ann Phys (N Y)
.
1
(
1
):
77
90
.

16

Irvine
 
RD
,
Stacey
 
FD
.
1975
.
Pressure dependence of the thermal Grüneisen parameter, with application to the Earth’s lower mantle and outer core
.
Phys Earth Planet Inter
.
11
(
2
):
157
165
.

17

Stacey
 
FD
.
1977
.
Applications of thermodynamics to fundamental earth physics
.
Geophys Surv
.
3
(
2
):
175
204
.

18

Stacey
 
F
.
1993
.
Properties of a harmonic lattice
.
Phys Earth Planet Inter
.
78
(
1–2
):
19
22
.

19

Toda-Caraballo
 
I
,
Wróbel
 
JS
,
Dudarev
 
SL
,
Nguyen-Manh
 
D
,
Rivera-Díaz-del Castillo
 
PEJ
.
2015
.
Interatomic spacing distribution in multicomponent alloys
.
Acta Mater
.
97
:
156
169
.

20

Liu
 
J
,
Nie
 
Y
,
Tong
 
H
,
Xu
 
N
.
2019
.
Realizing negative Poisson’s ratio in spring networks with close-packed lattice geometries
.
Phys Rev Mater
.
3
(
5
):
055607
.

21

Körmann
 
F
,
Ikeda
 
Y
,
Grabowski
 
B
,
Sluiter
 
MHF
.
2017
.
Phonon broadening in high entropy alloys
.
Npj Comput Mater
.
3
(
1
):
1
9
.

22

Troparevsky
 
MC
,
Morris
 
JR
,
Kent
 
PRC
,
Lupini
 
AR
,
Stocks
 
GM
.
2015
.
Criteria for predicting the formation of single-phase high-entropy alloys
.
Phys Rev X
.
5
(
1
):
011041
.

23

Otto
 
F
,
Yang
 
Y
,
Bei
 
H
,
George
 
EP
.
2013
.
Relative effects of enthalpy and entropy on the phase stability of equiatomic high-entropy alloys
.
Acta Mater
.
61
(
7
):
2628
2638
.

24

Ramsteiner
 
IB
, et al.
2008
.
Omega-like diffuse X-ray scattering in Ti-V caused by static lattice distortions
.
Acta Mater
.
56
(
6
):
1298
1305
.

25

Owen
 
LR
, et al.
2017
.
An assessment of the lattice strain in the CrMnFeCoNi high-entropy alloy
.
Acta Mater
.
122
:
11
18
.

26

Owen
 
LR
,
Jones
 
NG
.
2020
.
Quantifying local lattice distortions in alloys
.
Scr Mater
.
187
:
428
433
.

27

Chen
 
B
, et al.
2023
.
Correlating dislocation mobility with local lattice distortion in refractory multi-principal element alloys
.
Scr Mater
.
222
:
115048
.

28

Oh
 
HS
, et al.
2016
.
Lattice distortions in the FeCoNiCrMn high entropy alloy studied by theory and experiment
.
Entropy
.
18
(
9
):
321
.

29

Liu
 
S
,
Wei
 
Y
.
2017
.
The Gaussian distribution of lattice size and atomic level heterogeneity in high entropy alloys
.
Extreme Mech Lett
.
11
:
84
88
.

30

Dove
 
MT
,
Fang
 
H
.
2016
.
Negative thermal expansion and associated anomalous physical properties: review of the lattice dynamics theoretical foundation
.
Rep Prog Phys
.
79
(
6
):
066503
.

31

Wang
 
R
, et al.
2022
.
Classical and machine learning interatomic potentials for BCC vanadium
.
Phys Rev Mater
.
6
(
11
):
113603
.

32

Pei
 
Q-X
,
Jhon
 
MH
,
Quek
 
SS
,
Wu
 
Z
.
2021
.
A systematic study of interatomic potentials for mechanical behaviours of Ti-Al alloys
.
Comput Mater Sci
.
188
:
110239
.

33

Wen
 
T
, et al.
2021
.
Specialising neural network potentials for accurate properties and application to the mechanical response of titanium
.
Npj Comput Mater
.
7
(
1
):
206
.

34

Lakkad
 
SC
.
1971
.
Temperature dependence of the elastic constants
.
J Appl Phys
.
42
(
11
):
4277
4281
.

35

Wang
 
Y
, et al.
2010
.
A first-principles approach to finite temperature elastic constants
.
J Condens Matter Phys
.
22
(
22
):
225404
.

36

Ye
 
YF
, et al.
2018
.
Atomic-scale distorted lattice in chemically disordered equimolar complex alloys
.
Acta Mater
.
150
:
182
194
.

37

Otto
 
F
, et al.
2013
.
The influences of temperature and microstructure on the tensile properties of a CoCrFeMnNi high-entropy alloy
.
Acta Mater
.
61
(
15
):
5743
5755
.

38

Zhang
 
L
,
Xiang
 
Y
,
Han
 
J
,
Srolovitz
 
DJ
.
2019
.
The effect of randomness on the strength of high-entropy alloys
.
Acta Mater
.
166
:
424
434
.

39

He
 
QF
, et al.
2021
.
Understanding chemical short-range ordering/demixing coupled with lattice distortion in solid solution high entropy alloys
.
Acta Mater
.
216
:
117140
.

40

Jian
 
W-R
, et al.
2020
.
Effects of lattice distortion and chemical short-range order on the mechanisms of deformation in medium entropy alloy CoCrNi
.
Acta Mater
.
199
:
352
369
.

41

Chen
 
S-M
, et al.
2022
.
Phase decomposition and strengthening in HfNbTaTiZr high entropy alloy from first-principles calculations
.
Acta Mater
.
225
:
117582
.

42

Farkas
 
D
,
Caro
 
A
.
2018
.
Model interatomic potentials and lattice strain in a high-entropy alloy
.
J Mater Res
.
33
(
19
):
3218
3225
.

43

Tong
 
Y
, et al.
2018
.
A comparison study of local lattice distortion in Ni80Pd20 binary alloy and FeCoNiCrPd high-entropy alloy
.
Scr Mater
.
156
:
14
18
.

44

Hume-Rothery
 
W
,
Coles
 
BR
.
1954
.
The transition metals and their alloys
.
Adv Phys
.
3
(
10
):
149
242
.

45

Callister
 
WD
,
Rethwisch
 
DG
.
2018
.
Materials science and engineering: an introduction
.
Vol. 9
.
New York
:
Wiley
.

46

Wang
 
Z
,
Qiu
 
W
,
Yang
 
Y
,
Liu
 
CT
.
2015
.
Atomic-size and lattice-distortion effects in newly developed high-entropy alloys with multiple principal elements
.
Intermetallics
.
64
:
63
69
.

47

Thompson
 
AP
, et al.
2022
.
LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
.
Comput Phys Commun
.
271
:
108171
.

48

Zhang
 
Q
, et al.
2021
.
Deformation mechanisms and remarkable strain hardening in single-crystalline high-entropy-alloy micropillars/nanopillars
.
Nano Lett
.
21
(
8
):
3671
3679
.

Author notes

Competing Interest: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
Editor: Krzysztof Zur
Krzysztof Zur
Editor
Search for other works by this author on:

Supplementary data