Measuring the structure and equation of state of polyethylene terephthalate at megabar pressures

We present structure and equation of state (EOS) measurements of biaxially orientated polyethylene terephthalate (PET, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\hbox {C}}_{10} {\hbox {H}}_8 {\hbox {O}}_4)_n$$\end{document}(C10H8O4)n, also called mylar) shock-compressed to (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$155 \pm 20$$\end{document}155±20) GPa and (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6000 \pm 1000$$\end{document}6000±1000) K using in situ X-ray diffraction, Doppler velocimetry, and optical pyrometry. Comparing to density functional theory molecular dynamics (DFT-MD) simulations, we find a highly correlated liquid at conditions differing from predictions by some equations of state tables, which underlines the influence of complex chemical interactions in this regime. EOS calculations from ab initio DFT-MD simulations and shock Hugoniot measurements of density, pressure and temperature confirm the discrepancy to these tables and present an experimentally benchmarked correction to the description of PET as an exemplary material to represent the mixture of light elements at planetary interior conditions.


X-ray diffraction measurement
The X-ray diffraction (XRD) experiment was performed at the Matter in Extreme Conditions (MEC) endstation 14,15 of the Linac Coherent Light Source (LCLS) at the SLAC National Accelerator Laboratory. A 50 µm thick PET foil was irradiated with a laser pulse containing 32 J in 6 ns (full width at half maximum), focused to a spot size of 150 µm to 200 µm . A very thin layer of the target is turned into plasma and expands rapidly, driving a shock wave into the remaining material due to the ablation pressure. To avoid transmission of the driving pulse through the target before an absorbing corona has been formed, the sample was coated with 100 nm aluminium.
We probed the compressed sample with the LCLS free electron laser X-ray beam, collecting the diffraction pattern with area detectors capable of recording single-photon events 16 . This setup (a schematic is shown in Fig. 1) is similar to previous studies on hydrocarbons 17 but differs from experiments that found diamond formation in polystyrene 10,11 by the choice of a single-shock loading scheme, because the Hugoniot for PET (see "Comparison to different equations-of-state" section) is expected to realize colder states at similar pressures when compared to PS.
For the analysis of the X-ray diffraction data, the signal on the detector pixels is integrated azimuthally using DIOPTAS 18 . Moreover, corrections to account for the experiment geometry are applied to the resulting 1D datasets, since the X-ray absorption inside the sample and from an attenuator in front of the detector varies with the angle of incidence. Figure 2 shows the curves for increasing delays between the drive laser and the X-ray pulse, i.e. further transit of the shock front. The crystalline structure of the cold material results in clear peaks which recede over time as the plastic is shock-compressed and heated. By t = 6 ns these patterns have completely vanished as the shock breaks out at the rear surface, such that there is no remaining cold material. For even longer delays, the position of the signal shifts to smaller k which is due to the adiabatic expansion of the material after the shock breakout.
The XRD signal recorded close to the shock breakout has been compared to DFT-MD simulations of PET under different conditions, as this is the point at which the conditions inside the target are most homogeneous. For doing so the expected X-ray scattering intensity I has been calculated using (1) The resulting diffraction image has been recorded using a single-photon-sensitive area detector. www.nature.com/scientificreports/ where k is the length of the scattering vector, x a is the atomic fraction of species a (with a being O, C or H in the case of PET), f a the atomic form factor of the whole ion or atom a and f an the contribution to f a caused by the nth electron. The S ab denote the partial structure factors. The first line of Eq. (1) describes the elastic scattering and is obtained by assuming that the scattering occurs entirely from bound electrons, rather than screening clouds, as in the approach of Wünsch et al. 19 , which is justified by the low expected temperatures where only minor ionization should occur 20,21 . The second summand takes the inelastic scattering into account as derived by James 22 , supplemented by the summation over different species of ions or atoms. The form factors of the individual electrons in different orbitals were calculated by treating the low-Z atoms as hydrogen-like as it was performed by Pauling and Sherman 23 while all other quantities were obtained from DFT-MD simulations (see "Methods" section).
Scattering intensities for densities from 2.5 to 4.0 g/cm 3 and temperatures between 4000 and 10,000 K have been evaluated (see Fig 3). As expected, a higher density shifted the position of the peak in the intensity to higher scattering vectors while an increasing temperature results in an increased peak width. To fit the artificial XRD signal to the measurement, squared deviations of the DFT-MD predictions ( C i ) from the experimental data ( M i ) have been calculated and weighted with the variance ( σ 2 i ) of the data at the given angle (with N being the number of scattering vector points) for a wave number range ≈ 2.0 Å −1 < k < 5.0 Å −1 . At higher scattering vector lengths, i.e. large angles of incidence, the response function of the detector is uncertain while the underestimation of the signal for lower k is likely caused by rapid expansion of small fractions of the sample due to the onset of the shock release. The indicator ν 2 favours the prediction for T = 6000 K at P = 153 GPa and ρ = 3.5 g/cm 3 ( ν 2 = 1.8 ) over those for identical density and T = 7000 K at P = 162 GPa ( ν 2 = 2.5 ) and T = 5000 K at P = 143 GPa ( ν 2 = 5.1 ) as it can be seen in the subplot of Fig. 3. By setting a threshold of acceptance we determined the uncertainty of the comparison. Hence, the specified error-intervals are a result of a parameter variation and therefore estimations. By applying the described procedure, the conditions in the target achieved with one single shock wave have been estimated as Figure 2. Integrated XRD data at ambient conditions and different time delays between the pump and the probe laser. The cold material shows peaks from the crystalline PET structure as well as from the aluminium coating. With the propagation of the shock front, these features recede and a signal from the amorphous liquid plastic remains, showing a wide peak centred around k = 3.2 Å −1 . After shock breakout around t ≈ 6 ns the peak shifts to lower k, as it can be observed in the curves marked with "(b.o.)". Those later delays are offset for clarity. www.nature.com/scientificreports/ Using Eq. (1), it is possible to separate the total scattering intensity into the oxygen-containing and oxygenfree components, the black and red curves, respectively, in Fig. 4. The two contributions differ substantially for k < 4.5 Å −1 indicating that the considerable liquid correlations around k = 3.2 Å −1 are caused mainly by oxygen atoms. Being negative for small k, the partial structure factor S CO lowers the intensity in this regime compared to the pure carbon and hydrogen signal.

