Lattice stability of ordered Au-Cu alloys in the warm dense matter regime

We calculate the total energy and the phonon dispersions of Au-Cu alloys based on DFT and DFPT dfpt implemented in Quantum ESPRESSO (QE) code qe . In order to treat the exchange-correlation energy, we use the Perdew-Zunger (PZ) pz functional of the local density approximation and the Perdew-Burke-Ernzerhof (PBE) pbe functionals of the generalized gradient approximation. We use the ultrasoft pseudopotentials of Au.$xc$-n-rrkjus_psl.1.0.0.UPF and Cu.$xc$-dn-rrkjus_psl.1.0.0.UPF $(xc={\rm pz\ or\ pbe})$ that are provided in pslibrary.1.0.0 dalcorso . The cutoff energies for the wavefunction and the charge density are 80 Ry and 800 Ry, respectively. The self-consistent field and phonon dispersion calculations are performed by using 16$\times$16$\times$16 $k$ grid and 4$\times$4$\times$4 $q$ grid (3$\times$3$\times$3 $q$ grid for L1₂) MK , respectively, which are enough to study the dynamical stability of warm dense alloys. The total energy and forces are converged within $10^{-5}$ Ry and $10^{-4}$ a.u., respectively in geometry optimization. For phonon calculations, the threshold parameter for self-consistency (tr2_ph) is set to $10^{-14}$. The Marzari-Vanderbilt smearing smearing with a broadening of $\sigma=0.02$ Ry (0.27 eV) is used for the geometry optimization, while the Fermi-Dirac type smearing is used for the WDM regime (the electron temperature $T_{\rm e}$ in units of eV and the Boltzmann constant $k_{\rm B}=1$). The tetrahedron method is used to calculate the electron density-of-states (DOS) tetra . For the WDM regime, the number of electronic bands to be calculated is increased to more than three times larger the sizes in the ground state calculations in order to take into account the occupation of electronic states far above the Fermi level. Imaginary phonon energies $\hbar\omega$ ($\hbar$ the Planck constant and $\omega$ the phonon frequency) are represented by negative energies below.

Below is a brief description of the basic concepts for the lattice dynamics. We derive analytical expressions for the phonon frequencies at point R in order to study the instability of L1₀ AuCu in detail. The atomic motion in a crystal is regulated by maradudin

where $u_{\kappa\alpha}(\bm{R}_{i})$ is the displacement components for the atom $\kappa$ ($=1$ and 2 for Au and Cu, respectively) along the direction of $\alpha\ (=x,y,z)$ in the unit cell characterized by the lattice vector $\bm{R}_{i}=(R_{ix},R_{iy},R_{iz})$ expanded by the primitive vectors $\bm{a}_{1}=(a,0,0),\bm{a}_{2}=(0,a,0)$, and $\bm{a}_{3}=(0,0,c)$. $M_{\kappa}$ is the mass of the atom $\kappa$ and $D_{\alpha\beta}^{\kappa\kappa^{\prime}}(\bm{R}_{i},\bm{R}_{j})$ is the force constant matrix defined as

where $X_{\kappa\alpha}(\bm{R}_{i})=R_{i\alpha}+\tau_{\kappa\alpha}$ and $\tau_{\kappa\alpha}$ is the $\alpha$-component of the basis vector for $\kappa$: $(\tau_{1x},\tau_{1y},\tau_{1z})=(0,0,0)$ and $(\tau_{2x},\tau_{2y},\tau_{2z})=(a/2,a/2,c/2)$. Assuming the plane wave solution with the wavevector $\bm{q}$, one obtains the eigenvalue equation

where $\epsilon_{\kappa\alpha}$ is the $\alpha$-component of the polarization vector of the atom $\kappa$ and, using the translational symmetry of the crystal, the dynamical matrix $\tilde{D}$ is given by

where we used the notation of $D_{\alpha\beta}^{\kappa\kappa^{\prime}}(\bm{R}_{j})=D_{\alpha\beta}^{\kappa\kappa^{\prime}}(\bm{R}_{j},\bm{0})$. Also, due to the translational symmetry, the acoustic sum rule holds: $\sum_{\kappa^{\prime}j}D_{\beta\alpha}^{\kappa^{\prime}\kappa}(\bm{R}_{j})=0$.

We consider the force constants up to the fifth-order NN sites of L1₀ AuCu: for the Au atom at the origin, the position of the first, second, third, fourth, and fifth NNs are Cu$(a/2,a/2,c/2)$, Au$(a,0,0)$, Au$(0,0,c)$, Au$(a,0,c)$, and Cu$(3a/2,a/2,c/2)$, respectively, and these equivalent sites. The $6\times 6$ dynamical matrix is then expressed as

where the $3\times 3$ diagonal matrix is given by

and the $3\times 3$ off-diagonal matrix for $\kappa\neq\kappa^{\prime}$ is given by

