Negative linear compressibility in Se at ultra-high pressure above 120 GPa

Structural evolutions of selenium on strong compression are still not fully understood, therefore in situ X-ray diffraction measurements up to 210 GPa were carried out. An anomalous negative linear expansion was observed at pressures from 120 to 148 GPa, accompanied by a new isostructural phase transition.

To date, various pressure-induced phase-transition sequences have been reported on trigonal Se. Parthasarathy & Holzapfel (1988) reported the Se of the trigonal structure directly transformed to the monoclinic phase at 14 GPa, then to the tetragonal phase at 28 GPa. Akahama et al. (1993) reported that compressed trigonal Se experienced an intermediate phase from 14 to 23 GPa, and then transferred to a base-centered monoclinic phase at 23 GPa, i.e. trigonal ! intermediate phase ! base-centered monoclinic sequences. In the higher-pressure range, the phase transition from an incommensurately modulated monoclinic structure to a -Po type structure at 40-60 GPa was observed, whereas McMahon et al. (2004) suggested this transition occurs at 82 GPa by extrapolating the monoclinic structural parameters, and Degtyareva et al. (2005a) experimentally observed this phase transition at 80 GPa. So far, about seven selenium allotropes related to high-pressure and high-temperature conditions have been reported experimentally (Degtyareva et al., 2005b(Degtyareva et al., , 2007Akahama et al., 2021). Therefore, the pressure-induced phase transition sequence of trigonal Se should be as follows: trigonal phase (Se-I, space group P3 1 21) ! C-face-centered monoclinic phase (Se-II, space group P2 1 /c, 14 GPa) ! metallic base-centered monoclinic phase (Se-III, space group C2/m, 23 GPa) ! incommensurate monoclinic phase (Se-IV 28 GPa) ! rhombohedral phase with -Po structure (Se-V, $80 GPa) ! b.c.c. phase (Se-VI, 140 GPa). Phase Se-VI was reported to be stable up to 317 GPa, the highest recorded experimental pressure value for Se to the best of our knowledge (Akahama et al., 2021).
Theoretical studies based on density-functional calculations showed a structural phase transition for Se from the b.c.c. to the f.c.c. structure occurs around 260 GPa (Rudin et al., 2001). Other calculations indicated that Te transfers from a b.c.c. structure to an intermediate structure at 100 GPa, then to an f.c.c. structure at 255 GPa (Sugimoto et al., 2014). Given the previous calculations and the similarity between Te and Se in terms of pressure-induced phase transition sequences, it is important to explore the possible phase transitions of Se experimentally. In the current work, pressure-induced phasetransition sequences in selenium up to 210 GPa are explored.
Among the inconsistencies with pressure-induced phasetransition sequences, the specific phase transition from the Se-V phase to the Se-VI phase has drawn much theoretical and experimental attention. Though Akahama et al. (2021Akahama et al. ( , 1993Akahama et al. ( , 1992b reported the Se-V to Se-VI transitional pressure is around 140 GPa in experiment, many calculations indicate there might be a phase transition around 120, 110 GPa or even as low as 95 GPa. Rudin et al. (2001) predicted the phase transition occurrs at 120 GPa based on density-functional theory (DFT). This is consistent with calculations carried out by Otani & Suzuki (2001), which revealed the transition pressure at 110 GPa based on the first-principles calculations performed on electronic structures of Se-V and Se-VI phases. Similar results of phase transitional pressure at 110 GPa were reported by Nishikawa et al. (1993), who calculated the related band-structure and total energy. However, using the fullpotential linearized-augmented-plane-wave method, Geshi et al. (1998) revealed the transition was at a lower pressure of 95 GPa. Obviously, the controversies are not only between calculation and experiment, but also between various calculations themselves. In the current experiment, a phase transition around 120 GPa is observed and its interesting behavior at pressure above 120 GPa will be discussed.
Note that Akahama et al. (1992b) observed some discontinuity of the pressure dependence of bond length at about 120 GPa, which might propose a possible second-order phase transition. However, in their recent work, this issue was left unaddressed (Akahama et al., 2021). In the current work, a clear discontinuity of volume change with increasing pressure was observed via in situ X-ray diffraction (XRD) around 120 GPa, accompanied with a linear expansion in lattice parameter a. Compared with earlier works, where negative linear expansion (NLE) in materials was only observed under comparatively lower pressures, this work presents NLE in Se at pressures over 100 GPa, which is of great interest and significance in terms of new materials discovery.

Experimental procedure and calculation details
The trigonal Se sample (99.99% purity) was loaded in the hole of the rhenium gasket and compressed between the membrane-driven diamond anvil cell (DAC). Powder XRD data were collected at GSECARS and HPCAT, Advanced Photon Source, Argonne National Laboratory. Focused monochromatic beams of 0.4066 Å wavelength were used and the diffraction patterns were collected with a PILATUS 1M detector. Two-dimensional XRD patterns were integrated to create one-dimensional patterns using the Dioptas program (Prescher & Prakapenka, 2015). Gold and silicone oil were used as pressure marker and pressure transmitting media, respectively. The beveled diamond anvils with 70-300 mm culet diameter were used. Rhenium gaskets were pre-indented to a thickness of 20 mm and then drilled to make a sample chamber hole 30 mm in diameter. The sample of trigonal Se was ground, then loaded into the gasket hole along with a few several-micrometre-sized gold foils, and the silicone oil was added before sealing the chamber. The size of the X-ray spot was about 10 mm in diameter. The diffracted X-ray signal was detected with an exposure time of 0.5 s, while constant compression was applied by the membrane-controlled DAC system, with a compression rate at about 0.6 GPa s À1 in the chamber. Diffracted images were collected throughout the experiment, until the pressure in the sample chamber reached 210 GPa. Like the majority of high-pressure DAC experiments over 100 GPa, the DACs were broken in this run after the research papers 254 Shuhua Yuan et al. Negative linear compressibility in Se at ultra-high pressure pressure reached 210 GPa. Therefore, no diffraction data were available during the decompression process and the questions of possible hysteresis and phase metastability ranges on pressure release remain unanswered.
XRD from both the selenium and gold were measured simultaneously by selecting a spot on the sample where both were present. This approach minimized the possible systematic errors due to pressure gradients and ensured accurate in situ pressure measurements. Pressures were calculated using the equation of state (EoS) for Au under pressure (Anderson et al., 1989).
The first-principles calculation was performed using the plane wave DFT program VASP [Vienna ab initio simulation package (Kresse & Furthmü ller, 1996)], where the electronion interactions of Se atoms were represented by the projector augmented wave (PAW) pseudopotential (Blö chl et al., 1994). The exchange-correlation functional was described by the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) parameterization (Perdew et al., 1996). The calculations were computed with a 10 Â 10 Â 12 Monkhorst-Pack scheme (Monkhorst & Pack, 1976). The energy cut off was set at 550 eV. The geometry convergence criterion was set at 0.001 eV Å À1 for the maximal component of force and 0.01 GPa for stress.

Results and discussion
The in situ XRD patterns of Se were collected under constant compression and short-time exposure, and these high-resolution synchrotron diffraction patterns were able to show minute changes in the structure. More than 4300 XRD 2D images were collected while the sample was compressed from ambient pressure to 210 GPa. Fig. 1(a) representatively shows one of the typical XRD images, under 90 GPa. The 2D images were integrated into 1D XRD patterns using DIOPTAS (Prescher & Prakapenka, 2015), shown in Fig. 1(b). Note that the patterns from 90 to 140 GPa match the hexagonal Se-V (-  GPa at ambient temperature, evolving from Se-V to Se-V 0 , then to Se-VI. A zoomed-in view of (b) with more XRD patterns added between the plotted ones, from 90 to 160 GPa, of 2 (c) 11 to 13 ; (d) 13 to 21 .
Po-type rhombohedral) structure. With increasing pressure to 140 GPa, a new peak near 11.8 2 appears, indicating a phase transition to Se-VI (b.c.c. structure). This is consistent with previous works where the phase transition from Se-V to Se-VI was found at 140 GPa (Akahama et al., 2021(Akahama et al., , 1993Degtyareva et al., 2005a). On further compression, Se-VI remains stable with no other phase transitions observed from 140 to 210 GPa.
Interestingly, a closer look into the XRD patterns above 120 GPa reveals that the diffraction peak (110) of Se-V starts to shift towards a lower value of 2 with increasing pressure, while other diffraction peaks are continually moving towards the larger value, as one would expect for elevated pressures. This trend is more clearly shown in Fig. 1(c). The shift is observed during compression from 120 to 148 GPa. According to the Bragg equation, this indicates the (110) peak moves towards a larger d-spacing value, which implies that the interplanar distance of (110) planes increases with pressure. Given the hexagonal interplanar spacing relationship with lattice parameters at a given plane it can be concluded that lattice parameter a in hexagonal Se increases with pressure in this range. Fig. 1(d) shows interplanar spacing changes from other planes: peaks (201) and (102). From 90 to 120 GPa, both (201) and (102) peaks move towards larger 2 values with higher pressure at relatively high rates. However, once the pressure is increased to above 120 GPa, (201) moves (towards larger 2 values) at a decreased rate, while (102) moves at an increased rate. Considering the interplanar spacing relationship given above, one can deduce that once the pressure increases to above 120 GPa, the lattice along the c direction is more compressible whereas along a is the opposite case.
Noticeably, Akahama et al. (1993) also reported a relevant observation with Se at the same pressure region. They observed that when the Se-V phase was compressed to pressures above 120 GPa, the first-nearest neighbor bond length barely decreased with increasing pressure, while the secondnearest neighbor bond length decreased more drastically. However, there were not enough finely spaced data to obtain an accurate trend of the lattice parameter a change at this dedicated pressure range in their work. In the current work, with higher-resolution data as well as the fine control of exposure time and compression rate in this experiment, we were able to observe the abnormal increase of the distance between (110) planes, which implied that the first-nearest neighbor bond of a given Se atom increased with pressure during 120-148 GPa, rather than just 'an unusual hardening' as pointed by Akahama et al. (1993).
In order to take a closer look at the change in lattice parameters, Rietveld fitting was employed to extract lattice parameters using GSAS-II (Toby & Dreele, 2013). Some of the patterns with Rietveld refinement at different pressures are as shown in Fig. 2.
Lattice parameters a and c, the c/a ratio, and atomic volume versus pressure are plotted, as shown in Fig. 3. It is evident that below 120 GPa, a, c and c/a decrease with increasing pressure. However, an increase of the lattice parameter a occurs once the pressure exceeds 120 GPa, which corresponds with the shift of diffraction peak (110) to lower angles starting from 120 GPa. This increase lasts until the pressure reaches 148 GPa, during which a is calculated to increase by 2.22%. Such a linear increase in lattice parameters and the vast difference in bulk modulus indicate an intermediate phase existing between 120 and 148 GPa, referred to as Se-V 0 here.
Accompanied with the increase of lattice parameter a, lattice parameter c decreases at a higher rate, which guarantees the continuous volume decrease. It is obvious from the pressure dependence of atomic volume that the increase of a slows down the volume decrease to some extent. It is worth noting that the increase of a corresponds with the slower rate of volume decrease. The discontinuity of pressure dependence of atomic volume at 120 GPa and a different slope in the EoS fitting curve from 120 to 148 GPa indicate a phase transition happening with a lower compressibility of Se under the pressure from 120 to 148 GPa. This suggests an isostructural phase transition around 120 GPa, followed by another phase transformation to a b.c.c. structure (Se-VI) at 140 GPa. In Fig.  3, the shaded area in the graph represents the region where the Se-V 0 and Se-VI phases coexist.
In order to have a comprehensive understanding of such phenomenon, DFT calculations were carried out to study Se under high pressure. From the DFT calculations, a shows an increase above 120 GPa, which is consistent with the experiment, and then, a decrease with compression as shown in Fig.  3. Regardless of the absolute value difference between calculations and experimental data, the trend of c/a ratio is in fair agreement with the experimental data, which supports the fact that Se-V has experienced an isostructural phase transition at around 120 GPa. We further plotted the energy of Se-V 0 and Se-VI (relative to the Se-V structure) as a function of pressure in Fig. 4. The structure of Se-V 0 remained stable up to 148 GPa, while Se-VI was calculated to be more stable at pressures higher than 148 GPa. However, in the experiment, Pressure dependence of lattice parameters a and c, and atomic volume V from 90 to 210 GPa. The open symbols represent experimental data in this work, the solid symbols represent the calculated lattice constants a and c, and the c/a ratio, respectively. The b.c.c. structure remains stable up to 210 GPa, no new phase transition was observed from 140 to 210 GPa. The phases Se-V 0 and Se-VI coexist from 140 to 160 GPa.

Figure 4
Energy of Se-V 0 and Se-VI (relative to the Se-V structure) as a function of pressure from 110 to 160 GPa. discrepancy might be due to the limitation of the pseudopotential used in the VASP code, which is beyond the scope of our current work. Noticeably, the energy of the Se-VI phase relative to the Se-V structure is drastically reduced at 120 GPa, which indicates a huge adjustment of the lattice to accelerate the transition to the b.c.c. structure. Overall, the DFT calculation values are in adequate agreement with our experimental results, which have revealed the Se-V to Se-V 0 isostructural phase transition takes place around 120 GPa whereas the Se-V 0 to Se-VI structural phase transition takes place around 140 GPa.
The volume-pressure relationship was fitted according to the third-order Birch-Murnaghan EoS (Singh, 2005;Jeanloz, 1988;Birch, 1947), as shown on Fig. 3. The fitting results provide a value for the bulk modulus B 0 as 83 AE 2 GPa for Se-V, 321 AE 2 GPa for Se-V 0 and 266 AE 7 GPa for Se-VI. All the other fitting parameters are shown in Table 1. Recently, Akahama et al. (2021) reported the bulk modulus of Se-VI as 290 GPa using the diffraction data up to 317 GPa. Generally, more compression data covering a wider pressure region could yield a more precise bulk modulus. More P-V data, especially at a higher-pressure range between 210 and 317 GPa would achieve a better estimation of the bulk modulus. For easy comparison, the derivative of the bulk modulus B 0 was fixed at 3.2 during the fitting in our case, the same as the value of the Se-VI phase obtained in the recent report by Akahama et al. (2021). However, due to the flattening in the P-V curve near 300 GPa, their reported bulk modulus of Se-VI is slightly higher than the value from the current report. The other difference is the bulk modulus of Se-V, which they reported as 205 GPa, whereas in the current report, a smaller value of 83 AE 2 GPa for Se-V, and a larger value of 266 AE 7 GPa for Se-V 0 are reported. This could be due to the combined effect since P-V data from phase Se-V and Se-V 0 were not fitted separately in their report.
Interestingly, Se-V 0 has the biggest bulk modulus and the smallest compressibility among the three phases, which indicates Se-V 0 is the least compressible. As mentioned above, in Fig. 3 the increase of the lattice parameter a is accompanied by a decrease of volume at a reduced rate, which leads to a least compressible intermediate phase from 120 to 148 GPa. The increase of lattice constant a contributes to relaxation in the distortion from hexagonal Se-V, -Po-type rhombohedral structure, to Se-V 0 , which accelerates the structural evolution from Se-V to the b.c.c. structure Se-VI. Geshi et al. (1998) reported a similar phenomenon that the compressibility increased in the b.c.c. structure when compared with that in the previous phase, but no explanation was given. They also conducted calculations which suggested a phase transition around 95 GPa but wrongly assigned this phase transition as that observed experimentally at 140 GPa (Akahama et al., 1993).
Similar to the increase of lattice parameter a observed in Se-V 0 at high pressure, similar linear expansion has been observed in other materials under pressure, especially with structures like soft porous frameworks (Sobczak et al., 2020;Cai et al., 2015;Cai & Katrusiak, 2014), ferroelastics, tilting networks and helices (Cairns & Goodwin, 2015). This counterintuitive high-pressure behavior is usually referred to as negative linear compressibility (NLC). It arises whenever volume reduction can be coupled to linear expansion.
Notably, NLC has only been observed in materials under comparatively lower pressures. For example, NLC was found at pressures lower than 14 GPa in selenium with a helical structure, referred to as Se-I (Cairns & Goodwin, 2015;McCann et al., 1972;Keller et al., 1977;Akbarzadeh et al., 1993). The rare occurrence of NLC at much higher pressures (i.e. >100 GPa) is uncommon. As shown in this work, the NLC is observed in hexagonal selenium with a -Po-type rhombohedral structure between 120 and 148 GPa.
The marked four atoms in the -Po structural unit cell from Fig. 5(a) form a tetrahedron with lattice parameter a as sides of the base face; d 1 and d 2 are bond lengths from first-nearest and second-nearest neighbors, respectively; and is the angle formed by atoms 1, 2 and 3, as depicted in Fig. 5(a). From 120 to 148 GPa, with increasing pressure both d 1 and d 2 (corresponding to lattice parameter c) decrease, while increases, which results in a shorter and flatter tetrahedron. Such a geometrical change of structure leads to an increase of lattice parameter a under pressure. Note, in this process, d 1 decreases much slower compared with d 2 , and its decreasing rate reduces rapidly with compression. This indeed accelerated the transformation from a -Po-type rhombohedral structure to a cubic structure. Schematic of geometrical interpretation of negative linear compressibility in the structures of (a) Se-V 0 and (b) Se-I. Table 1 The EoS parameters of selected Se high-pressure phases.
Reference volume V 0 , bulk modulus B 0 and the derivative of the bulk modulus with respect to pressure B 0 0 , obtained from the least-squares fitting of the experimental data according to the third-order Birch-Murnaghan EoS. Values in brackets are estimated standard uncertainties and refer to the least significant digit; they reflect only the uncertainties derived from fitting V versus P points to the EoS. (2) 3.77 (6) Se-V 0 17.7 (1) 321 (2) 3.2 (Fixed) Se-VI 15.76 (2) 266 (7) 3.2 (Fixed) In terms of the NLC phenomenon with Se-I, which has a trigonal structure consisting of helices, c was observed to increase with elevating pressure, while a decreased. In this structure, the bonds between helices are much weaker than those within a specific helix. As a result of this anisotropy, the bond length in the a direction decreases while it increases in the c direction. The marked atoms in Fig. 5(b) form a parallelogram with two sides labeled l 1 and l 2 , between which the smaller angle is marked as '. Therefore, an increase of l 1 indicates an increase of c and vice versa, and l 2 corresponds to a. With increasing pressure, l 2 (parameter a) decreases rapidly, while both ' and l 1 (parameter c) increase. As a result, the pressure-induced densification of Se-I demonstrates preferential compression of a at the expense of some expansion in c. Table 2 shows previous works on NLC observations in Se, including Se-I, and the current work on the Se-V 0 structure. In the high-pressure structure, the absolute value of linear compressibility in Se-V 0 is about an order of magnitude smaller than that of Se-I. Cairns & Goodwin (2015) summarized some other materials with tilting network structures that exhibited NLC behavior under pressure, usually due to some correlated polyhedral tilts. This is similar to the mechanism of NLC in Se-V 0 discovered in this work.
Compared with the NLC occurring in Se-I, the NLC observation in Se-V 0 is more remarkable because the latter happened at pressures over 100 GPa. Considering the analogy of behavior and structure dependence among the elements in the VIA group under pressure, there is great promise in exploring possible NLC occurrences in S and Te under ultrahigh-pressure conditions. Interestingly, Hsueh et al. (2000) reported the NLC phenomenon in Te with c at pressures <5 GPa, which can be related to that of selenium in Se-I within 14 GPa. With this discovery of NLC in Se-V 0 at >100 GPa, there are more potential NLC phenomena in elements S and Te to explore under strong compression. This NLC behavior might be a universal feature in many elements under ultrahigh-pressure conditions.
Note that the pressure media used in this work, silicone oil, does not guarantee the hydrostatic pressure environment. Besides, the potential effects from varying compression rate would fine tune sample strain and stress state (Konô pková et al., 2015). A better pressure media like He or Ne is suggested for future works with various fast compression rates to substantially minimize the role of uniaxial stress and maximize the potential strain rate effect on this dedicated research topic.

Conclusions
In this work, selenium was compressed up to 210 GPa. An isostructural phase transition from Se-V to Se-V 0 was observed for the first time around 120 GPa, and the newly discovered Se-V 0 phase was estimated to have the largest bulk modulus of 321 AE 2 GPa, compared with that of Se-V (83 AE 2 GPa) and Se-VI (266 AE 7 GPa). The largest bulk modulus of Se-V 0 indicates its lowest compressibility among the three phases, which can be attributed to the occurrence of NLC during compression. Considering that NLC at higher pressure is always more significant in terms of fundamental mechanism and new materials discovery, and that it has barely been reported at such high pressures above 120 GPa, this discovery inspires future studies into NLC behavior in new materials at higher pressures. NLC behavior observed in selenium in this work is mainly attributed to the accuracy and fine steps of in situ XRD data collected at 0.5 s intervals. Such high-resolution measurements make it possible to map the evolution of sample structure at a continuously increasing pressure. The discoveries of the new phase Se-V 0 and its NLC behavior in selenium, a material considered one of the most studied elements in high-pressure research, illustrate a promising future where better understandings can be achieved by taking advantage of rapidly improving experimental techniques. Table 2 Linear compressibilities of selenium.

Funding information
Compressibilities of Se-I summarized from previous reported works and Se-V 0 from this work.