Shock Hugoniot measurements
For comparing the conditions indicated by the matching of the XRD signal to the DFT-MD predictions to established Hugoniot measurement techniques, additional experiments were performed at the Laboratoire pour l'Utilisation des Lasers Intenses (LULI). Multi-component targets (layers of ≈ 10 µm CH as an ablator, 48 µm Al (pusher), 20 µm to 61 µm quartz ( SiO 2 , standard) and 49 µm PET (sample)) were dynamically shock compressed by an optical drive laser with energies from ≈ 236 J to 945 J. The shock velocity was measured using two velocity interferometers for any reflector (VISAR) 24 while a calibrated streaked optical pyrometer (SOP) simultaneously provided temperature information (see Fig. 5).
To determine density and pressure in the PET, the shock velocities in quartz and PET were extracted from the VISAR fringe-shifts of the moving, reflective shock or transit times (see "Methods" section). By using a known Hugoniot for α-quartz and an analytical release model 25,26 , pressure, density and internal energy of the PET can be calculated from impedance matching and the Rankine-Hugoniot equations.
The temperature of the target was determined using a SOP system to measure the emission from the shocked material while assuming the target to be a grey-body with wavelength independent reflectivity 27 . The Hugoniot curve acquired by these VISAR and SOP measurements is shown in Fig. 6.

