Research Article

Ion dissolution mechanism and kinetics at kink sites on NaCl surfaces

Mark N. Joswiak, Michael F. Doherty, and Baron Peters
PNAS January 23, 2018 115 (4) 656-661; first published January 8, 2018; https://doi.org/10.1073/pnas.1713452115

  1. Edited by Donald G. Truhlar, University of Minnesota, Minneapolis, MN, and approved December 12, 2017 (received for review August 9, 2017)

Significance

Ion desolvation limits the rate of elementary ion-attachment steps during ionic crystal growth from aqueous solution. However, the desolvation mechanism and solvent coordinates are unknown, so crystal dissolution/growth rates remain difficult to compute. This work identifies solvent reaction coordinates for detachment/attachment of Cl− and Na+ from kink sites on NaCl crystals. We find that the larger Cl− limits the overall rate of NaCl dissolution in water because its size prevents water from accessing the kink site during detachment. Our findings show how ion and solvent characteristics influence the kinetics of NaCl growth and dissolution. This work also provides generalizable solvent coordinates and rate expressions for efficient calculations of crystal growth and dissolution rates.

Abstract

Desolvation barriers are present for solute–solvent exchange events, such as ligand binding to an enzyme active site, during protein folding, and at battery electrodes. For solution-grown crystals, desolvation at kink sites can be the rate-limiting step for growth. However, desolvation and the associated kinetic barriers are poorly understood. In this work, we use rare-event simulation techniques to investigate attachment/detachment events at kink sites of a NaCl crystal in water. We elucidate the desolvation mechanism and present an optimized reaction coordinate, which involves one solute collective variable and one solvent collective variable. The attachment/detachment pathways for Na+ and Cl− are qualitatively similar, with quantitative differences that we attribute to different ion sizes and solvent coordination. The attachment barriers primarily result from kink site desolvation, while detachment barriers largely result from breaking ion–crystal bonds. We compute ion detachment rates from kink sites and compare with results from an independent study. We anticipate that the reaction coordinate and desolvation mechanism identified in this work may be applicable to other alkali halides.

Crystal growth rates at mild supersaturations are determined by rates of solute attachment and detachment at kink sites (1, 2), which are special surface sites that resemble half-lattice positions. The barrier for solute detachment is believed to arise mostly from breaking solute–crystal interactions (3⇓–5), but it cannot entirely result from this; otherwise attachment would be a barrierless, diffusion-controlled process. Thus, the attachment barrier (and by microscopic reversibility, a part of the detachment barrier) must result from changes in solvation.

Desolvation is widely accepted as the rate-limiting process for growth of ionic crystals and biominerals; although a number of desolvation hypotheses have been proposed (1, 3, 4, 6⇓⇓⇓⇓⇓⇓–13), our actual understanding of the desolvation process remains largely incomplete. Only a few computational studies focus on attachment/detachment of a solute at a crystal surface (12, 14⇓⇓⇓⇓⇓⇓–21), and only recently have kink sites been examined (22, 23). None of these studies have examined the desolvation process during attachment at kink sites.

In this work, we investigate the attachment/detachment of Na+ and Cl− at kink sites of a NaCl crystal in water. Accurate rates can be computed using only ion positions without a complete understanding of the attachment/detachment mechanism (15, 22) [e.g., via trajectory-based rate calculations (24)], but these calculations provide only numerical rate estimates. They do not provide the kind of mechanistic insight that can be used to construct generalizable theories for other conditions and/or other ions. Therefore, we examine the role of desolvation by identifying important solvent collective variables (CVs) for attachment/detachment at kink sites (25). We elucidate the desolvation mechanism via transition path sampling, likelihood maximization, and committor tests, from which we find that kink attachment can be described with one solute collective variable and one solvent collective variable. By computing the respective free energy surfaces, we gain additional insight on similarities and differences between Na+ and Cl− detachment. Finally, we compute ion detachment rates from kink sites and make comparisons with other studies.

Methods

Simulation Setup.

