In the Kohn–Sham–Mermin scheme [54, 55], the occupation numbers of electrons follow the Fermi–Dirac distribution, and the electrons are considered to be in thermal equilibrium. When the electrons were excited by laser, the time scale of energy transfer between the electrons and the ions is much slower than that of the relaxation process of the electron subsystem. Thus, the electrons and the lattice can be considered to be in a transient state of nonequilibrium before the electron–phonon equilibration is completed. This justifies the application of the two-temperature model, and the electrons can be described by FTDFT with fixed ionic configuration according to Born–Oppenheimer approximation. There is yet another nonequilibrium process involved because the electron subsystem is also away from equilibrium. In fact, only electrons with specific energy can be photoionized during isochoric heating by the XFEL. To calculate such a complex nonequilibrium system, we use the plane-wave pseudopotential open-source package QUANTUM ESPRESSO (QE) [56, 57] with minor modifications to realize nonequilibrium distribution of electrons in the FTDFT calculations. The occupation number in the relevant inner shell Kohn–Sham orbitals is set as zero in every SCF step to create holes, and the same number of electrons as that of holes was added to the chemical potential to study the nonequilibrium transient state after laser excitation so that the total number of electrons remains unchanged in the system. This could be realized experimentally by carefully choosing the XFEL energy. The new chemical potential is decided self-consistently like in the original SCF process, which also determines the occupation numbers of all orbitals in return, except for those that were forced to be zero. In this way, we can go beyond the frozen core approximation and create holes directly in our FTDFT calculations for real materials, contrary to the common practice of creating holes in atomic calculation during the generation of pseudopotentials and treating the corresponding atom as an impurity in the FTDFT calculations [58, 59].

The DFPT and finite displacement method for calculating phonons are by now two different well-established tools for investigating vibrational properties. Here, we use both the DFPT and the finite displacement method to calculate the phonon dispersion relations of two metals, Al and Cu, at room temperature for benchmark. For nonequilibrium calculations, only the finite displacement method is used due to its straightforwardness.

For DFPT, the calculations are performed using the QE code; the pseudopotential takes the projector augmented wave (PAW) formalism [60] and the exchange-correlation functional takes the Perdew–Burke–Ernzerhof (PBE) parameterization [61]. For Al, we use a pseudopotential including both L-shell and M-shell electrons (i.e., 2s²2p⁶3s²3p¹) as valence electrons. For Cu, we use the pseudopotential which considers both M-shell and N-shell electrons (i.e., 3s²3p⁶3d¹⁰4s¹) as valence electrons. The plane wave cutoff energy is set to 170 Ry, a 16 × 16 × 16 Monkhorst-Pack K-point mesh [62] is used for the sampling of the Brillouin zone, and the grid of q-point of lattice dynamical matrices is 4 × 4 × 4. The lattice constant is relaxed to be 7.639 Bohr for Al and 6.848 for Cu. The convergence with respect to plane-wave cutoff energy, k-point, and q-point meshes is carefully checked for all our calculations.

For the finite displacement method, the calculations for Al and Cu are performed on a 3 × 3 × 3 supercell (108 atoms) using the ALAMODE [63, 64] as a postprocessing code for QE using the same pseudopotential as DFPT. Before isochoric heating with the XFEL, the samples are supposed to be at ambient conditions, which mean that the electron temperature of 300 K is used for the metallic system. The benchmark between the DFPT and finite displacement method for such a case is shown in Figures 3A,B (green triangles versus black curves) with perfect agreement. The advantage of the finite displacement method is that we can relate the calculations of electron excitation from multiatomic system to the excitation of electrons of metals in the experiments. As a result, we only use the finite displacement method to calculate phonon dispersions for photoexcitation cases rather than DFPT. The shifted Monkhorst-Pack k-point meshes are 6 × 6 × 6 for Al and 4 × 4 × 4 for Cu. The atomic displacement of 3 × 10⁻² Å and the plane-wave cutoff energy of 60 Ry are used for both Al and Cu. The convergence of plane-wave cutoff energy and the size of k-point mesh are also carefully checked in our calculations.

