Thermophysical properties for shock compressed polystyrene Cong Wang,1 Xian-Tu He,1,2 and Ping Zhang1,2,3,∗ 1LCP, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China 2Center for Applied Physics and Technology, Peking University, Beijing 100871, People’s Republic of China 3National Institute for Fusion Science, Toki-shi 509-5292, Japan Wehaveperformed quantummolecular dynamicsimulations for warm densepolystyreneat high pressures. The principal Hugoniot up to 790 GPa is derived from wide range equation of states, 1 wherecontributionsfrom atomic ionizationsaresemiclassically determined. Theoptical conductiv- 1 ity is calculated via the Kubo-Greenwood formula, from which the dc electrical conductivity and 0 opticalreflectivityaredetermined. Thenonmetal-to-metaltransitionisidentifiedbygradualdecom- 2 position of the polymer. Our results show good agreement with recent high precision laser-driven n experiments. a J PACSnumbers: 82.35.Lr,51.30.+i,52.27.Jt, 82.40.Fp,71.15.Pd 5 2 In high energy density (E/V≥1011J/m3) physics polystyrene, especially for high pressures at, or exceed- ] (HEDP), pressure-induced response of materials, which ing, 1 Mbar. Recently, the development and applica- i c can be probed through shock wave experiments [1–4], is tion of quantum molecular dynamic (QMD) simulation s of crucial interest for inertial confinement fusion (ICF). [16] techniques for WDM are available due to the enor- - l Fundamental experiments in typical ICF designs, such mous progress in computer capacity. QMD simulations, r as the investigation of hydrodynamic instabilities, have where quantum effects and correlations are systemati- t m introduced plastic compounds as basic candidates to callytreated,haveconductedhighlypredictiveresultsto . achieve high gain [5]. Furthermore, in recent fast ignitor describe WDM. t a experiments, the electrical properties of polymers have In the present work, QMD simulations are applied to m beenoftenimplicatedforthetransportofrelativisticelec- calculate a broadspectrum of thermophysicalproperties - tron beams in solid targets as compared to conducting for shock compressed polystyrene. The self-consistent d materials[6]. Understandingthebehaviorofpolymersat electronicstructurecalculationwithindensity functional n several megabar regime brings better insight into physi- theory (DFT) yields the charge density distribution in o cal processes in ICF. the simulation supercell at every time step. The ion-ion c [ Some of the ablators in the indirect-drive mode of pair correlation function (PCF), which is important for National Ignition Facility (NIF), will be made of glow- characterizing and identifying phase transitions, can be 1 v discharge polymer (GDP) with various levels of germa- given by molecular dynamics run. The Hugoniot curve, 3 niumdoping(Ge-GDP)[7]. Due tothe factthatnohigh whichisderivedfromEOSdataforawideregionofdensi- 9 pressuredataexistforGe-GDP,polystyrene(CH),which ties andtemperatures,is determinedandcomparedwith 7 is closest in structure, is considered as a coarse indica- availableexperimentalandtheoreticalresults. Asastart- 4 tor for shock timing simulations of NIF targets involv- ing point, we use the Kubo-Greenwood formula to eval- 1. ing such ablators [8]. The absolute equation of states uate the dynamic conductivity σ(ω), from which the dc 0 (EOS)ofpolystyrenealongthe principalHugoniotcurve conductivity,thedielectricfunctionǫ(ω),andtheoptical 1 has been studied by using gas gun up to 50 GPa [9– reflectivity can be settled. 1 12] and high energy laser-drivenshock waves up to 4000 We introduce ab initio plane wave code VASP [17, 18] : v GPa [13–15]. High energy dynamic compressions pro- to perform QMD simulations. The elements of our i duce the so-calledwarmdense matter(WDM), wheresi- calculations consist of a series of volume-fixed super- X multaneous dissociations, ionizations, and degenerations cells including N atoms, which are repeated periodically r a make the transition between condensed matter physics throughout the space. By involving Born-Oppenheimer to plasma physics. The EOSand the relative properties, approximation, electrons are quantum mechanically such as the Hugoniot curve, are important features in treated through plane-wave, finite-temperature DFT this context. Moreover, electrical properties such as the [19], where the electronic states are populated accord- dielectricfunctionarecloselyrelatedtothedynamiccon- ing to Fermi-Dirac distributions. Sufficient occupational ductivity, which determines the dc electrical conductiv- bands are included in the overall calculations (the occu- ity in the static limit. Then optical reflectance can also pational number down to 10−6 for electronic states are be extracted. All these quantities are indispensable in considered). The exchange-correlation functional is de- characterizing the unique behavior of shock compressed termined by generalized gradient approximation (GGA) with the parametrization of Perdew-Wang 91 [20]. The ion-electron interactions are represented by a projector augmented wave (PAW) pseudopotential [21]. Isokinetic ∗Correspondingauthor: [email protected] ensemble (NVT) is adopted in the present simulations, 2 10 50 TABLE I: Polystyrene Hugoniot results from QMD simula- 103 8 u (km/s)s234000 tions.ρ (g/cm3) P (GPa) T (eV) us (km/s) up (km/s) 10 e (GPa) ure (eV)6 00 10up (km/s)20 30 22..0140 2224..0661 00..2201 66..5895 33..1492 ur at Press102 PBGraaerssr ieognsu tne wt oarl.k Temper4 22..4539 5734..4584 00..3329 190..4983 56..3479 Cauble et al. 2.79 128.56 0.75 14.01 8.74 Ozaki et al. 2 Ozaki et al. QEOS 3.01 240.46 1.54 18.76 12.21 SESAME 3.08 293.98 1.85 20.60 13.59 0 2.0 2.5 3.0 3.5 4.0 4.5 0 200 400 600 800 1000 Density (g/cm3) Pressure (GPa) 3.18 352.98 2.14 22.40 15.01 3.29 429.62 2.53 24.51 16.69 FIG. 1: (Color online) The P-V (left panel) and T-P (right 3.39 554.96 3.11 27.67 19.10 panel) Hugoniot curves of shocked polystyrene. Inset is the (u ,u )diagram. Previousresultsarealsoshownforcompar- 3.49 789.07 4.24 32.79 22.92 s p ison. Experiments: gas-gun results (the gray filled squares); [9–12]absolutemeasurementsonNOVA(magentadiamonds); [13] high energy laser-driven measurements by Ozaki et al. siderable effect of atomic ionizations on EOS should be [14,15](blueupwardanddownwardtriangles)andBarrioset appropriately included beyond QMD; (ii) QMD energy al. (orange squares). [8] Theories: QEOS and SESAME [26] consists of kinetic and potential energy (for both ions results are denoted by green dashed line and purple dotted and electrons), while, the well known zero-point vibra- line, respectively. tionenergy(ZPVE),vanderWaalsenergy(VDWE),and so as the demanded energy for phase transitions, which are particularly significant in WDM, are excluded; (iii) wheretheionictemperatureT iscontrolledbyNo´sether- The pressure from QMD only describes the interaction i mostat [22], and the system is kept in local equilibrium contributions, and the ideal part for noninteracting par- bysetting the electron(T )andion(T )temperaturesto ticles are missed. To overcome these obstacles, in our e i be equal. presentcalculations for polystyrene,ZPVE,which is im- portantfor determining the low-pressureinternalenergy Thesimulationshavebeendonefor128atoms(namely, for the reference state along the Hugoniot curve, is sim- eightC H units)inasupercell. Aplane-wavecutoffen- 8 8 ply added onto the QMD output. Condensed matter to ergyof600.0eVisselected,sothataconvergenceofbet- plasma transition energy mainly comes from atomic ion- terthan3%issecured. The Brillouinzoneis sampledby Γpointand3×3×8Monkhorst-Packscheme[23]kpoints ization energy, which is also added onto QMD data, and the reference ionization degree is semiclassically deter- in molecular dynamic simulations and electronic struc- mined [25]. VDWE has been treated as small as neg- ture calculations, respectively, because EOS (conductiv- ligible. As for another important parameter—pressure, ity) can only be modified within 5% (15%) for the selec- tionofhighernumberof kpoints. Thedensitiesadopted contributions from both interacting and noninteracting in our simulations range from 1.05 g/cm3 to 3.50 g/cm3 parts(ionsandfreeelectrons,separately)areincludedin the present EOS data. andtemperaturesbetween300and50000K,whichhigh- High precision EOS data are essential for understand- light the regime of principal Hugoniot. All the dynamic ing the electrical and optical properties of polystyrene simulations are lasted 6000 steps, and the time steps for under extreme conditions. The present EOS have been the integrations of atomic motion are selected according examined theoretically through the Rankine-Hugoniot todifferentdensities(temperatures)[24]. Then,the sub- (RH) equations, whichfollow fromconservationof mass, sequent 1000 step simulations are used to calculate EOS momentum, and energy across the front of shock waves. as running averages. The locus of points in (E,P,V)-space described by RH The nature of high pressure behavior for warm dense equations satisfy the following relations: polystyrene is dominated by a two-stage transition, that is, dissociations and ionizations. QMD simulations have 1 E −E = (P +P )(V −V ), (1) been demonstrated to be powerful in describing as well 1 0 2 1 0 0 asunderstandingchemicaldecompositionsandelectronic excitations of shock compressed materials, namely, the P −P =ρ u u , (2) 1 0 0 s p shocked atomic and electronic structures can be effec- tively treated within QMD. Importantly, however, we V =V (1−u /u ), (3) 1 0 p s should address that the calculated EOS data from di- rect QMD simulations on shockedpolystyrene should be where subscripts 0 and 1 present the initial and shocked correctedduetothefollowingaspects: (i)Atlaser-driven state, and E, P, V denote internalenergy, pressure,vol- temperaturesashighasof5∼6eV(∼10Mbar),the con- ume,respectively. Theu andu correspondtotheshock s p 3 and mass velocities of the material behind the shock front. As the starting point along the Hugoniot, the ini- tial density is ρ =1.05 g/cm3, and the internal energy is 0 E =−87.44 kJ/g at a temperature of 300 K. Compared 0 tothehighpressureofshockedstatesalongtheHugoniot, the initial pressure P can be treated approximately as 0 zero. We use smooth functions to fit the internal energy andpressureintermsoftemperatureatsampleddensity, and derive Hugoniot point from Eq. 1. The calculated Hugoniot data are listed in Table I. The principal Hugoniot curve is plotted in Fig. 1, where previous theoretical and experimental results are also shown for comparison. The EOS of polystyrene has been previously probed by gas gun experiments, which have the advantage of high precision, but the pressure hardly exceeds 50 GPa. Recently, high energy laser- FIG.2: (Color online) Pair correlation functions for four dif- driven experiments have detected pressures up to 4000 ferent points along the principal Hugoniot curve. (a) ∼ (d) GPa. Meanwhile, however, large error bars were intro- correspond to the starting point, 53.48 GPa, 74.54GPa, and duced at above 100 GPa [13, 14], and thus the use of 128.56 GPa, respectively. C-H, C-C, and H-H bonds are de- noted by blue, black, and red lines, respectively. Insets are low-precision EOS of the polymer to predict the behav- thesampledatomicstructureandchargedensitydistribution ior of NIF Ge-GDP ablator materials provides an un- ateachstate(grayballsforcarbonandwhiteballsforhydro- acceptable uncertainty. Then, high precision EOS (up gen). to 1000 GPa) for polystyrene was obtained by using α- quartz as an impedance-matching (IM) standard [8]. As shown in Fig. 1, previous data by Ozaki et al. [14] (IM polystyrene under shock compressions by using PCF, with an aluminum standard) and Cauble et al. [13] (ab- which is evaluated at equilibrium during the molecular solute data) show clearly stiffer behavior compared to dynamics simulations. Along the Hugoniot, four points those results with quartz standard [8, 15] and SESAME ofPCFare shownin Fig. 2. At the initial state [see Fig. model [26]. The possible reason for this stiff behavior is 2(a)], the ideal condensed polymer phase is indicated by likelyduetox-raypreheatingofthesamples,ashasbeen one peak for C-H bond, meanwhile, two peaks (π-type stated by those authors. Then, thicker pushers and low- and sp3-type) for C-C bond. With the increase of pres- Z ablators were used in the newly measured data from sure (∼20 GPa), phenyl decomposes, which is indicated Ozaki et al. [15] to reduce pre-heating of the samples. by that the two C-C peaks merge together (a transverse The experiments, whereIMwith a quartzstandardwere from π bond to sp2, sp3 like bond). At this stage dia- also used, show results that are closer to those by Bar- mondlikecarbonnanoparticles(withdefects)areformed. rios et al.. In our QMD simulations, we find that direct The present phenomenon agrees well with the behavior calculatedEOSwithoutcorrectionsmentionedaboveare ofshockcompressedbenzene[28,29]. Furtherincreaseof onlyvalidinthelowpressure(P<100GPaandT<1eV) pressureproducesmolecularandatomichydrogen,which regime, beyond which no Hugoniot points can be found isindicatedbythepeakinH-HPCFaround0.75˚A[Fig. from the pure QMD data. For higher pressures, due to 2(b)]. Decomposition of hydrocarbons continues until our introduction of corrections to direct QMD results, 128 GPa, where PCF shows no obvious evidence for the the calculated principal Hugoniot curve are greatly soft- existence of C-H bond. Whereas, for higher pressures, ened, and the results show good agreements with those the PCFs do not show visible difference, and are thus experiments where quartz standard was used. not presented here. SESAME model suggested that the Temperature, as has been focused as one of the most softening behavior of Hugoniot is caused by the chem- important parameters in experiments, is difficult to be ical decomposition of the polymer at 200 ∼ 400 GPa. measured because of the uncertainty in determining the Whereas, our calculations show that the decomposition optical-intensity loss for ultravioletpart of the spectrum pressure lies between 20 and 128 GPa, which is much in adiabatic or isentropic shock compressions, especially lower than that of SESAME. Note that the chemicalde- forthetemperatureexceedingseveraleV[27]. QMDsim- compositions do not contribute much to the softening ulations can provide efficient predictions for shock tem- behavior, but the ionizations really do. peratures. AshasbeenshowninFig. 1(rightpanel),our Havingclarifiedthe EOS,letusturnnowto study the calculatedHugoniottemperaturesareaccordantwithex- electricaland optical properties for shockedpolystyrene. perimental results up to 2 eV, and discrepancy emerges Based on Kubo-Greenwood formula, the electrical and athighertemperatures. Beyond2eV,thepredictedtem- optical properties can be obtained as described in Ref. peraturesfrombothexperimentsandSESAMEmodelat [28]. In order to get converged results, 30 indepen- agivenpressurearehighercomparedtoourcalculations. dent snapshots, which are selected during one molecu- We also examine the structural transition of lar dynamics simulation at given conditions, are picked 4 0.8 ing shock temperature. Along the Hugoniot, we show ) inFig. 3(rightpanel)thevariationofopticalreflectance -1 cm104 0.6 at 532 nm as a function of shock velocity. Optical prop- (-1 onductivity 103 Reflectance00..24 GuebryePtsiaKeus)po,eoatnfnoidgp∼ose5ltyt0esaat%dyl.irlwey[ne3eir0ne]chorwaebviastesehirnbvugeesedrn.eoflfTeex1hcp1etie∼vrre1iitms6iueelskntmstraeb/alyslcyh(O8isnzt0gau∼kdv1iiae7eld0t- c Present work Dc KBoarernioigs eett aall.. al. indicated reflectivity from 16%to 42%in us rangeof 102 0.0 22∼27km/s(300∼500GPa)[15]. Recentexperimentsby 5 10 15 20 25 30 35 5 10 15 20 25 30 35 Barrioset al. [8]suggestedadrasticincreaseinreflectiv- us (km/s) us (km/s) ity around the Hugoniot pressure of 100 GPa, followed FIG. 3: (Color online) Dc conductivity (left panel) and op- by saturated value of 40% around 250∼300 GPa. Dis- tical reflectivity at the wavelength of 532 nm (right panel). crepancies in reflectivities among experiments have been Previous measurements by Koenig et al. (green squares) [30] claimed to be arisen from probe-beam stability or from and Barrios et al. (upward triangles) [8] are also plotted. differencesindiagnosticconfigurations. AlongtheHugo- niot, our QMD calculations provide the general feature for the optical reflectance—steep increase (from 15% to up to calculate dynamic conductivity as running aver- 60%)followedbysaturation(at∼128GPa). Thechange ages. The dc conductivity σ , which follows from the dc static limit ω→0 of σ (ω) is then evaluated and plot- inopticalreflectivitywith pressurecanbe interpretedas 1 the gradual insulator-conductor transition at above 20 ted in Fig. 3 (left panel) as a function of shock ve- GPa, where the atoms strongly fluctuate with neighbors locity. Initially, σ increases rapidly with shock ve- dc todense,partiallyionizedplasmaathighpressuresabove locity up to 14 km/s (128 GPa) towards the formation 128 GPa. of metallic state of polystyrene. Then, one can find from Fig. 3 that σ keeps almost invariant and the dc In summary, we have carried out QMD simulations shock compressed polystyrene maintains its metallic be- to study the thermophysical properties for warm dense havior. Theonsetofmetallization(σ >103Ω−1cm−1)is dc polystyrene. The Hugoniot EOS data of polystyrene up observedat around10 km/s (∼50 GPa), where dissocia- to ∼790 GPa, which is in agreement with dynamic ex- tionofthepolymerandformationsofhydrogen(molecu- periments in a wide range of shock conditions, has been larandatomicspecies)governthecharacteristicsofwarm evaluated through QMD calculations and corrected by dense polystyrene. We stress here that the nonmetal to taking into account the atomic ionization. A two-stage metal transition is induced by gradual chemical decom- (dissociationandionization)transitiongovernsthe char- positions and thermal activations of the electronic state, acteristic of polystyrene under extreme conditions. Con- instead of atomic ionization, which is not observed un- tributionfromchemicaldecompositiondemonstratesthe til 128 GPa with respect to the charge density distri- steep increase of σ and optical reflectance observed at dc bution in the QMD simulations. In the same pressure 20∼128 GPa. While, soften behavior of the Hugoniot is range (20∼200 GPa), we observe largerσ compared to dc dominated by atomic ionizations for higher pressure. experimental measurement [30]. This overestimation of σ values could be attributed to the use of DFT-based This work was supported by NSFC under Grants No. dc molecular dynamics, which is known to underestimate 11005012andNo. 51071032,bytheNationalBasicSecu- band gaps in many systems. rity Research Program of China, by the National High- Optical reflectance, from which emissivity can be de- Tech ICF Committee of China, and by the Core Univer- rived, is of great interest in experiment for determin- sity Program. [1] M. D. Knudson, and M. P. Desjarlais, Phys. Rev. Lett. versity of California Press, Berkeley, 1980). 103, 225501 (2009). [10] R. G. McQueen et al.,in High-Velocity Impact Phenom- [2] W. J. Nellis, Rep.Prog. Phys. 69, 1479 (2006). ena, edited by R. Kinslow (Academic, New York,1970), [3] D.G. Hicks et al.,Phys. Rev.B 78, 174102 (2008). Chap. VII,Sec. II,p. 293. [4] D.G. Hicks et al.,Phys. Rev.B 79, 014112 (2009). [11] M. Van Thiel, J. Shaner, and E. Salinas, Lawrence Liv- [5] J.D.Lindl,Inertial Confinement Fusion: The Quest for ermore National Laboratory Report No. UCRL-50108, IgnitionandEnergyGainUsingIndirectDrive(Springer- 1977. Verlag, New York,1998). [12] A. V.Bushman, et al.,JETP Lett. 82, 895 (1996). [6] H.Cai et al.,Phys. Rev.Lett. 102, 245001 (2009). [13] R. Cauble, et al.,Phys. Plasmas 4, 1857 (1997). [7] S.W. Haan et al.,Eur. Phys. J. D 44, 249 (2007). [14] N. Ozaki, et al.,Phys. Plasmas 12, 124503 (2005). [8] M. A.Barrios et al.,Phys. Plasmas 17, 056307 (2010). [15] N. Ozaki, et al.,Phys. Plasmas 16, 062702 (2009). [9] LASL Shock Hugoniot Data, Los Alamos Series on Dy- [16] A. Kietzmann, R. Redmer, M. P. Desjarlais, and T. R. namic Material Properties, edited by S. P. Marsh (Uni- Mattsson, Phys.Rev.Lett. 101, 070401 (2008). 5 [17] G.Kresse andJ. Hafner, Phys.Rev.B47, R558 (1993). ionic number density), k T presents the kinetic energy, B [18] G. Kresse and J. Furthmu¨ller, Phys. Rev. B 54, 11169 and m is the ionic mass. (1996). [25] C. Wang, X. T. He, and P. Zhang, J. Appl. Phys. 108, [19] V.Recoules,et al.,Phys.Rev.Lett.102,075002 (2009). 044909 (2010). [20] J. P. Perdew, Electronic Structure of Solids (Akademie [26] S.P.LyonandJ.D.Johnson,LosAlamosNationalLab- Verlag, Berlin, 1991). oratory Report No. LA-CP-98-100, 1998. [21] P.E. Bl¨ochl, Phys.Rev.B 50, 17953 (1994). [27] Q. F. Chen, privatecommunication. [22] S.Nos´e, J. Chem. Phys. 81, 511 (1984). [28] C. Wang and P. Zhang, J. Appl. Phys. 107, 083502 [23] H.J. Monkhorst and J. D. Pack, Phys.Rev. B 13, 5188 (2010). (1976). [29] W. J. Nellis, D. C. Hamilton, and A. C. Mitchell, J. [24] The time steps have been taken as ∆t=a/20pkBT/m, Chem. Phys.115, 1015 (2001). wherea=(3/4πn )1/3 istheionicsphereradius(n isthe [30] M. Koenig, et al.,Phys. Plasmas 10, 3026 (203). i i