We use SPC/E (extended simple point charge) (26) water model and the Joung–Cheatham (27) force field for crystalline NaCl as well as Na+ and Cl− ions in solution. As suggested in earlier studies, the Lennard-Jones interactions are cut off at 10 Å and shifted to zero, water geometry is held rigid with SHAKE (28), and long-range electrostatics are included via the particle–particle particle–mesh Ewald method. According to this force field, the solubility of NaCl in water is 3.6 mol/kg water (29⇓–31). A brief table of its other properties can be found in SI Appendix. All properties and trajectories were computed at 298 K in the NVT ensemble (i.e., constant number of particles, volume, and temperature) with a Langevin thermostat in LAMMPS (large-scale atomic/molecular massively parallel simulator) (32). Details on solution equilibration and NPT (constant number of particles, pressure, and temperature) simulations to determine appropriate box dimensions can also be found in SI Appendix.

The system consists of a NaCl crystal slab (854 formula units) with pure water (1,800 molecules) on both sides. The crystal is tilted [similar to that in Klimes̆ et al. (18) and Kolafa (33)] such that each exposed face contains two edges, each with a single kink site, resulting in an ever-kinked surface in periodic space (Fig. 1A). See SI Appendix for details on the crystal rotation. This system design eliminates kink–kink interaction effects because an ion pair can be removed without altering the distance between kink sites. Additionally, the slab thickness (ca. 36 Å) is sufficient to mitigate finite size effects, which can alter computed properties like solubility (34). Therefore, we expect this system to adequately mimic the bulk crystal in the slab interior without fixing atom locations in an inner crystal layer during simulations.

Fig. 1.
Fig. 1.

(A) Depiction of the NaCl crystal in the simulation box (black outline), showing the ever-kinked crystal face in periodic space. Cl− is teal, Na+ is yellow, and water is excluded for clarity. (B–E) Simulation snapshots of ion-attached states (B and D) and ion-detached/solvated kink states (C and E). The detaching ions remain docked in loosely bound states in C and E. The edge and kink are marked with black lines. Most water molecules are omitted for clarity. All simulation images are generated with VMD (Visual Molecular Dynamics) (35).

Reaction Coordinate Identification.

To identify important solvent reaction coordinates, the flexible-length (36) aimless shooting variety of transition path sampling (TPS) is used to separately harvest an ensemble of transition paths for Na+ and Cl− (NVT simulations with a Langevin thermostat). This variant of TPS has one adjustable parameter, δt (time between candidate shooting points), which was set to 30 fs and yielded acceptance rates of 41% for Na+ and 36% for Cl−. As shown in Fig. 1, the ion-attached basin corresponds to small ion–kink separation distances. The ion-detached basin is characterized by larger ion–kink separations and a water molecule in the kink site. See SI Appendix for details regarding basin definitions and determination of the ion–kink separation distance.

We use the inertial version of likelihood maximization (37) (iLmax) to rank the accuracy of an extensive list of trial reaction coordinates of the form q=c1rion,k+c2qsolvent−q‡, where rion,k is the ion–kink separation distance and qsolvent is a solvent CV. The adjustable parameters c1 and c2 determine the direction/gradient of the trial coordinate in (rion,k,qsolvent) space, while q‡ is the predicted value of c1rion,k+c2qsolvent at the optimal dividing surface location. The trial solvent CVs can be generally characterized as distances, local water densities, coordination numbers, joint coordinations, solvent structures, and solvation energy gaps. Parameters within the trial solvent CVs, e.g., cutoffs for computing coordination numbers, are also optimized to maximize reaction coordinate accuracy (38). These simple trial reaction coordinates help to identify a clear role for solvent degrees of freedom in the mechanism.

We chose rion,k as the solute CV not only because it is commonly used, but also because it outperformed (via iLmax) alternative choices of the solute CV (i.e., coordination number, height above surface, etc.). rion,k even outperforms pairs of other solute CVs, such as a coordination number and height above the surface (17). However, rion,k alone is not a good reaction coordinate, because it results in a hysteresis for Cl− (although not for Na+) as shown in SI Appendix. The hysteresis results from poor sampling of the water molecule entering and leaving the kink site.

In iLmax, the reaction probability is pRX=(1+erf[q+cVq˙])/2, where cV is an adjustable coefficient and q˙ is the reaction coordinate velocity, determined from a finite difference after a single NVE (constant number of particle, volume, and energy) time step from the shooting point. The reaction coordinates are scored relative to the baseline rion,k coordinate,Δln⁡L[q]=ln⁡L[q]−ln⁡L[rion,k]BIC[1]where ln⁡L[q] is the log-likelihood for reaction coordinate q and the Bayesian information criteria (BIC) helps determine whether differences between reaction coordinates are significant [BIC=ln(Ntraj)/2=4.95 in this work; Ntraj is the number of TPS trajectories]. The baseline log-likelihoods from iLmax for Cl− and Na+ are −12,402 and −13,204, respectively.