In the lattice dynamic calculations, we use the finite displacement method to calculate the second-order interatomic force constants and the phonon frequencies under quasi-harmonic approximation. To include the anharmonicity, we calculate the PES with excited electrons and illustrate the origin of lattice stability. The PES can be obtained by moving the atomic positions far away from their equilibrium positions along the eigenvector of different modes at specific high-symmetry points. The displacement of the j ^th atom in the l ^th unit cell along a phonon mode amplitude Q is denoted by d _jl and defined as d j l = Q N a m j R e exp ( i Φ ) e j ⁡ exp ( q ⋅ r jl ) , (1) where Q is the normal-mode coordinate, N _a is the total number of atoms in the supercell, and m _j and r _jl are the mass and position of the j ^th atom in the lth unit cell, respectively. e _j is the j ^th component of the unit normal coordinate eigenvector. Φ is an arbitrary phase factor and is set to zero by default. If the expansion of the supercell is commensurate with the wavevector q, the total energy should not change when a nonzero Φ is used because the displacement corresponds to a part of the wave period [65]. The unit normal eigenvector can be obtained from the QE code, and then, we can manage to move the atoms accordingly. The total energy can be calculated by FTDFT with new atomic positions, and then, we can obtain the PES by summarizing the data along different directions.

The absorption of photon by isochoric heating with the XFEL can change the electronic structure of metals to a transient exotic state. Experimentally, the XFEL which has tunable photon energy with pulse duration on the order of 10 fs and high intensity as much as 10¹⁷ W/cm² can heat and excite the electrons of metals and create a nonequilibrium electronic structure [66]. The L _I edge of Al is 117.8 eV and the M _I edge of Cu is 122.5 eV, both of which are in the soft X-ray range of an XFEL [67]. The specific holes can be created by predominant photoionization of 2s electrons of Al or 3s electrons of Cu to their Fermi surface, with other transitions being of minor importance. The electrons which are excited to higher energies are also found to relax to the chemical potential within a few femtoseconds [66]. Thus, this transient state could exist in a time scale longer than the time scale of the pulse and shorter than the time scale of Auger decay.

With the excitation of valence electrons in the inner shell, the reduced screening of core electrons [52] results in localization of the valence band. The DOS of electrons in such nonequilibrium transient states is plotted in Figures 1A,B, and we see similar change of the DOS as in [52], which shows small change for Al but considerable localization of valence bands for Cu. This is mainly due to the semilocal d band in Cu. The black curves show the calculated DOS without a hole, the red and blue ones, respectively, represent the DOS with a hole and two holes of every atom in the system. Due to the deeper Coulomb potential which the remaining electrons feel [52], it is clear that these atom-like features shift toward lower energy. The rise of chemical potential is also observed when the excited electrons are increased. In Figures 1C,D, the local density of states (LDOSs) is calculated by projecting the one-electron eigenstates onto the local atomic orbitals, which illustrates the contribution of different orbitals to the electronic structure of a specific atom.

To show the charge transfer process in real space, we plot the radial electron density for Al and Cu, as shown in Figures 2A,B, respectively. The black curve corresponds to the atoms without holes, and the red curve represents every atom with a hole in the system. For Al, the electrons of the 2s orbital are excited to the conduction band. As a result, the electron density is reduced near the nucleus and increased away from the nucleus, showing an increase in ionization degree. For Cu, the electrons of the 3s orbital are excited to the conduction band, and a similar trend is observed. This picture is in agreement with the shift of the DOS in Figure 1.

First, we calculate the phonon dispersion with the DFPT and finite displacement method using Al and Cu as illustrating examples. As shown in Figure 3, the result calculated at ambient condition is in perfect agreement with the experimental measurements by using the abovementioned two methods.

As shown in Figures 3A,C, the imaginary frequencies of phonon dispersions of Al can be observed along the direction of k-path Γ-X-K-Γ-L in the Brillouin zone of the FCC lattice [68] when there are more than four holes in the 2s orbital, which suggests that the lattice becomes unstable. In our calculation, we increase the number of excited electrons one by one until the electrons of the 2s orbital are all excited in the system. We found that the change of the system from the state of stable to unstable will happen with the increased excited electrons and the stability of lattice cannot recover when the lattice symmetry is satisfied again. For instance, the state of the system is still stable when the four electrons are excited; however, it becomes unstable when five or more electrons are excited. This suggests a process of lattice collapse for Al. It is to be noted that due to the incompleteness of our model, this picture is qualitatively right, but the number mentioned here is not quantitatively accurate.

