[Scientific journal] Science -.. July.20.2012

PERSPECTIVES certain critical fi eld strength, the Σ u state rep- beyond the hydrogen molecule, Lange et al. in the milli- or microsecond regime that are 3 resents the ground state of H 2 , which is then, have explored the He 2 spin-singlet ground already achieving fi elds of hundreds of teslas. 3 for even stronger fi elds, replaced by the Π u state that is weakly bound by dispersion Besides this, standard laboratory magnetic state ( 7, 8). forces if there is no external fi eld. In a strong fi elds of a few teslas are suffi cient to intro- 3 The Σ u state in fi eld-free space cannot fi eld on the order of a few 10 T, He 2 exhib- duce novel magnetic properties into highly 5 accommodate any vibrational states, nor can its a perpendicular orientation while its bond excited Rydberg states ( 3, 9), for which the it even for strong fi elds for molecular and distance decreases by a factor of 3 and its weakened Coulomb force is comparable to magnetic fi eld axes aligned parallel. How- dissociation energy increases several hun- the paramagnetic or diamagnetic forces, or ever, for a perpendicular magnetic fi eld, dra- dred times. Here, the same mechanism is at both. In this respect, the very weakly bound 3 matic changes were found by Lange et al. work as in the case of the Σ u state of H 2 , trig- Rydberg molecules with a huge bond length through extensive quantum chemical com- gered by the change of character of the anti- and a large electric dipole moment, originally putations. The antibonding orbitals respon- bonding orbitals. predicted by Greene et al. ( 10) and recently sible for the unbound character of the Σ u Atoms, molecules, and condensed matter found for a specifi c class experimentally by 3 state undergo a metamorphosis for a per- systems exposed to strong magnetic fi elds Pfau and co-workers ( 11), are of particular pendicular fi eld that stabilizes the molecule. represent a fascinating topic, and the work interest and might be an ideal system to probe The fi eld-induced bonding in this case is nei- by Lange et al. has added a key bonding strong magnetic fi eld effects in the laboratory. ther of covalent nor of ionic nature, but rep- mechanism. The competition between the resents a different bonding mechanism with anisotropy-introducing magnetic fi eld and References paramagnetic character. The latter is based the attractive and repulsive Coulomb forces 1. H. Ruder, G. Wunner, H. Herold, F. Geyer, Atoms in Strong Magnetic Fields (Springer-Verlag, Berlin, 1994). on the lowering of the kinetic energy and the is responsible for an enormous complex- 2. P. Schmelcher, W. Schweizer, Eds., Atoms and Molecules antiparallel orientation of the acquired spa- ity and diversity of the microscopic behav- in Strong External Fields (Plenum, New York, 1998). on July 19, 2012 tial angular momentum in the presence of a ior. Indeed, the existing investigations show 3. H. Friedrich, H. Wintgen, Phys. Rep. 183, 37 (1989). strong magnetic fi eld. A simple model for the that different excited states of a molecule, 4. K. K. Lange et al., Science 337, 327 (2012). 5. P. Schmelcher, L. S. Cederbaum, Phys. Rev. A 37, 672 fi eld-dependent molecular orbitals confi rms or ground states of similar molecules, can (1988). this basic mechanism. behave in a vastly different manner, and pro- 6. E. I. Tellgren et al., J. Chem. Phys. 129, 154114 (2008). The paramagnetic bonding mechanism vide a fi rst look at a largely unexplored area: 7. Yu. P. Kravchenko, M. A. Liberman, Phys. Rev. A 56, R2510 (1997). is not based on the correlations of the elec- the world of magnetized matter. 8. T. Detmer et al., Phys. Rev. A 57, 1767 (1998). tronic motion, although they slightly change It might seem that this world is hidden 9. T. F. Gallagher, Rydberg Atoms (Cambridge Univ. Press, it quantitatively. It already appears within from us by the requirement of enormous Cambridge, 1994). www.sciencemag.org 10. C. H. Greene et al., Phys. Rev. Lett. 85, 2458 (2000). Hartree-Fock theory, in which the electrons laboratory field strengths. However, avail- 11. V. Bendkowsky et al., Nature 458, 1005 (2009). move in their respective averaged fi eld. To able fi eld strengths are progressively increas- confi rm the appearance of this mechanism ing, in particular for pulsed magnetic fi elds 10.1126/science.1224869 CHEMISTRY Nanometer-scale polymeric materials are Nanomaterials for Drug Delivery increasingly used to surmount the barriers Downloaded from faced by drugs and vaccines on their way to their site of action. 1 Jeffrey A. Hubbell and Ashutosh Chilkoti 2 ll drugs face several transport bar- harsh acidic environment of endolysosomes, laries) ( 1). One challenge has been to atten- riers on their tortuous journey within which biomolecular drugs such as pro- uate the fi ltration rate of molecules from the Afrom their site of introduction to teins and oligonucleotides may be inactivated bloodstream, especially into the kidney, and to their molecular site of action. Critical barri- or degraded. Other barriers are the nuclear avoid premature clearance by the RES. ers include rapid fi ltration in the kidney and membrane and the multiple drug resistance To address this, the hydrodynamic clearance via the reticulo-endothelial sys- mechanisms that pathological cells can radius of protein drugs has been increased tem (RES)—particularly for drugs that spend develop. Recent studies illustrate some par- to ~10 nm by grafting the hydrophilic poly- a lot of time in the bloodstream—as well as ticularly promising ways in which nanoma- mer poly(ethylene glycol). This technology, transport from the bloodstream to target cells terials as drug or vaccine carriers can assist referred to as PEGylation, is based on the within tissues. At the tissue or cellular tar- in navigating these barriers, with a particular twin observations that larger hydrophilic mol- get, the drug must cross the plasma mem- focus on administration by injection. ecules are fi ltered by the kidney more slowly brane, and within the cell, it must escape the For drugs administered through injection, than smaller ones and that PEGylated proteins the vasculature provides both a road to the more successfully evade premature clearance destination—the site of disease—and many through the RES. Gao et al. recently created 1 Institute of Bioengineering, School of Life Sciences and School of Engineering, and Institute of Chemical Sciences detours for the drug to be lost in transit. Trans- site-specifi c, PEG-like conjugates by polym- and Engineering, School of Basic Sciences, Ecole Polytech- port in the vasculature and in the tissues that erizing very long polymer chains from one of nique Fédérale de Lausanne, CH-1015 Lausanne, Switzer- those vessels perfuse depends on convection the protein termini. This approach increased 2 land. Department of Biomedical Engineering and Center in the circulation as well as diffusion and con- the hydrodynamic radius from 3 nm for the for Biologically Inspired Materials and Materials Systems, Duke University, Durham, NC 27708, USA. E-mail: jeffrey. vection in the tissue interstitium (i.e., between native protein to >20 nm for the conjugate ( 2). hubbell@epfl .ch the blood capillaries and the lymphatic capil- In a similar, genetically encoded approach, www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 303 Published by AAAS

Nanomaterials for Drug Delivery Jeffrey A. Hubbell and Ashutosh Chilkoti Science 337 , 303 (2012); DOI: 10.1126/science.1219657 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/303.full.html This article cites 17 articles , 7 of which can be accessed free: http://www.sciencemag.org/content/337/6092/303.full.html#ref-list-1 www.sciencemag.org This article appears in the following subject collections: Chemistry http://www.sciencemag.org/cgi/collection/chemistry Medicine, Diseases Downloaded from http://www.sciencemag.org/cgi/collection/medicine Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

PERSPECTIVES certain critical fi eld strength, the Σ u state rep- beyond the hydrogen molecule, Lange et al. in the milli- or microsecond regime that are 3 resents the ground state of H 2 , which is then, have explored the He 2 spin-singlet ground already achieving fi elds of hundreds of teslas. 3 for even stronger fi elds, replaced by the Π u state that is weakly bound by dispersion Besides this, standard laboratory magnetic state ( 7, 8). forces if there is no external fi eld. In a strong fi elds of a few teslas are suffi cient to intro- 3 The Σ u state in fi eld-free space cannot fi eld on the order of a few 10 T, He 2 exhib- duce novel magnetic properties into highly 5 accommodate any vibrational states, nor can its a perpendicular orientation while its bond excited Rydberg states ( 3, 9), for which the it even for strong fi elds for molecular and distance decreases by a factor of 3 and its weakened Coulomb force is comparable to magnetic fi eld axes aligned parallel. How- dissociation energy increases several hun- the paramagnetic or diamagnetic forces, or ever, for a perpendicular magnetic fi eld, dra- dred times. Here, the same mechanism is at both. In this respect, the very weakly bound 3 matic changes were found by Lange et al. work as in the case of the Σ u state of H 2 , trig- Rydberg molecules with a huge bond length through extensive quantum chemical com- gered by the change of character of the anti- and a large electric dipole moment, originally putations. The antibonding orbitals respon- bonding orbitals. predicted by Greene et al. ( 10) and recently sible for the unbound character of the Σ u Atoms, molecules, and condensed matter found for a specifi c class experimentally by 3 state undergo a metamorphosis for a per- systems exposed to strong magnetic fi elds Pfau and co-workers ( 11), are of particular pendicular fi eld that stabilizes the molecule. represent a fascinating topic, and the work interest and might be an ideal system to probe The fi eld-induced bonding in this case is nei- by Lange et al. has added a key bonding strong magnetic fi eld effects in the laboratory. ther of covalent nor of ionic nature, but rep- mechanism. The competition between the resents a different bonding mechanism with anisotropy-introducing magnetic fi eld and References paramagnetic character. The latter is based the attractive and repulsive Coulomb forces 1. H. Ruder, G. Wunner, H. Herold, F. Geyer, Atoms in Strong Magnetic Fields (Springer-Verlag, Berlin, 1994). on the lowering of the kinetic energy and the is responsible for an enormous complex- 2. P. Schmelcher, W. Schweizer, Eds., Atoms and Molecules antiparallel orientation of the acquired spa- ity and diversity of the microscopic behav- in Strong External Fields (Plenum, New York, 1998). on July 19, 2012 tial angular momentum in the presence of a ior. Indeed, the existing investigations show 3. H. Friedrich, H. Wintgen, Phys. Rep. 183, 37 (1989). strong magnetic fi eld. A simple model for the that different excited states of a molecule, 4. K. K. Lange et al., Science 337, 327 (2012). 5. P. Schmelcher, L. S. Cederbaum, Phys. Rev. A 37, 672 fi eld-dependent molecular orbitals confi rms or ground states of similar molecules, can (1988). this basic mechanism. behave in a vastly different manner, and pro- 6. E. I. Tellgren et al., J. Chem. Phys. 129, 154114 (2008). The paramagnetic bonding mechanism vide a fi rst look at a largely unexplored area: 7. Yu. P. Kravchenko, M. A. Liberman, Phys. Rev. A 56, R2510 (1997). is not based on the correlations of the elec- the world of magnetized matter. 8. T. Detmer et al., Phys. Rev. A 57, 1767 (1998). tronic motion, although they slightly change It might seem that this world is hidden 9. T. F. Gallagher, Rydberg Atoms (Cambridge Univ. Press, it quantitatively. It already appears within from us by the requirement of enormous Cambridge, 1994). www.sciencemag.org 10. C. H. Greene et al., Phys. Rev. Lett. 85, 2458 (2000). Hartree-Fock theory, in which the electrons laboratory field strengths. However, avail- 11. V. Bendkowsky et al., Nature 458, 1005 (2009). move in their respective averaged fi eld. To able fi eld strengths are progressively increas- confi rm the appearance of this mechanism ing, in particular for pulsed magnetic fi elds 10.1126/science.1224869 CHEMISTRY Nanometer-scale polymeric materials are Nanomaterials for Drug Delivery increasingly used to surmount the barriers Downloaded from faced by drugs and vaccines on their way to their site of action. 1 Jeffrey A. Hubbell and Ashutosh Chilkoti 2 ll drugs face several transport bar- harsh acidic environment of endolysosomes, laries) ( 1). One challenge has been to atten- riers on their tortuous journey within which biomolecular drugs such as pro- uate the fi ltration rate of molecules from the Afrom their site of introduction to teins and oligonucleotides may be inactivated bloodstream, especially into the kidney, and to their molecular site of action. Critical barri- or degraded. Other barriers are the nuclear avoid premature clearance by the RES. ers include rapid fi ltration in the kidney and membrane and the multiple drug resistance To address this, the hydrodynamic clearance via the reticulo-endothelial sys- mechanisms that pathological cells can radius of protein drugs has been increased tem (RES)—particularly for drugs that spend develop. Recent studies illustrate some par- to ~10 nm by grafting the hydrophilic poly- a lot of time in the bloodstream—as well as ticularly promising ways in which nanoma- mer poly(ethylene glycol). This technology, transport from the bloodstream to target cells terials as drug or vaccine carriers can assist referred to as PEGylation, is based on the within tissues. At the tissue or cellular tar- in navigating these barriers, with a particular twin observations that larger hydrophilic mol- get, the drug must cross the plasma mem- focus on administration by injection. ecules are fi ltered by the kidney more slowly brane, and within the cell, it must escape the For drugs administered through injection, than smaller ones and that PEGylated proteins the vasculature provides both a road to the more successfully evade premature clearance destination—the site of disease—and many through the RES. Gao et al. recently created 1 Institute of Bioengineering, School of Life Sciences and School of Engineering, and Institute of Chemical Sciences detours for the drug to be lost in transit. Trans- site-specifi c, PEG-like conjugates by polym- and Engineering, School of Basic Sciences, Ecole Polytech- port in the vasculature and in the tissues that erizing very long polymer chains from one of nique Fédérale de Lausanne, CH-1015 Lausanne, Switzer- those vessels perfuse depends on convection the protein termini. This approach increased 2 land. Department of Biomedical Engineering and Center in the circulation as well as diffusion and con- the hydrodynamic radius from 3 nm for the for Biologically Inspired Materials and Materials Systems, Duke University, Durham, NC 27708, USA. E-mail: jeffrey. vection in the tissue interstitium (i.e., between native protein to >20 nm for the conjugate ( 2). hubbell@epfl .ch the blood capillaries and the lymphatic capil- In a similar, genetically encoded approach, www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 303 Published by AAAS

PERSPECTIVES Schellenberger et al. expressed therapeutic ABC proteins as fusions with a long, unstructured, hydrophilic polypeptide that seems to share many of the valuable attributes of PEG. Such fusions led to projected extension of circula- tion half-lives of exenatide, an important pep- tide drug for type 2 diabetes, from 2 hours to >100 hours ( 3). The tumor vasculature is much leakier to colloidal objects than the healthy vasculature, allowing micelles and associated drugs to accumulate there. For tumor-targeting protein drugs, the approaches described above can lead to enhancement in circulation lifetime by a factor of 10 and improved accumulation in An example of nanomaterials design. A complex structure helps siRNA to be carried to the tumor cell sur- face, internalized by endocytosis, and then released from the endolysosome into the cytoplasm. (A) Chains tumors relative to the native protein ( 2). of positively charged cyclodextrins (light blue) form electrostatic nanoparticles with the negatively charged Even larger self-assembled polymer siRNA (purple). These complexes are administered by injection. (B) In the bloodstream, nanoparticle clear- micelles are being developed to target specifi c ance is blocked by PEG chains terminated with hydrophobic moieties that bind as guests in the hydropho- tissues and even subcellular compartments in bic host center of the cyclodextrin. The nanoparticles are further functionalized with transferrin (green) as a tumors. In one such approach, MacKay et al. tumor-targeting moiety (many tumor cells overexpress the transferrin receptor), using the PEG as a linker. (C) ( 4) conjugated a hydrophobic cancer drug via In the tumor cell, the particle dissociates to release the siRNA into the cytoplasm. on July 19, 2012 an acid labile linker to an elastin-like poly- peptide. Driven by the hydrophobicity of the drugs such as paclitaxel ( 8). sosomal membranes and releases the drug conjugated drug, these polymers self-assem- Other features of tumors have been used into the cytoplasm ( 12). Additional features ble into micelles <50 nm in diameter. These for targeting to amplify a signal at the tumor to promote stability and targeting have been micelles then localize to tumors. In the acidic site. For example, von Maltzahn et al. have designed into this Lego-like nanomaterial (see environment of the endolysosomes of the targeted nanomaterials to molecular features the fi gure). Initial evidence for gene knock- tumor cells, the drug is cleaved to reach its unique to the tumor endothelium to induce down in tumors in patients supports the func- nuclear target, leading to a potent therapeutic coagulation there; the resulting coagulum tion of this nanomaterial design ( 12). www.sciencemag.org response in an animal model. was then targeted with drugs that bind bio- Vaccines are another area in which size In a remarkable feat of multifunctional chemically to the nascent clot ( 9). Because matters. In the interstitium of most tissues, polymer design, Murakami et al. have devel- coagulation is autocatalytic, substantial sig- fl uid fl ows between the blood capillaries and oped micelles that form a complex with a plat- nal amplifi cation occurs between the initial the draining lymphatics. After injection (e.g., inum-based DNA-targeting anticancer drug tumor targeting step and the fi nal binding of intradermally), nanomaterials that are suffi - – in a pH- and Cl -dependent manner, with drug the drug to the coagulum. Here, biomolecu- ciently small (sub-100 nm) to avoid entrap- release and activation favored at acidic pH and lar recognition was used in the fi rst step to ment in the tissue interstitium are effi ciently high Cl ( 5). The latter conditions are unique localize the procoagulant nanoparticles to transported by this interstitial fl ow into the Downloaded from – to late endosomes (just prior to endosomal the tumor endothelium through binding of draining lymphatics and lymph nodes. Here fusion with lysosomes) and lysosomes. Late a nanoparticle-conjugated domain of tissue they can be collected by lymph node–resident endosomes are localized near the cell nucleus, factor receptor, targeting the tissue factor that dendritic cells—the fi rst key cellular player in allowing directed delivery of the drug close to is naturally up-regulated on the pathological generating an immune response ( 13). its DNA target, thus escaping the cytoplasmic tumor endothelium ( 9). Particulate antigen from outside the cell resistance mechanisms that often inactivate Annexin 1 is also highly expressed on the is most frequently processed and presented the free drug. tumor endothelium. A short peptide sequence on class II major histocompatibility complex Other micelle approaches abound ( 6, 7). has recently been discovered that serves as (MHC) molecules, which results in humoral Notable among these is an optimization of a ligand for this marker, providing effi cient immunity. The resulting antibodies can be micelle shape for highly prolonged circula- drug localization to the tumor endothelium useful for preventing disease, for example by tion ( 8). Polymer micelles for use in drug ( 10). This ligand can presumably be used to binding and neutralizing a virus, but are less delivery are typically formed by self-assem- target micelles to tumor vasculature, as has useful for killing aberrant cells, such as cells bly of block copolymers with a hydrophobic been done with other peptides to target sites already infected by the virus. To achieve this, domain (which drives polymer association of atherosclerotic plaque after intravascular the antigen must be presented on MHC class in the presence of water) and a hydrophilic injection ( 11). I molecules, leading to the induction of lym- domain (which restricts this association to Biomacromolecular drugs typically can- phocytes that can kill virally infected cells and the nanometer dimension). The relative vol- not permeate through membranes to access even tumor cells. ume of the two polymer segments controls targets in the cytoplasm. An exciting exam- Antigen presented on MHC class I usually the shape of the resulting micelle. Chris- ple of nanomaterials to penetrate these bar- derives from protein within the cell, such as tian et al. have shown that polymer architec- riers is the delivery of small interfering viral proteins. However, under some circum- tures that lead to cylindrical micelles enable RNA (siRNA) using cationic nanocarriers. stances, antigen from extracellular sources highly prolonged circulation in the blood As the cell attempts to neutralize the basic can be shuttled to MHC class I. This “anti- CREDIT: P . HUEY/SCIENCE stream, thereby reducing off-target toxicity charge of the nanomaterials, the resulting gen cross-presentation” is particularly effec- and increasing effi cacy in tumor killing with osmotic imbalance destabilizes the endoly- tive when antigen is conjugated to nanopar- 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org 304 Published by AAAS

PERSPECTIVES ticles by a reducible disulfide bond ( 14), develop hollow nanocontainers with DNA without formation of untoward metabolites. which is cleaved in the reductive environment aptamer–based locks ( 17). These aptamers within the endosome, resulting in strong cel- also bind to particular biomolecular signals; References 1. J. W. Baish et al., Proc. Natl. Acad. Sci. U.S.A. 108, 1799 lular immunity that can kill aberrant cells. when the aptamer detects the binding antigen (2011). Such small nanoparticles can also very effec- key, the nanocontainer is opened and the pay- 2. W. Gao, W. Liu, T. Christensen, M. R. Zalutsky, A. Chilkoti, Proc. Natl. Acad. Sci. U.S.A. 107, 16432 (2010). tively target antigen-presenting cells after load is released. 3. V. Schellenberger et al., Nat. Biotechnol. 27, 1186 (2009). pulmonary administration to induce a potent These examples show how highly func- 4. J. A. MacKay et al., Nat. Mater. 8, 993 (2009). immune response in the lung and at other tional multicomponent polymeric nano- 5. M. Murakami et al., Sci. Transl. Med. 3, 64ra2 (2011). 6. R. A. Petros, J. M. DeSimone, Nat. Rev. Drug Discov. 9, mucosal surfaces, to protect from viral infec- structures can guide drugs to tissues, coax 615 (2010). tions such as infl uenza ( 15). Other nanoma- them through the biological barriers at the 7. J. W. Yoo et al., Nat. Rev. Drug Discov. 10, 521 (2011). terial designs have also been shown to lead surfaces of and within cells, and escape 8. D. A. Christian et al., Mol. Pharm. 6, 1343 (2009). 9. G. von Maltzahn et al., Nat. Mater. 10, 545 (2011). to strong cross-presentation, including mul- drug clearance and drug resistance. How- 10. S. Hatakeyama et al., Proc. Natl. Acad. Sci. U.S.A. 108, tilamellar vesicles, stabilized by interlamel- ever, in translation to human clinical trials, 19587 (2011). lar cross-linking, carrying antigen as well as the devil is often in the details: Materials 11. D. Peters et al., Proc. Natl. Acad. Sci. U.S.A. 106, 9815 (2009). adjuvant molecules ( 16). must be developed that maintain their sta- 12. M. E. Davis et al., Nature 464, 1067 (2010). The level of sophistication of designing bility, size distribution, and targeting speci- 13. S. T. Reddy et al., Nat. Biotechnol. 25, 1159 (2007). nanocarriers to be functionally responsive fi city in the complex and concentrated pro- 14. S. Hirosue et al., Vaccine 28, 7897 (2010). 15. C. Nembrini et al., Proc. Natl. Acad. Sci. U.S.A. 108, has recently been taken to new levels with tein environment of the body, that can be E989 (2011). compelling demonstrations in vitro. Doug- eliminated from the body at the same rate as 16. J. J. Moon et al., Nat. Mater. 10, 243 (2011). las et al. have used DNA origami—intri- planned administration to avoid accumula- 17. S. M. Douglas et al., Science 335, 831 (2012). cately folded nanoscale DNA structures—to tion, and that can be processed by the body 10.1126/science.1219657 on July 19, 2012 GEOCHEMISTRY Model studies point to a larger role of the The Marine Sulfur Cycle, Revisited sulfur cycle in maintaining atmospheric oxygen levels than previously recognized. Matthew T. Hurtgen www.sciencemag.org ulfur takes part in many biogeochemi- the Earth’s surface. As oxygen levels increase, However, over the past decade, scientists cal reactions that affect the global car- so, too, should sulfi de oxidation on land and have found that Phanerozoic sulfate concen- Sbon and oxygen cycles. On short time the resulting fl ux of sulfate via rivers to the trations varied considerably ( 3). Sulfate con- scales and in the absence of oxygen, many ocean. Because oxygen levels at the Earth’s centrations in the modern ocean are 28 mM, microbes use organic carbon to reduce sul- surface are thought to have been relatively but the composition of fl uid inclusions in fate to sulfi de, which may then react with iron stable since ~580 million years ago, marine halite suggests that Phanerozoic levels may to form pyrite. On much longer time scales, sulfate concentrations would be expected to have fl uctuated between ~5 and 25 mM. The Downloaded from the net addition of oxygen to the atmosphere have been largely unchanged through the Pha- time resolution of these data is quite low, but through organic carbon burial promotes sul- nerozoic (since ~542 million years ago). this work suggests that the fl ux of sulfur into fide oxidation on land and increases the amount of sulfate carried by rivers to the oceans. Two reports in this issue, by Wort- Sulfate weathering input via rivers; oxidative weathering of mann and Paytan ( 1) on page 334 and Halevy sulfides (flux 1) and dissolution of calcium sulfates (flux 2) et al. ( 2) on page 331, show that on time scales of millions of years, changes in the for- SO 2– 1 2 4 mation and dissolution of evaporite minerals can strongly impact marine chemistry and the 3 Sulfate outputs: calcium sulfate precipitation carbon cycle (see the fi gure). On longer time Evaporite (calcium sulfate) precipitation (flux 3) and microbial sulfate reduction and scales, the sulfur cycle played a much greater pyrite burial (flux 4) role in regulating atmospheric oxygen levels f = pyrite S burial than previously recognized. pyr pyrite S burial + evaporite S burial Earth scientists have long appreciated that marine sulfate concentrations have changed substantially through time as a result of vary- 4 ing sulfur fl uxes to and from the ocean. The story has generally been framed in terms of CREDIT: P . HUEY/SCIENCE changing inputs. The sulfate content of the Reassessing marine sulfur fl uxes. The main sulfur inputs to the ocean are the dissolution of evaporite ocean should track the oxidation history of minerals and oxidative weathering of pyrite; the main sulfur outputs are the formation and burial of evapo- rites and pyrite. During the Phanerozoic, pyrite deposition and weathering have dominated these fl uxes ( 2). However, episodic, massive changes in evaporite formation ( 1, 2) and dissolution ( 1) strongly impact the Earth and Planetary Sciences, Northwestern University, Evan- marine sulfur cycle on million-year time scales. ston, IL 60208, USA. E-mail: [email protected] www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 305 Published by AAAS

The Marine Sulfur Cycle, Revisited Matthew T. Hurtgen Science 337 , 305 (2012); DOI: 10.1126/science.1225461 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/305.full.html A list of selected additional articles on the Science Web sites related to this article can be found at: www.sciencemag.org http://www.sciencemag.org/content/337/6092/305.full.html#related This article cites 6 articles , 4 of which can be accessed free: http://www.sciencemag.org/content/337/6092/305.full.html#ref-list-1 Downloaded from This article appears in the following subject collections: Geochemistry, Geophysics http://www.sciencemag.org/cgi/collection/geochem_phys Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

PERSPECTIVES ticles by a reducible disulfide bond ( 14), develop hollow nanocontainers with DNA without formation of untoward metabolites. which is cleaved in the reductive environment aptamer–based locks ( 17). These aptamers within the endosome, resulting in strong cel- also bind to particular biomolecular signals; References 1. J. W. Baish et al., Proc. Natl. Acad. Sci. U.S.A. 108, 1799 lular immunity that can kill aberrant cells. when the aptamer detects the binding antigen (2011). Such small nanoparticles can also very effec- key, the nanocontainer is opened and the pay- 2. W. Gao, W. Liu, T. Christensen, M. R. Zalutsky, A. Chilkoti, Proc. Natl. Acad. Sci. U.S.A. 107, 16432 (2010). tively target antigen-presenting cells after load is released. 3. V. Schellenberger et al., Nat. Biotechnol. 27, 1186 (2009). pulmonary administration to induce a potent These examples show how highly func- 4. J. A. MacKay et al., Nat. Mater. 8, 993 (2009). immune response in the lung and at other tional multicomponent polymeric nano- 5. M. Murakami et al., Sci. Transl. Med. 3, 64ra2 (2011). 6. R. A. Petros, J. M. DeSimone, Nat. Rev. Drug Discov. 9, mucosal surfaces, to protect from viral infec- structures can guide drugs to tissues, coax 615 (2010). tions such as infl uenza ( 15). Other nanoma- them through the biological barriers at the 7. J. W. Yoo et al., Nat. Rev. Drug Discov. 10, 521 (2011). terial designs have also been shown to lead surfaces of and within cells, and escape 8. D. A. Christian et al., Mol. Pharm. 6, 1343 (2009). 9. G. von Maltzahn et al., Nat. Mater. 10, 545 (2011). to strong cross-presentation, including mul- drug clearance and drug resistance. How- 10. S. Hatakeyama et al., Proc. Natl. Acad. Sci. U.S.A. 108, tilamellar vesicles, stabilized by interlamel- ever, in translation to human clinical trials, 19587 (2011). lar cross-linking, carrying antigen as well as the devil is often in the details: Materials 11. D. Peters et al., Proc. Natl. Acad. Sci. U.S.A. 106, 9815 (2009). adjuvant molecules ( 16). must be developed that maintain their sta- 12. M. E. Davis et al., Nature 464, 1067 (2010). The level of sophistication of designing bility, size distribution, and targeting speci- 13. S. T. Reddy et al., Nat. Biotechnol. 25, 1159 (2007). nanocarriers to be functionally responsive fi city in the complex and concentrated pro- 14. S. Hirosue et al., Vaccine 28, 7897 (2010). 15. C. Nembrini et al., Proc. Natl. Acad. Sci. U.S.A. 108, has recently been taken to new levels with tein environment of the body, that can be E989 (2011). compelling demonstrations in vitro. Doug- eliminated from the body at the same rate as 16. J. J. Moon et al., Nat. Mater. 10, 243 (2011). las et al. have used DNA origami—intri- planned administration to avoid accumula- 17. S. M. Douglas et al., Science 335, 831 (2012). cately folded nanoscale DNA structures—to tion, and that can be processed by the body 10.1126/science.1219657 on July 19, 2012 GEOCHEMISTRY Model studies point to a larger role of the The Marine Sulfur Cycle, Revisited sulfur cycle in maintaining atmospheric oxygen levels than previously recognized. Matthew T. Hurtgen www.sciencemag.org ulfur takes part in many biogeochemi- the Earth’s surface. As oxygen levels increase, However, over the past decade, scientists cal reactions that affect the global car- so, too, should sulfi de oxidation on land and have found that Phanerozoic sulfate concen- Sbon and oxygen cycles. On short time the resulting fl ux of sulfate via rivers to the trations varied considerably ( 3). Sulfate con- scales and in the absence of oxygen, many ocean. Because oxygen levels at the Earth’s centrations in the modern ocean are 28 mM, microbes use organic carbon to reduce sul- surface are thought to have been relatively but the composition of fl uid inclusions in fate to sulfi de, which may then react with iron stable since ~580 million years ago, marine halite suggests that Phanerozoic levels may to form pyrite. On much longer time scales, sulfate concentrations would be expected to have fl uctuated between ~5 and 25 mM. The Downloaded from the net addition of oxygen to the atmosphere have been largely unchanged through the Pha- time resolution of these data is quite low, but through organic carbon burial promotes sul- nerozoic (since ~542 million years ago). this work suggests that the fl ux of sulfur into fide oxidation on land and increases the amount of sulfate carried by rivers to the oceans. Two reports in this issue, by Wort- Sulfate weathering input via rivers; oxidative weathering of mann and Paytan ( 1) on page 334 and Halevy sulfides (flux 1) and dissolution of calcium sulfates (flux 2) et al. ( 2) on page 331, show that on time scales of millions of years, changes in the for- SO 2– 1 2 4 mation and dissolution of evaporite minerals can strongly impact marine chemistry and the 3 Sulfate outputs: calcium sulfate precipitation carbon cycle (see the fi gure). On longer time Evaporite (calcium sulfate) precipitation (flux 3) and microbial sulfate reduction and scales, the sulfur cycle played a much greater pyrite burial (flux 4) role in regulating atmospheric oxygen levels f = pyrite S burial than previously recognized. pyr pyrite S burial + evaporite S burial Earth scientists have long appreciated that marine sulfate concentrations have changed substantially through time as a result of vary- 4 ing sulfur fl uxes to and from the ocean. The story has generally been framed in terms of CREDIT: P . HUEY/SCIENCE changing inputs. The sulfate content of the Reassessing marine sulfur fl uxes. The main sulfur inputs to the ocean are the dissolution of evaporite ocean should track the oxidation history of minerals and oxidative weathering of pyrite; the main sulfur outputs are the formation and burial of evapo- rites and pyrite. During the Phanerozoic, pyrite deposition and weathering have dominated these fl uxes ( 2). However, episodic, massive changes in evaporite formation ( 1, 2) and dissolution ( 1) strongly impact the Earth and Planetary Sciences, Northwestern University, Evan- marine sulfur cycle on million-year time scales. ston, IL 60208, USA. E-mail: [email protected] www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 305 Published by AAAS

PERSPECTIVES and out of the ocean has fl uctuated much reduced marine sulfate levels. The authors larger role in regulating Phanerozoic atmo- 34 more than previously recognized. attribute a 5‰ positive δ S sulfate shift about spheric oxygen levels than previously recog- To better constrain sulfur fl uxes to and 50 million years ago to an abrupt increase nized—perhaps by as much as 50%. How- from the ocean, Wortmann and Paytan use in marine sulfate concentrations as a result ever, both studies recognize the importance a model that tracks the mass and sulfur iso- of large-scale dissolution of freshly exposed of episodic evaporite burial on the sulfur tope composition of fl uxes into and out of the evaporites; they argue that the higher sulfate cycle, and Wortmann and Paytan show that ocean. Most sulfur in the ocean is brought concentrations led to more pyrite burial. large-scale deposition and dissolution of there by rivers as a result of oxidative weath- Halevy et al. likewise study past sulfur these evaporites over relatively short geo- ering of sulfi des and dissolution of calcium fl uxes to and from the ocean, but over lon- logic time scales can have an enormous sulfates. Sulfur is removed from the ocean ger time scales (the past 500 million years). impact on marine sulfate concentrations, through microbial sulfate reduction fol- They quantify sulfate evaporite burial rates pyrite burial rates, and the carbon cycle. lowed by pyrite burial, and by calcium sul- through time with the help of a database The studies by Wortmann and Paytan fate (evaporite) deposition. Evaporite depo- that contains surface and subsurface infor- and Halevy et al. highlight the dynamic sition does not have a large effect on the sul- mation on sedimentary strata from loca- nature of the Phanerozoic sulfur cycle and fur isotopic composition, but microbial sul- tions across North America and the Carib- its importance in regulating the chemistry of fate reduction does, because microbes pre- bean. They then scale these rates to obtain a the ocean-atmosphere system. Future work 32 fer to use S over S during the production global estimate. should focus on better constraining marine 34 of sulfi de ( 4). As a result, the sulfur isotope The results indicate that sulfate burial sulfate levels and global sulfate burial rates 34 composition (δ S sulfate ) of river inputs today rates were higher than previously estimated, through time, the relationship between is ~6 per mil (‰) ( 5), but that of seawater but also greatly variable. When Halevy et al. marine sulfate concentrations and pyrite sulfate is ~20‰. integrated these improved evaporite burial burial rates, and the role of sulfur in regulat- on July 19, 2012 34 Variations in δ S sulfate are traditionally fl uxes with seawater sulfate concentration ing ocean nutrient cycling, the carbon cycle, interpreted to reflect changes in the total estimates ( 3) and sulfur isotope constraints and ultimately climate. amount of sulfur buried in ocean sediments ( 6, 7), their calculations implied that Pha- as pyrite—a parameter referred to as f pyr . nerozoic f pyr values were ~100% higher on References Wortmann and Paytan conclude that a average than previously recognized. These 1. U. G. Wortmann, A. Paytan, Science 337, 334 (2012). 2. I. Halevy et al., Science 337, 331 (2012). 34 5‰ negative δ S sulfate shift in ~120-mil- surprisingly high and constant pyrite burial 3. T. K. Lowenstein et al., Geology 31, 857 (2003). lion-year-old rocks was caused by massive outputs must have been balanced by equally 4. D. E. Canfi eld, Geochim. Cosmochim. Acta 65, 1117 seawater sulfate removal that accompanied high and constant inputs of sulfate to the (2001). www.sciencemag.org 5. M. A. Arthur, Encyclopedia of Volcanoes (Academic Press, large-scale evaporite deposition during the ocean via sulfi de oxidation (weathering). San Diego, CA, 2000). opening of the South Atlantic Ocean. Evap- The relatively high and constant rates of 6. A. Paytan et al., Science 304, 1663 (2004). orite formation does not cause sulfur iso- pyrite weathering and burial over long geo- 7. A. Kampschulte, H. Strauss, Chem. Geol. 204, 255 tope fractionation; in the model, the nega- logic time scales, identifi ed by Halevy et al., (2004). 34 tive δ S sulfate shift is driven by lower pyrite suggest that the consumption and produc- burial rates resulting from substantially tion of oxygen via these processes played a 10.1126/science.1225461 ECOLOGY Downloaded from Ecological models have limited predictive The Art of Ecological Modeling power, but can provide insights into what makes an ecosystem vulnerable to disturbance. Ian L. Boyd cientists appear to have a very lim- issue—help to specify where the limits of evidence of the relative importance of dif- ited capacity to predict the behavior of prediction may lie. ferent species ( 7). Snonlinear dynamic systems, including Mougi and Kondoh used network models Forty years ago, May’s ( 8) suggestion that ecological systems ( 1). As understanding of to show that the specifi c mixture of agonistic complexity does not necessarily beget sta- these systems improves, uncertainty around and mutualistic interactions between species bility in ecological networks began a debate system behavior tends to increase rather than is likely to contribute to stabilizing popula- about a possible trade-off between complexity decline ( 2). The reason for this paradoxi- tions in ecological networks. Real ecosys- and stability in ecosystems. The recent results cal observation is that the uncertainty in the tems contain a rich mixture of these types of ( 3– 6) indicate that complexity does beget sta- models used for prediction has been under- interactions, meaning that community net- bility—but only in the presence of certain spe- estimated. Many of today’s ecological mod- work structure may affect the dynamics of cies and types of interaction between them. els have excessively optimistic ambitions to individual populations. Moreover, stability If only a small subset of all possible eco- predict ecosystem and population dynamics. and controllability of networks also depend logical networks with certain species abun- Some recent studies ( 3– 5)—including that by on how tightly the key nodes in these net- dances and types of interactions is stable, Mougi and Kondoh ( 6) on page 349 of this works are coupled ( 4, 5). These results sug- this result may explain the characteristic fre- gest that some species within ecosystems quency distribution of species abundance Scottish Oceans Institute, University of St Andrews, St are especially important to ecological sta- observed in ecological communities ( 9). It Andrews KY16 8LB, UK. E-mail: [email protected] bility, which is consistent with the empirical may also help to predict the consequences for 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org 306 Published by AAAS

The Art of Ecological Modeling Ian L. Boyd Science 337 , 306 (2012); DOI: 10.1126/science.1225049 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/306.full.html A list of selected additional articles on the Science Web sites related to this article can be found at: www.sciencemag.org http://www.sciencemag.org/content/337/6092/306.full.html#related This article cites 20 articles , 8 of which can be accessed free: http://www.sciencemag.org/content/337/6092/306.full.html#ref-list-1 Downloaded from This article appears in the following subject collections: Ecology http://www.sciencemag.org/cgi/collection/ecology Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

PERSPECTIVES and out of the ocean has fl uctuated much reduced marine sulfate levels. The authors larger role in regulating Phanerozoic atmo- 34 more than previously recognized. attribute a 5‰ positive δ S sulfate shift about spheric oxygen levels than previously recog- To better constrain sulfur fl uxes to and 50 million years ago to an abrupt increase nized—perhaps by as much as 50%. How- from the ocean, Wortmann and Paytan use in marine sulfate concentrations as a result ever, both studies recognize the importance a model that tracks the mass and sulfur iso- of large-scale dissolution of freshly exposed of episodic evaporite burial on the sulfur tope composition of fl uxes into and out of the evaporites; they argue that the higher sulfate cycle, and Wortmann and Paytan show that ocean. Most sulfur in the ocean is brought concentrations led to more pyrite burial. large-scale deposition and dissolution of there by rivers as a result of oxidative weath- Halevy et al. likewise study past sulfur these evaporites over relatively short geo- ering of sulfi des and dissolution of calcium fl uxes to and from the ocean, but over lon- logic time scales can have an enormous sulfates. Sulfur is removed from the ocean ger time scales (the past 500 million years). impact on marine sulfate concentrations, through microbial sulfate reduction fol- They quantify sulfate evaporite burial rates pyrite burial rates, and the carbon cycle. lowed by pyrite burial, and by calcium sul- through time with the help of a database The studies by Wortmann and Paytan fate (evaporite) deposition. Evaporite depo- that contains surface and subsurface infor- and Halevy et al. highlight the dynamic sition does not have a large effect on the sul- mation on sedimentary strata from loca- nature of the Phanerozoic sulfur cycle and fur isotopic composition, but microbial sul- tions across North America and the Carib- its importance in regulating the chemistry of fate reduction does, because microbes pre- bean. They then scale these rates to obtain a the ocean-atmosphere system. Future work 32 fer to use S over S during the production global estimate. should focus on better constraining marine 34 of sulfi de ( 4). As a result, the sulfur isotope The results indicate that sulfate burial sulfate levels and global sulfate burial rates 34 composition (δ S sulfate ) of river inputs today rates were higher than previously estimated, through time, the relationship between is ~6 per mil (‰) ( 5), but that of seawater but also greatly variable. When Halevy et al. marine sulfate concentrations and pyrite sulfate is ~20‰. integrated these improved evaporite burial burial rates, and the role of sulfur in regulat- on July 19, 2012 34 Variations in δ S sulfate are traditionally fl uxes with seawater sulfate concentration ing ocean nutrient cycling, the carbon cycle, interpreted to reflect changes in the total estimates ( 3) and sulfur isotope constraints and ultimately climate. amount of sulfur buried in ocean sediments ( 6, 7), their calculations implied that Pha- as pyrite—a parameter referred to as f pyr . nerozoic f pyr values were ~100% higher on References Wortmann and Paytan conclude that a average than previously recognized. These 1. U. G. Wortmann, A. Paytan, Science 337, 334 (2012). 2. I. Halevy et al., Science 337, 331 (2012). 34 5‰ negative δ S sulfate shift in ~120-mil- surprisingly high and constant pyrite burial 3. T. K. Lowenstein et al., Geology 31, 857 (2003). lion-year-old rocks was caused by massive outputs must have been balanced by equally 4. D. E. Canfi eld, Geochim. Cosmochim. Acta 65, 1117 seawater sulfate removal that accompanied high and constant inputs of sulfate to the (2001). www.sciencemag.org 5. M. A. Arthur, Encyclopedia of Volcanoes (Academic Press, large-scale evaporite deposition during the ocean via sulfi de oxidation (weathering). San Diego, CA, 2000). opening of the South Atlantic Ocean. Evap- The relatively high and constant rates of 6. A. Paytan et al., Science 304, 1663 (2004). orite formation does not cause sulfur iso- pyrite weathering and burial over long geo- 7. A. Kampschulte, H. Strauss, Chem. Geol. 204, 255 tope fractionation; in the model, the nega- logic time scales, identifi ed by Halevy et al., (2004). 34 tive δ S sulfate shift is driven by lower pyrite suggest that the consumption and produc- burial rates resulting from substantially tion of oxygen via these processes played a 10.1126/science.1225461 ECOLOGY Downloaded from Ecological models have limited predictive The Art of Ecological Modeling power, but can provide insights into what makes an ecosystem vulnerable to disturbance. Ian L. Boyd cientists appear to have a very lim- issue—help to specify where the limits of evidence of the relative importance of dif- ited capacity to predict the behavior of prediction may lie. ferent species ( 7). Snonlinear dynamic systems, including Mougi and Kondoh used network models Forty years ago, May’s ( 8) suggestion that ecological systems ( 1). As understanding of to show that the specifi c mixture of agonistic complexity does not necessarily beget sta- these systems improves, uncertainty around and mutualistic interactions between species bility in ecological networks began a debate system behavior tends to increase rather than is likely to contribute to stabilizing popula- about a possible trade-off between complexity decline ( 2). The reason for this paradoxi- tions in ecological networks. Real ecosys- and stability in ecosystems. The recent results cal observation is that the uncertainty in the tems contain a rich mixture of these types of ( 3– 6) indicate that complexity does beget sta- models used for prediction has been under- interactions, meaning that community net- bility—but only in the presence of certain spe- estimated. Many of today’s ecological mod- work structure may affect the dynamics of cies and types of interaction between them. els have excessively optimistic ambitions to individual populations. Moreover, stability If only a small subset of all possible eco- predict ecosystem and population dynamics. and controllability of networks also depend logical networks with certain species abun- Some recent studies ( 3– 5)—including that by on how tightly the key nodes in these net- dances and types of interactions is stable, Mougi and Kondoh ( 6) on page 349 of this works are coupled ( 4, 5). These results sug- this result may explain the characteristic fre- gest that some species within ecosystems quency distribution of species abundance Scottish Oceans Institute, University of St Andrews, St are especially important to ecological sta- observed in ecological communities ( 9). It Andrews KY16 8LB, UK. E-mail: [email protected] bility, which is consistent with the empirical may also help to predict the consequences for 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org 306 Published by AAAS

PERSPECTIVES ecosystem stability of weak- able to identify and measure the dynamics ening key nodes in these net- of infl uential network components, probably works through biodiversity in the form of key species ( 7). These species loss as a result of selective need not be the most common or those that exploitation or destruction of are exploited ( 7, 20). species ( 10), or because of the In future, the rational management of effects of climate change on exploited species will probably depend the relative strength of interac- increasingly on the status of these key species tions between species ( 11). in ecological networks. The legal designation However, if we are to have of protected species, together with an obli- any hope of managing eco- gation on authorities to monitor their status, logical systems, we require a presents an opportunity to establish ecolog- mechanism to identify which ically meaningful system indicators that are disturbances are likely to lead also relevant to public policy. to instability. Using cultured yeast populations, Dai et al. References ( 3) recently reported the pres- 1. S. A. Levin, Ecosystems 1, 431 (1998). 2. M. Maslin, P. Austin, Nature 486, 183 (2012). ence of telltale signs of tip- 3. L. Dai et al., Science 336, 1175 (2012). ping points, or bifurcations Ecosystem dynamics. Models show that the mixture of agonistic and 4. S. Allesina, S. Tang, Nature 483, 205 (2012). in the language of chaos the- mutualistic interactions in ecological networks is important for the sta- 5. Y. Y. Liu et al., Nature 473, 167 (2011). 6. A. Mougi, M. Kondoh, Science 337, 349 (2012). ory, when the system state can bility of communities. These complex interactions are typically not cap- 7. D. U. Hooper et al., Ecol. Monogr. 75, 3 (2005). on July 19, 2012 change rapidly and unpredict- tured by models used to manage ecological resources. 8. R. M. May, Nature 238, 413 (1972). ably. Predicting bifurcations in 9. F. W. Preston, Ecology 43, 185 (1962). 10. A. O. Shelton, M. Mangel, Proc. Natl. Acad. Sci. U.S.A. chaotic systems has many of the characteris- Predicting the dynamics of real ecosys- 108, 7075 (2011). tics of playing roulette, and some have rightly tems—or even of components of these eco- 11. P. L. Zarnetske et al., Science 336, 1516 (2012). cautioned against optimism about their pre- systems—will remain beyond the reach of 12. A. Hastings, D. B. Wysham, Ecol. Lett. 13, 464 (2010). dictability ( 12). However, it appears that per- even the best ecosystem models for the fore- 13. S. R. Carpenter et al., Science 332, 1079 (2011). 14. A. J. Veraart et al., Nature 481, 357 (2012). turbed populations take longer to return to seeable future (see the figure). However, 15. M. Scheffer et al., Nature 461, 53 (2009). their original state in advance of a bifurca- the emerging body of evidence ( 3– 6), sug- 16. R. Hilborn, Nat. Resour. Model. 25, 122 (2012). www.sciencemag.org tion. This supports the observation of chang- gests that ecological network models can 17. V. Christensen et al., Ecol. Modell. 220, 1984 (2009). 18. A. D. M. Smith et al., Science 333, 1147 (2011). ing dynamics within food webs up to a year in be used to describe ecosystem characteris- 19. G. Sugihara et al., Proc. Natl. Acad. Sci. U.S.A. 108, E1224 advance of an experimentally induced regime tics and general behavior, including provid- (2011). shift in a real ecosystem ( 13). In practice, ing indicators of when ecosystems are being 20. P. M. Cury et al., Science 334, 1703 (2011). this effect can be detected statistically ( 14, stressed. Successful management of eco- 15), suggesting the possibility of establishing logical processes will come down to being 10.1126/science.1225049 indicators of future system status. These results need to be translated into Downloaded from operational uses. Standard models used to APPLIED PHYSICS manage biological resources rarely have the associated with bifurcations. The main mod- Spin Twists in a Transistor capacity to model the dynamics typically els are general population models ( 16) and data-driven, heuristic, ecosystem models Igor Žutic´ and Jeongsu Lee ( 17, 18), which are rarely validated and often overparameterized. They are grounded in the By using electron spin in a transistor, researchers have developed a new approach for transferring concept of linear (predictable) rather than and processing information. nonlinear (unpredictable) systems and thus fail to capture the underlying complexity of ransistors are the centerpiece of con- as explored in the fi eld of spintronics ( 1, 2). ecological processes. ventional electronics, with two key Whereas harnessing spin for robust infor- These model limitations must be kept in Tfeatures of switching and amplifi ca- mation storage in computer hard drives and mind when reporting results, especially for tion. However, transistors relying on elec- magnetic random access memories has been use in decision-making. It may be tempting to tron charge are oblivious to another property very successful, delivering a spin transistor ignore model uncertainty in the name of prag- of electrons: their spin. In a simple picture, has been challenging ( 1). On page 324 of this CREDIT: RICH CAREY/SHUTTERSTOCK plication of ecological models may underlie “up” and “down,” spin lends itself to encod- ing that the pioneering work on another spin issue, Betthausen et al. ( 3) describe a newly these spins are compass needles, aligned by matism, but this will not lead to better deci- sions. Indeed, this type of systemic misap- discovered spin-transistor action. Consider- a magnetic fi eld. With different orientations, some cases of natural resource depletion. For transistor ( 4) waited two decades for experi- ing binary information as ones and zeroes, example, fi sheries models typically assume mental realization ( 5), it is a remarkable feat that the work of Betthausen et al. presents linear dynamics and rarely consider the impli- Department of Physics, University at Buffalo, The State cations of wider community interactions for the experiment and theory for their spin University at New York, Buffalo, NY 14260, USA. E-mail: transistor. the validity of their predictions ( 10, 19). [email protected]; [email protected] www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 307 Published by AAAS

Spin Twists in a Transistor Igor Zutic and Jeongsu Lee Science 337 , 307 (2012); DOI: 10.1126/science.1225219 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/307.full.html A list of selected additional articles on the Science Web sites related to this article can be found at: www.sciencemag.org http://www.sciencemag.org/content/337/6092/307.full.html#related This article cites 15 articles , 4 of which can be accessed free: http://www.sciencemag.org/content/337/6092/307.full.html#ref-list-1 Downloaded from This article appears in the following subject collections: Physics, Applied http://www.sciencemag.org/cgi/collection/app_physics Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

PERSPECTIVES ecosystem stability of weak- able to identify and measure the dynamics ening key nodes in these net- of infl uential network components, probably works through biodiversity in the form of key species ( 7). These species loss as a result of selective need not be the most common or those that exploitation or destruction of are exploited ( 7, 20). species ( 10), or because of the In future, the rational management of effects of climate change on exploited species will probably depend the relative strength of interac- increasingly on the status of these key species tions between species ( 11). in ecological networks. The legal designation However, if we are to have of protected species, together with an obli- any hope of managing eco- gation on authorities to monitor their status, logical systems, we require a presents an opportunity to establish ecolog- mechanism to identify which ically meaningful system indicators that are disturbances are likely to lead also relevant to public policy. to instability. Using cultured yeast populations, Dai et al. References ( 3) recently reported the pres- 1. S. A. Levin, Ecosystems 1, 431 (1998). 2. M. Maslin, P. Austin, Nature 486, 183 (2012). ence of telltale signs of tip- 3. L. Dai et al., Science 336, 1175 (2012). ping points, or bifurcations Ecosystem dynamics. Models show that the mixture of agonistic and 4. S. Allesina, S. Tang, Nature 483, 205 (2012). in the language of chaos the- mutualistic interactions in ecological networks is important for the sta- 5. Y. Y. Liu et al., Nature 473, 167 (2011). 6. A. Mougi, M. Kondoh, Science 337, 349 (2012). ory, when the system state can bility of communities. These complex interactions are typically not cap- 7. D. U. Hooper et al., Ecol. Monogr. 75, 3 (2005). on July 19, 2012 change rapidly and unpredict- tured by models used to manage ecological resources. 8. R. M. May, Nature 238, 413 (1972). ably. Predicting bifurcations in 9. F. W. Preston, Ecology 43, 185 (1962). 10. A. O. Shelton, M. Mangel, Proc. Natl. Acad. Sci. U.S.A. chaotic systems has many of the characteris- Predicting the dynamics of real ecosys- 108, 7075 (2011). tics of playing roulette, and some have rightly tems—or even of components of these eco- 11. P. L. Zarnetske et al., Science 336, 1516 (2012). cautioned against optimism about their pre- systems—will remain beyond the reach of 12. A. Hastings, D. B. Wysham, Ecol. Lett. 13, 464 (2010). dictability ( 12). However, it appears that per- even the best ecosystem models for the fore- 13. S. R. Carpenter et al., Science 332, 1079 (2011). 14. A. J. Veraart et al., Nature 481, 357 (2012). turbed populations take longer to return to seeable future (see the figure). However, 15. M. Scheffer et al., Nature 461, 53 (2009). their original state in advance of a bifurca- the emerging body of evidence ( 3– 6), sug- 16. R. Hilborn, Nat. Resour. Model. 25, 122 (2012). www.sciencemag.org tion. This supports the observation of chang- gests that ecological network models can 17. V. Christensen et al., Ecol. Modell. 220, 1984 (2009). 18. A. D. M. Smith et al., Science 333, 1147 (2011). ing dynamics within food webs up to a year in be used to describe ecosystem characteris- 19. G. Sugihara et al., Proc. Natl. Acad. Sci. U.S.A. 108, E1224 advance of an experimentally induced regime tics and general behavior, including provid- (2011). shift in a real ecosystem ( 13). In practice, ing indicators of when ecosystems are being 20. P. M. Cury et al., Science 334, 1703 (2011). this effect can be detected statistically ( 14, stressed. Successful management of eco- 15), suggesting the possibility of establishing logical processes will come down to being 10.1126/science.1225049 indicators of future system status. These results need to be translated into Downloaded from operational uses. Standard models used to APPLIED PHYSICS manage biological resources rarely have the associated with bifurcations. The main mod- Spin Twists in a Transistor capacity to model the dynamics typically els are general population models ( 16) and data-driven, heuristic, ecosystem models Igor Žutic´ and Jeongsu Lee ( 17, 18), which are rarely validated and often overparameterized. They are grounded in the By using electron spin in a transistor, researchers have developed a new approach for transferring concept of linear (predictable) rather than and processing information. nonlinear (unpredictable) systems and thus fail to capture the underlying complexity of ransistors are the centerpiece of con- as explored in the fi eld of spintronics ( 1, 2). ecological processes. ventional electronics, with two key Whereas harnessing spin for robust infor- These model limitations must be kept in Tfeatures of switching and amplifi ca- mation storage in computer hard drives and mind when reporting results, especially for tion. However, transistors relying on elec- magnetic random access memories has been use in decision-making. It may be tempting to tron charge are oblivious to another property very successful, delivering a spin transistor ignore model uncertainty in the name of prag- of electrons: their spin. In a simple picture, has been challenging ( 1). On page 324 of this CREDIT: RICH CAREY/SHUTTERSTOCK plication of ecological models may underlie “up” and “down,” spin lends itself to encod- ing that the pioneering work on another spin issue, Betthausen et al. ( 3) describe a newly these spins are compass needles, aligned by matism, but this will not lead to better deci- sions. Indeed, this type of systemic misap- discovered spin-transistor action. Consider- a magnetic fi eld. With different orientations, some cases of natural resource depletion. For transistor ( 4) waited two decades for experi- ing binary information as ones and zeroes, example, fi sheries models typically assume mental realization ( 5), it is a remarkable feat that the work of Betthausen et al. presents linear dynamics and rarely consider the impli- Department of Physics, University at Buffalo, The State cations of wider community interactions for the experiment and theory for their spin University at New York, Buffalo, NY 14260, USA. E-mail: transistor. the validity of their predictions ( 10, 19). [email protected]; [email protected] www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 307 Published by AAAS

PERSPECTIVES A BC Source Gate Drain Source Gate Drain Adiabatic process Channel Channel Ferromagnet Ferromagnet Channel with modulation of magnetic field Diabatic process ? Transistor switching. (A) In conventional transistors, the electron fl ow from the drain. (C) In adiabatic spin transistors, the “on” state is realized by the grad- source to drain is controlled by the gate voltage, V G . (B) In a spin transistor, V G ual (adiabatic) change of the magnetic fi eld (white arrows). The spins have time alters the spin-orbit coupling, which rotates the electron spin direction in the to respond to the direction dictated by the fi eld, and electrons fl ow smoothly channel. A ferromagnetic source initializes the spin alignment, and a ferromag- across the channel. For the “off” state, there is an abrupt (diabatic) change of the netic drain determines whether the fl ow is permitted (restricted), depending on magnetic fi eld. Electrons unable to align their spin along the fi eld are scattered the spin alignment (misalignment) between the electron exiting the channel and back and their fl ow is impeded. on July 19, 2012 Conventional transistors ( 6) are three- and spin-down orientation), and a high elec- be electrically controlled. This possibility terminal devices relying on classical charge tron mobility, while the strong magnetic fi eld would allow revisiting amplifi cation, as it is transfer (see the fi gure, panel A). Describing keeps their spins in check. Careful design of known in conventional transistors ( 6), where spin transistors (panels B and C), however, is the magnetic fi eld profi le combines a heli- the change in the gate voltage controls the more elusive because it requires quirky prob- cal pattern along the channel (see the fi gure, source-drain current (see the fi gure, panel abilistic laws of quantum mechanics. Simi- panel C), which results in a twisting of the A). Demonstrating the amplifi cation in adia- lar to tossing a coin, we cannot say whether spins, and a constant, but tunable, external batic spin transistors could also come with the tails or heads will be next, but we can fi eld that selects the transition mechanism. additional benefi ts of a versatile spin con- www.sciencemag.org only calculate the likelihood for either out- The new spin transistor gives the fi eld trol and a robust performance with respect to come to happen. In a spin transistor, the of spintronics a major boost, opening unex- slight variations, inevitable in fabrication of relevant outcomes are different spin orien- plored directions of not just storing but different devices. tations. By using spin twists, Betthausen et also transferring and processing spin-based At the heart of spintronics lies the quest al. have implemented transistor switching information ( 10– 12). At present, however, to elucidate physical phenomena and possi- action between the low- and high-resistance the spin transistor is unlikely to compete ble applications emerging from microscopic states (“on” and “off,” respectively). For the with its conventional counterpart. The adia- spins. Although spin twists are limited to compass-needle analogy, we can have two batic spin transistor is only realized at tem- small distances, they may offer a giant leap Downloaded from extreme cases for spin twists: For a magnet peratures below 1 K, and the reliance on an for effectively taming the spin transistor. slowly rotating above a compass, the needle external magnetic fi eld is not very practi- will follow the magnet, as expected for an cal. Furthermore, the amplifi cation in the References adiabatic transition. However, if we rotate transistor action has yet to be demonstrated. 1. I. Žutić, J. Fabian, S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004). the magnet too quickly, the compass needle However, none of these current limitations 2. J. Fabian et al., Acta Phys. Slov. 57, 565 (2007). will not have time to react and its orienta- are fundamental obstacles. 3. C. Betthausen et al., Science 337, 324 (2012). tion will stay put, a sign that the transition A higher operation temperature could 4. S. Datta, B. Das, Appl. Phys. Lett. 56, 665 (1990). 5. H. C. Koo et al., Science 325, 1515 (2009). is diabatic. be realized by adding more Mn impurities 6. P. Horowitz, W. Hill, The Art of Electronics (Cambridge, Betthausen et al. have shown how adia- to increase the spin splitting. Alternatively, New York, 1989). batic processes effectively shelter spin infor- one could consider other two-dimensional 7. J. Wunderlich et al., Science 330, 1801 (2010). mation, propagating undisturbed and with- magnetic channels—for example, mag- 8. N. Rangaraju, J. A. Peters, B. W. Wessels, Phys. Rev. Lett. 105, 117202 (2010). out scattering at distances 10 to 100 times netic semiconductors ( 13) or oxides ( 14)— 9. J. K. Furdyna, J. Appl. Phys. 64, R29 (1988). greater than the characteristic distance for showing a permanent magnetism with a 10. S. Sugahara, IEE Proc., Circ. Devices Syst. 152, 355 charge scattering. Although such a distance large spin splitting. Some caution is needed (2005). of spin-information transfer is still small, it because additional magnetic impurities typ- 11. H. Dery, Y. Song, P. Li, I. Žutić, Appl. Phys. Lett. 99, 082502 (2011). dwarfs what is currently available in other ically reduce the mobility and thus shorten 12. H. Dery et al., IEEE Trans. Electron. Dev. 59, 259 type of spin transistors ( 5, 7, 8) and even the the range of spin-information transfer. If (2012). characteristic dimensions of ever-shrinking an important functionality can be identifi ed 13. T. Jungwirth, J. Sinova, J. Mašek, J. Kučera, A. H. conventional transistors. These results were for the spin transistor, the low temperature MacDonald, Rev. Mod. Phys. 78, 809 (2006). 14. L. Li, C. Richter, J. Mannhart, R. C. Ashoori, Nat. Phys. 7, made possible by tailoring materials proper- alone would not be the deal breaker. With 762 (2011). ties and adding magnetic (Mn) impurities in either voltage-tunable permanent magnets 15. H. Ohno et al., Nature 408, 944 (2000). a two-dimensional channel material (Cd,Mn) ( 1, 15, 16) or spin-orbit coupling ( 1, 2, 5) 16. Y. Yamada et al., Science 332, 1065 (2011). Te ( 9). This choice enables both a large spin the need for an external magnetic fi eld could splitting (energy difference between spin-up be removed because the spin splitting would 10.1126/science.1225219 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org 308 Published by AAAS

Aaron Shatkin (1934 − 2012) Nahum Sonenberg and Witold Filipowicz Science 337 , 309 (2012); DOI: 10.1126/science.1226820 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online www.sciencemag.org version of this article at: http://www.sciencemag.org/content/337/6092/309.full.html This article appears in the following subject collections: Scientific Community http://www.sciencemag.org/cgi/collection/sci_commun Downloaded from Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

PERSPECTIVES RETROSPECTIVE Aaron Shatkin (1934–2012) A virologist discovered the eukaryotic mRNA 5´-terminal cap structure and revealed its role in gene expression. 2 1 Nahum Sonenberg and Witold Filipowicz aron Shatkin died of cancer at his plates for synthesis of progeny segments. The and collaborators remained in close contact home in Scotch Plains, New Jer- work with reovirus led Aaron and his associ- with Aaron, very often seeking his support A sey, on 4 June, at the age of 77. ate Yasuhiro Furuichi to discover the 5´-ter- and advice. Over 40 of us, from all over the When asked last year by the Bowdoin Col- minal m GpppN cap structure in 1975. world, gathered in September 2011 for the 7 lege (where he majored in chemistry) mag- The discovery that viral mRNAs have “Shatkin Reunion”—sadly the last one—to azine, what was the most rewarding part of their termini modifi ed with m GpppN caps show our appreciation for his mentorship and 7 his job, he responded: “Seeing students and (made together with Ber- friendship. The day of the others I’ve mentored at CABM (Center for nard Moss working on reunion was full of recol- Advanced Biotechnology and Medicine at vaccinia virus) had far- lections, short stories, and Rutgers University and the University of reaching consequences. news from the lives of our Medicine and Dentistry of New Jersey) and Caps are present on all families, which Aaron fol- elsewhere become independent, accom- nuclear-transcribed cellu- lowed very closely. What plished researchers and leaders in their fi elds lar mRNAs and mRNAs moved us most strongly on in many different countries.” His gracious- of many viruses, and are that day was Aaron’s engag- on July 19, 2012 ness and passion for science were among his involved in almost every ing presentation, deliv- fi nest attributes. Both of us, like all his train- aspect of mRNA metabo- ered already in frail health, ees, owe our progress in science largely to his lism. In 1975, Aaron and describing his current extraordinary mentorship. James Darnell showed research and reminiscing Aaron did his graduate studies with that long heterogeneous about the happy days with Edward Tatum, Nobel laureate, at the Rock- nuclear RNAs too harbor his students and colleagues. efeller University, investigating the morphol- m GpppN caps in addition Aaron was a recipi- 7 ogy of the fi lamentous fungus Neurospora to being polyadenylated ent of many awards and www.sciencemag.org crassa. However, the work coauthored with at the 3´ end. How both honors, including the Edward Reich and Richard Franklin, dem- these terminal structures Steel award from the U.S. onstrating that the antibiotic actinomycin D could be preserved in much National Academy of Sci- blocks cellular messenger RNA (mRNA) shorter mature mRNAs became evident only ences (1977) and the Association of Ameri- synthesis, had the most impact on his future upon the discovery of mRNA splicing. can Medical Colleges Award (2003). He was career and started his lifelong adventure with Aaron’s laboratory further led the field elected to the U.S. National Academy of Sci- animal viruses. The availability of actinomy- in conducting a comprehensive analysis of ences in 1981 and the American Academy of cin D allowed Aaron and others to investigate the mechanism by which the m GpppN cap Arts and Sciences in 1997. In addition to Downloaded from 7 virus gene expression in the absence of host facilitates translation initiation. This entailed his illustrious career as a researcher, Aaron mRNA synthesis, which greatly facilitated the discovery of new translation factors and greatly contributed to other aspects of sci- the study of virus-encoded proteins. Fol- promulgation of the scanning mechanism by ence. He was a founding editor-in-chief of lowing the move to the laboratory of Nor- which the mammalian 40S ribosome selects the journal Molecular and Cellular Biology. man Salzman at the National Institutes of the AUG initiation codon. His laboratory also In 1986, he became the founding director of Health in 1963, Aaron studied several human documented the importance of the cap for CABM in Piscataway, New Jersey. Under viruses including vaccinia virus and poliovi- mRNA stability. More recently, Aaron pur- Aaron’s leadership, CABM has grown into a rus, providing early insights into viral RNA sued his favorite topic of the capping mech- world-class research institution. In recogni- PHOTO CREDIT: UNIVERSITY OF MEDICINE AND DENTISTRY OF NEW JERSEY and protein synthesis. From the mid-1960s, anism and enzymes involved. His last paper tion of Aaron’s extraordinary leadership, an he focused almost entirely on the life cycle of in 2011 described the x-ray structure of the endowed annual lectureship was established reovirus, studied also independently by Wolf- carboxyl-terminal capping domain of the in his name at CABM, with Harold Var- gang Joklik and Bernard Fields. After moving human capping enzyme. Although focused mus, director of the National Cancer Insti- to the Roche Institute of Molecular Biology on mRNA capping, Aaron was always open tute, giving the inaugural lecture in April in Nutley, New Jersey, in 1968, Aaron showed to other ideas. His laboratory discovered the 2012. Aaron arrived a few minutes late for that the segmented genome of the reovirus pathway leading to RNA ligation and transfer the lecture, straight from Boston where he consists of 10 separate double-stranded RNA RNA (tRNA) splicing in animal cells in the was undergoing experimental therapy. He “chromosomes,” which are copied by a viral early 1980s. passed away a month later. He worked until a RNA polymerase, and that the resulting tran- Aaron was loved by his peers, students, few weeks before he died, a testament to his scripts function both as mRNAs and tem- and collaborators. He was a person of excep- devotion to science and his work. Aaron is tional integrity, and his warm character and survived by his son Greg. Joan, his beloved 1 Department of Biochemistry, McGill University, Goodman radiating optimism (and smile) were infec- wife for 51 years, passed away in 2009. Cancer Research Center, Montreal, Quebec H3A 1A3, Can- tious. He was an excellent mentor, always Till his last days, Aaron stopped daily for a ada. Friedrich Miescher Institute for Biomedical Research, 2 Maulbeerstrasse 66, 4058 Basel, Switzerland. E-mail: fostering openness, collaboration, and moment of silence in front of her grave. [email protected]; witold.fi [email protected] friendly dispute. Most of his former students 10.1126/science.1226820 www.sciencemag.org SCIENCE VOL 337 20 JULY 2012 309 Published by AAAS

The Exploration of Hot Nuclear Matter Barbara V. Jacak and Berndt Müller Science 337 , 310 (2012); DOI: 10.1126/science.1215901 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/310.full.html This article cites 52 articles , 3 of which can be accessed free: http://www.sciencemag.org/content/337/6092/310.full.html#ref-list-1 www.sciencemag.org This article appears in the following subject collections: Physics http://www.sciencemag.org/cgi/collection/physics Downloaded from Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

REVIEW can fill this gap by colliding heavy nuclei at high energies and generating the enormous temper- The Exploration of Hot Nuclear Matter atures required to produce QGP in the labora- tory, if only for a brief moment. 1 2 Barbara V. Jacak and Berndt Müller * Discoveries of the Past Decade The Relativistic Heavy Ion Collider (RHIC) at When nuclear matter is heated beyond 2 trillion degrees, it becomes a strongly coupled Brookhaven National Laboratory has explored plasma of quarks and gluons. Experiments using highly energetic collisions between heavy nuclei the QGP since 2000. RHIC collides two beams of have revealed that this new state of matter is a nearly ideal, highly opaque liquid. A description heavy ions, each with an energy up to 100 GeV based on string theory and black holes in five dimensions has made the quark-gluon plasma per nucleon. Proton-proton (p+p) and deuteron- an archetypical strongly coupled quantum system. Open questions about the structure and gold (d+Au) collisions provide control measure- theory of the quark-gluon plasma are under active investigation. Many of the insights are also ments without QGP formation. At top energy, relevant to ultracold fermionic atoms and strongly correlated condensed matter. the initial temperature reached in collisions be- tween two gold nuclei is inferred to uclear matter in today’s uni- lie between 300 and 600 MeV (5), verse hides inside atomic well above the QCD phase-transition Nnuclei and neutron stars. The temperature of ~150 MeV (6). RHIC nucleons (neutrons and protons) is a flexible, dedicated facility collid- are the building blocks of such mat- ing a wide range of nuclei at various ter and consist, in turn, of quarks. energies. This allows exploration of Quarks are bound together by the the phase diagram of QCD matter to on July 19, 2012 strong interaction, which is mediated experimentally pinpoint the conditions by the exchange of gluons. Unlike for the phase transition into QGP. the uncharged photons, which me- Today, two large experiments built diate electromagnetic interactions and maintained by international col- but do not interact with one anoth- laborations of scientists, PHENIX er, gluons have color, which is the and STAR (Fig. 2), continue to oper- strong interaction’sanalog of charge. ate whereas two smaller experiments, Colored gluons interact among them- BRAHMS and PHOBOS, have com- www.sciencemag.org selves, as well as with the quarks, pleted data taking. Each experiment making the theory of the strong in- was optimized for a different set of teraction, known as quantum chro- experimental observables, but com- modynamics (QCD), rich in structure mon capabilities allow crucial cross and at the same time extremely dif- checks. Together, PHENIX and STAR ficult to solve. use two kinds of plasma probes (7–11). Remarkably, the strong interac- Internal probes are particles emitted tion weakens at short distances—a from the plasma itself. “External” Downloaded from property known as “asymptotic free- probes are not external in the usual dom” (1, 2). Conversely, it is exceed- sense; they are energetic particles ingly strong at distances similar to generated in the first stage of the thesizeofa nucleon(10 −15 m), con- collision, which traverse the plasma fining quarks inside nucleons and and interact with it on their way to other quark-containing particles, Fig. 1. Phase diagram of QCD matter in the temperature–baryon density plane. the detectors. Baryons are hadrons containing three valence quarks; the most common are known as hadrons. Asymptotic free- protons and neutrons, shown at the lower left. Colored spheres indicate Most of the observed particles dom suggests that nucleons can be are hadrons. Their spectra are well individual quarks, which are not bound together in the quark-gluon plasma. “boiled” into a plasma of their con- RHIC (blue ovals) and LHC (green oval) explore matter with almost equal described by a thermal distribution stituent quarks and gluons when numbers of quarks and antiquarks. At lower beam energies, RHIC produces blue-shifted by radial expansion of the strong interaction among them matter with a surplus of quarks, corresponding to high net baryon density. the plasma. Particle correlations re- is weakened by increasing the den- There may be a critical point (yellow circle) in the phase diagram, at the end flect an anisotropic collective flow, sity or temperature of the matter. of a line indicating a first-order phase transition. [Credit: Brookhaven National known as “elliptic flow.” They exhibit Today, quarks are confined in nu- Laboratory] acos(2f) modulation in their azi- clei and neutron stars, which are muthal angular distribution with cold objects, but the early universe respect to the direction of the impact- was extremely hot (3). Its temperature exceeded conditions are sufficient for the quark-gluon plas- parameter vector between the two colliding ions 150 MeV (about 2 × 10 12 K) until about 10 ms ma (QGP), to exist (Fig. 1). (12). The amplitude of elliptic flow grows with after the Big Bang, and QCD predicts that such Understanding the evolution of our universe increasing impact parameter because the overlap thus requires knowledge of the structure and region of the incoming nuclei becomes more asym- dynamics of the QGP. Although numerical ab metric (Fig. 3, left). The dynamics within the plasma 1 Department of Physics and Astronomy, State University of initio simulations of the thermodynamic prop- as it expands translate the spatial asymmetry of 2 New York, Stony Brook, NY 11794, USA. Department of Phys- ics and Center for Theoretical and Mathematical Sciences, erties of hot QCD matter in equilibrium have the initial state into a final-state anisotropy in mo- Duke University, Durham, NC 27708, USA. made much progress over the past 30 years (4), mentum space. Higher Fourier components of the *To whom correspondence should be addressed. E-mail: the dynamical properties of the QGP remain out angular distribution are also observed in the cor- [email protected] of reach. Experimental study of hot QCD matter relation data; these arise primarily from fluctua- 310 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REVIEW tions in the initial positions of the nucleons within The opacity of the QGP is measured with ex- fits indicate nearly the same h/s ratio as at RHIC the nucleus (Fig. 3, right). ternal probes. During initial interpenetration (27); the hotter QGP produced at LHC is also The behavior of gases or liquids is often sim- of the two nuclei, quark and gluon constituents strongly coupled. Jet quenching measurements at ulated using hydrodynamics. Indeed, hydrody- (partons) can scatter with a large momentum LHC extend the kinematic range by a factor of 5. namics successfully reproduces the magnitude transfer, deflecting the struck quarks or gluons They are consistent with a linear, or slightly slower, and impact parameter dependence of elliptic by a large angle. These transit the plasma, losing gowth of the opacity with matter density. The flow (13, 14). Surprisingly, the most faithful match energy to it. As partons cannot exist in isolation, LHC’s higher energy produces higher-energy jets, to the data requires a nearly vanishing ratio of the they ultimately radiate multiple gluons, and the which simplifies reconstruction of complete jet ob- shear viscosity (the resistance to flow or the in- resulting parton cluster forms a spray of hadrons servables. Clusters of energy corresponding to back- ability of matter to transport momentum) to the known as a “jet.” In the absence of QGP, this pro- to-back jets are clearly visible in Fig. 4B (28). The entropy density, implying that the QGP is an al- cess is calculable and can be measured in p+p data reveal that even very energetic jets lose a siz- most ideal or “perfect” liquid. Including density collisions. RHIC experiments with Au+Au showed able fraction of their energy to the medium, where it fluctuations in the initial conditions of the hydro- that energetic hadrons in the jet are suppressed appears to thermalize rapidly. Analysis of jet shapes dynamical simulations also reproduces the high- relative to the production rate in p+p collisions and particle content help constrain the mechanism er harmonics with the same low shear viscosity (22, 23). Back-to-back jets of moderate energy of the parton-QGP interaction. The yields of D and (15–17). disappear entirely (24). Photons do not experi- B mesons, which contain heavy quarks, are also Quantum mechanics imposes a lower limit of ence the strong interaction (Fig. 4A) and are ob- much larger at LHC. Furthermore, the Z boson the shear viscosity h for a given particle density served to exit the plasma unscathed (25). The becomes available as a new electroweak probe by virtue of the uncertainty relation. For relativ- success of hydrodynamics with vanishingly small of the QGP. First, statistically limited, results for istic fluids like the QGP, which do not conserve h/s, together with the observed high opacity, show these “external” probes can be reproduced rea- particle number, the appropriate measure of den- that QGP cannot be the weakly coupled gas naïve- sonably well by extrapolation from RHIC. sity is the entropy density s. Constraining hydro- ly expected from asymptotic freedom. Kinetic the- dynamics calculations with the full suite of flow ory associates a small shear viscosity with a short Theoretical Tools on July 19, 2012 data (18)shows that h=s ¼ð1−2Þℏ=4pk B ,close mean free path, implying high opacity. A short Key theoretical tools to describe QGP properties to the quantum limit ℏ=4pk B (19–21). Low shear mean free path also requires strong coupling, and predict experimental observables are lattice viscosity per particle indicates correlations or co- because the scattering cross section is propor- gauge theory and transport theory. Lattice gauge www.sciencemag.org ordination within the QGP. Gases have very weak- tional to the coupling strength. theory is an ab initio formalism that simulates ly correlated constituents, whereas the molecules First results from Pb+Pb collisions at nearly the partition function of QCD on a space-time in crystals move in a highly coordinated manner. 14 times higher energy at the Large Hadron Col- lattice. Advances in algorithms and computer Liquids fall in between the two, exhibiting the lider (LHC) confirm the physics picture derived hardware now permit simulations with physical lowest shear viscosities and flowing freely, as from RHIC data (26). The initial temperature at quark masses on lattices that are simultaneously does the QGP. LHC is ~30% higher. Hydrodynamical model large and fine enough to be safely extrapolated Downloaded from Fig. 2. (A) The STAR detector has a time-projection chamber (TPC), which is es- has two spectrometers to measure photons, electrons, and hadrons at angles near sentially a three-dimensional digital camera to record trajectories of particles 90°; one is visible at left (Spec.). There are also two muon spectrometers in the produced in each collision. Surrounding detectors identify hadrons and tag high- beam direction; these detect decays of hadrons containing charm and bottom momentum electrons. STAR has large acceptance and is thus well suited to study quarks. A sample event display is shown on the right side of each detector. [Credit: multiparticle correlations and collisions at lower energies. (B) The PHENIX detector Brookhaven National Laboratory] SCIENCE 311 www.sciencemag.org VOL 337 20 JULY 2012

REVIEW monuclear fusion. In these plasmas, the ratio r of potential to kinetic energy is large, implying strong coupling; at sufficiently large r, such plasmas can even crystallize. The shear viscosity exhibits a min- imum at a certain value of r,where the dominant mechanism of mo- mentum transport changes from ballistic quasiparticle motion to some form of collective transport. An advantage of QCD matter over other strongly coupled sys- tems is that the interaction is well defined. The QGP thus offers a chance to understand how a strong- ly coupled fluid emerges from a Fig. 3. Elliptic (left) and triangular (right) flow patterns arise from the locations of individual nucleons at the microscopic theory that is pre- instant when two nuclei interpenetrate. The nucleons of one nucleus are shown in yellow and the other in orange. cisely known. The strongly cou- Red indicates those nucleons in the overlap region, which actually collide. (Left) Adapted with permission from figure 1 in (56) [Copyrighted by the American Physical Society]. (Right) Adapted with permission from figure 3 in pled QGP is also the only known (57). [Copyrighted by the American Physical Society] relativistic liquid. Its structure is not dominated by repulsive inter- actions, so it challenges the tradi- on July 19, 2012 to the thermodynamic and continuum limits (29). theory (33, 34) provides compelling evidence for tional concept of a liquid. However, the high The equation of state of hot QCD matter and cor- this conjecture. The gravity dual description offers temperature of the QGP, combined with the fun- relation functions, such as the screening distance an explanation of how a strongly coupled plasma damental nature of the QCD interaction, per- of the color force, are now within reach. However, of gauge fields can reach thermal equilibrium so mits ab initio techniques to address equilibrium reliable calculations in lattice QCD are still limited rapidly and why hydrodynamics furnishes a re- properties of hot QCD matter without any mod- to static properties, severely restricting our abil- liable description even at strong coupling, when el assumptions or approximations. The rapidly ity to address transport properties of the QGP. kinetic theory fails. Unfortunately, the coupling expanding capabilities to perform definitive cal- Transport theory describes the conversion of of true QCD is not as strong as would be required culations of this kind enable newly rigorous www.sciencemag.org the gluon fields in the incoming nuclei into thermal for rigorous application of the AdS/CFT duality. comparisons between the theory of strongly cou- QCD matter, the explosive expansion of the QGP, At the moment no gravity dual for true QCD is pled systems and experiment. and finally its disassembly into hadrons. A standard known, and it is unknown whether one exists. scenario of distinct reaction stages has emerged Open Questions and Challenges (30). In the first stage, gluons are liberated and form Interconnections The surprising experimental results present an en- a dense system of nonlinearly coupled fields, Understanding strongly coupled or strongly cor- tirely new set of questions about the QGP. Roughly known as the glasma. The second stage, the rapid related systems is at the intellectual forefront of following the time development of heavy ion col- expansion of the hot QGP, is effectively described multiple subfields of physics. One example is lisions, one must now ask: What is the nature of Downloaded from 6 by relativistic hydrodynamics with small viscous ultracold fermionic atoms, such as Li, where ap- QCD matter at low temperature but high density, effects. After the matter cools below 150 MeV, its plication of a magnetic field excites a strong res- and how does it affect plasma formation? How final expansion and freeze-out can be described by onance. When confined in an atomic trap, these can the plasma thermalize so rapidly? Does it ex- kinetic theory for hadrons. Whereas the experi- atoms form a degenerate Fermi liquid, which can hibit novel symmetry properties along the way? mental data provide solid evidence for the validity be manipulated and studied in detail (35). At The QCD plasma is strongly coupled, but at what of the description of stages 2 and 3, experimental temperatures below ~0.1 mK, the atoms inter- scales? Does it contain quasiparticles, or does the exploration of the glasma phase is just beginning. acting via the resonance form a superfluid (36). strong coupling completely wipe out long-lived Physicists were astounded to find that an en- The shear viscosity h and the entropy density s collective excitations? What impact does the cou- tirely different approach, using dualities that relate for this system can be measured separately, show- pling have on color screening? Is there a charac- QCD at strong coupling with weakly coupled ing that h/s falls with decreasing temperature. At teristic screening length, and if so, what is it? What gravitational theories, can yield insights into the very low temperatures, it approaches about four is the mechanism for parton-plasma interactions, dynamical properties of quantum inviscid liquids. times the universal quantum limit (37), only twice and how does the plasma respond to energy de- The duality of string theory in anti-de Sitter (AdS) as large as the value deduced for the QGP. posited in it? space with conformal quantum field theory (CFT) Strongly correlated electron systems in con- Gravity dual calculations show that thermal- provides an exact description of some strongly cou- densed matter provide an example of strong cou- ization propagates at the speed of light and all pled systems. The formalism, known as AdS/CFT pling where the elementary interaction is not anisotropies disappear quickly in the strong cou- correspondence (31, 32), holographically maps strong, but its role is amplified by the large num- pling limit (41, 42). Plasma instabilities may play the intractable strongly coupled quantum field ber of interacting particles and their ability to a role. The microscopic structure of the strongly theory onto a solvable classical gravity theory in dynamically correlate their quantum wave func- coupled QGP is still poorly understood; in the five dimensions. Thermalization of the quantum tions. Surprisingly, holographic gravity duals have gravity dual picture, no quasiparticles exist ex- field appears as formation of a black hole in the also helped to provide simple descriptions of cept phonons. Lattice simulations confirm that gravity dual theory. Although the formalism is ex- such complex systems (38, 39). the quantum numbers associated with quarks— act only in the limit of an infinite number of colors Strongly coupled systems in conventional plas- baryon number, electric charge, and flavor—are and at strong coupling, it is believed that for many ma physics include warm, dense matter and dusty carried by elementary, quarklike constituents at quantities of interest the three colors of QCD can plasmas (40) residing in astrophysical environ- temperatures above the critical temperature for be considered as a large number. Lattice gauge ments, such as the rings of Saturn, and in ther- QGP formation, T c . However, they are unable to 312 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REVIEW address the dynamic response of the plasma and contain bottom quarks. The small size of these discriminate between radiative and collisional provide no information about the presence or mesons enables their existence in QGP because energy loss. Because little difference has been absence of propagating quasiparticles. Energy screening occurs only at larger distance scales. found in the suppression of light and heavy loss of heavy quarks is sensitive to the spectrum Quarkonia have different excited states with vary- (mostly charm) quarks (52), separating the charm of excitations of the QGP and may provide a ing binding energies; loosely bound, larger states quark from the even heavier bottom quark is a handle on quasiparticles and their properties. are easier for the QGP to screen. Indeed, suppres- key experimental goal. Photons and leptons should preserve imprints sion of charm quarkonia (charmonia) in QGP, of the early stages of the collision. compared to p+p collisions, has already been ob- The Look Ahead at RHIC and LHC Shear viscosity and speed of sound are two served (47–51). It should be possible to infer the Recent RHIC upgrades have increased both the important indicators of the microscopic structure color screening scale with spectroscopy of differ- luminosity and the range of particle species avail- of any material. The viscosity probes how the ent quarkonium states as a function of beam en- able. Higher luminosity makes rare probes, such constituents of the material are coupled, whereas ergy, meson momentum, and the emission angle as jets and hadrons containing c and b quarks, the speed of sound is sensitive to both the mass with respect to the beam. Untangling initial state more accessible. PHENIX and STAR are being of the constituents and the strength of their in- effects on heavy quark production and final state upgraded with state-of-the-art silicon microvertex teraction. Because strongly coupled theories do effects, which can re-form bound states, requires detectors to enable precise measurements of heavy not allow particle-like excitations, the very na- control measurements in (p or d)+nucleus colli- quarks by tagging their decays. These will allow ture of the plasma constituents is a question. A sions, along with theoretical study of color screen- separation of the charm quark from the bottom promising way to measure both quantities in the inginanexpandingplasma. quark, which is three times heavier and should QGP is by systematic studies of the response Highly energetic partons created in the initial sail right through the QGP. A new ion source at on July 19, 2012 of the matter to initial density fluctuations. phase of the collision lose energy as they pass RHIC provides additional beam species, such as A B 2 R AA Au + Au S = 200GeV NN Neutral pion Direct photon 1.5 Single electron Jet 0, pt: 205.1 GeV ? Jet 1, pt: 70.0 GeV www.sciencemag.org 1 0.5 Downloaded from 0 0 51015 20 p [GeV/c] T Fig. 4. (A) Ratio of particle yield in Au+Au to p+p collisions by PHENIX at a green shape. The black arrows indicate the path of the energetic partons that RHIC, as a function of particle momentum transverse to the beam direction. create the two jets. The rest of the figure depicts a Pb+Pb event display at LHC Data from (25) and adapted with permission from (52, 58) [Copyrighted by from the Compact Muon Solenoid (CMS) (28) [Credit: CMS Collaboration; the American Physical Society]. The emission of pions and electrons from the reprinted with permission from https://twiki.cern.ch/twiki/bin/view/CMSPublic/ decay of heavy quarks is strongly suppressed in Au+Au, whereas photon HIEventDisplays]. The angular distribution of energy emitted has been unfolded emission is not suppressed. The suppression of hadrons is a measure of the onto a plane; the height of the peaks is proportional to the amount of energy color opacity of the QGP. (B) The cartoon illustrates energy measured in the jet observed. The Pb beams enter perpendicular to the page. This event was trig- cone on each side of a dijet, along with energy deposit into the QGP, shown as gered by the jet on the right, and a large energy loss by the jet on the left is seen. In the QGP, color is screened, akin to the elec- through by exciting modes of the medium (col- the highly deformed uranium nucleus. U+U colli- tromagnetic Debye screening observed in conven- lisional energy loss) or by radiating off gluons sions, along with asymmetric beam combinations, tional plasmas. Lattice QCD shows that above (radiative energy loss). The second mechanism, offer novel ways to control the nuclear geometry, T c , screening of color is incomplete (43–45); par- akin to bremsstrahlung of photons by electrons and thus the path length, for probes transiting the tial screening is characteristic of strongly coupled passing through matter, becomes less effective QGP. Excitation functions for rare probes can plasmas (46). Screening in the QGP can be probed as the mass of the parton increases. In a weakly now be used to study predicted features of the QCD experimentally by measuring the survival rate of coupled medium, thermalization of deposited en- matter phase diagram. These measurements, along heavy quark bound states. Charm or bottom quarks ergy occurs through a cascade of collisions among with RHIC’s polarized proton-proton collisions, are produced in pairs, which sometimes remain quarks and gluons in the plasma; in a strongly will systematically test theoretical models, pro- bound and are detected as heavy mesons called coupled medium, the energy is dissipated directly vide benchmarks to isolate the effects of hot QCD quarkonia. The mesons called J/y and y′ are com- into thermal excitations and sound waves. Mea- matter on observables, and map the parton struc- posed of charm quarks, the Upsilon mesons (U) surements of quarks with different mass should ture of nuclei in the relevant kinematic range. SCIENCE 313 www.sciencemag.org VOL 337 20 JULY 2012

REVIEW The wide range of beam energies available Challenging theoretical advances, including 20. P. K. Kovtun, D. T. Son, A. O. Starinets, Phys. Rev. Lett. at RHIC also makes it possible to explore the higher-order jet calculations and effective theo- 94, 111601 (2005). 21. The Kovtun-Son-Starinets bound (20) has been found phase diagram of QCD matter at higher baryon ries for heavy quarks that connect lattice simula- to be violated in certain strong coupled gauge theories densities, because nucleons are partially stopped tions with transport processes, are needed to extract [see (55)]. It would be interesting to determine whether in collisions at lower energies (53). According reliable values for the energy loss parameters and the quark-gluon plasma produced in experiments violates to some predictions, the transition between had- the color screening length in the plasma from this bound. 22. K. Adcox et al.; PHENIX Collaboration, Phys. Rev. Lett. 88, ron and quark matter becomes of first order high-precision data. Major numerical advances 022301 (2002). beyond a critical point in the phase diagram. The will be required to solve the transport equations 23. J. Adams et al.; STAR Collaboration, Phys. Rev. Lett. 91, results from an exploratory beam energy scan, describing rapid formation of an equilibrated 172302 (2003). 24. C. Adler et al., Phys. Rev. Lett. 90, 082302 (2003). which hopes to locate this critical point by search- QGP. Development of an exact gravity dual of 25. S. Afanasiev et al., PHENIX Collaboration, Measurement ing for the signatures of critical fluctuations, are QCD would enable realistic calculations of dy- of direct photons in Au + Au collisions at p S NN = 200 ffiffiffiffiffiffiffi expected soon. Higher beam luminosities will namical processes in the strong coupling limit. GeV, arXiv:1205.5759 [nucl-ex] greatly increase the sensitivity of future searches. So-called “stringy” corrections for finite num- 26. B. Müller, J. Schukraft, B. Wyslouch, First Results from Pb+Pb collisions at the LHC, arXiv:1202.3233 [hep-ex]. The upgrades already under way at RHIC bers of colors would allow bracketing the real- 27. H. Song, S. A. Bass, U. Heinz, Phys. Rev. C Nucl. Phys. will answer some of the questions posed in the world regime of intermediate coupling from both 83, 054912 (2011). preceding section. Fully exploiting the versatil- sides. Such advances will not only elucidate the 28. S. Chatrchyan et al.; CMS Collaboration, Phys. Rev. C ity and luminosity of RHIC requires further de- physics of the QGP but also address intellectual Nucl. Phys. 84, 024906 (2011). 29. S. Borsányi et al.; Wuppertal-Budapest Collaboration, tector capabilities. Accessing the early dynamics challenges of strong coupling in many areas of J. High Energy Phys. 2010, 73 (2010). calls for measuring photons over a larger range physics. The successes and limitations of the string 30. H. Petersen, J. Steinheimer, G. Burau, M. Bleicher, of angles, as well as the flow patterns of those theory–based approach suggest opportunities for H. Stöcker, Phys. Rev. C Nucl. Phys. 78, 044901 (2008). 31. J. M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998). photons. Quantifying the color screening length the development of novel mathematical tech- 32. O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri, entails precise measurements of heavy quarkonia niques applicable to strongly coupled systems. Y. Oz, Phys. Rep. 323, 183 (2000). over a large acceptance, as a function of the mass, Exploration of hot QCD matter has made 33. M. Panero, Phys. Rev. Lett. 103, 232001 (2009). on July 19, 2012 momentum, binding energy, path length through enormous progress during the past decade. Ex- 34. S. Datta, S. Gupta, Phys. Rev. D Part. Fields Gravit. Cosmol. 82, 114505 (2010). the plasma, and initial temperature of the sys- periments have discovered a new high-temperature 35. K. M. O’Hara, S. L. Hemmer, M. E. Gehm, S. R. Granade, tem. Disentangling initial- and final-state effects phase, the strongly coupled QGP, which persists J. E. Thomas, Science 298, 2179 (2002). on bound-state formation implies comparably to the highest temperature probed. Surprising 36. J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, sized sets of nucleon-nucleus collision data. To features of the QGP include near-perfect fluidity J. E. Thomas, Phys. Rev. Lett. 92, 150402 (2004). 37. C. Cao et al., Science 331, 58 (2011). pin down the role of gluons in the nuclear wave and extreme opaqueness to all colored probes. 38. S. Sachdev, Annu. Rev. Con. Mat. Phys. 3, 9 (2012). function (54), forward angle detectors for photons, The rapid development of theoretical and exper- 39. S. Hartnoll, Science 322, 1639 (2008). leptons, and hadrons are necessary. To determine imental tools promises quantitative insights into 40. C. L. Chan et al., Dusty plasma liquids, arXiv:physics/ www.sciencemag.org 0410042 [physics.plasm-ph]. the jet energy loss mechanism, full reconstruc- the still mysterious properties of the QGP during 41. V. Balasubramanian et al., Phys.Rev.Lett. 106, 191601 tion of moderate energy jets, which radiate par- the coming decade. These will also inform the (2011). tons near the QGP scale, will be indispensable. study of other strongly coupled systems in nature 42. P. M. Chesler, L. G. Yaffe, Phys. Rev. Lett. 106, 021601 Measurement of the energy redistribution in the and in the laboratory. (2011). 43. M. Asakawa, T. Hatsuda, Phys. Rev. Lett. 92, 012001 medium-modified jets coupled with modeling of (2004). energy and particle flow will help tease the dif- References and Notes 44. S. Datta, F. Karsch, P. Petreczky, I. Wetzorke, Phys. Rev. D ferent mechanisms apart. RHIC is the ideal fa- Part. Fields Gravit. Cosmol. 69, 094507 (2004). 1. D. J. Gross, F. Wilczek, Phys. Rev. Lett. 30, 1343 (1973). 45. H.-T. Ding et al., Proc. Sci. LATTICE 2010, 180 (2010); cility to do this, whereas the higher-energy jets 2. H. D. Politzer, Phys. Rev. Lett. 30, 1346 (1973). http://pos.sissa.it/archive/conferences/105/180/ Downloaded from at LHC yield the parton energy dependence of 3. E. W. Kolb, M. S. Turner, The Early Universe (Redwood Lattice%202010_180.pdf. the jet-QGP interaction. City, Addison-Wesley, 1988). 46. D. Sarmah, M. Tessarotto, M. Salimullah, Phys. Scr. 74, 4. P. Petreczky, Nucl. Phys. A. 830, 11c (2009). The LHC data will provide stringent tests of 288 (2006). 5. A. Adare et al.; PHENIX Collaboration, Phys. Rev. C Nucl. 47. B. Alessandro et al.; NA50 Collaboration, Eur. Phys. J. C jet quenching theory complementary to those at Phys. 81, 034911 (2010). 39, 335 (2005). RHIC—for example, via the momentum depen- 6. S. Borsányi et al., J. High Energy Phys. 2010,77 48. A. Adare et al.; PHENIX Collaboration, Phys. Rev. C dence of heavy quark energy loss, which is pre- (2010). Nucl. Phys. 84, 054912 (2011). 7. I. Arsene et al.; BRAHMS Collaboration, Nucl. Phys. A. 49. G. Aad et al.; ATLAS Collaboration, Phys. Lett. B 697, 294 dicted to be different in strongly and weakly 757, 1 (2005). (2011). coupled regimes of the QGP. The higher beam 8. K. Adcox et al.; PHENIX Collaboration, Nucl. Phys. A. 50. S. Chatrchyan et al., CMS Collaboration, J. High Energy energy at LHC makes the rate of rare probes 757, 184 (2005). Phys. 1205, 063 (2012). much higher than at RHIC. This opens a larger 9. B. B. Back et al.; PHOBOS Collaboration, Nucl. Phys. A. 51. B. Abelev et al., ALICE Collaboration, J/y production at p kinematic range for hadrons, photons, and b quarks 757, 28 (2005). low transverse momentum in Pb-Pb collisions at ffiffiffiffiffiffiffi S NN = 10. J. Adams et al.; STAR Collaboration, Nucl. Phys. A. 757, 2.76 TeV, arXiv:1202.1383 [hep-ex]. and accesses new probes such as the Z boson. 102 (2005). 52. A. Adare et al.; PHENIX Collaboration, Phys. Rev. Lett. 98, Study of hadrons inside jets energy is determined 11. B. Müller, J. L. Nagle, Annu. Rev. Nucl. Part. Sci. 56, 172301 (2007). by an opposing Z boson will provide new con- 93 (2006). 53. B. I. Abelev et al.; STAR Collaboration, Phys. Rev. C Nucl. straints on the parton-medium interaction. 12. K. H. Ackermann et al., Phys. Rev. Lett. 86, 402 Phys. 81, 024911 (2010). (2001). The higher energy density reached at LHC 54. L. McLerran, Nucl. Phys. A. 752, 355 (2005). 13. P. F. Kolb, P. Huovinen, U. Heinz, H. Heiselberg, 55. Y. Kats, P. Petrov, J. High Energy Phys. 2009, 044 should lead to stronger color screening and thus Phys. Lett. B 500, 232 (2001). (2009). to larger suppression of heavy quark bound states. 14. D. Teaney, J. Lauret, E. V. Shuryak, Phys. Rev. Lett. 86, 56. B. Alver et al., Phys. Rev. C Nucl. Phys. 77, 014906 This effect may be overwhelmed by coalescence 4783 (2001). (2008). 15. B. Schenke, S. Jeon, C. Gale, Phys. Rev. C Nucl. Phys. 85, of heavy quarks into such states when the QGP 57. B. Alver, G. Roland, Phys. Rev. C Nucl. Phys. 81, 054905 024901 (2012). (2010). converts back into hadrons; coalescence should 16. M. Luzum, J. Phys. G 38, 124026 (2011). 58. A. Adare et al.; PHENIX Collaboration, Phys. Rev. Lett. be more prominent at LHC than at RHIC be- 17. A. Adare et al.; PHENIX Collaboration, Phys. Rev. Lett. 101, 232301 (2008). cause more heavy quark pairs are produced in 107, 252301 (2011). 18. H. Song, S. A. Bass, U. Heinz, T. Hirano, C. Shen, each collision. High statistics spectroscopy of the Acknowledgments: This work was supported by grants from Phys. Rev. Lett. 106, 192301 (2011). U states will allow comparison of color screening 19. P. Danielewicz, M. Gyulassy, Phys. Rev. D Part. Fields 31, the U.S. Department of Energy. in QGP at different temperatures. 53 (1985). 10.1126/science.1215901 314 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

2.8 Million Years of Arctic Climate Change from Lake El'gygytgyn, NE Russia Martin Melles et al. Science 337 , 315 (2012); DOI: 10.1126/science.1222135 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/315.full.html Supporting Online Material can be found at: http://www.sciencemag.org/content/suppl/2012/06/21/science.1222135.DC1.html www.sciencemag.org This article cites 95 articles , 9 of which can be accessed free: http://www.sciencemag.org/content/337/6092/315.full.html#ref-list-1 This article appears in the following subject collections: Downloaded from Atmospheric Science http://www.sciencemag.org/cgi/collection/atmos Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

RESEARCH ARTICLE (7). The 170-m-deep lake has a bowl-shaped mor- phology with a diameter of ~12 km, a surface 2 2.8 Million Years of Arctic Climate area of 110 km , and a relatively small catchment 2 of 293 km (8). The modern continental Arctic climate produces herb-dominated tundra in the Change from Lake El’gygytgyn, catchment, 9 months per year of lake-ice cover, and oligotrophic to ultra-oligotrophic conditions NE Russia in the lake. Low productivity in combination with complete overturning of the water column during the ice-free period in summer leads to well- 2 1 3 Martin Melles, * Julie Brigham-Grette, Pavel S. Minyuk, Norbert R. Nowaczyk, 4 oxygenated bottom waters throughout the year. 2 5 1 Volker Wennrich, Robert M. DeConto, Patricia M. Anderson, Andrei A. Andreev, 1 Scientific deep drilling was performed in the 2 2 4 Anthony Coletti, Timothy L. Cook, † Eeva Haltia-Hovi, ‡ Maaret Kukkonen, 1 El’gygytgyn Crater in winter 2008/2009 (9). We 7 6 1 3 Anatoli V. Lozhkin, Peter Rosén, Pavel Tarasov, Hendrik Vogel, Bernd Wagner 1 used advanced high-resolution (logging/scanning) technologies and standard techniques to investi- The reliability of Arctic climate predictions is currently hampered by insufficient knowledge of gate the core composite from site 5011-1 (Fig. 1) natural climate variability in the past. A sediment core from Lake El’gygytgyn in northeastern of the International Continental Scientific Drill- (NE) Russia provides a continuous, high-resolution record from the Arctic, spanning the past ing Program (ICDP) for lithology as well as se- 2.8 million years. This core reveals numerous “super interglacials” during the Quaternary; for lected physical, chemical, and biological proxies marine benthic isotope stages (MIS) 11c and 31, maximum summer temperatures and annual (10). According to the age model (Fig. 2), which precipitation values are ~4° to 5°C and ~300 millimeters higher than those of MIS 1 and 5e. is based on magnetostratigraphy and tuning of Climate simulations show that these extreme warm conditions are difficult to explain with proxy data to the regional insolation and global on July 19, 2012 greenhouse gas and astronomical forcing alone, implying the importance of amplifying feedbacks marine isotope stratigraphy (fig. S1), the upper and far field influences. The timing of Arctic warming relative to West Antarctic Ice Sheet 135.2 m of the sediment record continuously retreats implies strong interhemispheric climate connectivity. represents the environmental history of the past 2.8 My. To display the obtained data versus time he effects of global warming are docu- from Lake El’gygytgyn, located ~100 km to the (Fig. 3) we removed volcanic ashes and other mented and predicted to be most pro- north of the Arctic Circle in northeastern Russia event layers caused by mass movement deposits Tnounced in the Arctic, a region that plays a (67.5°N, 172°E) (Fig. 1). The length, temporal con- (10). Though highly varied in nature, the result- crucial, but not yet well-understood role within tinuity, and centennial- to millennial-scale temporal ant record of pelagic sedimentation consists of www.sciencemag.org the global climate system (1). Reliable climate resolution (Fig. 2 and supplementary materials) three dominant lithofacies (see supplementary projections for high northern latitudes are, how- provide a detailed view of natural climatic and en- materials). Climate and environmental interpre- ever, hampered by the complexity of the under- vironmental variability in the terrestrial Arctic, a tation of the pelagic sedimentation record is based lying natural variability and feedback mechanisms better understanding of the representative nature on complementary biological and geochemical (2, 3). To date, information concerning the nat- of the last climate cycle for the Quaternary, and indicators, which show that distinct facies reflect ural climate variability in the Arctic is widely insight into how sensitively the terrestrial Arctic end-member glacial/interglacial climatic condi- restricted to the last glacial/interglacial cycle, the reacts to a range of forcing mechanisms. tions (9, 11). period covered by the longest ice-core records Lake setting, drilling, and core analyses. Glacial variability and proxies. Facies A is Downloaded from from the Greenland ice cap (4). A limited number Lake El’gygytgyn is located in a meteorite im- characterized by dark gray to black finely lam- of records extend deeper in time, from both the pact crater formed 3.58 million years ago (Ma) inated (<5 mm) silt and clay and may contain marine realm (5) and the Arctic borderland (6), but these records are restricted in terms of age control and temporal resolution. Fig. 1. Location of Lake Here, we present a time-continuous and high- El’gygytgyn in northeast- NW SE resolution record of environmental history from ern Russia (inserted map) the Arctic spanning the past 2.8 million years (My) andschematiccross-section of the El’gygytgyn basin stratigraphy showing the 5011-3 Lake El'gygytgyn 1 Institute of Geology and Mineralogy, University of Cologne, 55 2 Zuelpicher Strasse 49a, D-50674 Cologne, Germany. Depart- location of ICDP sites 141 5011-1 ment of Geosciences, University of Massachusetts, 611 North 5011-1 and 5011-3. At Permafrost sediments Lz1024 AB C 0 3 Pleasant Street, Amherst, MA 01003, USA. Far East Branch site 5011-1, three holes Quaternary 123 Russian Academy of Sciences, North-East Interdisciplinary (1A, 1B, and 1C) were sediments Lacustrine Scientific Research Institute, 16 Portovaya Street, 685000 drilled to replicate the Pliocene (m) 4 Magadan, Russia. Helmholtz Centre Potsdam, GFZ German Quaternary and upper- 318 Research Centre for Geosciences, Telegrafenberg C321, D-14473 Atlantic 420 5 Potsdam, Germany. Department of Earth and Space Sciences, most Pliocene sections. central rocks Impact University of Washington, Box 351310, Seattle, WA 98195– Hole 1C further pene- Green- uplift 517 0km2 land 6 1310, USA. Climate Impacts Research Centre, Umeå Univer- trated through the remain- Lacustrine sediments 7 sity,SE-98107Abisko,Sweden. InstituteofGeologicalSciences, ing lacustrine sequence Arctic Ocean Free University Berlin, Malteserstrasse 74-100, Haus D, D-12249 down to a 318-m depth Impact rock Berlin, Germany. Siberia Permafrost sediments and then ~200 m into Lake *To whom correspondence should be addressed. E-mail: the impact rock sequence El'gygytgyn Core loss/gap [email protected] underneath. Lz1024 is a Alaska supposed talik/permafrost border †Present address: Department of Physical and Earth Sciences, Worcester State University, Worcester, MA 01602, USA. 16-m-longpercussionpis- ‡Present address: Department of Geology, Lund University, ton core taken in 2003 Sölvegatan 12, S-223 62 Lund, Sweden. that fills the stratigraphic gap between the lake sediment surface and the top of drill cores 1A and 1B. SCIENCE 315 www.sciencemag.org VOL 337 20 JULY 2012

RESEARCH ARTICLE elongated sediment clasts of coarser grain sizes sonal ice cover promotes higher diatom produc- and high regional July insolation (Fig. 3, A and (fig. S2). This facies was deposited during times tivity, as indicated by high Si/Ti ratios (Fig. 3H B). The characteristics of facies C suggest that it of heavy global marine isotopic values (12)and and supplementary materials). In contrast, TOC represents particularly warm interglacials. High low regional July insolation (Fig. 3, A, B, and D) content is low (Fig. 3F), suggesting high organic Mn/Fe ratios (Fig. 3G) along with reddish-brown (13). Facies A represents peak glacial conditions, matter decomposition due to oxygenation of bot- sediment colors imply well-oxygenated bottom when perennial lake ice persisted, requiring mean tom waters as a consequence of wind- and density- waters. In contrast with facies B, however, the annual air temperatures at least 4 (T 0.5)°C lower driven mixing. Complete water-column ventilation sediments are distinctly laminated. This is traced than today (14). This resulted in a stagnant water is also indicated in the sediment colors, maxima back to a combination of factors, including a par- column with oxygen-depleted bottom waters, as in MS (Fig. 3E), and high Mn/Fe ratios (Fig. 3G). ticularly high primary production in spring and reflected by low Mn/Fe ratios (Fig. 3G) and mini- In addition, the lack of stratification indicates summer and anoxic bottom water conditions dur- ma in magnetic susceptibility (MS) (Fig. 3E), minor sediment homogenization by bioturbation. ing winter stratification under a seasonal ice cover, indicating reducing conditions with magnetite dis- Facies C consists of reddish-brown silt-sized which excludes bioturbation despite annual oxy- solution (see supplementary materials). Dark lami- sediment with distinct fine laminations (<5 mm) genation. High primary production, presumably nations along with maxima in the content of total (fig. S2). This facies is irregularly distributed in caused by a longer ice-free season and enhanced organic carbon (TOC) (Fig. 3F) reflect the absence the record compared with facies A and B (Fig. nutrient supply from the catchment relative to of bioturbation and enhanced preservation of or- 3D). Facies C coincides with some periods of other interglacials, is indicated by exceptionally ganic matter. Low Si/Ti ratios (Fig. 3H) and a ro- light values in the global marine isotope record high Si/Ti ratios (Fig. 3H). Anoxic bottom water bust correlation between Si/Ti ratios and biogenic silica contents (see supplementary materials), how- ever, suggest relatively low primary production. Brunhes Matuyama Gauss Gilbert Facies A first appears 2.602 to 2.598 Ma, dur- ing marine isotope stage (MIS) 104 (Fig. 3D), Mammoth corresponding with pollen assemblages that indi- Jaramillo Cobb Mountain Olduvai Réunion Kaena on July 19, 2012 cate substantial cooling at the Pliocene/Pleistocene 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 boundary (see supplementary materials). This cooling coincides with distinct climatic deterio- ration at Lake Baikal (15) and may be associated 0 0 with the poorly dated Okanaanean Glaciation in eastern Chukotka at the beginning of the Pleis- tocene (16). On the other hand, the first occur- rence of facies A at Lake El’gygytgyn clearly www.sciencemag.org 50 50 postdates the onset of stratification across the western subarctic Pacific Ocean at 2.73 Ma, an event believed to have triggered the intensifi- cation of Northern Hemispheric glaciation (17). Hence, the onset of full glacial cycles in central 100 100 Chukotka cannot directly be linked to changes in Quaternary thermohaline circulation in the Pacific. From the long-term succession of facies A Pliocene Downloaded from (Fig. 3D) and Mn/Fe ratios (Fig. 3G), pervasive 150 150 glacial episodes at Lake El’gygytgyn gradually Composite Depth (mblf) increase in frequency from ~2.3 to ~1.8 Ma, even- tually concurring with all glacials and several stadials reflected globally in stacked marine iso- 200 200 tope records (12). The full establishment of glacial/ Sedimentation rates: interglacial cycles by ~1.8 Ma at Lake El’gygytgyn coincides well with enhanced glacial erosion in British Columbia (18) and the onset of subpolar 4 cm kyr -1 cooling in both hemispheres with an average 250 5 cm kyr -1 250 bipolar temperature drop of 4° to 5°C due to the emergence of the tropical Pacific cold tongue (19). Nevertheless, this event clearly predates the 50 cm kyr -1 mid-Pleistocene transition, when the dominance 300 300 of 41 thousand years (ky) of obliquity was glob- ally replaced by the 100-ky cycle between 1.25 and 0.7 Ma (20). Interglacial variability and proxies. Facies B 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 is characterized by massive to faintly banded silt Age (Ma) that is olive gray to brownish in color (fig. S2). Fig. 2. Age/depth model with resulting sedimentation rates for the ICDP 5011-1 core composite based on This facies comprises the majority of sediment magnetostratigraphy and correlation between sediment proxy data, the LR04 marine isotope stack (12), that accumulated in Lake El’gygytgyn during the and regional spring and summer insolation (13). Initial first-order tie points are indicated by black past 2.8 My, representing 79% of the Quaternary diamonds; second- and third-order tie points are denoted by the blue curve. The red star marks the time of 39 40 history (Fig. 3D). Facies B reflects a wide range the impact inferred from Ar/ Ar dating (7)at3.58(T0.04) Ma. Black and white bars denote normal and of glacial to interglacial settings and includes the reversed polarity, respectively. Mass movement deposits and core gaps greater than 50 cm in thickness style of modern sedimentation. As such, a sea- areindicated on theright y axis in gray and blue, respectively. mblf, meters below lake floor. 316 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

RESEARCH ARTICLE conditions during winter are implied by high TOC 91, and 93 (red bars in Fig. 3), suggesting that tional biological proxies and pollen-based climate contents (Fig. 3F), reflecting high primary pro- these interglacials represent unusual “super in- reconstructions (Fig. 3, I to L). duction and incomplete decomposition compared terglacials” in the Arctic throughout the Quater- Sediments formed in Lake El’gygytgyn dur- with facies B, and variable MS values (Fig. 3E), nary. The exceptional character of these interglacial ing MIS 1 and 5e have Si/Ti ratios only slight- reflecting partial dissolution of magnetite. conditions becomes evident based on a compar- ly higher than those formed during glacial and The described characteristics of facies C are ison of MIS 1 and 5e (facies B) with MIS 11c and stadial conditions of MIS 2, 5d, and 6 (Fig. 3K). most pronounced for MIS 11c, 31, 49, 55, 77, 87, 31 (super interglacials of facies C), using addi- Pollen data show distinct increases in tree and Fig. 3. (A to H)(A) LR04 global marine isotope stack (12) and (B) mean July in- solation for 67.5°N (13)for the past 2.8 My compared with(C)magnetostratigraphy, (D) facies, (E) magnetic sus- ceptibility, (F) TOC contents, (G) Mn/Fe ratios, and (H) Si/Ti ratios in the sediment re- cordfrom Lake El’gygytgyn (magnetic susceptibility and x-ray fluorescence data are smoothed using a 500-year weighted running mean to on July 19, 2012 improve the signal-to-noise ratio). Super interglacials at Lake El’gygytgyn are high- lighted with red bars. (I to L)Expanded views into the interglacials MIS 1, 5e, 11c, and31andadjoiningglacials/ stadials. (I) Reconstructed MTWM and (J) PANN based www.sciencemag.org on the pollen spectra and bestmodernanalogapproach [modern values from (56)]. (K) Mean July insolation for 67.5°N (13) compared with El’gygytgyn Si/Ti ra- tios, smoothed by five-point weighted running mean. (L) Tree and shrub pollen per- Downloaded from centagescomparedwithspruce pollencontent.SimulatedJuly surface air temperatures (red andgreen dots)atthe lo- cation of the lake are shown for comparison. The location of the dots relative to the x axis corresponds with the GHGandorbitalforcingused in each interglacial simula- tion (see supplementary ma- terials). Simulated modern and preindustrial tempera- tures are close to observed values, so model tempera- tures are not corrected for bias. The green dot indicates the results derived with a deglaciated Greenland and increasedheatfluxunderArc- −2 tic Ocean sea ice by 8 W m . SCIENCE 317 www.sciencemag.org VOL 337 20 JULY 2012

RESEARCH ARTICLE shrub pollen (Fig. 3L) and suggest that birch and have occurred ~1.2 Ma (37) and thus may be cor- mer. The net effect of this orbital configuration alder shrubs dominated the vegetation (fig. S4). relative with MIS 31. Another possibly correla- produces high-intensity summer insolation at the Pollen-based climate reconstructions (see supple- tive site is at Fosheim Dome on Ellesmere Island. lake, >50 Wm −2 greater than today (Fig. 3K). mentary materials) suggest that the mean temper- This site includes terrestrial deposits dated to Similarly, peak warmth during MIS 1 and 11c ature of the warmest month (MTWM; i.e., July) ~1.1 Ma, which enclose fossil beetle (Coleoptera) also coincides with perihelion during boreal sum- and the annual precipitation (PANN) during the assemblages, suggesting temperatures 8° to 14°C mer, but lower eccentricity (and lower obliquity peak of MIS 1 and 5e were only ~1° to 2°C and, above modern values (38). at MIS 11c) attenuates the effect of precession with one exception, ~50 mm higher than today, Other Arctic sites potentially correlative with relative to MIS 5e and 31, making summer inso- respectively (Fig. 3, I and J). one or more of the older Early Pleistocene super lation forcing less intense, albeit longer in duration. This observation is consistent with temper- interglacials recorded in Lake El’gygytgyn (Fig. GHG radiative forcing from a combination of ature reconstructions for the Holocene thermal 3) include the Kap København Formation in CO 2 ,CH 4 ,and N 2 O atmospheric mixing ratios maximum, which indicate +1.6 (T0.8)°C warm- northern Greenland, currently dated to ~2.4 Ma. determined from ice cores (see supplementary inginthe westernArctic (21)and +1.7 (T0.8)°C At this time, sea ice was strongly reduced and materials) is similar during MIS 5e and 11c (+0.16 across the entire Arctic (3) relative to modern, con- forests reached the Arctic Ocean about 1000 km and +0.19 Wm −2 relative to preindustrial GHG firming that Lake El’gygytgyn records regional further to the north than today (39). Another can- concentrations, respectively). Early MIS 1 is clear- rather than just local climate change (14). In con- didate is the balmy Bigbendian Transgression of ly an exception, with substantially lower CO 2 trast, temperature reconstructions for the MIS 5e the Gubik Formation dated to ~2.6 Ma (36). levels [~260 parts per million by volume (ppmv)] thermal maximum are more variable, indicating Interglacial forcings and feedbacks. Compar- around the time of peak Holocene warmth [~9 +5(T1)°C across the entire Arctic, albeit with ing the relative warmth of the Pleistocene inter- thousand years ago (ka)] producing –0.44 Wm −2 smaller anomalies reconstructed for the Pacific glacials recorded at Lake El’gygytgyn (Fig. 3I) in less radiative forcing relative to preindustrial lev- sector (3, 22). The warmer climate across the the context of orbital and greenhouse gas (GHG) els. MIS 31 (~1.072 Ma) lies beyond the oldest Arctic during MIS 5e compared with MIS 1 is forcing (40), we find that peak summer warmth ice cores, so no direct information on atmospher- thought to have caused a size reduction of the during MIS 5e and 31 corresponds to the con- ic composition is available. However, a proxy- on July 19, 2012 Greenland Ice Sheet equivalent to 1.6 to 2.2 m in gruence of high obliquity, high eccentricity, and based reconstruction of mid-Pleistocene partial global sea-level rise (23). precession aligning perihelion with boreal sum- pressure of CO 2 based on boron isotopes in plank- Strongly enhanced primary productivity dur- ing the super interglacials MIS 11c and 31 com- 1 pared with MIS 1 and 5e, as inferred from higher AB 1 Si/Ti ratios (Fig. 3K), is associated with compa- 3 3 1 1 −1 rable maxima in tree and shrub pollen but is 3 1 5 3 3 3 5 1 7 5 1 3 5 1 1 marked by distinct differences in pollen com- 5 1 3 3 3 www.sciencemag.org 5 3 1 5 1 3 7 1 1 5 1 3 position (Fig. 3L and supplementary materials). 5 1 3 3 3 5 5 3 5 5 3 1 3 1 3 ΔT °C 1 For instance, substantial spruce pollen is present 3 1 3 1 5 1 1 3 7 9 3 1 1 3 5 1 3 5 3 16 9 1 1 3 1 5 1 5 3 7 1 3 3 1 3 7 during MIS 11c and 31 but is missing during 7 1 5 1 3 3 1 3 5 1 1 1 1 5 5 1 3 1 5 1 5 5 5 15 shrub-dominated MIS 1 and 5e interglacials. Ac- 7 7 1 1 1 3 3 1 1 1 3 1 5 5 5 3 5 3 3 cording to the best modern analog (BMA; see 7 3 3 3 9 9 5 5 5 3 1 5 5 14 3 3 7 1 supplementary materials) climate reconstruction, 5 5 3 5 1 3 3 3 13 5 3 9 7 maximum MTWM and PANN were up to 4° to 1 1 3 9 1 3 1 12 1 7 3 1 5°C and ~300 mm higher than those of MIS 1 and 3 5 5 11 Downloaded from 5e, respectively (Fig. 3, I and J). 1 3 1 10 Sediment records of MIS 11 are rare in the 1 3 1 9 Arctic, and their reconstructed temperature sig- 8 nals are inconclusive (22). However, there are C D 7 indications that the Greenland Ice Sheet was 3 2 2 much smalleroreven absent (24, 25), with forests 8 1 9 6 4 4 3 5 11 covering at least South Greenland (26). Relative 4 4 12 10 6 2 3 3 5 5 3 7 5 4 5 13 5 sea level may have been significantly higher than 4 2 6 2 2 7 3 5 5 3 4 4 2 2 2 16 2 5 3 5 17 5 1 5 3 today (25, 27). Particularly warm conditions are 4 6 12 14 0 0 2 3 15 1 3 4 2 4 5 3 7 13 1 1 2 2 6 10 2 2 7 1 3 5 also suggested by records from Lake Biwa (28), 2 6 8 2 4 3 1 11 9 5 1 3 5 2 2 0 2 2 2 1 3 5 5 1 3 2 4 Lake Baikal (29), the mid-latitude Atlantic (30), 4 2 4 0 0 0 2 0 4 7 5 1 1 1 3 1 1 3 4 0 0 2 4 2 7 3 3 3 1 4 and the Belize Reef (31) and may have been 0 0 1 1 2 5 5 associated with megadroughts in the southwest- 3 1 3 3 0 6 2 2 1 ern United States (32). 6 2 2 0 0 0 0 2 7 1 3 −1 As yet, MIS 31 is not unambiguously re- 4 4 2 5 5 3 1 3 2 5 3 2 corded in the Arctic, but it is known for substan- 2 3 tial changes in and around Antarctica, including a 2 3 4 southward shift of the subtropical front and warm- 4 2 3 3 er waters in the Southern Ocean (33, 34)and the Fig. 4. Simulated interglacial warming (2-m surface temperature in degrees Celsius) relative to pre- collapse of the West Antarctic Ice Sheet (WAIS) industrial levels. (A) MIS 1 (9-ky orbit and GHGs). (B) MIS 5e (127-ky orbit and GHGs). (C) MIS 11c (409-ky (35). In the Northern Hemisphere, the Plio/ orbit, GHGs, no Greenland Ice Sheet, and 8 W m −2 enhanced oceanic heat convergence under Arctic sea Pleistocene Gubik Formation of northern Alaska ice). (D) MIS 31 (1072-ky orbit, GHGs, and no Greenland Ice Sheet). Orbital and GHG forcing for MIS 5e includes at least five sea-level high stands as- and 11c follow that used by Yin and Berger (40). Forcing for MIS 31 follows that used by DeConto et al. sociated with episodes of warm climate and re- (42). ThelocationofLakeEl’gygytgyn is shown with a star near the bottom-center of each panel. Areas of duced sea ice (36). One of these episodes, the no shading (white) roughly correspond to statistically insignificant anomalies at the 95% confidence Fishcreekian transgression, is now thought to interval. 318 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

RESEARCH ARTICLE tonic foraminifera (41) indicates that the highest ocean overturning (ignored in our simulations) ditional warming (<0.7°C) (Figs. 3I and 4C) in mid-Pleistocene CO 2 levels (~325 ppmv) occurred generally have a cooling effect on the Northern the Beringian interior. around 1 Ma, roughly coinciding with the excep- Hemisphere, adding to the difficulty in explain- Fully testing these ideas will require additional tional warmth of MIS 31. Though uncertain, these ing the exceptional warmth at MIS 11c relative to climate-ocean modeling, explicitly accounting for reconstructed CO 2 levels at MIS 31 would have MIS 1 and 5e. glacial/interglacial changes in regional sea level added ~0.84 Wm −2 of radiative forcing, even if The super interglacials at Lake El’gygytgyn (paleobathymetry and gateways), changes in land- CH 4 and N 2 O mixing levels remained close to coincide remarkably with diatomite layers in ice distributions, and melt-water inputs in both po- preindustrial values, which is unlikely consider- the Antarctic ANDRILL 1B record (see supple- lar regions, as well as contemporaneous sediment ing the ubiquitous correlation of elevated CH 4 mentary materials), which reflect periods of a records from the Arctic and North Pacific Oceans. and N 2 O during late Pleistocene interglacials. In diminished WAIS and open water in the Ross ThepaleoclimaticrecordfromLakeEl’gygytgyn sum, much of the warmth during MIS 31 can be Embayment (35, 45). The higher number of events provides a benchmark of Arctic change from explained by elevated greenhouse gas levels (42). at Lake El’gygytgyn does not necessarily reflect an area that has otherwise been a data desert for To investigate potential reasons for the super a higher frequency, but it may also reflect the time-continuous terrestrial records of the Pliocene interglacials at Lake El’gygytgyn, we tested the discontinuity of the ANDRILL 1B record (46). and Pleistocene. The sediments provide a fresh equilibrated response of a Global Climate Model Linkages between extraordinary warmth at window into the environmental dynamics of the (GCM) with an interactive vegetation component Lake El’gygytgyn and Antarctic ice volume im- Arctic from a terrestrial high-latitude site for com- (see supplementary materials) to the orbital and ply strong intrahemispheric climate coupling that parison with other Arctic records. Marine cores GHG forcing corresponding to the timing of peak could be related to reductions in Antarctic Bot- from the Arctic basin, such as those from the summer warmth at MIS 1, 5e, 11c, and 31. Com- tom Water (AABW) formation (47) during times ACEX/Lomonosov Ridge or HOTRAX expedi- parisons with a preindustrial control simulation of ice sheet/shelf retreat and elevated fresh water tions (55), still lack the comparable resolution (Fig. 4) show that differences in MTWM maxima input into the Southern Ocean. This is supported and length to test for perennial versus seasonal at Lake El’gygytgyn during MIS 1 and 5e (+2.1° by distinct minima in AABW inflows into the sea-ice conditions during interglacials over the and +4.2°C) were in the same range as those southwest Pacific during MIS 11 and 31 (48). As past 2.8 My. The attenuated response of Arctic on July 19, 2012 during MIS 11c and 31 (+2.2° and +3.5°C) (Fig. a consequence, changes in thermohaline circula- SSTs in model simulations of the interglacials 3I, red dots, Fig. 4, and supplementary materials). tion during MIS 11 and 31 might have reduced (Fig. 4) (43) relative to surrounding continents The same holds true for the modeled differences upwelling in the northern North Pacific (49), as hints that deep Arctic Ocean cores might not in PANN (0 and –37 mm/a and +38 and 0 mm/a, indicated by distinctly lower BSi concentrations provide a complete perspective of the pacing or respectively). The results are similar to previous compared with other interglacials at Ocean Drill- magnitude of climate change in the Arctic border- interglacial simulations using an intermediate ing Program site 882 (50, 51). A stratified water lands. The observed response of the region’s complexity model (40), with the combined effect column during the super interglacials would have climate and terrestrial ecosystems to a range of of orbital and GHG forcing at MIS 5e producing resulted in higher SSTs in the northern North interglacial forcing provides a challenge for mod- www.sciencemag.org the greatest summer warming among the four Pacific, with the potential to raise air temperatures eling and important constraints on climate sen- interglacials modeled here. Our simulated sum- and precipitation rates over adjacent land masses sitivity and polar amplification. The marked mer warming (4.2°C) over the Beringian interior via effects on the dominant pressure patterns coherence of interglacial warmth across the west- at MIS 5e also closely matches the warming sim- (Siberian high and Aleutian low) that dominate ern Arctic with repeated deglaciation events in ulated by a coupled atmosphere-ocean GCM (43). the modern climatology at the lake (52). West Antarctica supports the notion of strong Consequently, the distinctly higher observed val- An alternative mechanism linking Lake teleconnections between the polar regions over ues of MTWM and PANN at MIS 11c cannot El’gygytgyn with Antarctica could be related to the past 2.8 My. readily be explained by the local summer orbital higher relative sea level due to the combined Downloaded from forcing or GHG concentrations alone and sug- retreats of the WAIS (44) and the Greenland Ice References and Notes gest that other processes and feedbacks contrib- Sheet (24), resulting in enhanced warm-water 1. ACIA, Arctic Climate Impact Assessment (Cambridge Univ. Press, Cambridge, 2005). uted to the extraordinary warmth at this interglacial intrusion into the Arctic Ocean. Potential gate- 2. J. Christensen et al., in Fourth Assessment Report of the and the relatively muted response to the strongest ways are the Denmark Strait and Barents Sea Intergovernmental Panel on Climate Change (Cambridge forcing at MIS 5e. from the Atlantic Ocean and the Bering Strait Univ. Press, Cambridge, 2007), pp. 847–940. Vegetation-land surface feedbacks are ac- from the Pacific Ocean. In the northeastern At- 3. G. H. Miller et al., Quat.Sci.Rev. 29, 1779 (2010). 4. North Greenland Ice Core Project members, Nature 431, counted for in our model, and the simulated lantic, however, SSTs, at least during MIS 11, 147 (2004). poleward advance of evergreen needle-leaf forest were lower than during MIS 9, 5e, and 1 (53). 5. M. O’Regan, C. J. Williams, K. E. Frey, M. Jakobsson, during the interglacials provides a good match Today, Bering Strait throughflow is restricted by Oceanography 24, 66 (2011). with our reconstructions (see supplementary ma- shallow waters of only ~50 m in depth, result- 6. CAPE Last Interglacial Project Members, Quat. Sci. Rev. 25, 1383 (2006). terials), yet the warming effect of boreal forest ing in an average northward transport of ~0.8 7. P. Layer, Meteorit. Planet. Sci. 35, 591 (2000). 6 −1 3 expansion does not provide a satisfactory expla- sverdrups (1 sverdrup = 10 m s )(54). Sub- 8. M. Nolan, J. Brigham-Grette, J. Paleolimnol. 37, 17 (2007). nation for the warmth of MIS 11c. A deglaciated stantial interannual variability in flow rate can 9. M. Melles et al., Sci. Drill. 11, 29 (2011). 20 Greenland has been shown to have important re- produce elevated heat fluxes (5 to 6 × 10 J/year 10. Materials and methods are available as supplementary gional effects on surrounding sea surface temper- in 2007), which can be amplified in the Arctic by materials on Science Online. 11. M. Melles et al., J. Paleolimnol. 37, 89 (2007). atures (SSTs) and sea-ice conditions,but widespread internal-feedback mechanisms (3). No evidence 12. L. E. Lisiecki, M. E. Raymo, Paleoceanography 20, warming in the circum-Arctic (and Beringia in as yet exists for substantial changes in tem- PA1003 (2005). particular) has been shown to be minimal (43, 44). perature or flow rates during super intergla- 13. J. Laskar et al., Astron. Astrophys. 428, 261 (2004). This observation is supported by our simulations, cials; however, as a first exploration of this 14. M. Nolan, Clim. Past Disc. 8, 1443 (2012). 15. D. Demske, B. Mohr, H. Oberhänsli, Palaeogeogr. showing that the loss of the Greenland Ice Sheet idea, we increased the heat flux convergence Palaeoclimatol. Palaeoecol. 184, 107 (2002). at MIS 11c raises summer temperatures at Lake under Arctic sea ice in our interglacial climate 16. A. F. Fradkina et al., Geol. Soc. Am. Spec. Pap. El’gygytgyn by only 0.3°C. Furthermore, Green- model simulations by 8 W m −2 (reflecting an 382 (2005). land was likely reduced in size during MIS 5e extreme ~fourfold increase in warmer Bering 17. G. H. Haug et al., Nature 433, 821 (2005). 18. D. L. Shuster, T. A. Ehlers, M. E. Rusmoren, K. A. Farley, and perhaps other interglacials, offering little help Strait throughflow). The additional heat flux re- Science 310, 1668 (2005). in differentiating Beringia’s response from one sults in substantial reductions in seasonal sea ice 19. A. Martínez-Garcia, A. Rosell-Melé, E. L. McClymont, interglacial to the next. Meltwater impacts on and warmer Arctic SSTs but contributes little ad- R. Gersonde, G. H. Haug, Science 328, 1550 (2010). SCIENCE 319 www.sciencemag.org VOL 337 20 JULY 2012

20. P. U. Clark et al., Quat. Sci. Rev. 25, 3150 (2006). 43. B. L. Otto-Bliesner et al., Science 311, 1751 (2006). Russian Global Lake Drilling 800 drilling system was developed 21. D. S. Kaufman, J. Brigham-Grette, Quat. Sci. Rev. 12, 44. S. J. Koenig, R. M. DeConto, D. Pollard, Clim. Dyn. 37, and operated by DOSECC. We thank all participants of the 21 (1993). 1247 (2011). expedition for their engagement during recovery of the ICDP 22. G. H. Miller et al., Quat. Sci. Rev. 29, 1679 (2010). 45. D. Pollard, R. M. DeConto, Nature 458, 329 (2009). 5011-1 cores. Funding of core analyses was provided by BMBF 23. E. J. Colville et al., Science 333, 620 (2011). 46. R. McKay et al., Quat. Sci. Rev. 34, 93 (2012). (grant 03G0642), German Research Foundation (Deutsche 24. E. Willerslev et al., Science 317, 111 (2007). 47. A. Foldvik, J. Geophys. Res. 109, C02015 (2004). Forschungsgemeinschaft, grants ME 1169/21 and ME 1169/24), 25. M. E. Raymo, J. X. Mitrovica, Nature 483, 453 (2012). 48. I. R. Hall, I. N. McCave, N. J. Shackleton, G. P. Weedon, the NSF (grant 0602471), Vetenskapsrådet, The Swedish Research 26. A. de Vernal, C. Hillaire-Marcel, Science 320, 1622 (2008). S. E. Harris, Nature 412, 809 (2001). Council, the Kempe Foundation, and the Civilian Research and 27. G. A. Milne, J. X. Mitrovica, Quat.Sci.Rev. 27, 2292 (2008). 49. E. D. Galbraith et al., Nature 449, 890 (2007). Development Foundation (grant RUG1-2987-MA-10). We thank 28. P. E. Tarasov et al., Earth Sci. Rev. 108, 64 (2011). 50. G. H. Haug, M. A. Maslin, M. Sarnthein, R. Stax, N. Mantke, A. Shahnazarian, and numerous students (Univ. 29. A. A. Prokopenko et al., Clim. Past 6, 31 (2010). R. Tiedemann, Proc. ODP Sci. Results 145, 293 of Cologne) for their competent help in core processing; 30. R. Stein, J. Hefter, J. Grützner, A. Voelker, B. D. A. Naafs, (1995). T. Matrosova for contributing modern surface pollen data; and Paleoceanography 24, PA2203 (2009). 51. S. L. Jaccard, E. D. Galbraith, D. M. Sigman, G. H. Haug, R. McKay for providing a modified lithological log of the ANDRILL 31. E. Gischler, R. N. Ginsburg, J. O. Herrle, S. Prasad, Quat. Sci. Rev. 29, 206 (2010). 1B record. We also thank four anonymous reviewers for Sedimentology 57, 1049 (2010). 52. C. J. Mock, P. J. Bartlein, P. A. Anderson, Int. J. Climatol. supportive and constructive comments that greatly improved 32. P. J. Fawcett et al., Nature 470, 518 (2011). 18, 1085 (1998). the manuscript. The data reported in this paper are available in 33. R. P. Scherer et al., Geophys. Res. Lett. 35, L03505 (2008). 53. J. P. Helmke, H. A. Bauch, Paleoceanography 18, 1036 the databases of PANGAEA via www.pangaea.de (DOI: 10.1594/ 34. P. Maiorano,M.Marino, J.-A.Flores, Mar. Micropaleontol. (2003). PANGAEA.783305) and of the U.S. National Climatic Data 71, 166 (2009). 54. R. A. Woodgate, T. Weingartner, R. Lindsay, Geophys. Center via www.ncdc.noaa.gov/paleo/paleolim/paleolim_data.html. 35. T. Naish et al., Nature 458, 322 (2009). Res. Lett. 37, L01602 (2010). 36. J. Brigham-Grette, L. D. Carter, Arctic 43, 74 (1992). 55. L. Polyak, M. Jakobsson, Oceanography 24, 52 (2011) 37. G. A. Goodfriend, J. Brigham-Grette, G. H. Miller, and references therein. Supplementary Materials Quat. Res. 45, 176 (1996). 56. M. New, D. Lister, M. Hulme, I. Makin, Clim. Res. 21, www.sciencemag.org/cgi/content/full/science.1222135/DC1 Materials and Methods 38. S. A. Elias, J. V. Matthews Jr., Can. J. Earth Sci. 39, 1 (2002). 911 (2002). Supplementary Text Figs. S1 to S6 39. S. Funder et al., Bull. Geol. Soc. Den. 48, 177 (2001). Acknowledgments: Drilling operations were funded by the Tables S1 to S5 40. Q. Z. Yin, A. Berger, Clim. Dyn. 36, 1 (2011). ICDP, the NSF, the German Federal Ministry of Education and 41. B. Hönisch, N. G. Hemming, D. Archer, M. Siddall, Research (BMBF), Alfred Wegener Institute and Helmholtz References (57–106) on July 19, 2012 J. F. McManus, Science 324, 1551 (2009). Centre Potsdam (GFZ), the Russian Academy of Sciences Far 19 March 2012; accepted 1 June 2012 42. R. M. DeConto, D. Pollard, D. Kowalewski, Global Planet. East Branch, the Russian Foundation for Basic Research, and Published online 21 June 2012; Change 88-89, 45 (2012). the Austrian Federal Ministry of Science and Research. The 10.1126/science.1222135 REPORTS range to image the full set of relevant dopants. www.sciencemag.org Imaging the Impact of Single We present atomically resolved STM spectros- copy on Bi 2+y Sr 2–y CaCu 2 O 8+x (BSCCO) that Oxygen Atoms on Superconducting doubles the energy range of previous work (8) and demonstrates the impact of single dopant atoms on electronic states. Bi 2+y Sr 2–y CaCu O of states shows energy variation of ~100%, on In BSCCO, a prominent gap in the density 2 8+x a 2- to 3-nm length scale, across a wide range of 3 2 1 2 2 1 Ilija Zeljkovic, Zhijun Xu, Jinsheng Wen, Genda Gu, Robert S. Markiewicz, Jennifer E. Hoffman * Downloaded from doping (9). Recent studies (10, 11) suggest that this spectral inhomogeneity results primarily from High-temperature cuprate superconductors display unexpected nanoscale inhomogeneity in variations in the pseudogap (PG), a depression in essential properties such as pseudogap energy, Fermi surface, and even superconducting critical the density of states near the Fermi level e F that temperature. Theoretical explanations for this inhomogeneity have ranged from chemical disorder persists far above the superconducting T c , most to spontaneous electronic phase separation. We extend the energy range of scanning tunneling dominantly on the underdoped side of the phase spectroscopy on Bi 2+y Sr 2–y CaCu 2 O 8+x , allowing a complete mapping of two types of interstitial diagram. Nanoscale variation in the PG offers the oxygen dopants and vacancies at the apical oxygen site. We show that the nanoscale spatial opportunity to uncover the variable(s) determining variations in the pseudogap states are correlated with disorder in these dopant concentrations, its local strength, thus setting the stage for con- particularly that of apical oxygen vacancies. trol of this mysterious phase, whose energy scale anticorrelates with superconductivity (12, 13) any of today’s prominent materials, such chemical inhomogeneity. To fully understand and which many believe is a competitor to su- as high–transition temperature (T c )su- and harness these materials, and to drive the dis- perconductivity (11, 14). Mperconductors, doped semiconductors, covery of new materials, it is crucial to understand Several previous studies searched for chem- and colossal magnetoresistance materials, are the impact of single atoms on these inhomoge- ical origins of spectral gap variation in BSCCO. nonstoichiometric and electronically inhomo- neous electronic states. Kinoda et al. observed atomic defects in Pb-doped geneous at the nanoscale (1). Both desirable and Cuprate superconductors are quasi–two- BSCCO with resonance energy +1.7 eV, identi- undesirable electronic properties may arise from dimensional (2D) materials that arise from off- fied as Bi 3+ substitutions at the Sr 2+ site (15). stoichiometric doping of a Mott insulator by These defects may correlate weakly with regions 1 Department of Physics, Harvard University, 17 Oxford Street, oxygen intercalation and/or cation substitu- of large PG, but no quantitative analysis was 2 Cambridge, MA 02138, USA. Condensed Matter Physics and tion. Nanoscale electronic inhomogeneity in the presented (16). McElroy et al. observed localized Materials Science Department, Brookhaven National Labora- cuprates has been predicted (2, 3) and detected conductance signatures in BSCCO at –0.96V, 3 tory, Upton, NY 11973–5000, USA. Department of Physics, Northeastern University, 360 Huntington Avenue, Boston, MA (4–7), but the roles of spontaneous electronic identified as interstitial oxygen dopants (8), but 02115, USA. phase separation (2) and chemical disorder (3) several mysteries were left unresolved. First, the *To whom correspondence should be addressed: jhoffman@ are unresolved. Scanning tunneling microscopy counted oxygen fell ~50% short of the expected physics.harvard.edu (STM) has thus far lacked the necessary energy density for each T c (17). Second, the correlation 320 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

Imaging the Impact of Single Oxygen Atoms on Superconducting Bi 2+y Sr 2 − y CaCu 2 O 8+x Ilija Zeljkovic et al. Science 337 , 320 (2012); DOI: 10.1126/science.1218648 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/320.full.html Supporting Online Material can be found at: http://www.sciencemag.org/content/suppl/2012/07/19/337.6092.320.DC1.html www.sciencemag.org This article cites 51 articles , 8 of which can be accessed free: http://www.sciencemag.org/content/337/6092/320.full.html#ref-list-1 This article appears in the following subject collections: Downloaded from Physics http://www.sciencemag.org/cgi/collection/physics Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

20. P. U. Clark et al., Quat. Sci. Rev. 25, 3150 (2006). 43. B. L. Otto-Bliesner et al., Science 311, 1751 (2006). Russian Global Lake Drilling 800 drilling system was developed 21. D. S. Kaufman, J. Brigham-Grette, Quat. Sci. Rev. 12, 44. S. J. Koenig, R. M. DeConto, D. Pollard, Clim. Dyn. 37, and operated by DOSECC. We thank all participants of the 21 (1993). 1247 (2011). expedition for their engagement during recovery of the ICDP 22. G. H. Miller et al., Quat. Sci. Rev. 29, 1679 (2010). 45. D. Pollard, R. M. DeConto, Nature 458, 329 (2009). 5011-1 cores. Funding of core analyses was provided by BMBF 23. E. J. Colville et al., Science 333, 620 (2011). 46. R. McKay et al., Quat. Sci. Rev. 34, 93 (2012). (grant 03G0642), German Research Foundation (Deutsche 24. E. Willerslev et al., Science 317, 111 (2007). 47. A. Foldvik, J. Geophys. Res. 109, C02015 (2004). Forschungsgemeinschaft, grants ME 1169/21 and ME 1169/24), 25. M. E. Raymo, J. X. Mitrovica, Nature 483, 453 (2012). 48. I. R. Hall, I. N. McCave, N. J. Shackleton, G. P. Weedon, the NSF (grant 0602471), Vetenskapsrådet, The Swedish Research 26. A. de Vernal, C. Hillaire-Marcel, Science 320, 1622 (2008). S. E. Harris, Nature 412, 809 (2001). Council, the Kempe Foundation, and the Civilian Research and 27. G. A. Milne, J. X. Mitrovica, Quat.Sci.Rev. 27, 2292 (2008). 49. E. D. Galbraith et al., Nature 449, 890 (2007). Development Foundation (grant RUG1-2987-MA-10). We thank 28. P. E. Tarasov et al., Earth Sci. Rev. 108, 64 (2011). 50. G. H. Haug, M. A. Maslin, M. Sarnthein, R. Stax, N. Mantke, A. Shahnazarian, and numerous students (Univ. 29. A. A. Prokopenko et al., Clim. Past 6, 31 (2010). R. Tiedemann, Proc. ODP Sci. Results 145, 293 of Cologne) for their competent help in core processing; 30. R. Stein, J. Hefter, J. Grützner, A. Voelker, B. D. A. Naafs, (1995). T. Matrosova for contributing modern surface pollen data; and Paleoceanography 24, PA2203 (2009). 51. S. L. Jaccard, E. D. Galbraith, D. M. Sigman, G. H. Haug, R. McKay for providing a modified lithological log of the ANDRILL 31. E. Gischler, R. N. Ginsburg, J. O. Herrle, S. Prasad, Quat. Sci. Rev. 29, 206 (2010). 1B record. We also thank four anonymous reviewers for Sedimentology 57, 1049 (2010). 52. C. J. Mock, P. J. Bartlein, P. A. Anderson, Int. J. Climatol. supportive and constructive comments that greatly improved 32. P. J. Fawcett et al., Nature 470, 518 (2011). 18, 1085 (1998). the manuscript. The data reported in this paper are available in 33. R. P. Scherer et al., Geophys. Res. Lett. 35, L03505 (2008). 53. J. P. Helmke, H. A. Bauch, Paleoceanography 18, 1036 the databases of PANGAEA via www.pangaea.de (DOI: 10.1594/ 34. P. Maiorano,M.Marino, J.-A.Flores, Mar. Micropaleontol. (2003). PANGAEA.783305) and of the U.S. National Climatic Data 71, 166 (2009). 54. R. A. Woodgate, T. Weingartner, R. Lindsay, Geophys. Center via www.ncdc.noaa.gov/paleo/paleolim/paleolim_data.html. 35. T. Naish et al., Nature 458, 322 (2009). Res. Lett. 37, L01602 (2010). 36. J. Brigham-Grette, L. D. Carter, Arctic 43, 74 (1992). 55. L. Polyak, M. Jakobsson, Oceanography 24, 52 (2011) 37. G. A. Goodfriend, J. Brigham-Grette, G. H. Miller, and references therein. Supplementary Materials Quat. Res. 45, 176 (1996). 56. M. New, D. Lister, M. Hulme, I. Makin, Clim. Res. 21, www.sciencemag.org/cgi/content/full/science.1222135/DC1 Materials and Methods 38. S. A. Elias, J. V. Matthews Jr., Can. J. Earth Sci. 39, 1 (2002). 911 (2002). Supplementary Text Figs. S1 to S6 39. S. Funder et al., Bull. Geol. Soc. Den. 48, 177 (2001). Acknowledgments: Drilling operations were funded by the Tables S1 to S5 40. Q. Z. Yin, A. Berger, Clim. Dyn. 36, 1 (2011). ICDP, the NSF, the German Federal Ministry of Education and 41. B. Hönisch, N. G. Hemming, D. Archer, M. Siddall, Research (BMBF), Alfred Wegener Institute and Helmholtz References (57–106) on July 19, 2012 J. F. McManus, Science 324, 1551 (2009). Centre Potsdam (GFZ), the Russian Academy of Sciences Far 19 March 2012; accepted 1 June 2012 42. R. M. DeConto, D. Pollard, D. Kowalewski, Global Planet. East Branch, the Russian Foundation for Basic Research, and Published online 21 June 2012; Change 88-89, 45 (2012). the Austrian Federal Ministry of Science and Research. The 10.1126/science.1222135 REPORTS range to image the full set of relevant dopants. www.sciencemag.org Imaging the Impact of Single We present atomically resolved STM spectros- copy on Bi 2+y Sr 2–y CaCu 2 O 8+x (BSCCO) that Oxygen Atoms on Superconducting doubles the energy range of previous work (8) and demonstrates the impact of single dopant atoms on electronic states. Bi 2+y Sr 2–y CaCu O of states shows energy variation of ~100%, on In BSCCO, a prominent gap in the density 2 8+x a 2- to 3-nm length scale, across a wide range of 3 2 1 2 2 1 Ilija Zeljkovic, Zhijun Xu, Jinsheng Wen, Genda Gu, Robert S. Markiewicz, Jennifer E. Hoffman * Downloaded from doping (9). Recent studies (10, 11) suggest that this spectral inhomogeneity results primarily from High-temperature cuprate superconductors display unexpected nanoscale inhomogeneity in variations in the pseudogap (PG), a depression in essential properties such as pseudogap energy, Fermi surface, and even superconducting critical the density of states near the Fermi level e F that temperature. Theoretical explanations for this inhomogeneity have ranged from chemical disorder persists far above the superconducting T c , most to spontaneous electronic phase separation. We extend the energy range of scanning tunneling dominantly on the underdoped side of the phase spectroscopy on Bi 2+y Sr 2–y CaCu 2 O 8+x , allowing a complete mapping of two types of interstitial diagram. Nanoscale variation in the PG offers the oxygen dopants and vacancies at the apical oxygen site. We show that the nanoscale spatial opportunity to uncover the variable(s) determining variations in the pseudogap states are correlated with disorder in these dopant concentrations, its local strength, thus setting the stage for con- particularly that of apical oxygen vacancies. trol of this mysterious phase, whose energy scale anticorrelates with superconductivity (12, 13) any of today’s prominent materials, such chemical inhomogeneity. To fully understand and which many believe is a competitor to su- as high–transition temperature (T c )su- and harness these materials, and to drive the dis- perconductivity (11, 14). Mperconductors, doped semiconductors, covery of new materials, it is crucial to understand Several previous studies searched for chem- and colossal magnetoresistance materials, are the impact of single atoms on these inhomoge- ical origins of spectral gap variation in BSCCO. nonstoichiometric and electronically inhomo- neous electronic states. Kinoda et al. observed atomic defects in Pb-doped geneous at the nanoscale (1). Both desirable and Cuprate superconductors are quasi–two- BSCCO with resonance energy +1.7 eV, identi- undesirable electronic properties may arise from dimensional (2D) materials that arise from off- fied as Bi 3+ substitutions at the Sr 2+ site (15). stoichiometric doping of a Mott insulator by These defects may correlate weakly with regions 1 Department of Physics, Harvard University, 17 Oxford Street, oxygen intercalation and/or cation substitu- of large PG, but no quantitative analysis was 2 Cambridge, MA 02138, USA. Condensed Matter Physics and tion. Nanoscale electronic inhomogeneity in the presented (16). McElroy et al. observed localized Materials Science Department, Brookhaven National Labora- cuprates has been predicted (2, 3) and detected conductance signatures in BSCCO at –0.96V, 3 tory, Upton, NY 11973–5000, USA. Department of Physics, Northeastern University, 360 Huntington Avenue, Boston, MA (4–7), but the roles of spontaneous electronic identified as interstitial oxygen dopants (8), but 02115, USA. phase separation (2) and chemical disorder (3) several mysteries were left unresolved. First, the *To whom correspondence should be addressed: jhoffman@ are unresolved. Scanning tunneling microscopy counted oxygen fell ~50% short of the expected physics.harvard.edu (STM) has thus far lacked the necessary energy density for each T c (17). Second, the correlation 320 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS between oxygen locations and PG was weak and section I.) However, the type-A oxygens were cal density of states at –1 eV, acquired over the same of unexpected sign: Although it is well established predicted to have resonances even farther below area as Fig. 1A, containing atomic-scale features that globally increasing oxygen content corre- e F than the –0.96-V type Bs. This experimental of similar form and concentration to the previously lates with decreasing PG (12, 13), McElroy found challenge has prevented their observation to date. observed type-B oxygen dopants (8). Figure 1C → that the oxygen dopant locations were weakly We set out to observe the high-energy dopants extends the energy range down to show gðr; V ¼ correlated with regions of increased gap. and elucidate their role in cuprate inhomogeneity. 1:5VÞ, resolving a second set of atomic-scale To resolve the latter issue, one carefully tuned BSCCO was grown with a floating zone tech- features presumed to be the predicted type-A inter- theory suggested that the dominant local effect nique. As is common to facilitate the growth stitial oxygen dopants (20). Figure 1D extends the → of interstitial oxygens may be strain, which in- process, the starting composition for all crystals energy range up to show gðr; V ¼þ1VÞ, reveal- creases the local pairing strength (18, 19). A studied here contained Bi:Sr in a ratio 2.1:1.9. ing a third set of atomic-scale features. We observe → second proposal by Zhou et al.(20) explained By some combination of cation substitution and no other distinct atomic-scale features in gðr; VÞ both the missing oxygen and the counterintuitive oxygen intercalation, as-grown crystals tended images at biases between –2Vand +1.6 V. correlation by postulating the existence of two to be optimally doped with T c ~91K. Some To clarify the identity of these three features, types of interstitial oxygen dopants: “type-B” ox- samples were then annealed at 550°C in vacuum we repeat the same measurements on four dif- ygens observed by McElroy, which live around to remove oxygen and reduce T c toward the un- ferent samples and show the density of each type the BiO plane and contribute only delocalized derdoped side of the phase diagram. We present of dopant versus T c in Fig. 1E. As expected, we charge, and “type-A” oxygens, which live around STM data from four samples with decreasing observe a monotonic decrease in the number the SrO plane and have an immediate electro- oxygen content and T c = 91, 82, 68, and 55 K. The of both types of interstitial oxygen dopants with static effect, locally hole-doping the adjacent samples were cleaved in cryogenic ultra-high vac- falling T c .However,the +1 V features increase CuO 2 layer. Therefore, Zhou et al. predicted that uum and immediately inserted into our home- from virtually zero for optimally doped BSCCO the type-A oxygens would follow the expected built STM, where they were imaged at T =6 K. to almost 1% per CuO 2 plaquette for the T c =55K global trend: a strong anticorrelation with the Figure 1A shows a topography of the BiO sample. Previous studies detected no change in PG. The weak correlation observed between surface of BSCCO with T c =55K, demonstrat- cation concentration on annealing up to 840°C on July 19, 2012 McElroy’s type-B oxygens and the PG could be ing atomic resolution at +1 V bias. Figure 1B (21), so our observed increasing concentration explained as a side effect of slight repulsion be- shows the local differential tunneling conductance after 550°C vacuum annealing suggests oxygen- → tween type-B and type-A oxygens. (See addi- gðr; V ¼1VÞ (where r is the local position and site vacancies instead. There are three distinct tional discussion in supplementary materials, V is voltage), approximately proportional to the lo- oxygen lattice sites, but because the +1 V features A B E 5 -1V (type-B oxygen) www.sciencemag.org -1.5V (type-A oxygen) 4 +1V (apical O vacancy) % per CuO 2 3 2 1 Downloaded from 0 50 55 60 65 70 75 80 85 90 95 5 nm 5 nm T [K] C Low High Low High C D F type-B oxygen BiO 80 dI/dV [arb. units.] 60 background Ca 2 SrO type-A oxygen CuO 2 apical O vacancy CuO 40 SrO BiO 20 0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 5 nm 5 nm Bias [V] Fig. 1. (A) Atomically resolved topographic image of BSCCO with T c = the four samples studied. (F) dI/dV spectra over each dopant type, from 55 K, acquired at V sample = +1 V and setpoint current I = 150 pA over a the T c = 82 K sample. Each colored spectrum is the average of all imaged 35-nm area (inset, 3× magnification). (B to D) Differential conductance dopants of its type within a 30-nm field of view; the black trace shows (dI/dV) images in the same area as (A), acquired at V sample = –1, –1.5, the average spectrum far from dopants. (Inset) Top half of the BSCCO and +1 V, respectively. (E) Density of each type of dopant imaged within unit cell. SCIENCE 321 www.sciencemag.org VOL 337 20 JULY 2012

REPORTS are laterally aligned with the Bi site, we identify ground, are shown in Fig. 1F. The defining fea- cancy, respectively, are robust for different set- them as missing apical oxygen atoms in the SrO tures used to identify the impurities, such as the up conditions, locations, and samples (fig. S4). layer (see supplementary materials, section II). peak around –1 V for a type-B oxygen, and the We investigated the relationship of these dop- Typical differential conductance spectra at sudden changes of the slope around –1.2 V and ants to the inhomogeneous PG. Figure 2, A to C each of the three types of dopants, and the back- +0.8 V for type-A oxygen and apical oxygen va- shows the three types of dopants (two interstitials and one vacancy), superimposed on a PG map acquired for the T c = 55 K sample (see sup- plementary materials, section IV for details of PG determination). Any correlation between type-B oxygens and PG is too weak to be seen by the naked eye. In contrast to Zhou’s pre- diction, the type-A oxygens are correlated with regions of large PG. One immediately sees the strongest correlation between apical oxygen va- cancies and PG. Figure 2D shows the cross- correlations between the positions of the dopant atoms and the corresponding gap map. Figure 2E shows the average PG as a function of distance from the nearest dopant. This behavior is consist- ent across the three underdoped samples studied (see supplementary materials, section V). We examined more carefully the surprising departure from Zhou’spredictioninFig.3,which on July 19, 2012 plots the local PG energy versus the local con- centration of each dopant type. For the type-B interstitial oxygens in Fig. 3A, the local trend within each sample shows correlation between increased local oxygen concentration and in- creased local PG; the global trend between samples shows anticorrelation between global oxygen concentration and global PG. In con- www.sciencemag.org trast to Zhou’s prediction, Fig. 3B shows that a qualitatively similar juxtaposition of local and global trends holds for the type-A interstitial ox- ygens in the three underdoped samples. How- ever, Fig. 3C shows excellent alignment of local and global trends for the correlation between apical oxygen vacancy density and PG. This sug- gests that, in underdoped samples, variations in Downloaded from the local hole concentration, and thus the local PG, are governed primarily by the local removal Fig. 2. Comparison of dopant locations and PG map. Locations of (A) type-B interstitial oxygens (green circles), of holes by apical oxygen vacancies, rather than (B) type-A interstitial oxygens (red circles), and (C) apical oxygen vacancies (blue circles) superimposed on top of the strain or the donation of localized holes the PG map of the T c =55Ksample. (D) Cross-correlation coefficient relating the PG to the distance to the from the interstitial oxygens. nearest dopant. (E) Average PG versus distance from the nearest dopant of each type. In (D) and (E), green, red, The interstitial oxygens, previously predicted and blue lines represent type-B oxygens, type-A oxygens, and apical oxygen vacancies, respectively. (3) to determine the local gap, do contribute de- A 150 type-B oxygen T =55K B 150 type-A oxygen T c =55K C 150 apical O vacancy T =55K c c T =68K T =68K T c =68K c c T =82K c T =82K 100 T c =82K 100 c Local gap [meV] 50 Local gap [meV] 50 Local gap [meV] 50 T c =55K T =91K T c =91K T =91K c c 100 T c =68K T c =82K 2.5 3.0 3.5 4.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 Local % per CuO 2 Local % per CuO 2 Local % per CuO 2 Fig. 3. Local PG versus local dopant density. For each sample, we create a local of type-B interstitial oxygen dopants, from the four different samples used in density map of the number of a given dopant type per region of radiusx ~2nm. this study. (B) Local PG versus local density of type-A interstitial oxygen dopants, The values in the density map are then binned based on the PG value at the from the four different samples used in this study. (C)Local PG versus local corresponding pixel. The average PG value in each bin is then plotted versus the density of apical oxygen vacancies, from the three underdoped samples used in average dopant density value within each bin. (A)Local PG versus localdensity this study. Trend lines in (A) to (C) are guides to the eye. 322 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS localized holes, but in underdoped samples even samples by reducing the local hole concentration, Our high-bias atomically resolved STM re- the type-A interstitial oxygen positions remain which locally strengthens the PG, tying up anti- sults suggest a possible route to increase T c in correlated to regions of decreased local hole con- nodal states (11, 14), and locally decreasing the BSCCO: underdope to increase the pairing po- centration, likely due to a tendency to remain close Fermi level carriers that would otherwise be tential (31) but explore different annealing re- to the apical oxygen vacancies. In the optimally available for coherent pairing. Calculations also cipes to allow interstitial oxygen removal without doped sample, the apical oxygen vacancy concen- showed a correlation between T c and apical creating apical oxygen vacancies, the defects most tration is negligible, so the local hole concentration oxygen height and emphasized the importance favorable to the PG and possibly competitive to of the axial orbital for the hopping and phase superconductivity. Recently, a ~15% increase in may be set by the next-closest dopant to the CuO 2 plane, the type-A interstitial oxygen, as predicted coherence necessary for superconductivity (24, 25). maximum T c has been predicted from alternate by Zhou et al.(20). Indeed, Fig. 3B shows that the Thus, apical oxygen vacancies may lower the lo- dopant arrangements (32). Our spectroscopy local PG decreases with increasing local type-A cal superconducting T c and/or critical current J c . methods may also prove useful in future studies concentration in the T c =91 K sample. We there- Another manifestation of electronic inhomo- of the effects of single-atom impurities in other fore reconcile the apparent contradiction between geneity is a disordered “checkerboard” (CB) mod- fragile quasi-2D materials. the global (12, 13) and local (7, 8) doping trends ulation of the spectral weight that is static but for both types of interstitial oxygens. mostnoticeableatenergieswithinandnearthePG References and Notes Does dopant inhomogeneity cause PG in- energy. Field (26), doping (9), and temperature 1. E. Dagotto, Science 309, 257 (2005). homogeneity or merely pin intrinsic inhomogene- (27) dependence suggest that the CB is in fact the 2. V. J. Emery, S. A. Kivelson, O. Zachar, Phys. Rev. B 56, 6120 (1997). ity, caused perhaps by strong correlations? Because electronic ordered phase associated with the PG. 3. I. Martin, A. Balatsky, Physica C 357-360, 46 (2001). the 2- to 3-nm length scale of PG inhomogene- The CB wave vector tracks the antinodal nesting 4. S.-H. Pan et al., Nature 413, 282 (2001). ity remains the same across a wide range of dop- wave vector across a wide range of doping (28). A 5. K. M. Lang et al., Nature 415, 412 (2002). 6. K. K. Gomes et al., Nature 447, 569 (2007). ant concentrations (fig. S7A, inset), intrinsic PG closely related disordered periodic inhomogeneity 7. W. D. Wise et al., Nat. Phys. 5, 213 (2009). inhomogeneity seems plausible. However, our arises from elastic scattering between degenerate 8. K. McElroy et al., Science 309, 1048 (2005). experiment shows clearly that the dopant loca- states; this dispersing quasiparticle interference 9. K. McElroy et al., Phys. Rev. Lett. 94, 197005 (2005). on July 19, 2012 tions that are fixed at high T, particularly the ap- (QPI) may exist at similar wave vectors to the 10. M. C. Boyer et al., Nat. Phys. 3, 802 (2007). 11. A. Pushp et al., Science 324, 1689 (2009). ical oxygen vacancies, do govern the local PG static CB, but only at a limited range of energies 12. G. Deutscher, Nature 397, 410 (1999). strength that is subsequently determined on cool- within the superconducting gap (29). 13. S. Hüfner, M. A. Hossain, A. Damascelli, G. A. Sawatzky, ing through the PG transition temperature T*. It was previously claimed that type-B oxygen Rep. Prog. Phys. 71, 062501 (2008). 14. T. Kondo, R. Khasanov, T. Takeuchi, J. Schmalian, Our experiment directly measures the strong dopants are found in the minima of the QPI pat- A. Kaminski, Nature 457, 296 (2009). correlation between apical oxygen vacancies and terns (8), at both positive and negative energies. 15. G. Kinoda et al., Phys. Rev. B 71, 020502 (2005). the PG, but is there a relationship between apical However, QPI has opposite spatial phase for 16. G. Kinoda, T. Hasegawa, Phys.Rev.B 67, 224509 (2003). oxygens and superconductivity itself? In fact, filled and empty states (30), suggesting that the 17. M. Presland, J. Tallon, R. Buckley, R. Liu, N. Flower, www.sciencemag.org Physica C 176, 95 (1991). apical oxygen effects may be credited for the dis- previously observed correlation (8) relates instead 18. T. S. Nunner, B. M. Andersen, A. Melikyan, P. J. Hirschfeld, covery of high-T c superconductivity, as Müller to the CB. Figure 4A shows the distribution of Phys. Rev. Lett. 95, 177003 (2005). was originally driven to explore the LaBaCuO sys- the three types of observed atomic dopants on 19. W. Chen, M. Gabay, P. J. Hirschfeld, New J. Phys. 14, → tem by the expectation of strong electron-phonon top of a filtered g(r; V ¼ 36 mV) image (figs. 033004 (2012). 20. S. Zhou, H. Ding, Z. Wang, Phys. Rev. Lett. 98, 076401 S10 to S12). A histogram of the distance from coupling due to the Jahn-Teller effect in the CuO 6 (2007). octahedral environment, specifically the displace- each dopant to the center of the nearest “checker” 21. D. B. Mitzi, L. W. Lombardo, A. Kapitulnik, S. S. Laderman, ment of the Cu-apical-O bond (22). However, sub- (Fig. 4B) demonstrates the strong tendency of R. D. Jacowitz, Phys. Rev. B Condens. Matter 41, 6564 (1990). 22. J. G. Bednorz, K. A. Müller, Zeitschrift für Physik B 64, sequent isotope effect measurements suggested apical oxygen vacancies to lie in the peaks of the 189 (1986). Downloaded from that any phonon contribution to superconductivity imaged CB and weaker tendency for interstitial 23. D. Zech et al., Nature 371, 681 (1994). is dominated by CuO 2 plane oxygen (23). We con- oxygen dopants to lie in the troughs. We there- 24. Y. Ohta, T. Tohyama, S. Maekawa, Phys. Rev. B Condens. jecture that the apical oxygen vacancies influence fore conclude that the apical oxygen vacancies Matter 43, 2968 (1991). 25. E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen, the superconductivity indirectly in underdoped play the primary role in pinning the CB. O. K. Andersen, Phys. Rev. Lett. 87, 047003 (2001). 26. J. E. Hoffman et al., Science 295, 466 (2002). 27. C. V. Parker et al., Nature 468, 677 (2010). 28. Y. Kohsaka et al., Nature 454, 1072 (2008). 29. J. E. Hoffman et al., Science 297, 1148 (2002). 30. T. Hanaguri et al., Nat. Phys. 3, 865 (2007). 31. E. Berg, D. Orgad, S. Kivelson, Phys. Rev. B 78, 094509 (2008). 32. L. Goren, E. Altman, Phys. Rev. B 84, 094508 (2011). Acknowledgments: We thank A. Bansil, I. Bozovic, S. Davis, E. Hudson, P. Hirschfeld, A. Kapitulnik, D.-H. Lee, and J. Nieminen for useful discussions. J.E.H. acknowledges support from the NSF CAREER grant DMR-0847433 and AFOSR PECASE grant FA9550-06-1-0531. G.D.G. acknowledges support from the U.S. Department of Energy (DOE) under contract DE-AC02-98CH10886. R.S.M. acknowledges support from the DOE under contract DE-FG02-07ER46352. Supplementary Materials www.sciencemag.org/cgi/content/full/337/6092/320/DC1 Materials and Methods Supplementary Text Fig. 4. Comparison of dopant locations and the checkerboard. (A) Fourier-filtered dI/dV image of the T c = Figs. S1 to S12 55 K sample at +36 mV, showing a clear CB (setup: V sample = –150 mV; I = 800 pA). Type-B oxygens, type- References (33–52) A oxygens, and apical oxygen vacancies are superimposed as green, red, and blue circles, respectively. (B) 3 January 2012; accepted 29 May 2012 Dopant density of each type versus distance from the center of the nearest CB maximum. 10.1126/science.1218648 SCIENCE 323 www.sciencemag.org VOL 337 20 JULY 2012

Spin-Transistor Action via Tunable Landau-Zener Transitions C. Betthausen et al. Science 337 , 324 (2012); DOI: 10.1126/science.1221350 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/324.full.html Supporting Online Material can be found at: http://www.sciencemag.org/content/suppl/2012/07/18/337.6092.324.DC1.html www.sciencemag.org A list of selected additional articles on the Science Web sites related to this article can be found at: http://www.sciencemag.org/content/337/6092/324.full.html#related This article cites 36 articles , 5 of which can be accessed free: http://www.sciencemag.org/content/337/6092/324.full.html#ref-list-1 This article has been cited by 1 articles hosted by HighWire Press; see: Downloaded from http://www.sciencemag.org/content/337/6092/324.full.html#related-urls This article appears in the following subject collections: Physics, Applied http://www.sciencemag.org/cgi/collection/app_physics Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

REPORTS Spin-Transistor Action via Tunable The stripes are premagnetized in a field of 8 T, which gives rise to a remnant magnetization of 6 0.7 × 10 A/m, determined by superconducting Landau-Zener Transitions quantum interference device (SQUID) measure- ments. In the plane of the 2DEG, the stray field 3 1 2 2 C. Betthausen, T. Dollinger, H. Saarikoski, V. Kolkovsky, G. Karczewski, 3 of the periodic grating of magnets is approxi- 1 2 3 T. Wojtowicz, K. Richter, D. Weiss * mately helical with an amplitude of ~50 mT (Fig. 2B). The magnetic field texture translates to a helical spin polarization via the giant Zeeman Spin-transistor designs relying on spin-orbit interaction suffer from low signal levels resulting coupling. In our devices, we use the beating from low spin-injection efficiency and fast spin decay. Here, we present an alternative approach in which spin information is protected by propagating this information adiabatically. We pattern of the Shubnikov–de Haas oscillations at low fields (Fig. 2E) to estimate a Zeeman split- demonstrate the validity of our approach in a cadmium manganese telluride diluted magnetic semiconductor quantum well structure in which efficient spin transport is observed over device ting of 1 meV for a 50-mT stray field at 25 mK (18, 19). The Zeeman energy E Z is the largest spin- distances of 50 micrometers. The device is turned “off” by introducing diabatic Landau-Zener transitions that lead to a backscattering of spins, which are controlled by a combination of a dependent energy scale; Rashba and Dresselhaus spin-orbit splittings are more than one order of helical and a homogeneous magnetic field. In contrast to other spin-transistor designs, we find magnitude smaller (20). Spin polarization p = that our concept is tolerant against disorder. (n ↓ – n ↑ )/(n ↑ + n ↓ ) ≈ E Z /2E F in the 2DEG ranges from 8 to 15% in a 50-mT stray field for our he use of electron spin to store and pro- genstates. This causes diabatic Landau-Zener samples (Fig. 2, C and D). Here, n ↑,↓ denotes cess information calls for the ability to transitions between the spin eigenstates (11, 12), the density of spin-up/down electrons; E F is the Tinject, propagate, and manipulate spin with which can be used to selectively backscatter spin- Fermi energy. high efficiency (1, 2). Much of the recent research polarized charge carriers, thus giving rise to tran- To measure the longitudinal resistance r xx (B) on July 19, 2012 in the field has been motivated by attempts to sistor action. in a Hall bar geometry, the magnetic field is first realize the seminal spin transistor proposed by We demonstrate the validity of our approach ramped up to 8 T under an angle q (Fig. 2A) to Datta and Das (3). Their original concept re- in a (Cd,Mn)Te diluted magnetic semiconductor magnetize the stripes, and then data are taken on quires spin injection from ferromagnetic contacts quantum-well (QW) structure in which the giant the sweep back to –0.5 T. In the presence of a into a narrow channel of a two-dimensional elec- Zeeman splitting, attributed to s-d exchange inter- helical stray field, r xx (B) shows distinct peaks at tron gas (2DEG), where transport is ballistic and action between electronic states and the localized about T50 mT (Fig. 2, E and F) for the studied spins precess in a gate-controlled spin-orbit field Mn spins (13), gives rise to intrinsic spin polariza- devices A, B, and C (fig. S1). The samples differ into “on” or “off” orientations. Even though such tion. We modulate the diabatic spin-backscattering in QW thickness, carrier density, electron mo- www.sciencemag.org spin manipulation in a spin-orbit field has been probability by a combination of a spatially rotat- bility, and Mn concentration (table S1), but the demonstrated in a nonlocal measurement (4), ing magnetic field B s and a homogenous mag- results are independent of such details. The peaks signal levels remain small as a result of low spin- netic field B. B s is created by premagnetized do not appear in unmodulated reference sam- injection efficiency and limited spin lifetime. The ferromagnetic stripes placed above the device ples. Apart from a small dip near B =0,rem- latter can be enhanced if spin polarization is pro- (Fig. 2, A and D). In the absence of the ho- iniscent of weak antilocalization (black traces in tected against scattering processes by an SU(2) mogeneous field, eigenfunctions change slowly, Fig. 2F), these samples do not show any notable symmetry that emerges when the Rashba spin- leading to adiabatic spin transport and a heli- orbit interaction strength is equal to the Dressel- cal spin wave resulting from U(1) symmetry. If Downloaded from haus term (5, 6). This creates a persistent spin-helix scattering is spin-independent, spin orientation e 2 state (6) that has been observed in optical measure- depends only on the position along the transport ments (7). Aside from 2D systems, spin-transistor direction and not on the specific path of the action was also explored in bulk bipolar semi- spin: the device is in its on state. However, if E F conductors (8, 9). the helical field component is equal in ampli- Energy diabatic Here, we present an alternative efficient spin- tude to the homogeneous component, diabatic transistor design that uses adiabatic spin propa- transitions and backscattering of spin occur. gation to protect spin information and tunable Consequently, the device is in the off state. As adiabatic e 1 diabatic Landau-Zener transitions between spin in the original Datta-Das concept, we measure eigenstates for spin-transmission control (Fig. 1). the source-drain resistance, which intensifies as x According to the adiabatic theorem of quantum the spin-backscattering rate is increased. The mechanics (10), a spin that is initially in an ei- degree of adiabaticity is therefore monitored Fig. 1. Spin-transistor action via tunable diabatic genstate will remain in its instantaneous eigen- electrically in the channel resistance, in contrast transitions in a system of two orthogonal spin- state during transport if external perturbations act to interference measurements (14, 15)of geomet- eigenstates e 1 and e 2 (blue arrows) with energy sep- on it slowly. Spin information is hence protected rical (Berry) phases (16), associated with adia- aration G at the closest approach. A spin that is transported along a spatial coordinate x through against decay, allowing spin to propagate over batic spin evolution only. the system, starting from e 1 ,willremaininthis device distances of several micrometers. For Our setup consists of three building blocks: instantaneous eigenstate if evolution (here defined rapidly changing perturbations, spin cannot adapt (i) a high-mobility 2DEG with a giant Zeeman as rotation of eigenstate spin orientation) is slow in its state, which becomes a superposition of ei- splitting, (ii) a grating of ferromagnetic stripes the reference frame of the spin (small dx/dt). In a on top of the sample (Fig. 2, A and D), and (iii) fast evolution (large dx/dt), spin cannot adapt to 1 Department of Experimental and Applied Physics, Regensburg a homogeneous B field generated by a super- 2 University, 93040 Regensburg, Germany. Department of Theo- conducting coil. A modulation doped (Cd,Mn)Te the changing environment, leading to a diabatic retical Physics, Regensburg University, 93040 Regensburg, evolution (red arrows) corresponding to a Landau- 3 Germany. Institute of Physics, Polish Academy of Sciences, quantum-well structure (17, 18) is used to confine Zener transition to e 2 .The e 2 state lies above E F , 02668 Warsaw, Poland. electrons to two dimensions. The ferromagnetic causing wave-function decay and spin backscatter- *To whom correspondence should be addressed. E-mail: grating is made of 75-nm-thick dysprosium stripes ing. Spin-transistor action involves either tuning [email protected] patterned to periods a ranging from 0.5 to 8 mm. G or dx/dt. 324 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS features at low B.The r xx peaks in the modu- to the high mobility of our samples, let us assume For one pair of spin-polarized states (that is, for lated samples vanish rapidly with increasing that the spin transport is ballistic on the scale of fixed wave vector k y ), the total magnetic field, temperature T and are no longer recognizable the period of the Dy grating. Figure 2C shows the Zeeman-split energy levels, spin directions, and above ~0.6 K (Fig. 3A). This behavior, con- Zeeman-split subbands of a 2DEG.At E F there local energy bands are depicted in Fig. 3E (case firmed by full numerical transport calculations is an imbalance of spin states moving in the x B = B s /2) and Fig. 3F (B = B s ). Because of the in Fig. 3B (18), is expected for the spin po- direction because some of the occupied spin-split large Zeeman splitting associated with the stray larization of (Cd,Mn)Te, as the dominant s-d modes exhibit unoccupied counterparts above E F . field at B =0, the Larmor frequency w L = exchange part in the g factor contributes appre- For a Zeeman-split state pair close to E F , the upper g eff m B B s /ħ (where ħ is Planck’sconstant h divided ciably only at low T. Because of their orbital na- state (solid blue line in Fig. 2C) is unoccupied, by 2p) is higher than the frequency of stripe ture, magnetic commensurability effects persist whereas the lower state (solid red line) is occupied modulation in a ballistic flight through the device to much higher T (~40 K) (21) and can therefore and carries spin. Assuming a helical stray field w f =2pv F,x /a (where v F,x = ħk x /m* is the com- be excluded. Hence, the T dependence of the (22, 23) with amplitude B s , the total magnetic ponent of the Fermi velocity parallel to the trans- peaks confirms a spin-related effect. field, B+B s (x), gives rise to a Zeeman splitting port direction, m* is the effective mass of CdTe, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p 1 2 The r xx peaks reflect reduced spin transmission E Z,TðxÞ ¼ T g eff m B s 1 þ g þ 2gcosð2px=aÞ, and k x is the wave vector in x direction). This 2 B through the 2DEG channel, which can be intui- where g = B/B s , g eff is the (giant) effective g keeps spin transport predominantly in the adia- tively explained within the following model: Due factor of (Cd,Mn)Te, and m B is Bohr’s magneton. batic regime, as expressed by the condition Q = w L /w f >1,where Q is ameasure of thedegreeof adiabaticity in ballistic systems (16, 24, 25). For A Drain B B = B s , a degeneracy point emerges in the energy θ bands midway between adjacent stripes, because B the total magnetic field vanishes at those points (Figs. 2B and 3F). In the reference frame of the 20 μm 20 μm x * ** tates continuously as it approaches the degen- V spin, the direction of the magnetic field first ro- on July 19, 2012 eracy point adiabatically, but then the field reverses direction. This violates the conditions for adiabaticity, because an instant 180° spin- Source y z x (μm) rotation would be needed to stay on the lower C D B s band. Hence, a transition to the upper Zeeman Energy band occurs close to the degeneracy. On the up- www.sciencemag.org per band the spin is parallel to the magnetic field, Dy giving rise to a tunneling barrier because the po- tential energy rises above E F ;thisleadstoback- scattering, and spin-transport is blocked, increasing the resistance of the sample. Only spin-polarized 2DEG modes are affected by this blocking, because both spin bands of spin-compensated modes remain everywhere below E F . The probability of a diabatic transition to the Downloaded from E F upper Zeeman band is finite, even if there are no degeneracy points. Using Landau-Zener formal- ism, we show (18) that the spin-backscattering probability associated with a diabatic transition becomes increasingly probable as the distance between the Zeeman-split levels at the closest approach decreases or transport velocity v F,x in- creases (11, 12, 26). The periodic stripe grating forms a spatially repeating sequence of potential B(T) B(T) transition points that enhances the probability of a diabatic transition. Fig. 2. (A) Electron micrograph of a Hall-bar segment covered with Dy stripes. (Inset) Orientations of The actual stray-field profile of the magne- the external homogeneous B field and the magnetic moments of the premagnetized stripes (magenta tized stripes is anisotropic and only approximate- arrows), defined by the tilt angle q with respect to the z axis in the xz plane. (B) Calculated magnetic ly helical (Fig. 2B). This provides an additional stray-field components B s,x (red) and B s,z (blue) in the plane of the 2DEG for sample C with a grating of means to confirm that the r xx peaks stem from a =1 mm. B = B C,+ and B = B C,– denote the external field needed to generate points of vanishing total spin-transmission blocking. To do that, we studied magnetic field (marked by asterisks) if the external field is applied in positive or negative z direction, samples with different stripe periodicities a ranging respectively. (C) Schematic band structure of the (Cd,Mn)Te QW. Black solid lines, Zeeman-split spin-up and spin-down subbands for wave vectors in the y direction, k y = 0; dashed lines, spin-split subbands from 0.5 to 8 mm (shown for q =45° in Fig. 3C) for k y ≠ 0, but with spin-up and spin-down states occupied at E F ; colored solid lines, spin-split subbands and measured r xx (B) at different tilt angles q for for k y ≠ 0, but with an empty upper band. Only Landau-Zener transitions between filled and empty agiven value of a (shown for a =1 mminFig. subbands contribute to spin backscattering. Spin-up and spin-down densities n ↑ and n ↓ are represented 3D). We calculated the stray field of the stripes by blue and red areas, respectively. (D) In the adiabatic transport regime at B =0,electronspins in the within the dipole approximation for different 2DEG (spheres with arrows) keep anti-aligned with the stray-field B s of the stripes, and a spin-helix values of q and a, which gives the positions of forms. (E and F) r xx with and without modulation at T =25mKand q = 0° (upper curves, shifted for the spin-blocking peaks as field values where clarity). Mn ions in the QW cause a distinct beating pattern with nodes (open triangles) at higher B (E). zeros form in the total field. Our theoretical data Resistance peaks (black triangles) are associated with blocking of spin transmission. are compared with the corresponding values SCIENCE 325 www.sciencemag.org VOL 337 20 JULY 2012

REPORTS A BFE Adiabatic regime at γ = 0.5 Diabatic regime, γ = 1 2DEG Energy Energy blocking spin spin transport 0101 x (μm) x (μm) C D Energy Energy Fig. 3. Measured (A)and simulated(B) magnetoresistance traces at various T in sample C with a =1 mmand q = 0°. Calculations were done with a ballistic transport model [top part of (B)] and a disorder model with a 4-mm mean free path [bottom part of (B)]. Magnetoresistance traces at various periods a (C) and tilt angles q (D)at25mK, q = 45° in (C), and a =1 mm in (D). All curves, except the lowest ones, in (A) to (D) are shifted for clarity. (E and F)Total magnetic fieldactingonanelectron(toppanels) andZeeman-split 1 energy-bands of a polarized state (middle panels) with spin-orientations (arrows) in the (E) adiabatic transport regime at g = B/B s = / 2 and in the (F) diabatic regime g =1.The x component of magnetic field reverses direction in the center at g = 1, which leads to spin backscattering. (Bottom panels) Corresponding on July 19, 2012 local Zeeman-split energy bands. extracted for sample C in Fig. 4A (and for sam- AB ple B in fig. S3) and show the expected an- isotropy. Note that no fit parameters were used in the calculations. In the ballistic model, all spin-polarized trans- www.sciencemag.org port modes are blocked at B = B s , increasing the resistance of the sample. A single degeneracy point in the Zeeman bands is sufficient to block all spin transmission; hence, the peak height does not depend on the number of modulations, which is in line with the measurements of dif- ferent modulation periods at a fixed device length of 50 mm (Fig. 3C). Only the 8-mm modulation Fig. 4. (A) Measured and calculated anisotropy of the stray field, as reflected by the resistance peak separation Δ = B C,+ – B C,– .(B) Relative resistance peak height Δr xx /r xx as a function of temperature in Downloaded from shows a reduced peak height, which is probably due to the highly anisotropic stray field. sample C (solid and open squares for B C,+ and B C,– , respectively), in the ballistic theoretical model (solid line), and in a model that includes disorder (triangles). In the latter case, the mean free path is 4 mm, The ballistic model gives an estimate E Z (T,2B s )/ (4E F ) for the relative magnetoresistance peak and the device length is 18 mm. The stray-field modulation period is 1 mm. The corresponding reference values of the spin-compensated limit are estimated from the resistances at 1 K. Error bars denote the height Δr xx /r xx (black line in Fig. 4B). In the uncertainty in the disorder-averaged values. presence of disorder, diabatic transitions are perturbed and, hence, peak heights are reduced. However, we found that disorder does not strong- and off states. In contrast to the Datta-Das con- teraction can be tuned by means of an electric ly affect the measured peak heights in sample C. cept, our approach is tolerant against disorder. field involve, for instance, magnetic semicon- The experimental values in Fig. 4B are compa- Our proof of concept device shows source-drain ductors (27). rable to our disorder model at an electron mean resistance modulation of up to ~10%, which is free path of l e =4 mm (blue triangles in Fig. 4B). expected to be considerably enhanced in mate- References and Notes The peaks disappear in the calculations at l e = rials with high spin polarization (18). Further- 1. D. D. Awschalom, D. Loss, N. Samarth, Eds., Semiconductor 0.65 mm, which equals the measured mean free more, the homogeneous B field, provided here Spintronics and Quantum Computation (Springer, Berlin, path in sample C. Orbital effects may enhance spin by coils, might be replaced with a magnetic gate 2009). 2. I. Žutić, J. Fabian, S. Das Sarma, Rev. Mod. Phys. 76, 323 effects and contribute to a better-than-expected that is switched by spin torque, discussed in (2004). device performance at low T (18). Nevertheless, (18). This would enable the device to be controlled 3. S. Datta, B. Das, Appl. Phys. Lett. 56, 665 (1990). experiments indicate that adiabatically transported purely by electric fields. As we use a paramag- 4. H. C. Koo et al., Science 325, 1515 (2009). spin is stable over device distances of 50 mm, netic material, our device operates only at low 5. J. Schliemann, J. C. Egues, D. Loss, Phys. Rev. Lett. 90, 146801 (2003). which are much longer than electron mean free T. However, we stress that the design concepts 6. B. A. Bernevig, J. Orenstein, S.-C. Zhang, Phys. Rev. Lett. paths in the samples (between 0.47 and 1.39 mm). are not restricted to a particular choice of ma- 97, 236601 (2006). Our results demonstrate the concept of an terials, temperature, methods of spin injection, 7. J. D. Koralek et al., Nature 458, 610 (2009). adiabatic spin transistor: Adiabatic transport pro- manipulation, or detection. Transferring our con- 8. J. Fabian, I. Žutić, S. Das Sarma, Appl. Phys. Lett. 84,85 (2004). tects spin information, whereas diabatic transitions cept to higher T requires employing exchange or 9. N. Rangaraju, J. A. Peters, B. W. Wessels, Phys. Rev. Lett. lead to spin backscattering (thus depolarizing spin-orbit splitting rather than Zeeman splitting. 105, 117202 (2010). the current) and allow for switching between on Possible material systems in which exchange in- 10. M. Born, V. Fock, Z. Phys. A 51, 165 (1928). 326 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS 11. L. Landau, Phys. Sov. Union 2, 46 (1932). 23. M. Calvo, Phys. Rev. B 18, 5073 (1978). Forschungsgemeinschaft through SFB 689, WE 247618, 12. C. Zener, Proc. R. Soc. London Ser. A 137, 696 (1932). 24. O. Zaitsev, D. Frustaglia, K. Richter, Phys. Rev. B 72, FOR 1483, and Elitennetzwerk Bayern. Our research in Poland 13. J. K. Furdyna, J. Appl. Phys. 64, R29 (1988). 155325 (2005). (V.K., G.K., T.W.) was partially supported by the European 14. M. König et al., Phys. Rev. Lett. 96, 076804 (2006). 25. M. Popp, D. Frustaglia, K. Richter, Phys. Rev. B 68, Union within the European Regional Development Fund, 15. T. Koga, Y. Sekine, J. Nitta, Phys. Rev. B 74, 041302(R) 041303(R) (2003). through Innovative Economy grant POIG.01.01.02-00-008/08. (2006). 26. J. P. Davis, P. Pechukas, J. Chem. Phys. 64, 3129 (1976). 16. M. V. Berry, Proc. R. Soc. London Ser. A 392, 45 (1984). 27. H. Ohno et al., Nature 408, 944 (2000). Supplementary Materials 17. B. A. Piot et al., Phys. Rev. B 82, 081307 (2010). www.sciencemag.org/cgi/content/full/337/6092/324/DC1 18. Materials and methods are available as supplementary Acknowledgments: We thank M. Wimmer for useful Materials and Methods materials on Science Online. discussions and providing the code for the transport equation Supplementary Text 19. F. J. Teran et al., Phys. Rev. Lett. 88, 186803 (2002). solver; M. Kiessling for SQUID measurements; M. Wiater Figs. S1 to S8 20. M. Cardona, N. E. Christensen, G. Fasol, Phys. Rev. B 38, for technical assistance in molecular beam epitaxy growth; Table S1 1806 (1988). and C. Back, G. E. W. Bauer, J. Fabian, V. I. Falko, References (28–39) 21. P. H. Beton et al., Phys. Rev. B 42, 9689 (1990). C. Strunk, and G. Woltersdorf for fruitful discussions. 1 March 2012; accepted 14 June 2012 22. C. Jia, J. Berakdar, Phys. Rev. B 81, 052406 (2010). We acknowledge financial support from the Deutsche 10.1126/science.1221350 A Paramagnetic Bonding of our recently developed LONDON code, which is capable of treating molecular systems accu- Mechanism for Diatomics in rately in all field orientations. Our studies not only confirm the bonding of triplet H 2 but also provide an elementary molecular orbital (MO) Strong Magnetic Fields explanation that involves neither charge displace- ment nor dispersion: Nonbonding molecular elec- tronic states are stabilized by the reduction of the on July 19, 2012 1 1 1 Kai K. Lange, E. I. Tellgren, M. R. Hoffmann, 1,2 T. Helgaker * paramagnetic kinetic energy of antibonding MOs when these are oriented perpendicular to the mag- Elementary chemistry distinguishes two kinds of strong bonds between atoms in molecules: netic field. The generality of the proposed bonding the covalent bond, where bonding arises from valence electron pairs shared between neighboring mechanism is confirmed by calculations on He 2 , atoms, and the ionic bond, where transfer of electrons from one atom to another leads to previously not studied in strong magnetic fields. Coulombic attraction between the resulting ions. We present a third, distinct bonding mechanism: To represent the molecular electronic states in perpendicular paramagnetic bonding, generated by the stabilization of antibonding orbitals in magnetic fields, we use the full configuration- their perpendicular orientation relative to an external magnetic field. In strong fields such as interaction (FCI) method [implemented using www.sciencemag.org 5 those present in the atmospheres of white dwarfs (on the order of 10 teslas) and other stellar string-based techniques (10–12)], where the objects, our calculations suggest that this mechanism underlies the strong bonding of H 2 in the N-electron wave function is expanded linearly in S (1s g 1s*) triplet state and of He 2 in the S (1s 1s* ) singlet state, as well as their 3 þ 1 þ 2 2 Slater determinants, u u g g u preferred perpendicular orientation in the external field. jFCI〉 ¼ ∑ C n detjf , f , ::: f jð1Þ p 1n p 2n p Nn n hemical bonding mechanisms are not only the behavior of molecules under Earth-like con- well understood phenomenologically and ditions. In the absence of direct measurements whose coefficients C n are determined by the Ctheoretically, but are also accurately de- and observations, ab initio (as opposed to semi- Rayleigh-Ritz variation principle (13). Each de- Downloaded from scribed by the methods of modern quantum chem- empirical) quantum mechanical simulations play terminant is an antisymmetrized product of N or- istry. Molecular atomization energies, for example, a crucial role in unraveling the behavior of mol- thonormal spin MOs f p ; the summation is over all are today routinely calculated to an accuracy of ecules in strong magnetic fields and may be useful determinants that may be generated from a given a few kilojoules per mole—the “chemical ac- in the interpretation of white dwarf spectra (3, 4). set of MOs. The exact solution to the Schrödinger curacy” characteristic of modern measurements Over the years, many quantum chemical studies equation is reached in the limit of a complete set (1). However, nearly all our knowledge about have been performed on one- and two-electron mol- of MOs, making it possible to approach this so- chemical bonding pertains to Earth-like condi- ecules in strong magnetic fields (5). Some of these lution in a systematic manner. tions, where magnetic interactions are weak rel- demonstrate how certain otherwise unbound one- The FCI model makes no assumptions about ative to the Coulomb interactions responsible for electron molecules become bound in strong fields. the structure of the electronic system; in par- bonding. By contrast, in the atmospheres of rap- Intriguingly, Hartree-Fock calculations by Žaucer ticular, it makes no assumptions regarding the idly rotating compact stellar objects, magnetic and Ažman in 1978 (6) and by Kubo in 2007 (7) dominance of one Slater determinant (assumed fields are orders of magnitude stronger than those suggest that the otherwise dissociative lowest triplet in Hartree-Fock and coupled-cluster theories). 3 þ that can be generated in laboratories. In particu- state S ð1s g 1s*Þ of H 2 becomes bound in the This model is therefore capable of describing all u u lar, some white dwarfs have fields as strong as perpendicular orientation of the molecule relative bonding situations and dissociation processes in 10 5 10 T, and fields up to 10 Texist on neutron stars to the field. The binding has also been noted in an unbiased manner, which is essential when un- and magnetars. Under these conditions, magnet- simple model calculations and rationalized in terms familiar phenomena are studied. Equally impor- ism strongly affects the chemistry and physics of of van der Waals binding (dispersion) (8)and a tant, the FCI model provides a uniform description molecules, playing a role as important as that of shift of electronic charge density toward the mo- of different electronic states and is therefore able Coulomb interactions (2). To understand this un- lecular center (9). Bearing in mind that the un- to describe the complicated evolution of such familiar chemistry, we cannot be guided solely by correlated Hartree-Fock model often strongly states that occurs with increasing field strength. overestimates the binding energy in the absence The FCI method is a standard technique of 1 Centre for Theoretical and Computational Chemistry, Depart- of magnetic fields, these findings must be con- quantum chemistry, often used to benchmark less ment of Chemistry, University of Oslo, N-0315 Oslo, Norway. 2 firmed by more advanced quantum chemical expensive and less accurate methods, and was Chemistry Department, University of North Dakota, Grand Forks, ND 58202, USA. simulations. previously used by Schmelcher and Cederbaum *To whom correspondence should be addressed. E-mail: Here, we report highly accurate calculations in their study of H 2 in strong parallel magnetic [email protected] on H 2 in strong magnetic fields, taking advantage fields (14). Our FCI implementation differs from SCIENCE 327 www.sciencemag.org VOL 337 20 JULY 2012

A Paramagnetic Bonding Mechanism for Diatomics in Strong Magnetic Fields Kai K. Lange et al. Science 337 , 327 (2012); DOI: 10.1126/science.1219703 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/327.full.html A list of selected additional articles on the Science Web sites related to this article can be found at: www.sciencemag.org http://www.sciencemag.org/content/337/6092/327.full.html#related This article has been cited by 1 articles hosted by HighWire Press; see: http://www.sciencemag.org/content/337/6092/327.full.html#related-urls Downloaded from This article appears in the following subject collections: Chemistry http://www.sciencemag.org/cgi/collection/chemistry Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

REPORTS 11. L. Landau, Phys. Sov. Union 2, 46 (1932). 23. M. Calvo, Phys. Rev. B 18, 5073 (1978). Forschungsgemeinschaft through SFB 689, WE 247618, 12. C. Zener, Proc. R. Soc. London Ser. A 137, 696 (1932). 24. O. Zaitsev, D. Frustaglia, K. Richter, Phys. Rev. B 72, FOR 1483, and Elitennetzwerk Bayern. Our research in Poland 13. J. K. Furdyna, J. Appl. Phys. 64, R29 (1988). 155325 (2005). (V.K., G.K., T.W.) was partially supported by the European 14. M. König et al., Phys. Rev. Lett. 96, 076804 (2006). 25. M. Popp, D. Frustaglia, K. Richter, Phys. Rev. B 68, Union within the European Regional Development Fund, 15. T. Koga, Y. Sekine, J. Nitta, Phys. Rev. B 74, 041302(R) 041303(R) (2003). through Innovative Economy grant POIG.01.01.02-00-008/08. (2006). 26. J. P. Davis, P. Pechukas, J. Chem. Phys. 64, 3129 (1976). 16. M. V. Berry, Proc. R. Soc. London Ser. A 392, 45 (1984). 27. H. Ohno et al., Nature 408, 944 (2000). Supplementary Materials 17. B. A. Piot et al., Phys. Rev. B 82, 081307 (2010). www.sciencemag.org/cgi/content/full/337/6092/324/DC1 18. Materials and methods are available as supplementary Acknowledgments: We thank M. Wimmer for useful Materials and Methods materials on Science Online. discussions and providing the code for the transport equation Supplementary Text 19. F. J. Teran et al., Phys. Rev. Lett. 88, 186803 (2002). solver; M. Kiessling for SQUID measurements; M. Wiater Figs. S1 to S8 20. M. Cardona, N. E. Christensen, G. Fasol, Phys. Rev. B 38, for technical assistance in molecular beam epitaxy growth; Table S1 1806 (1988). and C. Back, G. E. W. Bauer, J. Fabian, V. I. Falko, References (28–39) 21. P. H. Beton et al., Phys. Rev. B 42, 9689 (1990). C. Strunk, and G. Woltersdorf for fruitful discussions. 1 March 2012; accepted 14 June 2012 22. C. Jia, J. Berakdar, Phys. Rev. B 81, 052406 (2010). We acknowledge financial support from the Deutsche 10.1126/science.1221350 A Paramagnetic Bonding of our recently developed LONDON code, which is capable of treating molecular systems accu- Mechanism for Diatomics in rately in all field orientations. Our studies not only confirm the bonding of triplet H 2 but also provide an elementary molecular orbital (MO) Strong Magnetic Fields explanation that involves neither charge displace- ment nor dispersion: Nonbonding molecular elec- tronic states are stabilized by the reduction of the on July 19, 2012 1 1 1 Kai K. Lange, E. I. Tellgren, M. R. Hoffmann, 1,2 T. Helgaker * paramagnetic kinetic energy of antibonding MOs when these are oriented perpendicular to the mag- Elementary chemistry distinguishes two kinds of strong bonds between atoms in molecules: netic field. The generality of the proposed bonding the covalent bond, where bonding arises from valence electron pairs shared between neighboring mechanism is confirmed by calculations on He 2 , atoms, and the ionic bond, where transfer of electrons from one atom to another leads to previously not studied in strong magnetic fields. Coulombic attraction between the resulting ions. We present a third, distinct bonding mechanism: To represent the molecular electronic states in perpendicular paramagnetic bonding, generated by the stabilization of antibonding orbitals in magnetic fields, we use the full configuration- their perpendicular orientation relative to an external magnetic field. In strong fields such as interaction (FCI) method [implemented using www.sciencemag.org 5 those present in the atmospheres of white dwarfs (on the order of 10 teslas) and other stellar string-based techniques (10–12)], where the objects, our calculations suggest that this mechanism underlies the strong bonding of H 2 in the N-electron wave function is expanded linearly in S (1s g 1s*) triplet state and of He 2 in the S (1s 1s* ) singlet state, as well as their 3 þ 1 þ 2 2 Slater determinants, u u g g u preferred perpendicular orientation in the external field. jFCI〉 ¼ ∑ C n detjf , f , ::: f jð1Þ p 1n p 2n p Nn n hemical bonding mechanisms are not only the behavior of molecules under Earth-like con- well understood phenomenologically and ditions. In the absence of direct measurements whose coefficients C n are determined by the Ctheoretically, but are also accurately de- and observations, ab initio (as opposed to semi- Rayleigh-Ritz variation principle (13). Each de- Downloaded from scribed by the methods of modern quantum chem- empirical) quantum mechanical simulations play terminant is an antisymmetrized product of N or- istry. Molecular atomization energies, for example, a crucial role in unraveling the behavior of mol- thonormal spin MOs f p ; the summation is over all are today routinely calculated to an accuracy of ecules in strong magnetic fields and may be useful determinants that may be generated from a given a few kilojoules per mole—the “chemical ac- in the interpretation of white dwarf spectra (3, 4). set of MOs. The exact solution to the Schrödinger curacy” characteristic of modern measurements Over the years, many quantum chemical studies equation is reached in the limit of a complete set (1). However, nearly all our knowledge about have been performed on one- and two-electron mol- of MOs, making it possible to approach this so- chemical bonding pertains to Earth-like condi- ecules in strong magnetic fields (5). Some of these lution in a systematic manner. tions, where magnetic interactions are weak rel- demonstrate how certain otherwise unbound one- The FCI model makes no assumptions about ative to the Coulomb interactions responsible for electron molecules become bound in strong fields. the structure of the electronic system; in par- bonding. By contrast, in the atmospheres of rap- Intriguingly, Hartree-Fock calculations by Žaucer ticular, it makes no assumptions regarding the idly rotating compact stellar objects, magnetic and Ažman in 1978 (6) and by Kubo in 2007 (7) dominance of one Slater determinant (assumed fields are orders of magnitude stronger than those suggest that the otherwise dissociative lowest triplet in Hartree-Fock and coupled-cluster theories). 3 þ that can be generated in laboratories. In particu- state S ð1s g 1s*Þ of H 2 becomes bound in the This model is therefore capable of describing all u u lar, some white dwarfs have fields as strong as perpendicular orientation of the molecule relative bonding situations and dissociation processes in 10 5 10 T, and fields up to 10 Texist on neutron stars to the field. The binding has also been noted in an unbiased manner, which is essential when un- and magnetars. Under these conditions, magnet- simple model calculations and rationalized in terms familiar phenomena are studied. Equally impor- ism strongly affects the chemistry and physics of of van der Waals binding (dispersion) (8)and a tant, the FCI model provides a uniform description molecules, playing a role as important as that of shift of electronic charge density toward the mo- of different electronic states and is therefore able Coulomb interactions (2). To understand this un- lecular center (9). Bearing in mind that the un- to describe the complicated evolution of such familiar chemistry, we cannot be guided solely by correlated Hartree-Fock model often strongly states that occurs with increasing field strength. overestimates the binding energy in the absence The FCI method is a standard technique of 1 Centre for Theoretical and Computational Chemistry, Depart- of magnetic fields, these findings must be con- quantum chemistry, often used to benchmark less ment of Chemistry, University of Oslo, N-0315 Oslo, Norway. 2 firmed by more advanced quantum chemical expensive and less accurate methods, and was Chemistry Department, University of North Dakota, Grand Forks, ND 58202, USA. simulations. previously used by Schmelcher and Cederbaum *To whom correspondence should be addressed. E-mail: Here, we report highly accurate calculations in their study of H 2 in strong parallel magnetic [email protected] on H 2 in strong magnetic fields, taking advantage fields (14). Our FCI implementation differs from SCIENCE 327 www.sciencemag.org VOL 337 20 JULY 2012

REPORTS theirs in being invariant with respect to gauge invariance is carefully imposed. For the results to position relative to the center of the Gaussian K origin and hence capable of describing all ori- be reliable, it is essential that the calculations be (here an atomic center), a > 0 is the Gaussian entations of the molecule in a field equally well. rigorously invariant with respect to gauge origin exponent, and N ijk is the normalization constant. To understand this point, we recall that the ki- in all molecular orientations. The kinetic energy These AOs depend on the field B andonthe netic energy operator (including the spin-Zeeman operator in Eq. 2 depends quadratically on p i gauge origin O in a physically reasonable man- term) in the magnetic field B is given, in atomic and therefore both linearly and quadratically ner, being correct to first order in the external units, by on B. For the field strengths considered here, magnetic field and ensuring gauge origin in- 5 on the order of B 0 =2.35 × 10 T (one atomic variance of all computed expectation values. The 1 2 T ¼ ∑ p þ B ⋅∑ s i , unit), the linear and quadratic field contributions use of such field-dependent orbitals, introduced i 2 i i ð2Þ to the Hamiltonian are equally important, result- by London in 1937 (15), is a standard technique 1 p i ¼ −i∇ i þ B ðr i − OÞ ing in a complicated chemistry in this regime. in perturbative treatments of molecular magnetic 2 To ensure gauge origin invariance, we ex- phenomena (16–19). The use of London orbitals where i is the imaginary unit, and p i and s i are the pand the MOs linearly in a set of field-dependent in nonperturbative studies is technically more kinetic momentum and spin operators of electron atom-fixed Cartesian Gaussian atomic orbitals complicated and therefore uncommon. Indeed, i, respectively. The kinetic energy operator de- (AOs)ofthe form London orbitals have previously been used only + pends parametrically on the gauge origin O,an in calculations on the one-electron H 2 molecule j k i arbitrary point in space where the field contri- c ðr,K,B,OÞ¼ N ijk x y z (20–22), on the two-electron H 2 molecule (6, 7), ijk K K K bution to the operator vanishes. In exact theory, and on larger molecules by our group (23, 24), all choices of O yield the same energy and other 1 2 all at the uncorrelated Hartree-Fock level of theory  exp iB ðO − KÞ ⋅ r expð−ar Þð3Þ properties of the system; in approximate calcu- 2 K [see also the Heitler-London model of Basile lations (except in parallel orientations with the et al.(9)]. Our gauge origin–invariant FCI code gauge origin on the molecular axis), the calcu- where r is the position of the electron relative allows us to study molecules in different elec- lated results depend on O unless gauge origin to the origin of the coordinate system, r K is its tronic states and arbitrary orientations in a reliable on July 19, 2012 AB C 5 5 www.sciencemag.org 4 4 3 3 2 2 1 1 Downloaded from DE 0 1000 B = 2.25B 0 –1000 R 50 100 150 200 250 (pm) –2000 B = 0.00 E (kJ mol -1 ) –1000 B = 1.50B 0 –3000 –4000 B = 0.75B 0 –2000 B = 0.75B 0 B = 1.50B 0 –5000 B = 0.00 B = 2.25B 0 –3000 –6000 R (pm) 0 50 100 150 200 250 300 350 Fig. 1. The H 2 molecule in an external magnetic field. (A) Schematic il- axis relative to the field (high values in red, low values in blue). (D and E) 2 1 þ 3 þ lustration of the chemical bonding in parallel orientation (left) and perpen- Potential energy curves E(R, q)ofthe S g (1s g ) state (D) and the S u (1s g 1s u ) dicular orientation (right) relative to the magnetic field, represented by red state (E) calculated at different field strengths in parallel orientation (q =0°, 2 1 þ arrows. (B and C) Potential energy surfaces E(R, q)ofH 2 in the S g (1s g ) state solid lines) and perpendicular orientation (q = 90°, dashed lines). The areas 3 þ (B) and the S u (1s g 1s u ) state (C) calculated at the FCI/un-aug-cc-pVTZ level between the full and dashed lines represent the energy for intermediate of theory in a field of strength B =1.0B 0 , using a polar coordinate system (skew) orientations in the field. The minimum of each curve is marked with where R is the internuclear separation and q is the angle of the molecular a black dot. 328 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS and unbiased manner. A more flexible but com- their orientation in a magnetic field, as illustrated ponent in Fig. 1E (the ground state for all finite putationally demanding scheme has been pro- for the parallel and perpendicular orientations in fields considered here) is lowered paramagneti- posed by Kennedy and Kobe, who equipped Fig. 1A. In the covalently bound singlet state, the cally and the system becomes more compact, as the wave function with a variationally optimized electronic energy of H 2 increases diamagnetically expected from the shrinking atomic size. The mol- phase factor (25). with increasing field strength (Fig. 1D). More- ecule thereby becomes bound, with a preferred In all calculations reported here, we use the over, the molecule becomes shorter and more perpendicular orientation in the field. At B = correlation-consistent aug-cc-pVTZ basis set strongly bound, reflecting a more pronounced 2.25B 0 , for instance, the bond distance is 92 pm –1 (26, 27) in uncontracted form (denoted un-aug- diamagnetic behavior in the dissociation limit and the dissociation energy 38 kJ mol .The bar- –1 cc-pVTZ). For B on the order of or greater than than in the united-atom limit. As expected for a rier to rotation is only slightly lower, 34 kJ mol , B 0 , the anisotropic distortion of the electronic covalently bound state, the energy is lower in the with saddle points at a large internuclear sep- distribution by the magnetic field is most eco- parallel orientation than in the perpendicular aration of 241 pm. The other two triplet com- nomically described using anisotropic AOs, orientation (Fig. 1A), owing to the greater in- ponents behave in the same manner but with the with different Gaussian exponents in the parallel teratomic overlap in this orientation. In the polar energy shifted upward by the Zeeman interaction. and perpendicular directions relative to the field plot for B = B 0 (Fig. 1B), we therefore observe The different shapes of the H 2 singlet and triplet (28, 29). The isotropic un-aug-cc-pVTZ basis global minima at inclination angles q =0°, 180° surfaces are immediately apparent from the polar that we use is sufficient for an accurate descrip- connected by saddle points at q =90°, 270°.As plots in Fig. 1, B and C. Although both surfaces are tion of electronic systems in fields up to B 0 but the field increases from 0 to 2.25B 0 (the strongest distinctively prolate in the field direction, the plots becomes progressively less suited in stronger fields, field considered here), the bond distance R e de- reveal two very different states: a fairly isotropic, as the systems become more compact and an- creases by 24% from 74 to 56 pm, while the compact singlet state with parallel global minima isotropic. Selected FCI calculations carried out dissociation energy D e (without the zero-point connected by perpendicular saddle points, con- in the larger un-aug-cc-pVQZ basis and in the vibrational contribution) increases by 83%, from trasted with a more anisotropic, diffuse triplet –1 un-aug-cc-pVTZ basis with extra orbitals added 455 kJ mol –1 to 834 kJ mol . Because of the state with a larger inaccessible inner region and confirm that the un-aug-cc-pVTZ basis provides prolate shape of the atoms, the saddle points perpendicular minima connected by parallel sad- on July 19, 2012 a qualitatively correct description of the systems for rotation occur at a slightly shorter distance dle points. –1 studied here. of 50 pm, with a barrier to rotation of 239 kJ mol . To understand the mechanism responsible for Because of the shrinking prolate shape of the The triplet state S ð1s g 1s*Þ behaves very the field-induced perpendicular bonding, we ex- 3 þ u u www.sciencemag.org constituent atoms, diatomic molecules become differently from the singlet state. With increasing amine the behavior of the bonding 1s g and an- smaller and develop an energy dependence on field strength, the energy of the bb triplet com- tibonding 1s* orbitals in an external magnetic u AB C (kJ mol -1 ) (pm) (pm) Downloaded from Fig. 2. Orbital energies of the bonding 1s g (blue) and antibonding 1s u (red) exponent a = 1 for different internuclear separations R.(B) The field- H 2 orbitals in parallel orientation (solid lines) and perpendicular orientation induced change in the orbital bonding energy DE p (R, q, B 0 ) calculated (dashed lines) relative to an external magnetic field. (A) The field-induced with an optimized exponent for different internuclear separations R.(C) change in the orbital bonding energy DE p (R, q, B 0 ) calculated with a fixed The MO energy level diagram in a magnetic field. AB 2 2 Fig. 3. (A) Potential energy curve of He 2 in the S g (1s g1s u Þ state calculated using FCI/un-aug-cc-pVTZ theory in parallel orientation (solid lines) and 1 þ 2 perpendicular orientation (dashed lines) for 0 ≤ B ≤ 2.5B 0 .(B)Sameas(A) forHe 2 in the S u (1s g1s u2s g) state. The energy minimum of each curve is marked 3 þ with a black dot. SCIENCE 329 www.sciencemag.org VOL 337 20 JULY 2012

REPORTS field. In a minimal basis consisting of two 1s the molecular configuration (R, q) relative to preferred parallel field orientation, as observed Gaussian orbitals of the form given in Eq. 3 the change observed in the dissociation limit. by Kubo using Hartree-Fock theory (7). with i = j = k = 0, the normalized bonding and With a fixed orbital exponent a =1,weobservethe To demonstrate that the behavior observed antibonding MOs of H 2 located on the z axis in expected stabilization of 1s* in the perpendicular for H 2 is a general phenomenon, we have calcu- u an external field B of arbitrary orientation are orientation (and a smaller 1s g destabilization in lated the potential energy curves of He 2 in its S ð1s 1s* Þ singlet (Fig. 3A) and given by the bonding region), and neither stabilization lowest 1 þ 2 g u 2 g S ð1s 1s*2s g Þ triplet (Fig. 3B) states. With- nor destabilization in hio the parallel orientation 3 þ 2 a n −1=2 u g u 2 2 1s g=u ¼ 2 þ 2 exp − ð1 þ ˜ B þ ˜ B ÞR 2 (Fig. 2A). However, when the exponent a is var- out a field, the He 2 singlet ground state is weakly − y x 2 iationally optimized for each (R, q) (Fig. 2B), 1s g bound by dispersion, with R e = 297 pm and D e = − is stabilized in the united atom (where the ori- 0.092 kJ mol (30). In the field, the energy in- ð1s A þ 1s B Þð4Þ –1 entation no longer matters for this orbital), where- creases diamagnetically. Moreover, He 2 assumes where 1s A and 1s B are 1s orbitals with exponent as 1s* (for which the orientation in the united a perpendicular orientation, becoming smaller u a on the two atoms, R is the internuclear sep- atom matters) is destabilized in the parallel ori- and more strongly bound, with a bond distance aration, and ˜ B x ¼ B x =4a and ˜ B y ¼ B y =4a are entation but further stabilized in the perpendicular of 94 pm and a dissociation energy of 31 kJ –1 scaled perpendicular field components. There is orientation. These changes lead to the modified mol at B =2.5B 0 . The nonmonotonic variation no contribution from B z to the MOs. When R MO energy level diagram shown in Fig. 2C and of the saddle point (in the parallel orientation) tends to zero, these MOs transform smoothly into to the following energy ordering of the lowest may be an artifact arising from basis set incom- helium AOs of the same exponent: H 2 states in a magnetic field: pleteness. In the triplet state, the covalently bound He 2 molecule behaves similarly to H 2 in the sin- lim 1s g ¼ 1s ð5Þ 2 2   glet state. As the field increases to 2.5B 0 ,the u u g R→0 E ∥ ð1s Þ≤ E ⊥ ð1s Þ≤ E ⊥ ð1s g 1s Þ≤ E ∥ ð1s g 1s Þ molecule aligns with the field and shortens from g ð8Þ 104 pm to 80 pm while the dissociation energy 2 2 −1=2 –1 –1 lim 1s* ¼ð1 þ ˜ B þ ˜ B Þ increases from 178 kJ mol to 655 kJ mol .In on July 19, 2012 R→0 u x y The field-induced bonding of H 2 in the per- the bb component of the triplet, He 2 becomes ð2p þ i ˜ B x 2p − i ˜ B y 2p Þð6Þ x y z pendicular orientation is thus not covalent in diamagnetic at B ≈ 2.2B 0 . Whereas 1s g transforms into an 1s orbital, 1s* u nature, nor does it depend on dispersion. Instead, We have presented advanced FCI calculations transforms into a combination of 2p orbitals. In it arises from a lowering (relative to the atomic on H 2 and He 2 in strong magnetic fields and particular, for ˜ B ¼ 1, 1s* becomes 2p 0 in the limit) of the kinetic energy associated with the explained their behavior in terms of elementary u parallel orientation and 2p –1 in the perpendicular induced paramagnetic rotation of the electron in concepts of MO theory. To examine the role of orientation. In a magnetic field, the 2p –1 orbital the antibonding orbital. The pivotal role of the electron correlation, we compare the FCI and has a lower energy than the 2p 0 orbital (by the kinetic energy is confirmed by FCI calculations Hartree-Fock potential energy curves of triplet www.sciencemag.org orbital-Zeeman interaction); therefore, 1s* favors on H 2 ; although the field-induced changes in the H 2 (Fig. 4A) and singlet He 2 (Fig.4B).These u a perpendicular orientation relative to the mag- FCI kinetic and electrostatic energies at a given plots demonstrate that the paramagnetic perpen- netic field. No such orientational preference is (R, q) are of the same order of magnitude, the ap- dicular bonding discussed above does not require observed for 1s g , which transforms into the same pearance of a minimum in the dissociation curve electron correlation for its qualitative description; 1s orbital in all orientations. By this argument, is almost entirely due to the lowering of the ki- Hartree-Fock theory, in which electronic inter- H 2 adopts a perpendicular orientation in the netic energy. actions are described in an averaged, mean-field triplet state, with its singly occupied bonding As shown in Fig. 2, the bonding 1s g orbital manner, recovers all the main effects of para- and antibonding orbitals. favors a perpendicular orientation in a magnetic magnetic bonding, underestimating the bond dis- Downloaded from Let E p (R, q, B) be the orbital energy of f p at field for intermediate bond distances. In general, tance slightly and the dissociation energy more configuration (R, q) and in the field B, and therefore, the orientation of H 2 in the triplet state strongly. The bonding is clearly not van der Waals consider the quantity depends on a balance between the preference of in nature, although electron correlation is neces- 1s g for a parallel orientation and the preference sary for its quantitative description. In short, we DE ðR,q,BÞ¼½E ðR,q,BÞ − E ðR,q,0ފ − p p p of1s* for a perpendicular orientation. Indeed, for have identified a distinct mechanism for chemical u B =2.25B 0 ,the H 2 minimum shifts slightly away bonding in strong magnetic fields, arising from ½E ð∞,q,BÞ − E ð∞,q,0ފ ð7Þ p p from the perpendicular orientation. Eventually, the stabilization of antibonding orbitals in a which represents the field-induced change in the H 2 ground state changes from S ð1s g 1s*Þ perpendicular orientation relative to the magnetic 3 þ u u 3 the orbital energy E p (R, q, B) – E p (R, q,0)at to P u (1s g 1p u ), which is covalently bound with a field. This stabilization leads to the bonding of Fig. 4. (A) The FCI and AB unrestricted Hartree-Fock (HF) dissociation curves 3 þ of H 2 in the S (1s g1s u) u state for B ⊥ = 2.25B 0 calculated in the un-aug- cc-pVTZ basis set. (B) Same as (A) for He 2 in 2 2 1 þ the S g (1s g1s u ) statefor B ⊥ =2.5B 0 . The Hartree- Fock model overestimates the bond distances of H 2 and He 2 by 1.5 and 4.1%,respectively,whereas the dissociation energies are underestimated by 25 and 49%. The counterpoise corrections for the basis-set superposition error (not added to the plotted curves) –1 are 4 kJ mol –1 and 2 kJ mol , respectively, for H 2 and He 2 at the equilibrium distances. 330 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS species of zero bond order, which are either un- 10. P. J. Knowles, N. C. Handy, Chem. Phys. Lett. 111, 24. E. I. Tellgren, T. Helgaker, A. Soncini, Phys. Chem. Chem. bound or bound by dispersion in the absence of a 315 (1984). Phys. 11, 5489 (2009). 11. W. Duch, GRMS or Graphical Representation of Model 25. P. K. Kennedy, D. H. Kobe, Phys. Rev. A 30, 51 (1984). magnetic field. This bonding is sufficiently strong Spaces I: Basics (Springer-Verlag, New York, 1986). 26. T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). to affect the chemistry of molecules in strong 12. J. Olsen, B. O. Roos, P. Jørgensen, H. J. A. Jensen, 27. D. E. Woon,T.H.Dunning, J. Chem. Phys. 100, 2975 (1994). magnetic fields. J. Chem. Phys. 89, 2185 (1988). 28. P. Schmelcher, L. S. Cederbaum, Phys. Rev. A 37, 672 13. T. Helgaker, P. Jørgensen, J. Olsen, Molecular (1988). Electronic-Structure Theory (Wiley, Chichester, UK, 29. U. Kappes, P. Schmelcher, J. Chem. Phys. 100, 2878 References and Notes 2000). 1. T. Helgaker, W. Klopper, D. P. Tew, Mol. Phys. 106, (1994). 14. P. Schmelcher, L. S. Cederbaum, Phys. Rev. A 41, 4936 30. M. Jeziorska, W. Cencek, K. Patkowski, B. Jeziorski, 2107 (2008). (1990). 2. R. H. Garstang, Rep. Prog. Phys. 40, 105 (1977). K. Szalewicz, J. Chem. Phys. 127, 124303 (2007). 15. F. London, J. Phys. Radium 8, 397 (1937). 3. S. Jordan, P. Schmelcher, W. Becken, W. Schweizer, 16. R. Ditchfield, J. Chem. Phys. 56, 5688 (1972). Acknowledgments: Supported by the Norwegian Research Astron. Astrophys. 336, L33 (1998). 17. K. Wolinski, J. F. Hinton, P. Pulay, J. Am. Chem. Council through Centre for Theoretical and Computational 4. S. Jordan, P. Schmelcher, W. Becken, Astron. Astrophys. Phys. Soc. 112, 8251 (1990). Chemistry (CTCC) grant 179568/V30 and through grant 376, 614 (2001). 18. W. Kutzelnigg, Isr. J. Chem. 19, 193 (1980). 197446/V30 and by the European Research Council (ERC) 5. D. Lai, Rev. Mod. Phys. 73, 629 (2001). 19. M. Schindler, W. Kutzelnigg, J. Chem. Phys. 76, 1919 under the European Union’s Seventh Framework Program 6. M. Žaucer, A. Ažman, Phys. Rev. A 18, 1320 (1978). (1982). through the advanced grant ABACUS, ERC grant agreement 7. A. Kubo, J. Phys. Chem. A 111, 5572 (2007). 20. U. Kappes, P. Schmelcher, Phys. Lett. A 210, 409 (1996). 267683. M.R.H. acknowledges support from the CTCC during a 8. Y. E. Lozovik, A. V. Klyuchnik, Phys. Lett. A 66, 282 21. U. Kappes, P. Schmelcher, Phys. Rev. A 53, 3869 (1996). sabbatical stay at the University of Oslo in 2010. (1978). 22. U. Kappes, P. Schmelcher, Phys. Rev. A 54, 1313 (1996). 9. S. Basile, F. Trombetta, G. Ferrante, Nuovo Cim. 9, 23. E. I. Tellgren, A. Soncini, T. Helgaker, J. Chem. Phys. 129, 26 January 2012; accepted 25 May 2012 457 (1987). 154114 (2008). 10.1126/science.1219703 Sulfate Burial Constraints on the evaporites from seawater and obtain estimates of on July 19, 2012 It is possible to measure the sink of sulfate 34 the influx magnitude and d S by mass balance, Phanerozoic Sulfur Cycle though previous volume estimates of Phanero- zoic evaporites (mostly halite, but some sulfate) 3 Itay Halevy, 1,2 * Shanan E. Peters, Woodward W. Fischer 2 have been considered too coarse or uncertain to accurately constrain past rates of sulfate burial (16–18). We quantified sulfate burial over Phan- The sulfur cycle influences the respiration of sedimentary organic matter, the oxidation state of erozoic time, using a comprehensive macrostrati- the atmosphere and oceans, and the composition of seawater. However, the factors governing the graphic database (19, 20), which includes 23,843 www.sciencemag.org major sulfur fluxes between seawater and sedimentary reservoirs remain incompletely understood. lithostratigraphic rock units in 949 geographic Using macrostratigraphic data, we quantified sulfate evaporite burial fluxes through Phanerozoic locations across North America and the Carib- time. Approximately half of the modern riverine sulfate flux comes from weathering of recently bean (NAC). Data were binned by age, and sulfate deposited evaporites. Rates of sulfate burial are unsteady and linked to changes in the area of burial rates were obtained by dividing evapo- marine environments suitable for evaporite formation and preservation. By contrast, rates of rite volume by bin duration. Macrostratigraphy- pyrite burial and weathering are higher, less variable, and largely balanced, highlighting a based estimates of sulfate burial rates are higher greater role of the sulfur cycle in regulating atmospheric oxygen. than those derived from other compilations. This is due to the improved spatial and lithological 2– Downloaded from ulfate (SO 4 ) is the fourth most abundant The chemical composition of fluid inclusions in resolution of this data set, which includes sedi- ion in modern seawater and a major com- halite constrains the concentration of major ions mentary rocks in the surface and subsurface, and Sponent of the alkalinity budget, which gov- in seawater, including sulfate (9, 10). many comparatively thin but widespread depos- 34 erns the pH of seawater (1). Bacterial sulfate Variability in the d S records of seawater its not included in previous compilations. Nota- reduction accounts for ~50% of sedimentary sulfate and sedimentary pyrite is typically inter- bly, the NAC burial rates are highly variable, with organic matter respiration (2), and precipitation preted to reflect changes in the fraction of sulfur values 2 to 14 times the average occurring mainly of pyrite (FeS 2 ) is one of the major exit channels removed from the oceans as pyrite, f pyr . Because in Paleozoic intervals (Fig. 1C). of sulfur from the ocean (3). Because reduction pyrite is depleted in 34 S by several percent rel- The macrostratigraphic database currently pro- of riverine sulfate and burial of the sulfide leave ative to the sulfate reservoir from which it formed, vides comprehensive coverage only in NAC, but 34 oxidized products in the ocean-atmosphere sys- times of high seawater sulfate d S are interpreted can be scaled globally (Fig. 1D) by using mech- tem, pyrite burial is considered a major indirect as times of high rates of pyrite burial. By as- anistic relationships between the observations source of oxygen to the atmosphere (4, 5). suming a steady state and constant input mag- and environmental controls on sulfate evaporite 34 Several time series data sets constrain aspects nitude and d S, or by scaling inputs and outputs deposition (20). The volume-weighted average of the Phanerozoic sulfur cycle (Fig. 1A). The to modern values, models of the Phanerozoic sul- ratio of global to North American sulfate deposit 34 34 sulfur isotope composition, d S, of carbonate- fur cycle explain long-term trends in d Svalues volumesis~8 (16, 17). In comparison, the area- associated sulfate, sulfate evaporites, and barite by changes in f pyr between ~0.2 and ~0.6 (4, 11–13). weighted ratio of global to NAC submerged 34 34 (BaSO 4 ) records the d S of seawater sulfate, Recognizing that the magnitude and d Softhe continental area in latitudes of net evaporation, 34 whereas the d S of sedimentary pyrite captures influxes to the ocean have likely varied in time, estimated from paleogeographic reconstructions the products of microbial sulfate reduction (6–8). thereby influencing the isotopic record, some mod- (20, 21), is ~7. This close agreement reflects a els included parameterized influxes and solved primary requirement for massive sulfate evapo- 1 Environmental Sciences and Energy Research, Weizmann In- mass balance equations for the outfluxes and the rite deposition—hydrographic isolation of large, 2 stitute of Science, Rehovot 76100, Israel. Geological and Plan- value of f pyr (4, 13). The parameterizations are marine-fed basins at latitudes of net evaporation etary Sciences, California Institute of Technology, Pasadena, CA 3 91125, USA. Geoscience, University of Wisconsin-Madison, uncertain, however, because they are largely based (22). Such basins are created by rifting, small Madison, WI 53706, USA. on a scaling of modern influxes by debated fac- changes in sea level or the development of a *To whom correspondence should be addressed. E-mail: tors, such as the relative rates of seafloor spread- barrier to circulation (22), often at the shoreward [email protected] ing and continental runoff (14, 15). edge of submerged continental shelves. Indeed, SCIENCE 331 www.sciencemag.org VOL 337 20 JULY 2012

Sulfate Burial Constraints on the Phanerozoic Sulfur Cycle Itay Halevy et al. Science 337 , 331 (2012); DOI: 10.1126/science.1220224 This copy is for your personal, non-commercial use only. If you wish to distribute this article to others , you can order high-quality copies for your colleagues, clients, or customers by clicking here. Permission to republish or repurpose articles or portions of articles can be obtained by following the guidelines here. The following resources related to this article are available online at on July 19, 2012 www.sciencemag.org (this information is current as of July 19, 2012 ): Updated information and services, including high-resolution figures, can be found in the online version of this article at: http://www.sciencemag.org/content/337/6092/331.full.html Supporting Online Material can be found at: http://www.sciencemag.org/content/suppl/2012/07/18/337.6092.331.DC1.html www.sciencemag.org A list of selected additional articles on the Science Web sites related to this article can be found at: http://www.sciencemag.org/content/337/6092/331.full.html#related This article cites 37 articles , 14 of which can be accessed free: http://www.sciencemag.org/content/337/6092/331.full.html#ref-list-1 This article has been cited by 1 articles hosted by HighWire Press; see: Downloaded from http://www.sciencemag.org/content/337/6092/331.full.html#related-urls This article appears in the following subject collections: Geochemistry, Geophysics http://www.sciencemag.org/cgi/collection/geochem_phys Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2012 by the American Association for the Advancement of Science; all rights reserved. The title Science is a registered trademark of AAAS.

REPORTS species of zero bond order, which are either un- 10. P. J. Knowles, N. C. Handy, Chem. Phys. Lett. 111, 24. E. I. Tellgren, T. Helgaker, A. Soncini, Phys. Chem. Chem. bound or bound by dispersion in the absence of a 315 (1984). Phys. 11, 5489 (2009). 11. W. Duch, GRMS or Graphical Representation of Model 25. P. K. Kennedy, D. H. Kobe, Phys. Rev. A 30, 51 (1984). magnetic field. This bonding is sufficiently strong Spaces I: Basics (Springer-Verlag, New York, 1986). 26. T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). to affect the chemistry of molecules in strong 12. J. Olsen, B. O. Roos, P. Jørgensen, H. J. A. Jensen, 27. D. E. Woon,T.H.Dunning, J. Chem. Phys. 100, 2975 (1994). magnetic fields. J. Chem. Phys. 89, 2185 (1988). 28. P. Schmelcher, L. S. Cederbaum, Phys. Rev. A 37, 672 13. T. Helgaker, P. Jørgensen, J. Olsen, Molecular (1988). Electronic-Structure Theory (Wiley, Chichester, UK, 29. U. Kappes, P. Schmelcher, J. Chem. Phys. 100, 2878 References and Notes 2000). 1. T. Helgaker, W. Klopper, D. P. Tew, Mol. Phys. 106, (1994). 14. P. Schmelcher, L. S. Cederbaum, Phys. Rev. A 41, 4936 30. M. Jeziorska, W. Cencek, K. Patkowski, B. Jeziorski, 2107 (2008). (1990). 2. R. H. Garstang, Rep. Prog. Phys. 40, 105 (1977). K. Szalewicz, J. Chem. Phys. 127, 124303 (2007). 15. F. London, J. Phys. Radium 8, 397 (1937). 3. S. Jordan, P. Schmelcher, W. Becken, W. Schweizer, 16. R. Ditchfield, J. Chem. Phys. 56, 5688 (1972). Acknowledgments: Supported by the Norwegian Research Astron. Astrophys. 336, L33 (1998). 17. K. Wolinski, J. F. Hinton, P. Pulay, J. Am. Chem. Council through Centre for Theoretical and Computational 4. S. Jordan, P. Schmelcher, W. Becken, Astron. Astrophys. Phys. Soc. 112, 8251 (1990). Chemistry (CTCC) grant 179568/V30 and through grant 376, 614 (2001). 18. W. Kutzelnigg, Isr. J. Chem. 19, 193 (1980). 197446/V30 and by the European Research Council (ERC) 5. D. Lai, Rev. Mod. Phys. 73, 629 (2001). 19. M. Schindler, W. Kutzelnigg, J. Chem. Phys. 76, 1919 under the European Union’s Seventh Framework Program 6. M. Žaucer, A. Ažman, Phys. Rev. A 18, 1320 (1978). (1982). through the advanced grant ABACUS, ERC grant agreement 7. A. Kubo, J. Phys. Chem. A 111, 5572 (2007). 20. U. Kappes, P. Schmelcher, Phys. Lett. A 210, 409 (1996). 267683. M.R.H. acknowledges support from the CTCC during a 8. Y. E. Lozovik, A. V. Klyuchnik, Phys. Lett. A 66, 282 21. U. Kappes, P. Schmelcher, Phys. Rev. A 53, 3869 (1996). sabbatical stay at the University of Oslo in 2010. (1978). 22. U. Kappes, P. Schmelcher, Phys. Rev. A 54, 1313 (1996). 9. S. Basile, F. Trombetta, G. Ferrante, Nuovo Cim. 9, 23. E. I. Tellgren, A. Soncini, T. Helgaker, J. Chem. Phys. 129, 26 January 2012; accepted 25 May 2012 457 (1987). 154114 (2008). 10.1126/science.1219703 Sulfate Burial Constraints on the evaporites from seawater and obtain estimates of on July 19, 2012 It is possible to measure the sink of sulfate 34 the influx magnitude and d S by mass balance, Phanerozoic Sulfur Cycle though previous volume estimates of Phanero- zoic evaporites (mostly halite, but some sulfate) have been considered too coarse or uncertain to 3 Itay Halevy, 1,2 * Shanan E. Peters, Woodward W. Fischer 2 accurately constrain past rates of sulfate burial (16–18). We quantified sulfate burial over Phan- The sulfur cycle influences the respiration of sedimentary organic matter, the oxidation state of erozoic time, using a comprehensive macrostrati- the atmosphere and oceans, and the composition of seawater. However, the factors governing the graphic database (19, 20), which includes 23,843 www.sciencemag.org major sulfur fluxes between seawater and sedimentary reservoirs remain incompletely understood. lithostratigraphic rock units in 949 geographic Using macrostratigraphic data, we quantified sulfate evaporite burial fluxes through Phanerozoic locations across North America and the Carib- time. Approximately half of the modern riverine sulfate flux comes from weathering of recently bean (NAC). Data were binned by age, and sulfate deposited evaporites. Rates of sulfate burial are unsteady and linked to changes in the area of burial rates were obtained by dividing evapo- marine environments suitable for evaporite formation and preservation. By contrast, rates of rite volume by bin duration. Macrostratigraphy- pyrite burial and weathering are higher, less variable, and largely balanced, highlighting a based estimates of sulfate burial rates are higher greater role of the sulfur cycle in regulating atmospheric oxygen. than those derived from other compilations. This is due to the improved spatial and lithological 2– Downloaded from ulfate (SO 4 ) is the fourth most abundant The chemical composition of fluid inclusions in resolution of this data set, which includes sedi- ion in modern seawater and a major com- halite constrains the concentration of major ions mentary rocks in the surface and subsurface, and Sponent of the alkalinity budget, which gov- in seawater, including sulfate (9, 10). many comparatively thin but widespread depos- 34 erns the pH of seawater (1). Bacterial sulfate Variability in the d S records of seawater its not included in previous compilations. Nota- reduction accounts for ~50% of sedimentary sulfate and sedimentary pyrite is typically inter- bly, the NAC burial rates are highly variable, with organic matter respiration (2), and precipitation preted to reflect changes in the fraction of sulfur values 2 to 14 times the average occurring mainly of pyrite (FeS 2 ) is one of the major exit channels removed from the oceans as pyrite, f pyr . Because in Paleozoic intervals (Fig. 1C). of sulfur from the ocean (3). Because reduction pyrite is depleted in 34 S by several percent rel- The macrostratigraphic database currently pro- of riverine sulfate and burial of the sulfide leave ative to the sulfate reservoir from which it formed, vides comprehensive coverage only in NAC, but 34 oxidized products in the ocean-atmosphere sys- times of high seawater sulfate d S are interpreted can be scaled globally (Fig. 1D) by using mech- tem, pyrite burial is considered a major indirect as times of high rates of pyrite burial. By as- anistic relationships between the observations source of oxygen to the atmosphere (4, 5). suming a steady state and constant input mag- and environmental controls on sulfate evaporite 34 Several time series data sets constrain aspects nitude and d S, or by scaling inputs and outputs deposition (20). The volume-weighted average of the Phanerozoic sulfur cycle (Fig. 1A). The to modern values, models of the Phanerozoic sul- ratio of global to North American sulfate deposit 34 34 sulfur isotope composition, d S, of carbonate- fur cycle explain long-term trends in d Svalues volumesis~8 (16, 17). In comparison, the area- associated sulfate, sulfate evaporites, and barite by changes in f pyr between ~0.2 and ~0.6 (4, 11–13). weighted ratio of global to NAC submerged 34 34 (BaSO 4 ) records the d S of seawater sulfate, Recognizing that the magnitude and d Softhe continental area in latitudes of net evaporation, 34 whereas the d S of sedimentary pyrite captures influxes to the ocean have likely varied in time, estimated from paleogeographic reconstructions the products of microbial sulfate reduction (6–8). thereby influencing the isotopic record, some mod- (20, 21), is ~7. This close agreement reflects a els included parameterized influxes and solved primary requirement for massive sulfate evapo- 1 Environmental Sciences and Energy Research, Weizmann In- mass balance equations for the outfluxes and the rite deposition—hydrographic isolation of large, 2 stitute of Science, Rehovot 76100, Israel. Geological and Plan- value of f pyr (4, 13). The parameterizations are marine-fed basins at latitudes of net evaporation etary Sciences, California Institute of Technology, Pasadena, CA 3 91125, USA. Geoscience, University of Wisconsin-Madison, uncertain, however, because they are largely based (22). Such basins are created by rifting, small Madison, WI 53706, USA. on a scaling of modern influxes by debated fac- changes in sea level or the development of a *To whom correspondence should be addressed. E-mail: tors, such as the relative rates of seafloor spread- barrier to circulation (22), often at the shoreward [email protected] ing and continental runoff (14, 15). edge of submerged continental shelves. Indeed, SCIENCE 331 www.sciencemag.org VOL 337 20 JULY 2012