Comparison to different equations of state
To validate the results obtained in the "X-ray diffraction measurement" and "Shock Hugoniot measurements" sections, three Hugoniot curves for PET in the WDM regime have been calculated using different EOS. Subsequently two tables, SESAME 07550 Mylar 12 and PrOpacEOS 4.0.0 C5H4O2 13 , will be compared to an EOS obtained from ab initio DFT-MD. Ambient conditions ( ρ 0 = 1.38 g/cm 3 at room temperature and normal pressure) were set as the starting points for the Hugoniots. The resulting graphs are illustrated in Fig. 6 and allow us to compare the EOS tables with the experimental observations.  Comparison of the X-ray diffraction line-outs at the shock breakout to predictions made using density functional theory and molecular dynamics calculations. The gray region marks the k interval in that the simulations are fitted to the measurement. In the small plot, squared-deviations ν 2 (see Eq. (2)) are plotted for different simulation parameters. The coloured lines (representing T = 7000 K at P = 162 GPa , T = 6000 K at P = 153 GPa and T = 5000 K at P = 143 GPa ) fit the experimental data best. www.nature.com/scientificreports/ The temperatures that SESAME predicts are higher than the results of calculations based on the PrOpacEOS equation of state and DFT-MD for similar pressures. SESAME and DFT-MD agree on the pressure-density relation except in an interval between 2.5 and 3.0 g/cm 3 , where pressure-induced dissociation of the molecules is expected, while PrOpacEOS predicts higher required pressures to reach the same density. The discrepancies between these graphs indicate the poor understanding of PET under WDM conditions. The blurred regions in Fig. 6 depict the conditions that gave best agreement with the XRD pattern in the "X-ray diffraction measurement" section. The ab initio Hugoniot is compatible with this area while differences in density or temperature are found for PrOpacEOS and SESAME, respectively.
When comparing Hugoniot curves from different tables to the VISAR and SOP measurements described in the previous section, good agreement with the ab initio EOS can be found. Both SESAME and PrOpacEOS diverge considerably from both the simulated and measured data across the range of conditions probed, and provide therefore no adequate description of PET, or elemental mixtures of similar composition, at intra-planetary conditions.
To check the influence of a potential pre-heating of the sample, hydrodynamic simulations including radiation transport have been performed with the software package HELIOS from PRISM computational sciences, Inc. 13 for SESAME and PrOpacEOS. We assumed the temperature of electrons and ions to be equal and estimated the thermal conductivity with the Spitzer-model 28 . The results of the hydro-simulations agree with the calculations from the Rankine-Hugoniot equations within the uncertainties and therefore suggest only a minor change caused by pre-heating while validating the starting conditions and interpolations used for the calculation of the Hugoniot from the EOS tables.

Conclusions and outlook
We presented new data of polyethylene terephthalate from DFT-MD simulations and laser-shock experiments and compared it to tabulated EOS values to infer the macroscopic thermodynamic properties of PET at intraplanetary conditions. Furthermore, in-situ X-ray diffraction measurements were performed and analyzed by fitting DFT-MD based scattering image predictions to the experiment, allowing us to gain insight in both the microscopic structure and the temperature, pressure and density conditions of the sample, simultaneously.
The different methods of XRD, classical shock-state measurements and DFT-MD EOS calculations determine the state of our target consistently, complementing each other. Given this accordance, it might be worth exploring if XRD, in combination with DFT, could be utilized as a diagnostic not only for density but also temperature, e.g. for states off the Hugoniot curve.
We stress that the time scales relevant in DFT-MD simulations (ps), laser shock-experiments (ns) and planetary live-times (millions of years) differ vastly. Effects like superheating-that can occur during rapid shock compression 27,29 but are not likely to be relevant for astrophysical objects-impede a direct analogy. Nevertheless, nano-second scale experiments allow experimental access to the WDM regime that can be used to gain valuable insight into the physics under extreme conditions.
Using the described method, we constrain the conditions inside the region probed by the X-rays to T = (6000 ± 1000) K , P = (155 ± 20) GPa and ρ = (3.5 ± 0.4)g/cm 3 . This result disagrees with the Hugoniots of both SESAME 07550 Mylar 12 (which predicts higher temperatures) and PrOpacEOS 4.0.0 C5H4O2 13 (that www.nature.com/scientificreports/ underestimates the density). However, the conditions do fall on the Hugoniot of the EOS derived from our own DFT-MD simulations. This dataset has been compared to SOP and VISAR measurements of shock compressed PET where we found agreement within the error bars. The regime characterized by the T, P and ρ stated above is believed to be realized inside the icy planets 2,30 and is close to the conditions where diamond formation-which might be supported by the presence of oxygen 31 -has been observed in comparable experiments with polystyrene samples 10,11 . While no conclusive statement about the crystallization of diamond from PET can be given here based on the analyzed dataset due to the notable liquid background, further experiments for obtaining insight into the possibility of this process, e.g. with thicker samples (and therefore more time for diamonds to form) or colder, weaker shocks (since the temperature we found lies only little under the diamond melting line 32 ), seem promising-even though correlations of oxygen with carbon and hydrogen may impede the detection. PET seems to be a more suitable model system for the icy giants' interior layers than PS, since it contains oxygen, but also due to the technical simplicity of generating relevant conditions using just a single shock wave.  www.nature.com/scientificreports/ Methods XRD. We probed the sample using a ≈ 50 fs long X-ray beam with ≈ 3 mJ pulse energy focussed to a spot size of 20 µm (full width at half maximum). The photon energy amounted E ph = 8.1 keV ± 0.3 % . The diffraction measurement was performed in an angular range from 2θ = 20 • to 95 • (corresponding to scattering vector lengths from k = 1.4 Å −1 to 6.0 Å −1 ) with a 8 cm × 8 cm Cornell-Stanford Pixel Array detector 16 . Gaps in the detector and malfunctioning pixels have been masked out before the signal was integrated azimuthally. In the latter step, corrections for the polarization of the X-ray beam were applied.
The electronic density in the simulation box with periodic boundary conditions was represented by a plane wave expansion with a cut-off energy of E cut = 1000 eV . We used the Mermin formulation of DFT to optimize the Helmholtz free energy at a given temperature 37 . The electron-ion interaction was modeled using the projector augmented wave (PAW) approach, specifically the hard PAW pseudopotentials for hydrogen (H_h, 06Feb2004), carbon (four valence electrons, C_h Feb2004), and oxygen (six valence electrons, O_h Feb2004) as provided with VASP 38,39 . The exchange-correlation potential was taken in generalized gradient approximation in Perdev-Burke-Ernzerhof parametrisation (GGA-PBE) 40 . We generally sampled the Brillouin zone of the supercell at the Ŵ-point only. The electronic bands where populated using a Fermi distribution at the chosen temperature. We performed simulations with two different system sizes: 32 (64) protons, 40 (80) atoms of carbon, and 16 (32) atoms of oxygen (this corresponds to four (eight) PET subunits), and their movements were calculated using the Hellman-Feynman forces derived from the electron densities of DFT under the Born-Oppenheimer approximation. The DFT-MD run covered a time span of 2 ps to 4 ps, the time step was t = 0.2 fs , and the ion temperature was controlled by a Nosé-Hoover thermostat 41 .

