Edinburgh Research Explorer Electronically driven 1D cooperative diffusion in a simple cubic crystal

Atomic diffusion is a spontaneous process and significantly influences properties of materials, such as fracture toughness, creep-fatigue properties, thermal conductivity, thermoelectric properties, etc. Here, using extensive molecular dynamics simulations based on both ab initio and machine-learning potentials, we demonstrate that an atomic one dimensional cooperative diffusion exists in the simple cubic high-pressure finite-temperature phase of calcium in the premelting regime, where some atoms diffuse cooperatively as chains or even rings, while others remain in the solid state. This intermediate regime is triggered by anharmonicity of the system at high temperature and is stabilized by the competition between the internal energy minimization and the entropy maximization, and has close connections with the unique electronic structures of simple cubic Ca as an electride with a pseudogap. This cooperative diffusion regime explains the abnormal enhancement of the melting line of Ca under high pressure and suggests that the cooperative chain melting is a much more common high-temperature feature among metals under extreme conditions than hitherto thought. The microscopic electronic investigations of these systems combining ab initio and machine-learning data point out the direction for further understanding of other metallic systems such as the glass transition, liquid metals, etc. DOI:


I. INTRODUCTION
Melting as a macroscopic phenomenon is ruled by temperature and pressure, with competition between the internal energy and entropy to minimize the free energy of the system. Several unusual states induced by this competition show intriguing behaviors. For instance, in a superionic state [1][2][3][4][5][6][7][8][9], part of the atoms form a solidlike fixed sublattice, while the other part loses long-range order and diffuses freely as in liquids. Multiple superionic and plastic states in helium compounds have also been found within a large pressure range in recent studies [10][11][12]. Another example is the "chain-melting" behavior in highpressure phases of some elemental metals such as potassium and rubidium and in some compounds at ambient conditions, all of which feature incommensurate host-guest (HG) structures [13][14][15]. These HG structures consist of two sublattices with 1D guest chains located in the interstitial space of the host structure. In experiments, Bragg peaks of the guest sublattice disappear on heating, which is a signature of chain melting [13][14][15][16]. Theoretical simulations revealed the loss of intrachain and interchain correlations during the chain-melting process, establishing the "chain melt" as a true thermodynamic state [17].
In bcc Cu, a synchronized cyclic motion of four atoms (four-atom ring mechanism) was first proposed by Zener in 1950 [18] and collective atomic motion was seen in classical potential (embedded-atom method) simulations, as well as in Lennard-Jones [19,20]. In bcc Ta, a shearinduced anisotropic plastic flow was studied using classical molecular dynamics (MD) [21]. Recently, ab initio MD has been applied to a low-viscosity phase of bcc Fe at Earth's inner core conditions [22,23] and bcc Ti at high temperatures [24]. Similar cooperative diffusion was also seen in these simulations. In all cases, unusual yet very specific quasisolid or partially molten states are the result of an interplay of crystalline symmetry and a shifting balance between internal energy and entropy at elevated temperatures. However, it remains unclear how general this phenomenon is. We show here that even the arguably simplest possible crystal structure, a simple cubic (sc) system, can support such a chainmelting regime. We tie the presence of this regime to the peculiar electronic structures of sc Ca.
The group II alkaline-earth-element calcium (Ca) plays an important role in a wide variety of applications, especially in medical and chemical industries. Under extreme conditions, Ca exhibits a variety of complex and interesting behaviors. Previous high-pressure experimental studies have investigated a series of solid phase transitions in Ca [25][26][27][28][29] [30,31].
In recent years there has been an increasing interest in the simple cubic phase Ca III. Phase transitions from Ca I to Ca II then to Ca III are accompanied by a decrease in coordination, which is counterintuitive from a close-packed sphere concept. Even though the sc phase has a lower enthalpy than the bcc and fcc phases in its stable pressure range, several theoretical calculations showed the zerotemperature phonon instability of sc Ca and suggested a number of structures (e.g., Cmmm and I4 1 =amd) with a lower zero-temperature enthalpy than that of the sc phase [32][33][34][35][36]. Experimental attempts have been made to observe these low-temperature phases, and anharmonic and/or entropy effects have been suggested to be the reason for the stability of the sc phase at room temperature [37][38][39].
Theoretical calculations proposed two post-bcc phases with orthorhombic Cmcm and Pnma structures instead of the sc structure, and the solid-liquid phase transitions at high temperatures have also been explored in simulations [40]. A few experiments have been performed to study the high-pressure and high-temperature behavior of Ca up to 80 GPa and 3000 K [41][42][43]. An unusual melting curve of Ca has been reported, which stays flat at T m ∼ 1500 K between 5 and 32 GPa. Very interestingly, above 32 GPa, the melting temperature increases dramatically until it reaches another plateau at about 45 GPa and around 2750 K. A recent experiment that combined heated diamond-anvil cells and synchrotron x-ray diffraction has established phase boundaries between the Ca I, Ca II, and Ca III phases and the stability of the sc phase Ca III up to 50 GPa and 800 K [44].
Over the past decades, most studies in Ca focused on the phase transitions at low temperature and room temperature while little attention has been paid to the high-temperature behavior of Ca. In particular, the abnormal increase in the melting temperature of sc Ca around 30-45 GPa still needs to be elucidated, and it is highly needed to know whether there is a high-temperature phase above the sc Ca. In this paper, we report extensive molecular dynamics simulations with both ab initio and machine-learning (ML) potentials to explore these issues and reveal the existence of a premelting region with 1D cooperative diffusion along chains and rings between solid simple cubic calcium and the liquid state. This unexpected behavior is reported in a high symmetry simple cubic system for the first time and suggests that chain melting is a widespread feature among elemental solids.
We use the Vienna ab initio simulation package (VASP) [45] to perform ab initio molecular dynamics simulations and electronic structure calculations. The machine-learning molecular dynamics was performed by employing LAMMPS [46] together with a machine-learning neural network potential for calcium trained by the deep potential molecular dynamics package (DeePMD) method [47]. The intrinsic properties including free energy, electronic structures, and barrier of cooperative motion related to chain melting are discussed using several advanced methods. More extensive details of calculations and the validity of parameters can be found in Sec. VI.