Committor tests are performed at the dividing surface for Na+ and Cl− detachment. Configurations within the transition-state ensemble are collected from umbrella-sampling simulations, with each configuration separated by at least 2 ps. Two hundred configurations are tested, each with 20 trajectories initiated with different Boltzmann distributed velocities (accounting for rigid water molecules), and integrated under NVE conditions until the trajectory enters the ion-attached or ion-detached basin. The committor, pB, for each configuration is the fraction of trajectories which enter the ion-detached (product) state.

Rate Calculations.

Free energy calculations are performed via umbrella sampling (39) with parabolic bias potentials on both rion,k and the solvent CV ρw (discussed in Reaction Coordinate). This technique consists of separate simulations (i.e., umbrella windows) constraining the system to different parts of the transition between the ion-attached and -detached states via an added bias potential. Within each window, data were collected from at least 3 ns of thermostated molecular dynamics simulation. The overlapping distributions from the collection of biased windows and the weighted histogram analysis method (40, 41) yield the overall free energy surface. Note that free energies are normalized by kBT=β−1 in this work, where kB is Boltzmann’s constant and T is absolute temperature.

Rate constants were obtained from generalized transition-state theory (42) using the optimized reaction coordinate q via k=κ[q]kTST[q], where κ is the transmission coefficient that accounts for dynamical barrier recrossing events and kTST is the transition-state theory rate constant. Note that the overall rate constant is independent of q. The effective positive flux method (43) was used to calculate κ. Briefly, NVE trajectories are initiated from a Boltzmann distribution of configurations at the dividing surface. Trajectories are tested to determine whether they cross in the positive (P) direction, whether they are the first (F) crossing, and whether they are effective (E) at reaching the product (detached) state. A PFE trajectory has χ=1 and otherwise χ=0. The transmission coefficient isκ=⟨q˙0χ⟩‡12⟨|q˙|⟩‡,[2]where ⟨…⟩‡ indicates a conditional equilibrium average over the dividing surface. Additional details on the rate constant calculation are provided in SI Appendix.

Desolvation and Kink Attachment Mechanism

Reaction Coordinate.

According to iLmax, the mechanism for solute attachment/detachment at kinks involves a concerted exchange of solute and solvent for the kink site. The optimal reaction coordinate for Cl− and Na+ combines rion,k (solute CV) and the local water density near the kink site (solvent CV),ρw=1(2πσ2)3/2∑wexp(−|𝐫w−𝐫probe|22σ2),[3]where the sum is over water molecules, σ is an effective probe radius (given in Table 1), 𝐫w is the location of a water molecule (best characterized by the oxygen atom), and 𝐫probe is the probe location. We tested different probe locations in the vicinity of the kink site; for Na+, the kink site is the best probe location. For Cl−, the optimal probe location is slightly offset from/deeper into the kink site, which we anticipate to more closely reflect the preferred oxygen location in the kink site (SI Appendix). These differences are likely a function of different water coordination environments and different sizes of the kink sites (due to different ion sizes).

View this table:
Table 1.

Coefficients in the reaction coordinate for Na+ and Cl− attachment/detachment at kink sites as well as the effective probe radii, σ (units of Å)

By first comparing pure-solute and pure-solvent reaction coordinates, we can assess the importance of the solute and solvent (recall Δln⁡L[rion,k]=0). For each ion, the pure-solvent reaction coordinate outperforms the pure-solute reaction coordinate (Δln⁡L[ρw]=11 for Cl− and 34 for Na+). The synergistic benefit of accounting for the solute and solvent results in an increase of Δln⁡L by over 100, as shown in Table 1 along with reaction coordinate coefficients. Consistent with the observed hysteresis for Cl−, but not for Na+ in F(rion,k), the results in Table 1 show that the solvent CV is more important for Cl− than for Na+. Additionally, the significantly nonzero cV coefficients indicate that the dynamics of attachment/detachment of Na+ and Cl− are partially inertial, like the previously reported inertial association of these ions with each other (44).

Committor Tests.