By solving Eq. (3), one obtains the phonon energies at point R in ascending order

with

Based on these expressions for $\omega_{k}\ (k={\rm I}\cdots{\rm VI})$, we interpret the lattice vibrations as follows: the vibration of the Au (Cu) atoms polarized to the $x$ direction propagates along the $\bm{q}=(0,\pi/a,\pi/c)$ direction with the frequency $\omega_{\rm II}$ ($\omega_{\rm IV}$); the vibration of the Au atoms polarized to the $z$ direction is coupled with the vibration of the Cu atoms polarized to the $y$ direction, yielding $\omega_{\rm I}$ and $\omega_{\rm V}$; and the vibration of the Au atoms polarized to the $y$ direction is coupled with the vibration of the Cu atoms polarized to the $z$ direction, yielding $\omega_{\rm III}$ and $\omega_{\rm VI}$. Due to the mode mixing, a frequency gap must be present between $\omega_{\rm I}$ and $\omega_{\rm V}$ ($\omega_{\rm III}$ and $\omega_{\rm VI}$). Negative value in the expression of $\omega_{\rm I}^{2}$ (i.e., $-\tilde{C}_{zy}(\bm{q})$) will be responsible for the instability of L1₀ AuCu in the WDM regime. The expressions of $\tilde{D}_{\alpha\beta}^{\kappa\kappa^{\prime}}(\bm{q})$ as a function of $D_{\alpha\beta}^{\kappa\kappa^{\prime}}(\bm{R}_{j})$ are provided in the Appendix A. The values of $D_{\alpha\beta}^{\kappa\kappa^{\prime}}(\bm{R}_{j})$ are extracted from DFPT calculations described in Sec. II.1.

It has been known that the cohesive energies and bulk moduli of alloys may be sensitive to the exchange-correlation functionals used zhang ; isaacs ; ruzsinszky2019 ; ruzsinszky2020 . In general, the PBE functional underestimates those values for weakly bonded systems such as Au-Cu alloys. Such a discrepancy can be improved by using the hybrid exchange-correlation functional zhang or the random phase approximation ruzsinszky2019 . For more accuracy, it may be better to use the latter approaches zhang ; ruzsinszky2019 . However, we use the standard functionals described above to reduce the computational costs in phonon dispersion calculations.

For studying numerical accuracy in the present calculations, we list the optimized lattice parameters $a$ and $c$ (only for L1₀) predicted by using PZ and PBE functionals in Table 1. These values obtained at $T_{\rm e}=0$ agree with the experimental values okamoto : The mean absolute error (MAE) is estimated to be less than 0.1 Å for both cases. For Au-rich phases the PZ results are in good agreement with the experimental lattice constants, while for Cu-rich phases the PBE results are better.

Figure 2 shows the phonon dispersions along the symmetry lines for Au-Cu systems: (a) fcc Au, (b) fcc Cu, (c) L1₀ AuCu, (d) L1₂ CuAu₃, and (e) L1₂ AuCu₃. Overall, the values of $\omega$s calculated by using the PBE functional seem to be smaller than those calculated by using the PZ. For simple metals, the PZ and PBE results agree with the experiments of Au lynn and Cu nicklow , respectively. This tendency is similar to the previous calculations tang . For AuCu the lowest $\omega$ at point R is sensitive to the functional used: $\hbar\omega=-2.3$ and $1.8$ meV for PZ and PBE, respectively. For CuAu₃ the phonon softening at the point M is observed. For AuCu₃, the agreement between the calculated and experimental curves katano is good when the PBE functional is used. Below we will use the PBE functional to study the dynamical stability of AuCu in detail.

Figure 3 shows the phonon dispersion curves of (a) Au, (b) Cu, (c) AuCu, (d) CuAu₃, and (e) AuCu₃ in the WDM regime, where $T_{\rm e}$ is set to 0, 4, and 6 eV. For both simple metals (Au and Cu), the phonon energy increases monotonically (i.e., phonon hardening) with $T_{\rm e}$, which is consistent with the previous calculations recoules ; yan . For CuAu₃ and AuCu₃ in the L1₂ structure, the phonon hardening behaviors are also observed when $T_{\rm e}$ is increased, while nonmonotonic increases are observed along the line of M-R for the transverse acoustic phonon branch. As shown in Fig. 3(c), L1₀ AuCu becomes unstable against $T_{\rm e}$. While the optical phonons are hardened significantly, the transverse acoustic mode is unstable because around point R the phonon energy is imaginary.