II. RESULTS
The unique and complex properties of Ca usually demands a first-principles description. Therefore, we have carried out ab initio molecular dynamics (AIMD) in a cubic supercell with 125 atoms to investigate the atomistic behavior of sc Ca at high temperature. Results from AIMD simulations at 50 GPa along different temperatures (900, 1300, 1600, 2100, and 3000 K) are presented in Fig. 1. Mean-square displacements (MSD) of the Ca atoms reveal three states of sc Ca with different diffusion coefficients D: solid at 900 K (D 1 ≤ 2 × 10 −12 m 2 =s), intermediate at 1300 K, 1600 K (D 2 ¼ 0.137 × 10 −8 m 2 =s, D 3 ¼ 0.305 × 10 −8 m 2 =s), and liquid at 2100 K, 3000 K (D 4 ¼ 1.308 × 10 −8 m 2 =s, Fig. 1(a)], where the diffusion coefficients of liquid are one order of magnitude larger than in the intermediate states.
The vibrational density of states (VDOS) FðνÞ obtained from the Fourier transform of the velocity autocorrelation function (VACF) of AIMD trajectories shows the difference between the VDOS and the vibrational states at zero frequency Fð0Þ of the three different states, where ðkT=mÞFð0Þ is the selfdiffusion coefficient. As shown in Fig. 1(b), the shapes of the VDOS and the values of ðkT=mÞFð0Þ immediately divide the five representative trajectories into three different groups. The diffusion coefficients calculated by this method are very similar to the numbers from the slope of the MSD Interestingly, the trajectories of simulations at 1300-1700 K show a difference in the atomic motion behavior along the three directions ðx; y; zÞ. MSD from the simulation at 1300 K is shown in Fig. 1(c). More specifically, the total nonzero diffusion coefficient D 2 is mostly contributed from atomic motion in the x direction, while there is no diffusion in the y and z direction. Since the three xyz directions in sc Ca are equivalent, we found that chains in sc Ca can diffuse randomly along h100i in different simulations run at the same conditions. The VDOS projected onto each direction in sc Ca at 1300 K illustrates there is no significant difference in the frequency peak positions along the x direction compared to y and z [see Fig. S1(a) in Supplemental Material [48]]. However, the zero (nonzero) vibrational DOS at zero frequency confirms the fixed (diffusive) behavior along the different directions. These results indicate that an intriguing diffusion along the cubic axis directions appears at finite temperature in sc Ca without changing its eigenfrequency modes. The average structure extracted from MD trajectories over a 2 ps window also shows the persistence of the simple cubic framework [see Fig. S1(b) in Supplemental Material [48]]. This 1D chainlike diffusion regime can also be seen in AIMD simulations at 50 GPa and 1400 K of a supercell initiated with a Frenkel defect.
In order to reveal atomic details of this phenomenon, we examined the movement of atoms in different directions at 1300 K and 50 GPa. Further analyses showed atomic motion patterns can be divided into two types: one group of atoms vibrating around their equilibrium positions (solid state), and the other one diffusing along the h100i directions (diffusive state). Focusing on a [100] case as before, the diffusing atoms form a chain that extends across the entire supercell. We distinguish atoms of different diffusive behavior with different colors in Fig. 1(d). The MSD along the chain's direction (cMSD) and the VDOS in the x direction also show obvious differences [Figs. 1(e) and 1(f)]. From the x displacement of atoms over time, it can be seen that the atoms in the diffusive chain hop in a concerted manner between equilibrium lattice sites on a picosecond timescale, while the solidlike atoms only oscillate around their original equilibrium positions [see Fig. S2(a) in Supplemental Material [48]]. The sliding of different atoms on the chain does not occur simultaneously, but the time difference between them is on the femtosecond scale. The migration of the first atom from a lattice site to an interstitial site caused by thermal fluctuation promptly induces a Frenkel defect which triggers the surrounding atom to form a chainlike motion. As the temperature increases, chains of different diffusion directions coexist and can even exchange atoms. For example, at 50 GPa and 1500 K, a chain in the y direction slides into a new equilibrium position at t ¼ 2.5 ps, and then another chain containing an atom moving from the adjacent y position slides to the z direction within 1 ps [ Fig. 1(g)]. The specific change of the atoms' displacement over time in these two directions can be found in Fig. S2(b) in Supplemental Material [48].
The number of melted chains and the crossing between the chains increase with temperature. The chain-melting onset inside the sc solid phase is gradual, unlike the firstorder solid-liquid phase transition. We quantified the proportion of atoms in the diffusive state at different temperatures by using a two-phase thermodynamics with memory function (2PT MF) model [49,50]. This model separates the system into a solidlike (vibratory) components ð1 − f g ÞF s ðνÞ for which the VDOS goes smoothly to zero at zero frequency and a highly gaslike (diffusive) f g F g ðνÞ state where diffusion makes migration calculations is ∼0.8 eV, which is evaluated via Arrhenius linear regression, where D(T) is the diffusion constant in a defect-free lattice at temperature T and D 0 is the exponential prefactor.
In addition, as shown in Fig. 3, we found that the chainmelting regime can be distinguished from the solid and liquid by pair distribution function gðrÞ. When heating from 900 to 1500 K, as sc Ca evolves from solid into the chain-melting regime, the third and fourth peaks of gðrÞ are combined and the sixth peak of the solid vanishes, but it is apparent that the chain-melting regime still retains some of the solid properties, especially the existence of the second peak in gðrÞ. Upon further heating to 2700 K, the system completely turns into liquid, where the second peak disappears and gðrÞ approaches 1.0 at large r.
Since the chains are commensurate with the surrounding lattice and diffusion of infinite straight chains is impossible in reality, periodic boundary conditions might spuriously introduce it in a small size simulation. To rule that out, simulations with larger lattice scale are required. AIMD simulations based on density functional theory (DFT) have high accuracy but are computationally too costly for largescale simulations. On the other hand, machine-learning force fields developed rapidly in recent years due to a nice balance between accuracy and computational cost [47,[51][52][53]. ML potentials have been successfully applied to many different types of systems [54,55].
Here we build a machine-learning potential surface of Ca using DeePMD [47], which has been reported can reach the accuracy of DFT in some systems [47,55]. We use AIMD data as training and test dataset. Details of ML potential training are described in Sec. VI and the Supplemental Material [48]. The root mean square errors are 1.68 meV=atom for the energies, 146 meV=Å for the forces on the testing set at 50 GPa, and 1.34 meV=atom, 123 meV=Å for the dataset at 40 GPa. Comparisons between the pair distribution function at some representative temperatures computed by the ML potential and AIMD show very good agreement (see Fig. S3 in Supplemental Material [48]).
To gauge the potential finite size effects in the small supercells of the AIMD simulations, we performed MD with the ML potential using 20 × 20 × 20 (8000 atoms) and 25 × 25 × 25 (15625 atoms) supercells, respectively, and a simulation length reaching the nanosecond scale. To describe the movement of atoms more obviously, we define a flux function F x ðy; zÞ to qualify the chainlike melting (see Sec. VI for details). F x ðy; zÞ represents the net directional atom flux density through all planes parallel to the [100] direction accumulated over a certain time range. Similarly, we can define F y ðx; zÞ and F z ðx; yÞ to investigate the chain melting along the other two directions.
Typical results of flux functions are presented in Figs. 4(a) and 4(b) accumulated over t ¼ 100 ps. Chain melting along the cubic axes can still be seen to exist in large supercells, symbolized by "hot spots" in F and atom flux in both positive [100] and negative [100] direction. We found that the length of the longest chain is shorter than the diameter of these large supercells, which suggests that chain diffusions with infinite length indeed do not exist. Instead, some adjacent chains with the same or opposite diffusion directions are observed. This represents chains that will make an L turn and form rings (not necessary in a 2D plane) on a long timescale [ Fig. 4(c)]. These behaviors were not observed in shorter AIMD runs with smaller supercells. More details can be found in the movies in the Supplemental Material [48].

III. DISCUSSIONS
Using DFT and ML molecular dynamics simulations, we have found that sc Ca should become diffusive at high temperature, where it gradually crosses over to a cooperative diffusive regime prior to full melting. The reason why nature should choose such an unusual state is still missing. In the following, we provide more insight into this question.

A. Energy barriers
The migration barriers of concerted migrations of Ca atoms along the h100i edge, the h110i face diagonal and the h111i body diagonal directions were calculated employing the nudged-elastic-band (NEB) method [56] based on ab initio computation. As shown in Fig. 5(a), the barrier along the h100i direction is only 0.18 eV per chain atom, while the other two directions h110i; h111i are 3.47 and 3.6 eV, respectively. Interestingly, moving atoms by 0.1 unit cell lattice offset along the h100i direction will even reduce the energy of the system, by about −1.2 meV per chain atom. This is consistent with what is discussed about the total energy per atom when the atoms are displaced along the phonon mode in the Ref. [35] and can be connected to the anharmonic effects of sc Ca at low temperature [33], which leads to the known lowtemperature I4 1 =amd and Cmmm phases [36][37][38][39]. Meanwhile, the vacancy formation enthalpy in sc Ca at 50 GPa is about 2.67 eV, which is larger than the migration barrier of cooperative motion (0.18 eV) and suggests the preference for atomic cooperative motion. We also find that the collective motion along other directions (such as h111i) triggered by a vacancy defect and assisted by atomic movement along h100i and its potential barrier can be the same level with directly diffusing along h100i. More details can be seen in Figs. S5 and S6 [48].

B. Lattice dynamics at high temperature
A previous study discussed that lattice dynamics at high temperature may play a role in this kind of atomic cooperative diffusion [24]. The softening of the phonon frequency below the melting triggers the atomic cooperative motion. To demonstrate this link in sc Ca, we calculated the phonon spectrum at different temperatures (0, 600, 1200 K) [ Fig. 5(b)]. Our results show that at 600 K the lattice instability of sc Ca is fully eliminated by anharmonicity. This temperature-induced stability is similar to the results in the previous studies [34,35]. However, as temperature increases to 1200 K, there is a resoftening of the phonon before the atomic cooperative motion. Therefore, the anharmonicity at high temperature may trigger the atomic cooperative motion in sc Ca. Such an intriguing behavior is important to help us to understand the relationship between the microscopic atomic behaviors and its lattice dynamics at high temperature.

C. Electronic structures: Electride and pseudogap
The above calculations highlight the low barrier of cooperative motion and its trigger of high-temperature lattice dynamics. However, how the electronic structures influence the energy surface and the lattice dynamics at high temperature is elusive. We calculated the temperature-dependent electronic distribution and density of states to verify that this atomic cooperation motion is electronically driven.
A recent study revealed two qualitatively different contributions to the electride character of the HG alkali metal phases, which are (i) fully localized basins and (ii) 1D partially localized electrons along channels in the host structure [57]. Bader charge calculations [58] of sc Ca at 50 GPa show that some electrons localize in the interstitial space with an integrated charge value of 0.61e (Table I). In addition, our electronic localization function (ELF) calculation of sc Ca at 50 GPa (shown in Supplemental Material Fig. S7 [48]) also confirmed the electride picture. Therefore, sc Ca at 50 GPa is actually an electride. These localized electrons form a kind of anionic sublattice. And for sc Ca with a displaced chain of atoms, the ELF maxima near the chain indicate there is a tendency to form 1D localized electrons which least affects the other pseudoanions in the body diagonal direction. [Fig. 5(c)]. This gives the 1D cooperative motion an energetic preference and the low barrier of motion.
To corroborate the former, we analyze in detail the electronic structures of the different regimes. As clearly shown in Fig. 5(d)  obviously important for its stability. However, there is a filling of the pseudogap at the Fermi level as temperature increases and the pseudogap exists in both the solid and chain-melting regime, but not in the liquid. The orbital nature of the electronic states near the Fermi level is d orbital, for all states including solid, chain-melting state, and liquid. The occupied states of solid sc Ca are mainly contributed from the d xy , d yz , d xz orbitals, and the unoccupied states are formed by the d x 2 −y 2 and d z 2 orbitals (see Fig. S8 in the Supplemental Materials [48]). In the liquid state, all orbitals are mixed together and contribute equally to both the occupied and unoccupied states. And the situation in the chain-melting regime is just in between those extrema. Therefore, the persistence of the pseudogap demonstrates that diffusive chains do not raise the electronic energy of the system significantly. The resulting internal energy stabilization relative to the liquid suggests itself as the main mechanism for the enhanced melting temperature, thus allowing in principle an entropy gain to be realized in the intermediate chain-melting regime. Meanwhile, the increase of electron density at the Fermi level NðE f Þ screens the atomic forces and induces the resoftening of the lattice dynamics [59]. Here we show a clear interplay between the electronic structures, anharmonicity at high temperature and 1D cooperative motion in sc Ca.

D. Entropy
As the temperature increases in a system with constant volume, the vibrational amplitude of atoms and displacements from their equilibrium positions become larger, which increases the internal energy U. Meanwhile, phonon excitations and diffusion in the system at finite temperature increase the disorder and thus the entropy S. The competition between U and S leads to a state with a minimum of the Gibbs free energy G. Their relationship is written as G ¼ U þ PV − TS. In our MD simulations with NVT ensemble, for a sc Ca structure optimized to the specific pressure, the fluctuation of P in its equilibrium state at different temperatures is less than 1 GPa. Detailed data can  be found in the Supplemental Material [48]. This allows us to determine the entropies and Gibbs free energies of the system from 600 to 2700 K by applying the 2PT MF models [49,50,60,61] along an isobar of 50 GPa. From the above calculations of the diffusion energy barriers in sc Ca, we discover that the 1D chain diffusion only contributes a tiny amount of energy to U, while at the same time the entropy change induced by this disorder reduces the free energy by TS. However, as already discussed above, the persistence of the pseudogap and electride also demonstrates that diffusive chains do not raise the electronic energy of the system significantly. Therefore, entropy gains, see

E. Elastic constant
There are very few elements in nature possessing the simple cubic crystal structure. As far as we know, only calcium, polonium, and phosphorus belong to this class. A previous study investigated the strong elastic anisotropy of sc Po by evaluating the elastic constants [62]. In their calculations, the elastic anisotropy factor A ¼ C 44 =C 0 [where C 0 ¼ 1=2ðC 11 − C 12 Þ of sc Po is no larger than 0.17, which is a very small value among solids. Therefore, we calculated the elastic constants and A of sc Ca as well as some selected structures for comparison. As listed in Table II, the elastic anisotropy factor A for sc Ca at 40 GPa is only 0.046. To the best of our knowledge, there is no other solid that has such a low value of the elastic anisotropy factor. Such low elastic anisotropy indicates that the angular forces in the first coordination shell, which consists of six atoms in simple cubic Ca, are very small. This allows the shear along the h100i direction to be very easy.

F. Phase diagram
We studied the melting behavior of Ca from ambient pressure to near 60 GPa starting from different phases using extensive AIMD simulations. The results are summarized in Fig. 6. Our simulations successfully reproduced the experimental melting line of bcc Ca, which has unique features including a small increase at low pressure and subsequent flat line at high pressures [42,43]. We also successfully located the bcc-fcc phase boundary close to experiments. The positive outcomes of these tests gave us more confidence in our method. We then applied the same method to sc Ca at higher pressures. Excitingly, we revealed an intermediate regime with cooperative diffusion through chains of atoms between the solid sc Ca and the fully liquid state. In Fig. 6, the PT range of this chainmelting regime is shown with gradient color. The crossover boundary between the cooperative diffusion state and solid sc Ca has a negative dT=dP slope. The melting line of sc Ca, on the other hand, has a positive slope so that a crossover region exists between hard solid sc Ca and the liquid state. The cooperative diffusion regime fills that region with a positive dT=dP slope to the liquid phase, which agrees with the experimental melting line. Note that the slope of the theoretical melting line is somewhat less steep than that of the experimental melting line. It should be pointed out that we also saw cooperative diffusion through chains along the h111i direction at around 1400 K and 10 GPa in bcc Ca. But its temperature range is not very large, only 100 K below the melting temperature.

IV. CONCLUSIONS
Using extensive ab initio and machine-learning molecular dynamics simulations, we revealed a one-dimensional cooperative diffusion state interposed between the solid simple cubic Ca state and its liquid state in a wedgelike PT range (pressure from around 35-60 GPa and temperature from around 1200-2200 K). In this cooperative diffusion state, it was found that some Ca atoms experience atomic diffusion and exchanges, while all other atoms remain in the solid state. More interestingly, large supercell simulations using the machine-learning potential indicated that some of the atomic chains can form closed rings on a large timescale. This phenomenon looks very similar to the "cooperative ring exchange," which was originally proposed to understand the quantum Hall effect of 2D electrons [63], and the Grotthuss mechanism in water [64]. This cooperative 1D diffusion regime in calcium has a superioniclike dynamic behavior. But generally, superionicity exists in systems with at least two different ionic species which occupy different effective charges, while sc Ca is elemental and metallic.
Collective atomic motion has been seen in fcc metals in the superheated regime above homogeneous melting and in bcc iron, tantalum and titanium; most but not all of these previous studies are at atmospheric pressure [19][20][21][22][23][24]. Our results in sc Ca show that the electron DOS at the Fermi level increases with temperature and screens the interatomic forces, which induces a resoftening of phonon in the lattice dynamics at high temperature. This anharmonicity triggers the movement of atoms and the electride properties stabilize the cooperative diffusion in a simple cubic metal, over a wide temperature range below the melting point. This kind of interplay between electronic structure, anharmonicity, and microscopic atomic novel motion may explain the sudden rise of the melting line in many systems after pressureinduced phase transitions, as seen for instance in potassium due to the chain melt. [15,17]. Furthermore, by analyzing the free-energy properties of this system, we have found that the appearance of this cooperative diffusion state is permitted by the competition between an increase of the internal energy and the entropy gain provided by finite-length self-diffusing Ca chains. The Gibbs free energy of this regime is the lowest at certain pressures and temperatures when the chain entropy increase overcomes the internal energy change caused by reduction albeit without demise of crystallinity. The specific crystal structure of sc Ca therefore plays a major role to enable this unexpected behavior. The particularly small elastic anisotropy factor A of simple cubic calcium also reflects its unique properties. We constructed the PT phase diagram based on extensive AIMD simulations, which shows the trends in the melting curve consistent with previous experiments. The existence of the cooperative 1D diffusion regime emerging from sc Ca, which gains entropy without excessive energy cost, very well explains the abnormal increase in the melting line of calcium at around 32-45 GPa. In this manner the thermal activation of additional degrees of freedom associated with diffusive motion reduces the free energy of sc Ca and pushes its complete melting to higher temperature. This effect could be more pronounced in experiment due to the presence of defects that will facilitate the cooperative diffusion (see Fig. S5 [48]).

V. PERSPECTIVE
We suggest that the concept of atomic cooperative diffusion as a precursor to the proper melting first-order phase transition is quite a general phenomenon that can shape the phase diagrams of materials under combined highpressure and high-temperature conditions and the influence of the electrons on the melting mechanism, especially in some novel melting systems, cannot be neglected. We have identified the following indications that would support it in other materials: (a) low-energy barriers for collective atomic motion, associated with anharmonicity at high temperature, also helped by persistence of an initial electride electronic structure as temperature increases, and (b) a steep increase of the melt line without a distinct high-temperature phase. In sc Ca, this kind of premelting regime shares the feature of cooperative 1D chain diffusion also seen in chain-melting HG phases of rubidium, potassium, and scandium, and the bcc phases of iron, titanium, and tantalum [21][22][23][24].
Their experimental detection, beyond steep increases of the melting line, brings challenges. Note that this regime in sc Ca might already have been observed in experiment but not identified as such; additional impetus for experimental groups to detect it might be needed. The HG chain-melting states can be observed in experiments through the disappearance of a subset of certain Bragg peaks on heating [5,16,17]. According to our calculations, it is possible that diffraction peaks of sc Ca should gradually broaden upon heating but still retain a fraction of their crystalline Bragg singularities because in the chain diffusive regime there is always, like in a superionic state, a finite proportion of solid lattice left. In addition, a recent experiment that performed in situ x-ray absorption, x-ray diffraction, and Raman scattering measurements has confirmed the liquidliquid transition (LLT) in high-pressure sulfur [65], which was also investigated by ab initio molecular dynamics [66]. One of the valid experimental quantities to suggest the existence of a LLT is changes in the pair distribution function gðrÞ, which corresponds to the Fourier transform of the structure factor SðQÞ. In our study, gðrÞ of different states (solid, chain-melting, liquid) shows different behaviors. We believe this is a possible observable quantity for experiments. Moreover, dynamic properties of the atomic collective motion in this regime will also influence experimentally measured quantities such as electrical and thermal conductivity, mechanical modulus, sound velocity, thermoelectric properties, etc. We believe this will be the case not only in calcium, but also other metals under high pressure.
It has been reported that the diffusion barrier for cooperative diffusion is much less than that for diffusion of an isolated ion in lithium batteries [67]. Therefore, cooperative motion may appear more commonly than previously thought, especially in some industrial processes, such as under strain or at elevated temperature. The idea of lattice-dynamics-triggered and electronically driven (influenced by high-T electronic behavior) 1D cooperative motion can be applied to study atomic diffusion and related properties (including fracture toughness, creep-fatigue properties, thermal conductivity, and thermoelectric properties) in other metallic systems such as alloys and metal glasses at ambient conditions. This also has potential implications for materials' responses to dynamic compression, e.g., in impact or detonation events, and for rocky mineral properties under planetary mantle or core conditions. Although pure Ca is not a main element of Earth's inner core, our results show a 1D cooperative diffusion induced by high pressure. This behavior could exist in some other materials, such as iron at extreme conditions [22]. Dynamic properties of the atomic collective motion will also influence quantities which govern the interiors of planets, such as electrical and thermal conductivity, mechanical modulus, magnetic field, sound velocity, viscosity, thermoelectric properties, etc. The electronic properties, with a chainmelting mediated switching between a pseudogapped electridelike solid and a frankly metallic liquid, would be revealing if observable. A deep understanding of it may also be helpful to understand the thermodynamic and kinetic properties of other metallic systems such as the glass transition, liquid metals [68], and certain astronomical phenomena such as calcium-rich supernovae [69].

VI. METHOD
A. First-principles calculations DFT calculations of Ca were performed using the Vienna ab initio simulation package (VASP) [45], together with the projector-augmented wave (PAW) method [70] retaining 10 electrons (3s 2 3p 6 4s 2 ) in the valence with small inner-core radius (1.22 Å), and the generalized gradient approximation exchange correlation functional with the Perdew-Burke-Ernzerhof [71] parametrization. A plane wave cutoff of 380 eV was employed. The cutoff energy was tested to an energy convergence. The optimization of lattice constant using these parameters provides a good agreement with the experimental equation of state (Fig. S9 [48]).
Ab initio molecular dynamics simulations were performed with a NVT ensemble using cubic supercells and periodic boundary conditions. For three room-temperature phases in the pressure range between 0 and 50 GPa with pressure step 10 GPa, the simulation cells consisted of 128, 108, and 125 atoms for fcc, bcc, and sc Ca solid phases, respectively. The equations of motion were integrated with ionic time steps of 0.5 fs, the Brillouin zone was sampled at the Γ point, and the ionic temperature was controlled with a Nosé-Hoover thermostat [72,73]. The same "chain melting" appeared if a 2 × 2 × 2 Monkhorst-Pack mesh k-point sampling was used on a 64 atom supercell at 50 GPa and 1500 K (Fig. S10 [48]). Simulations ran for at least 5 ps, and some of the trajectories were extended up to 10 ps to check the stability of the simulations. The initial 2 ps were used for thermalization and the last 6000 steps were used to extract the statistical quantities. These parameters ensure a stable canonical ensemble sampling and guarantee us to extract an accurate vibrational density of state ( Fig. S11 [48]).

B. Machine-learning neural network potential
We trained a machine-learning neural network potential for calcium using the DeePMD method [47], which directly uses relative atomic coordinates to describe local atomic environments to obtain atomic energies, while forces are obtained by taking the derivatives of the energies. It builds for every atom a local coordinate frame in order to keep invariance to translation and rotation. Our dataset was built by choosing 7000 configurations in each trajectory from AIMD of sc Ca in the canonical (NVT) ensemble at 40 and 50 GPa with a temperature step of 200 K from 900 to 2700 K including solid, cooperative diffusion state, and liquid, with 80% of it as a training set and 20% as a test set. Further details of DeePMD can be found in Ref. [47].

C. ML MD
MD with machine-learning neural network potential has been performed using LAMMPS [46] with periodic boundary conditions and a time step of 1 fs. The same range of temperature and pressure in the AIMD simulations of sc Ca has been used. To reproduce and verify the cooperative diffusion state in sc Ca, simulations were conducted in 15 × 15 × 15, 20 × 20 × 20, and 25 × 25 × 25 cubic supercells with 3375, 8000, and 15 625 atoms, respectively. The total time of ML MD was at least 100 ps, and some of the trajectories were extended up to 2 ns to investigate the stability of this state.

D. Mean-square displacement and velocity autocorrelation function
The averaged mean-square displacement is defined by where r i ðtÞ is the position of the atom i at the time t, and hi represent an average on the time steps and/or the particles. And the velocity autocorrelation function is defined by where v i ðtÞ ¼ _ r i and N is the number of atoms. The diffusion coefficient D can be calculated using one of the below two formulas:

E. Flux functions
Cutting a number of planes perpendicular to the h100i edge along the z direction in the middle of two atoms, z 0 is a set of coordinates representing their positions. x ði;tÞ , y ði;tÞ , z ði;tÞ are coordinates of the ith atom at time of t after trajectory smoothing of 1 ps. N is the total number of atoms. F z ðx; yÞ represents how many atoms pass through the given planes parallel to the xy plane. Hðz i;t − z 0 Þ − Hðz i;tþΔt − z 0 Þ only has a value when the i atom at time t and t þ Δt spans a given plane, where Δt equals 1 ps. The positive and negative value of F z ðx; yÞ depends on the direction of the atom's diffusion. R 1 and R 2 are parameters of the cutoff function f c ðrÞ in the XY plane. We take where d is the nearest neighbor atom spacing. Similarly, we can define F x ðy; zÞ and F y ðx; zÞ to investigate the chain melting in the other two directions: HðxÞ ¼ 0 F. Quantization of diffused chain atoms and thermodynamic properties The first-principles calculations of entropy and free energy were conducted in each phase utilizing the twophase thermodynamic model [49] and its enhanced version based on the memory function formalism [50]. It has been applied successfully to study the melting curve of aluminum [60] and the superionic water [61]. The key quantities in the 2PT MF method are the velocity autocorrelation function and the phonon DOS FðvÞ (normalized to 3) derived from the Fourier transform, FðvÞ ¼ 12m kT where ⃗ v i ðtÞ is the velocity of atom i at time t, which can be obtained from AIMD, and k is Boltzmann's constant, T is the temperature, and m is the mass of the atom. In this study, for chain melting and liquid, the total correlation function and the phonon DOS can be decomposed into solidlike and gaslike components. Since the entropy S in these systems includes an electronic contribution, S e , that is part of the free-energy functional minimized in AIMD, and a contribution due to the nuclei motion, and S i , we evaluated the latter contribution explicitly in this study. The details of this method can be found in Refs. [49,50].
In the 2PT MF method, f g represents the proportion of atoms with diffusing properties in the liquid, so we calculated values of f g in the chain-melting sc Ca at different temperatures and compared them with the value in pure liquid state to obtain the tendency of the proportion of diffused chains. The specific values of f g at different temperatures can be found in Table S1 [48].

G. Movies
The movies in Supplemental Material are produced using Visual Molecular Dynamics (VMD) [74].

H. Energy landscape, averaged electronic DOS and lattice dynamics, elastic constant
The nudged-elastic-band (NEB) calculations were performed in 3 × 3 × 3 supercell models using a Γ-centered 6 × 6 × 6 k-point grid with a tolerance on forces of 30 meV/Å. 5 × 5 × 5 supercells were used to calculate the electronic DOS for solid, cooperative diffusion state, and liquid with 8 × 8 × 8 Monkhorst-Pack mesh k-point sampling. The results for chain melting and liquid were averaged on five uncorrelated snapshots extracted from AIMD simulations.
Lattice dynamics were calculated in a 2 × 2 × 2 supercell with 13 × 13 × 13 Monkhorst-Pack mesh k-point sampling. The harmonic 0 K phonon sprectrum was calculated using the density-functional perturbation theory method as implemented in the PHONOPY package [75] and finitetemperature renormalized phonon spectra were calculated using phonon mode decomposition and AIMDs' projection technique implemented in the DYNAPHOPY package [76].
The total elastic moduli calculations for different structures used a grid density of 0.016 × 2π Å −1 for k-points sampling. We used a five-electron pseudopotential with 1.0 Å inner-core radius and 350 eV plane wave cutoff for sc P.

I. Bader charge analysis
Based on charge densities calculated by DFT, the Bader charge [58] analysis code is used to divide the Ca cations and localized electron "anions." ACKNOWLEDGMENTS J. S. gratefully acknowledges financial support from the National Key R&D Program of China (Grant