A preliminary analysis of committors on the iLmax dividing surface revealed higher pB values on one end of the dividing surface than on the other. We used the solvent CV identified by iLmax, but adjusted the coefficients c1 and c2 [i.e., tilted the dividing surface in (rion,k,ρw) space] to correct the iLmax prediction (see SI Appendix for details). The committor distributions at the optimized dividing surface for kink attachment/detachment are shown in Fig. 2A for Cl− and Fig. 2B for Na+. Both distributions appear approximately uniform. For an ideal reaction coordinate, the expected distribution is peaked at pB=1/2, with no contributions in the limits of 0 and 1. Therefore, our reaction coordinate is not ideal. However, it is a significant improvement on rion,k (Fig. 2 C and D), indicating that ρw is mechanistically important. It also greatly improves upon the only previously reported committor distribution for a similar process, where few realizations near pB=1/2 were observed (16).

Fig. 2.
Fig. 2.

(A–D) Committor distributions computed at the q (A and B) and rion,k (C and D) dividing surfaces for Cl− (A and C) and Na+ (B and D).

Upon ion detachment and solvation of the kink, the ion is in an adsorbed state near the kink site (Fig. 1 C and E); i.e., the ion is not a fully hydrated “aqua ion” (45). Similar “docked” or partially hydrated ion configurations near CaCO3 step edges have been predicted for CO32− in water (9). We have computed the probability, psol, that the docked ion tends to go into bulk solution rather than reattach to the kink site. We find that psol=1 for Cl− and psol=0.86 for Na+, so detachment is essentially complete for both ions once they reach the docked state. This suggests that diffusion to/from the surface (via surface or from bulk solution) would be a negligible resistance to ion attachment/detachment. Our results confirm that the rate-limiting step is exchange of the ion between the adsorbed and the kink site and that solvent CV is important for the rate-limiting step.

Free Energy Surfaces.

We gain additional mechanistic insight on desolvation and kink attachment from the computed free energy surfaces, shown in Fig. 3A for Cl− and Fig. 3B for Na+ (note they have different color scales). The surfaces are plotted in terms of ln⁡ρw to clearly visualize the impact of ρw. The basins are labeled according to whether the ion or a water molecule occupies the kink site. A high free energy ridge on the free energy surface for Cl− separates the solvated and desolvated kink sites (ca. rCl,k=1.8 Å, lnρw=−5.5 in Fig. 3A). The ridge, which is not present on the free energy surface of Na+, explains why sampling along rion,k gave a hysteresis for Cl− and not for Na+. Committor estimates for configurations in the predicted transition-state ensemble are grayscale points to represent the committor (pB=0 is black and pB=1 is white). The intermingled black and white points on the dividing surface shows that some third component of the reaction coordinate remains unknown. The dotted lines indicate the minimum free energy path (MFEP) between the attached and detached states, as discussed below.

Fig. 3.
Fig. 3.

(A and B) Free energy surfaces for kink attachment/detachment of (A) Cl− and (B) Na+. The minima (0kBT) correspond to the ion in the kink. Points along the dividing surface (black line) are colored in grayscale based on their pB values. The dotted line indicates the minimum free energy path.

Overall, the free energy surfaces show that the kink attachment/detachment processes are similar between the two ions. As the ion detaches (increasing rion,k), a void grows at the kink site (i.e., water cannot immediately enter the kink). Similarly, the ion attachment process is initiated by void growth at the kink, resulting in kink desolvation (decreasing ln⁡ρw). The work required to grow these voids will be reflected in the free energy barriers for ion detachment and ion attachment/kink desolvation.

As equivalent ion–crystal bonds are broken during detachment, the ion-attached basins for Cl− and Na+ appear quite similar. However, the free energy surfaces begin to differ near the dividing surface as solvent motion begins to become important. Additionally, the different dividing surface locations in rion,k show that a larger void must be grown for Cl− detachment than for Na+ detachment. We hypothesize that this occurs because the larger Cl− must move farther away before a water molecule can access the kink site.

The kink desolvation work (i.e., moving to more negative ln⁡ρw in Fig. 3 at fixed rion,k) appears to be markedly different between the two ions, with the Cl− kink being much more difficult to desolvate. We anticipate that this occurs due to the different water coordination environments in the vacant kinks, which consist of a “nest” of three counterions. As water binds more strongly to Na+ than to Cl−, we expect, and our results (Fig. 3) show, that water binds more strongly to the vacant Cl− kink.