For Cu, however, as shown in Figures 3B,D, we found that the lattice becomes unstable when the symmetry is broken, and the lattice becomes stable again when the symmetry is satisfied. We also found that the change of the system from the state of stable to unstable will happen when we increase the number of excited electrons of the 3s orbital gradually. For instance, the state of the system is stable when the two electrons are excited; however, it becomes unstable when three or more electrons are excited. It is possibly related to the symmetry breaking of the perfect lattice since the holes are randomly located within the atom in the sample. However, it is surprising that the lattice stability recovers when each atom excites one electron or two electrons from the 3s orbital. It may suggest that the Cu lattice becomes more stable when laser fluence is increased, which corresponds to the increased phonon frequencies. We suspect that the lattice might undergo a process of phase transition in this case. As for Cu, the phenomenon of phonon hardening in all three acoustic branches is observed significantly and which means the perfect lattice becomes harder to break when the inner shell electrons of every atom are excited. It is very different from aluminum.

Phonon dispersion with quasi-harmonic approximation might not be persuasive, especially for cases far away from equilibrium. Focusing on the crystal instability and phonon hardening [46, 51, 65, 69–71] cases, we plotted the potential energy surface (PES) of Al and Cu in the normal coordinates for the transverse acoustic (TA) and longitudinal acoustic (LA) modes at the high-symmetry point L ( 1 2 , 1 2 , 1 2 ) in the first Brillouin zone to include anharmonicity effect. The scaling factor Q ₁ of the lowest TA mode and the scaling factor Q ₂ of the LA mode are applied in the calculation, and the PES is a function of scaled coordinates E = E(Q ₁, Q ₂); we do not show the second-lowest TA mode because it is degenerate with the lowest TA mode. The finite displacements are imposed along the corresponding mode eigenvectors so that we can obtain the shape of the PES around the equilibrium position. To compare the curve of the potential surface, we both set the potential energy as zero in the equilibrium position.

The PES is plotted in Figures 4A,B, and the PES of Al without a hole and with a hole in the 2s orbital is represented, respectively. In Figure 4A, the PES of atoms which vibrate along the corresponding eigenvector in the system is plotted, and the potential energy curve of Al along the Q ₁ direction and Q ₂ direction is plotted using the black lines in Figures 4C,D. In Figure 4B, the PES of atoms with a hole is also plotted, and the potential energy curve of Al along the Q ₁ direction and Q ₂ direction is plotted using the red lines in Figures 4C,D. We can see that the slope of the TA mode is lower than that of the LA mode in the two cases by associating the PES and potential energy curve along different directions with the phonon dispersions, which is shown in Figure 3C. The transition from a single-well to double-well potential along the Q ₁ direction and the phonon softening along the Q ₂ direction shown in Figures 4C,D indicates probable softening of chemical bonds or displacive phase transition. So, we can know that the dynamic instability of aluminum originates from the fact that it is the local maximum of the double-well, which results in the negative curvature ∂ 2 E ∂ Q 1 2 < 0 and the eigenvalue of the TA mode ω ² is also negative with the excitation of electrons in the inner shell of aluminum. Meanwhile, the potential energy curve of the LA mode along the vibration direction becomes much smoother when electrons are excited to the conduction band, while for Cu, as shown in Figure 5, the phonon hardening originates from the steeper potential energy curve along both Q ₁ and Q ₂ direction. From the PES of Cu, we can intuitively observe the increased trends of energy of both the LA and TA mode in three dimensions, corresponding to the increased frequencies in phonon dispersions shown in Figure 3D. The evidence indicates the reduced screening of orbitals and an increase of the effective ion–ion potential, which results in the localized electron of the d orbital being closer to the nucleus in real space and phonon hardening as a whole.