REPORTS the long time-scale variability in the burial rate ability of suitable environments in which evap- ing constraints from sulfur isotope measurements data is well explained by the estimated NAC sub- oration of seawater could lead to saturation, and sulfate concentration data from fluid inclu- merged continental area at latitudes T10° to 50° precipitation, accumulation, and long-term pres- sions. In any isotope mass balance model, the 34 (linear product moment correlation coefficient of ervation of sulfate evaporites. Unsteady sulfate time-dependent concentration and d S of sea- 0.47 at the temporal resolution of the paleogeo- burial rates, governed by the interactions between water sulfate can be expressed by two equations graphic reconstructions), and we use this relation- sea level, tectonics, and paleogeography, suggest dM ship to derive average global burial rates. We add that the critical statistic derived from isotope ¼ J in − J e − J p ð1Þ shorter time scale variability to this average using mass balance studies, f pyr , convolves information dt correlations between NAC fluxes and the rate of regarding the relative activity of sulfate-reducing dMd change in eustatic sea level (20, 23). Additional microbiota with the availability of environments ¼ J in d in − J e d − J p ðd − DÞð1Þ ð2Þ factors that influence evaporite deposition, such suitable for sulfate burial. For example, a de- dt as geodynamic controls on basin subsidence and crease in the area of shallow seas due to a drop in Here, M is the concentration of seawater sul- 34 climate patterns at the regional-to-local scale (22), sea level or continental migrations could decrease fate and d is its d Svalue. J in is the total influx of are not explicitly represented by this scaling meth- sulfate burial rates and increase f pyr . The rate of sulfur to the oceans, with contributions from odology. They are, however, implicitly included pyrite burial itself need not change. Scaled glob- evaporite weathering, oxidative weathering of in the scaling because they have contributed to ally and corrected for the decay of surviving rock sedimentary and igneous sulfide minerals, and 34 the observed NAC sulfate burial rates (20). with time (20), the macrostratigraphic data reveal volcanic outgassing of sulfur volatiles. The d S Intervals of rapid sulfate evaporite burial (NAC an average Phanerozoic sulfate evaporite burial rate of this influx, d in , depends on the relative con- 11 –1 11 and global) occur with equal frequency during of ~3.3 × 10 to 4.5 × 10 mol year , depending tributions of these three components. J e and J p times of high and low marine sulfate concen- on bin duration. This is only ~10 to 30% of the are the burial fluxes of sulfate evaporites and py- trations (Fig. 1). This reveals that sulfate burial estimated riverine influx of sulfate to the oceans rite, respectively, and D is the average difference 12 –1 12 34 rates were disconnected from changes in the ac- [~1.5 × 10 to 3.5 × 10 mol year ;(12, 24)], (in permil) between the d S of contemporaneous 2– 2+ tivity product of calcium and sulfate [aCa aSO 4 ; implying that the value of f pyr has been ~0.7 to sulfate evaporite and pyrite sedimentary sinks. on July 19, 2012 e.g., (13)], which have varied by up to ~15% 0.9 if the sulfur cycle operated close to steady state. The values of M, d,and D can be constrained by 34 (9, 10). We suggest, instead, that the first-order The macrostratigraphic data motivated us to d S and fluid inclusion data, and J e , by macro- control on sulfate burial was the episodic avail- reevaluate models of the sulfur cycle by integrat- stratigraphic data. www.sciencemag.org 34 δ S NAC shelf area NAC J e AB C (10 mol yr ) D 6 −1 12 2 (10 km ) (‰ VCDT) −40 −20 0 20 40 0 5 10 15 0.0 0.5 1.0 1.5 0 Δ 100 sulfate Downloaded from 200 pyrite Age (Ma) 300 400 500 0 10 20 30 −2 0 2 4 0 50 100 150 0 1 2 −2 [SO ] Change in eustatic sea level NAC J e : shelf area Global J 4 e −1 12 −1 3 −1 −2 (mM) (m Myr ) (10 mol yr km ) (10 mol yr ) Fig. 1. Observational constraints on the Phanerozoic sulfur cycle. (A) merged continental area estimates and the rate of change of global mean Sulfur isotope composition of seawater sulfate and sedimentary pyrite (8), sea level (23). (C) Estimates of NAC total and per-area sulfate evaporite and seawater sulfate concentration from fluid inclusions in halite (9, 10). burial rates. (D) NAC sulfate evaporite burial rates scaled to global fluxes Isotope compositions are reported relative to the standard VCDT (Vienna and corrected for decay of the surviving record (20), including uncertainty Canyon Diablo Troilite). (B) North America and the Caribbean (NAC) sub- (shaded envelope). 332 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org