By analyzing a sample of TPS trajectories, we find that multiple attachment/detachment pathways exist; i.e., water can vacate the kink in different directions. However, a projection into the (rion,k,ln⁡ρw) space gives a free energy surface that appears to have just one pathway, i.e., along the minimum free energy path. The MFEPs for Cl− and Na+ indicate distinct pathways; however, the MFEPs do not yield dynamical insights. In SI Appendix, we show the evolution of select trajectories on the free energy surfaces. Along a typical detachment pathway, rCl,k increases until a sudden “jump” in ρw as the water enters the void left behind. The rapid change in ρw suggests a nonadiabatic rate formulation like that of Schenter and coworkers (46) for future attempts to develop a general kink attachment/detachment model. The two-step nature of the Cl− attachment/detachment pathway is likely due to water being weakly coordinated to Cl−. Conversely, the pathways for Na+ indicate a more simultaneous motion of ion and water, which occurs because Na+ drags water molecules with it due to its strong coordination to oxygen.

Kink Detachment Rates

The free energy profiles for kink attachment/detachment in Fig. 4 are obtained by projecting the free energy surfaces from (rion,k,ρw) onto q. The transition state is located at q=0 with the ion-attached state at q<0 and the ion-detached/solvated kink state at q>0. We note that projecting the free energy surfaces solely onto rion,k decreases the free energy maximum by 2 kBT for Cl− and 1 kBT for Na+ compared with F‡=F(q=0), as expected from variational transition-state theory (47). The barrier for Cl− detachment is about 3.5 kBT greater than that for Na+ detachment, which we attribute to the previously discussed larger void that must be formed before water enters the kink site.

Fig. 4.
Fig. 4.

Free energy profiles for kink attachment/detachment of Cl− (blue) and Na+ (red). The minimum free energy is set to 0kBT, corresponding to the ion in the kink. Inset shows the committor over the course of an ion-detaching trajectory for Cl− (blue; 4.5-ps duration) and Na+ (red; 35-ps duration). The dashed line indicates the pB=1/2 dividing surface, which is recrossed. Shaded regions represent error bars of pB(1−pB)/n, where n is the number of trajectories used to determine pB.

Detachment rates are computed viak−=12⟨|q˙|⟩‡κe−βΔF‡,[4]where ΔF‡ is the free energy difference between the transition state and the ion-attached stateΔF‡=F‡−[−kBTln(ℓq−1∫−∞0e−βF(q)dq)],[5]where ℓq is a unit length of q. The average reaction coordinate velocities at the dividing surface are ⟨|q˙|⟩‡=0.0024 fs−1 and 0.0026 fs−1 for Na+ and Cl−, respectively. For more details on arriving at Eq. 4, see SI Appendix. We note that k− is formally a detachment rate constant. It is equivalent to the detachment rate because the detachment process is zeroth order in solute concentration (48).

From effective positive flux we find that κNa=0.10±0.01 and κCl=0.25±0.02. These lower transmission coefficients may result from intrinsic friction in the dynamics (i.e., a rigid solvent structure) and/or reaction coordinate error. Alternative reaction coordinates may improve the committor distributions, but we anticipate little improvement upon κ because Fig. 4, Inset shows that even the ideal dividing surface (pB=1/2) is crossed several times along typical detachment trajectories. Thus, the maximum possible transmission coefficient is less than unity.

We find that the Cl− detachment rate (kCl−=6.3×106 s−1) is over an order of magnitude slower than the Na+ detachment rate (kNa−=1.6×108 s−1). Thus, Cl− detachment limits dissolution kinetics. From the direct simulation of NaCl dissolution in water by Lanaro and Patey (49) (using identical force fields), we estimate the average detachment rate of ions from a kink site is 4.5×107 s−1 (deconvolution of rates for Na+ and Cl− is not possible with the reported data). This is in good agreement with the (harmonic) average detachment rate in this work, 1.2×107 s−1. Note that Lanaro and Patey (49) examined dissolution of an approximately spherical NaCl nanoparticle with a rougher surface than our stepped surface. Therefore, a somewhat faster detachment rate is expected from their study.