In general, the phonon hardening is attributed to an increase in the internal pressure caused by weakened electron screening. When the screening effect is weakened, bared ions are created, leading to a peak shift and a shrinkage of the width in the electron energy spectrum. This has been illustrated by the electron DOS as a function of $T_{\rm e}$ recoules . Figure 4 shows the DOS of Au-Cu alloys for several $T_{\rm e}$s. At $T_{\rm e}=0$, due to the presence of $d$ bands, the high DOS is observed below the Fermi level from $-2$ to $-7$ eV, while the band width is narrow in Cu (from $-2$ to $-5$ eV). When $T_{\rm e}$ is increased, the peak position and the width become deep and small, respectively, as described above. Although the variation of DOS with $T_{\rm e}$ may be correlated with the phonon hardening of Au, Cu, CuAu₃, and AuCu₃, it never explains the phonon softening observed in L1₀ AuCu.

To understand the anomalous softening of AuCu in the WDM regime, let us remind that the ion-ion interactions are decomposed into the direct and indirect parts: the former is the repulsive Coulomb interaction that is independent of the magnitude of $T_{\rm e}$, while the latter is the electron-mediated attractive interaction that is significantly modified in the WDM regime. At ambient condition, the former and the latter may be partly canceled with each other, producing the repulsive and attractive interactions for the short- and long-range parts, respectively. Given small indirect contribution, the repulsive Coulomb forces are dominant for the short-range part. On the other hand, for the long-range part, the attractive forces compete with the repulsive forces, creating relatively small interaction forces. This situation is exactly realized in WDMs. We thus hypothesize that the lowest energy mode at point R of L1₀ AuCu is stabilized by the long-range interaction between atoms at $T_{\rm e}=0$ eV. With this assumption, it is reasonable to observe the softening of the phonon mode when L1₀ AuCu is in the WDM regime having weakened long-range interatomic interactions. Similarly, we can expect that when a phonon mode is stabilized by only the short-range interactions, such a mode will show hardening behavior in the WDM regime.

To demonstrate this concept, we list the values of $\omega_{k}$ in Eq. (16) for several models accounting for up to the $p$-th NN interactions ($p$NN models) in Table 2. As $p$ increases, the phonons at point R are stabilized and their energies approach the DFPT results. It is noteworthy that the sizes of $\omega_{\rm V}$ and $\omega_{\rm VI}$ within the 1NN model are almost comparable to those within DFPT. This means that the optical phonon energies are well determined by only the short-range interaction between atoms. When the interatomic interactions are considered up to the second order NN, $\omega_{\rm III}$ becomes positive. However, $\omega_{\rm I}$ is positive only within 5NN and DFPT. This implies that the long-range forces between Au and Cu atoms are important to yield positive $\omega_{\rm I}$, i.e., the dynamical stability of L1₀ AuCu. Table 3 also lists the values of $\omega_{k}$ but for $T_{\rm e}=4$ eV. As in Table 2, $\omega_{k}$ approaches the DFPT results with the inclusion of higher order NN interactions. Again, the sizes of $\omega_{\rm V}$ and $\omega_{\rm VI}$ within the 1NN model almost agree with those within DFPT, indicating that the optical phonon hardening is due to an enhanced short-range forces. The most important fact is that the convergence of $\omega_{\rm I}$ is fast compared to the case of $T_{\rm e}=0$ eV: $\omega_{\rm I}=-5.8$ meV within 3NN is already similar to $-5.1$ meV within DFPT. This reflects the short-range nature of WDM, leading to the instability of the $\omega_{\rm I}$ mode.

The stability of AuCu in the L1₀ structure can also be understood by calculating the free-energy as a function of the ratio $c/a$ with the equilibrium volume $V=V_{0}$ fixed (i.e., along the Bain path grimvall ): $c/a=1$ and $\sqrt{2}$ correspond to the bcc (or B2) and fcc structures, respectively. Figure 5(a) shows the $c/a$-dependence of the free-energy of AuCu for several $T_{\rm e}$s, where $V=V_{0}$ is the equilibrium volume of L1₀ at $T_{\rm e}=0$ eV. At $T_{\rm e}=0$ eV, the free energy takes minimum values at $c/a=1.26$ (L1₀) and $0.9$, whereas that takes a maximum value at $c/a=1$ (B2). As expected, AuCu in the B2 structure is dynamically unstable (see Fig. 5(b) (blue)). When $T_{\rm e}$ is increased, the profile of such a double-well-like curve is modified, indicating a non-thermal solid-to-solid transformation, as discussed in warm dense tungsten giret . At $T_{\rm e}=2.5$ eV, the free-energy minimum at smaller $c/a$ vanishes, while that at larger $c/a$ shifts to $c/a\simeq 1.22$; at $T_{\rm e}=3$ eV, the free-energy curve becomes relatively flat in the region of $c/a\in[1.0,1.2]$ and takes a minimum only at $c/a=1$; and at higher $T_{\rm e}$s, the curvature around $c/a=1$ becomes larger, stabilizing the B2 phase. Such a B2 structure is found to be dynamically stable, as shown in Fig. 5(b) (orange), where $T_{\rm e}=4$ eV is assumed. We hope that our prediction for the phase transformation in AuCu, before melting occurs, can be confirmed by future experiments.