REPORTS It is possible to solve Eqs. 1 and 2 for the latter is related to the former; if J in is small, then hypothesize that this is because parameterizations outputs using the inputs (“parameterized input”) the change in d in required to drive an observed for the sulfate influx into the ocean are calibrated 34 or for the inputs using the outputs (“param- change in seawater sulfate d S must be large. to modern riverine fluxes, which are higher than eterized output”), and the solutions can be eval- The models of submerged continental area- the expected long-term average due to the ex- uated with independent physical understanding dependent or constant J p yield reasonable values posure and weathering of recently deposited evap- of the sulfur cycle. We examined three variants of of J in and low values of d in (near the average orites (supplementary text). This notion is supported 34 the parameterized input model (20). The first as- sedimentary pyrite d S) that become higher dur- by anomalously high per-area sulfate deposition sumed constant values of J in , d in ,and D (11, 12, 24) ing times of increased J e . This is consistent with rates obtained from the NAC macrostratigraphic and either a steady state or a dynamic mass the idea of rapid sediment recycling (4), where compilation during the last 10 million years (Fig. balance. In the second and third models, we used low values of d in closely follow high relative 1C); by widespread massive Neogene-age (~23 influx parameterizations similar to those of (18) rates of pyrite burial due to oxidative weathering to 2.6 Ma) evaporites in Europe, Asia, and Africa and (13). We examined three variants of the of recently deposited pyrite; high d in values result (22, 27–31); and by the recognition that deposi- parameterized output model (20). The first is a from weathering of newly deposited sulfates. No- tion rates decrease with increasing observation model of constant f pyr . The second, motivated by tably, dynamic and steady-state solutions diverge time scale to constant, long-term values (32). This the observation that pyrite burial is common in only during short intervals of rapid change in implies that the high abundance of Neogene-age shallow, organic-rich sediments (25), is a model seawater sulfate concentrations, implying that the sulfate evaporites represents gross deposition. Ul- in which J p scales with normalized global sub- system operates much of the time close to steady timately much of this material will not be pre- merged continental area. This scaling may be weak state. All three parameterized output models, served, and net deposition rates will converge to because although high submerged areas may lead when constrained by the macrostratigraphic data, long-term Phanerozoic rates (33). to high bacterial sulfate reduction rates, inefficient indicate high average values of f pyr (~0.7to0.9). The results presented here describe a Phaner- delivery of reactive iron across the shelf could The parameterized input models yield values ozoic sulfur cycle in which the majority of net limit pyrite burial to near-shore environments of f pyr similar to those in previous modeling inputs and outputs are oxidative weathering and (26). We therefore considered a third model of studies of the Phanerozoic sulfur cycle (Fig. 2, D burial of sedimentary pyrite, respectively. Over on July 19, 2012 constant J p . Using the macrostratigraphy-based and F). However, the model of constant d in at the time intervals resolved by our macrostrati- sulfate burial rates as additional constraints, we steady state occasionally requires unrealistically graphic data, the long-term value we estimate for 34 solved Eqs. 1 and 2 for the time-dependent val- low values of J in to reproduce the d S records f pyr is generally high, and large downward excur- ues of J in and d in . [Fig. 2E; e.g., 200 Ma (millions of years ago)], sions in its value are associated with high rates of 34 Of the parameterized output models (Fig. 2, supporting the notion that the d S value of sulfate evaporite burial, rather than times of less A to C), the models of constant f pyr yield un- inputs has varied with time. The values of J e pyrite burial. Although spatial heterogeneity and realistically low values of J in (often lower than required for mass balance under the other param- short–time scale variability in the value of J p have estimates of the volcanic influx and occasionally eterizations are ~three times as large as our burial likely occurred (34), global models of constant or www.sciencemag.org negative) and wildly varying values of d in .The rate estimates (Fig. 2, E and G, and fig. S1). We slowly varying J p (in response to changes in sub- 4 1.0 A B C 40 34 0.8 3 20 δ S e 0.6 J in [10 12 mol yr −1 ] 2 ƒ pyr Constant J p δ 34 S in [‰ VCDT] 0 0.4 Constant J p , steady state 1 δ S Downloaded from 34 Shelf area-dependent J p −20 p 0.2 f pyr = 0.6 0 −40 f pyr = 0.8 J 0.0 v 500 400 300 200 100 0 500 400 300 200 100 0 500 400 300 200 100 0 Age [Myr] 1.0 D Constant δ in , steady state E 4 4 , Constant J in δ in, dynamic mass balance Fig. 2. Model results. (A) f pyr in the param- 0.8 Macrostrat J e 3 3 eterized output models. (B and C) J in and d in 0.6 required for mass balance constrained by time ƒ pyr J e [10 12 mol yr −1 ] 2 2 J in [10 12 mol yr −1 ] series in Fig. 1. (D) f pyr in the constant-input 0.4 1 J 1 models. (E) J e calculated with the dynamic 0.2 v model (mass balance) compared to the macro- stratigraphic data, and J in (=J out )requiredto 0 0 34 reproduce seawater sulfate d Sinthe steady- 0.0 500 400 300 200 100 0 500 400 300 200 100 0 state model. (F) f pyr in the parameterized input models. (G) J e calculated by mass balance com- 1.0 F Parameterization 1 G pared to the macrostratigraphic data. Shaded 3 Parameterization 2 envelopes reflect uncertainty in J e estimates 0.8 Macrostrat J e (see Fig. 1). 0.6 J e [10 12 mol yr −1 ] 2 ƒ pyr 0.4 1 J v 0.2 0 0.0 500 400 300 200 100 0 500 400 300 200 100 0 Age [Myr] Age [Myr] SCIENCE 333 www.sciencemag.org VOL 337 20 JULY 2012

REPORTS merged continental area) yield results that are con- 12. D. E. Canfield, Am. J. Sci. 304, 839 (2004). 32. P. M. Sadler, GeoResearch Forum 5, 15 (1999). sistent both internally and with existing observations 13. A. Kampschulte, H. Strauss, Chem. Geol. 204, 255 33. S. E. Peters, J. Geol. 114, 391 (2006). (2004). 34. B. C. Gill et al., Nature 469, 80 (2011). 34 of seawater sulfate concentration and d S. Large 14. D. B. Rowley, Geol. Soc. Am. Bull. 114, 927 (2002). and stable pyrite weathering and burial fluxes 15. J. M. Edmond, Y. Huh, in Tectonic Uplift and Climate Acknowledgments: We thank D. Canfield and J. Adkins highlight the importance of oxidation-reduction Change, W. F. Ruddiman, W. Prell, Eds. (Plenum, New for helpful discussion, and C. Scotese for help with the feedbacks between carbon, iron, and sulfur (24) York, 1997), pp. 329–351. paleogeographic reconstructions. I.H. acknowledges support 16. A. B. Ronov, Int. Geol. Rev. 24, 1313 (1982). from a Texaco Postdoctoral Fellowship in Geological and and imply a greater role for the sulfur cycle in 17. M. A. Zharkov, History of Paleozoic Salt Accumulation, Planetary Sciences at the California Institute of Technology regulating Phanerozoic atmospheric oxygen. A. L. Yanshin, Ed. (Springer, New York, 1981). and a Sir Charles Clore Prize for Outstanding Appointment in 18. R. A. Berner, Am. J. Sci. 304, 438 (2004). the Experimental Sciences at the Weizmann Institute of References and Notes 19. S. E. Peters, N. A. Heim, Paleobiology 36, 61 (2010). Science. S.E.P. was funded by NSF grant EAR-0819931. 1. F. M. M. Morel, J. G. Hering, Principles and Applications 20. Methods are available on Science Online. W.W.F. acknowledges support from the Agouron Institute of Aquatic Chemistry (Wiley, New York, 1993). 21. C. R. Scotese, Atlas of Earth History (PALEOMAP Project, and a David and Lucile Packard Foundation Fellowship 2. B. B. Jørgensen, Nature 296, 643 (1982). Arlington, TX, 2001). for Science and Engineering. The binned macrostratigraphic 3. R. A. Berner, Geochim. Cosmochim. Acta 48, 605 (1984). 22. J. K. Warren, Earth Sci. Rev. 98, 217 (2010). data are available as a supplementary table on 4. R. A. Berner, Am. J. Sci. 287, 177 (1987). 23. B. U. Haq, A. M. Al-Qahtani, GeoArabia 10, 127 (2005). Science Online. 5. D. E. Canfield, Annu. Rev. Earth Planet. Sci. 33, 24. L. R. Kump, R. M. Garrels, Am. J. Sci. 286, 337 (1986). 1 (2005). 25. D. E. Canfield, R. Raiswell, Am. J. Sci. 299, 697 (1999). Supplementary Materials 34 34 34 6. d S=( R sam / R ref – 1) × 1000, where 34 R is the ratio 26. R. Raiswell, Elements 7, 101 (2011). www.sciencemag.org/cgi/content/full/337/6092/331/DC1 of 34 Sto 32 S. 27. J. M. Rouchy, D. Noel, A. M. A. Wali, M. A. M. Aref, Materials and Methods 7. H. Strauss, Palaeogeogr. Palaeocl. 132, 97 (1997). Sediment. Geol. 94, 277 (1995). Supplementary Text 8. N. P. Wu, J. Farquhar, H. Strauss, S. T. Kim, D. E. Canfield, 28. T. M. Peryt, Sediment. Geol. 188–189, 379 (2006). Figs. S1 to S5 Geochim. Cosmochim. Acta 74, 2053 (2010). 29. W. Krijgsman,F. J. Hilgen, I. Raffi,F.J.Sierro,D. S.Wilson, References (35–46) 9. T. K. Lowenstein, M. N. Timofeeff, S. T. Brennan, Nature 400, 652 (1999). Database S1 L. A. Hardie, R. V. Demicco, Science 294, 1086 (2001). 30. J. M. Rouchy, A. Caruso, Sediment. Geol. 188–189,35 10. J. Horita, H. Zimmermann, H. D. Holland, Geochim. (2006). Cosmochim. Acta 66, 3733 (2002). 31. H. Rahimpour-Bonab, Z. Shariatinia, M. G. Siemann, 7 February 2012; accepted 5 June 2012 on July 19, 2012 11. R. M. Garrels, A. Lerman, Am. J. Sci. 284, 989 (1984). Geol. J. 42, 37 (2007). 10.1126/science.1220224 Rapid Variability of Seawater Chemistry sporadically in the geologic record and are not always visible, and the chemical composition of seawater does not depend on the mass of ex- Over the Past 130 Million Years isting evaporites (Fig. 1A) but on the amount of salts that have already been eroded (Fig. 1B) www.sciencemag.org 1 2 Ulrich G. Wortmann * and Adina Paytan * and/or were originally extracted (Fig. 1C). The past 130 million years saw only two BSE events (Fig. 1C): one caused by the desiccation of the Fluid inclusion data suggest that the composition of major elements in seawater changes Mediterranean during the Messinian (12), and slowly over geological time scales. This view contrasts with high-resolution isotope data the second related to the Early Cretaceous opening that imply more rapid fluctuations of seawater chemistry. We used a non–steady-state of the South Atlantic (11, 12). 34 box model of the global sulfur cycle to show that the global d S record can be explained Until recently, all available data supported the by variable marine sulfate concentrations triggered by basin-scale evaporite precipitation assumption of slow secular changes to major and dissolution. The record is characterized by long phases of stasis, punctuated by seawater constituents and their related biogeo- Downloaded from short intervals of rapid change. Sulfate concentrations affect several important biological chemical fluxes. However, as the temporal reso- processes, including carbonate mineralogy, microbially mediated organic matter lution of proxy data increases, it becomes evident remineralization, sedimentary phosphorous regeneration, nitrogen fixation, and sulfate that the rate of change is often faster than ex- aerosol formation. These changes are likely to affect ocean productivity, the global pected from the residence time of major species. 34 carbon cycle, and climate. High-resolution data sets of seawater d S record two major events, at 130 to 120 million years he chemical composition and mineralogy Evaporites play an important role in the lat- ago (Ma) and 55 to 45 Ma, that require large of skeletal limestones and biogenic car- ter process because their precipitation/dissolution changes to the S-fluxes and/or their isotopic com- 34 Tbonates vary systematically through time, rates exceed those of other sediments by three position. Previous interpretations of the d S indicating that the Mg/Ca ratio as well as other orders of magnitude (7, 8). Halite is the domi- record called for changes in the planetary de- constituents of seawater have also changed (1). nant evaporite phase, but the effect of halite gassing flux and/or S burial rates (6). Clearly, Fluid inclusion studies are also consistent with precipitation/dissolution on seawater chemistry volcanic activity and pyrite burial are major pa- − + variable magnesium, calcium, sodium, and sul- is limited because the marine Na and Cl reser- rameters, but we propose to broaden the discourse 20 20 fate concentrations through time (2, 3). Several voirs are large (≈ 6.47 × 10 and 7.5 × 10 mol by including the effects of evaporite precipitation hypotheses have been proposed to explain these respectively). Sulfur-bearing evaporites (such as and dissolution. If these processes occur on a secular trends, including changes in global weath- CaSO 4 ) comprise on average 20% of an evap- basin-wide scale, they will modify the flux and 34 ering patterns (4), sea-floor spreading rates (5), or orite sequence (9) but have a 20 times smaller d S of the sulfur input/output and change the − burial rates of these elements (6). marine reservoir size than that of Cl . marine sulfate concentration, which in turn af- The role of sulfur-bearing salts in controlling fects pyrite burial in a nonlinear way (11), further seawater chemistry is often overlooked. In the affecting seawater isotopic composition. Under 1 Geobiology Isotope Laboratory, Department of Geology, Uni- modern ocean, pyrite burial is not limited by modern conditions, this effect is negligible. How- 2 versity of Toronto, Toronto, ON M5S 3B1, Canada. Institute of Marine Sciences University of California Santa Cruz, Santa sulfate availability (10), and the strong link be- ever, with decreasing sulfate concentrations, the Cruz, CA 95064, USA. tween pyrite burial rates and sulfate concentration importance of sulfate availability increases until *To whom correspondence should be addressed. E-mail: uli. has just recently been recognized (11). Further- it becomes the dominant parameter controlling [email protected] (U.G.W.); [email protected] (A.P.) more, basin-scale evaporites (BSEs) occur only pyrite burial (11). 334 20 JULY 2012 VOL 337 SCIENCE www.sciencemag.org