VISAR.
To determine the shock velocity without ambiguity, two VISAR systems with different sensitivities were included in the setup. For each, the velocity per VISAR fringe VPS 0 was calculated using 42,43 where L is the wavelength of the applied laser (1064 nm / 532 nm), c the speed of light, n = 1.4497 or 1.4607 the wavelength dependent index of refraction of the etalon material, d its thickness (8.06 mm/10.04 mm) and δ = 0.0163 or 0.0318 a correction for the dispersion in the etalon 43,44 . We found for the two systems VPF 0 = 12.81 km/s/fringe and VPF 0 = 4.98 km/s/fringe , respectively. To calculate the velocity per fringe in a medium, VPF 0 has to be divided by its index of refraction which is given as n Qz = 1.54 and n PET = 1.63 for our targets. For the two lowest-energy shots, where this procedure was not applicable in the quartz-region due to the low reflectivity, transit times were used to calculate the shock velocities instead. To measure this quantity, the time of shock-entry and -breakout was determined up to ±0.15 ns , each, and the uncertainty in the target's thickness amounted ±1 µm.

SOP.
For analyzing the SOP data, we treated the target as a grey body with wavelength-independent reflectivity R, measured using the VISAR system and normed to the known reflectivity in quartz 45 . With these assumptions, the spectral radiance L of the target is given by In this equation, c denotes the speed of light, h the Planck's constant, the radiation's wavelength and k b Boltzmann's constant 46,47 . Using a narrow bandwidth filter that can be described by a δ function around a given wavelength 0 , Eq. (4) can be inverted and the temperature inside the sample calculated using where I is the number of counts accumulated on the CCD chip. T 0 = hc 0 k b only depends on the filter used 45 while A is a constant containing various machine-parameters 46 and was obtained using the known temperature to shock velocity relation of α-quartz 27,48 as a reference 49 . Errors where estimated from the uncertainty of this calibration (by fitting upper and lower bounds to the reference data, resulting in relative derivations of T 0 : +2.5% −15% and A : −10% −14% or A : +39% −22% at a fixed T 0 , depending on the filter 49 ), the reflectivity measurement and the signal's variance.

Data availability
The datasets generated during the experiments are available from the corresponding author on reasonable request.