The dissolution of NaCl has been the focus of many studies (17⇓–19, 49⇓⇓⇓⇓⇓⇓⇓–57), in particular for thin water films that occur at moderate relative humidities. While some works focus on direct dissolution simulations, a few use ab initio molecular dynamics to compute detachment barriers for ions at edge corner sites (i.e., the corner of a 2D nucleus) (17, 19). While Liu et al. (17) found a similar Cl− detachment barrier compared with our results (although the barriers for Na+ differ), this most likely results from a fortuitous cancellation due to the different surface sites, force fields, and reaction coordinates.

Stack et al. (15) determined the Ba2+ detachment rate from a barite step site was approximately 104 s−1. Darkins (23) computed similar detachment rates of CO32− from two different kink sites on calcite. We anticipate that the three orders of magnitude difference from our computed rates largely results from different ion charges, which impacts the bonds that must be broken to detach. In other works, Stack et al. (3) and Stack and coworkers (10) inferred even slower kink detachment rates by fitting models to experimental measurements of CaCO3 and MgCO3 crystal growth. Dogan et al. (22) found kink detachment rates on the order of 105 s−1 via molecular simulations for aspirin and alpha-lactose monohydrate crystals. While electrostatics cannot explain these slower rates, the differences may result from slower rotational motion of these larger molecules and/or the more collective nature of breaking many atom–atom interactions. It is not yet clear how different solute sizes, charges, shapes, and solvent characteristics impact attachment/detachment rates.

Conclusion

Ionic crystal growth and dissolution rates near saturation in aqueous solution are limited by desolvation and attachment of solute growth units at kink sites. We have examined the ion attachment mechanism in NaCl crystal growth/dissolution with path sampling and reaction coordinate identification techniques. The mechanism involves a competition between solute and solvent for the kink site. Likelihood analyses and committor distributions show that the local density of water molecules in the kink site is a mechanistically important solvent CV. Free energy surfaces, transition-state theory rates, and transmission coefficient corrections provide detachment rates that agree with independent estimates from direct simulations of dissolution (49). While the absolute dissolution rates may be sensitive to force-field details, mechanistic conclusions are expected to be more robust. For example, Cl− detachment limits the overall kinetics because, compared with Na+, the larger Cl− must move farther away from the kink site before water can solvate the kink. The optimized solvent CV identified in this work is similar to a solvent CV that was accurate for ion-pair dissociation across several alkali halides (44, 58). Thus, the reaction coordinates identified here may also be useful for understanding other alkali halide crystals. More generally, with the strategy used in this work it should be possible to fully understand how ion sizes, charges, shapes, and interactions with the solvent influence the rates of attachment and detachment at kink sites in ionic crystal growth and dissolution.

Acknowledgments

We thank David Quigley for stimulating discussions. We are grateful to our funding sources, the US National Science Foundation (NSF) (Awards CBET-1335694 and CHE-1465289) and Eli Lilly and Company. We acknowledge support from the Center for Scientific Computing from the California NanoSystems Institute (CSNI) and Materials Research Laboratory (MRL): NSF Materials Research Science and Engineering Center (Grant DMR-1121053) and NSF Grant CNS-0960316.

Footnotes

  • ↵1To whom correspondence should be addressed. Email: baronp{at}ucsb.edu.

  • Author contributions: M.N.J., M.F.D., and B.P. designed research; M.N.J. performed research; M.N.J. analyzed data; and M.N.J., M.F.D., and B.P. wrote the paper.

  • The authors declare no conflict of interest.

  • This article is a PNAS Direct Submission.

  • This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1713452115/-/DCSupplemental.

Published under the PNAS license.

View Abstract
Back to top
Article Alerts
Email Article
Citation Tools
Request Permissions
Dissolution kinetics at NaCl kink sites
Mark N. Joswiak, Michael F. Doherty, Baron Peters
Proceedings of the National Academy of Sciences Jan 2018, 115 (4) 656-661; DOI: 10.1073/pnas.1713452115
Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo
Proceedings of the National Academy of Sciences: 115 (4)
Table of Contents

Submit

Article Classifications

Jump to section

You May Also be Interested in

Large piece of gold
News Feature: Tracing gold's cosmic origins
Astronomers thought they’d finally figured out where gold and other heavy elements in the universe came from. In light of recent results, they’re not so sure.
Image credit: Science Source/Tom McHugh.

Similar Articles