University of Birmingham Sodiation and Desodiation via Helical Phosphorus Intermediates in High-Capacity Anodes for Sodium-Ion Batteries

.


■ INTRODUCTION
Na-ion batteries (NIBs) are promising alternatives for longterm sustainability in terms of both cost and natural abundance compared to Li-ion systems. 1−5 Na is widely available and evenly distributed worldwide, lessening the political tensions that may arise from continued Li use.While numerous viable cathode materials for NIBs have been identified from Li analogues, 6,7 many intercalation and alloying anode materials that work well for Li-ion batteries (LIBs) fail for Na chemistries.For example, the larger atomic size of Na does not allow intercalation into graphite 8 and the safety issues associated with microstructural formation are exacerbated when Na metal is used. 9−13 Due to the high natural abundance of P, the use of Na−P chemistries has the potential to further drive down the cost of NIBs while achieving performances comparable to those of LIBs.Such high performing phosphorus-containing batteries have potential use in stationary systems, such as grid storage; we note, however, that the use of phosphorus anodes is associated with a safety riskreduced phosphides can react with water, for example, to produce phosphine gas.However, with an emphasis on battery safety, monitoring tools, and quality assurance, the likelihood that these risks can be mitigated increases.
−21 Further improvements were achieved by using single-layer black P, i.e. phosphorene, in combination with graphene, which displayed a reversible capacity of 2440 mA h g −1 in the first cycle. 22−33 Interpretation of electrochemical signatures relies on a thorough understanding of the chemistry that occurs on sodiation and desodiation.Thus, tracking the (de)sodiation pathways provides a platform with which to elucidate the structure−property relationships responsible for high capacities as well as to devise strategies to mitigate the performance degradation.Amine and co-workers 19 used ex situ and in situ synchrotron high-energy X-ray diffraction (XRD) to study electrochemical cycling of black P, and the only crystalline phases that were observed were the pristine material (orthorhombic black P) and crystalline Na 3 P (c-Na 3 P).No Bragg reflections were observed for any intermediate Na x P (0 < x < 3) compositions studied, indicating the lack of long-range ordering, consistent with the formation of amorphous/ metastable phases (a-Na x P; 0 < x < 3).Unlike other alloying anodes used in NIBs (e.g., Sn or Sb), P is a relatively light element and thus total scattering techniques such as pair distribution function (PDF) analysis of X-ray data, (which have been successfully used to aid structural assignment of both Na x Sn 34 and Na x Sb, 35 both systems where disordered phases are observed along with crystalline phases) are difficult to model and interpret due to the relatively poor signal-to-noise ratio. 36Conversely, NMR is well-suited to study the shortrange order in the amorphous/metastable structures in a-Na x P due to the good sensitivities and 100% natural abundances of the NMR-active 31 P and 23 Na nuclei.
The Na−P composition space has previously been explored by some of the authors by using density functional theory (DFT) calculations to evaluate the formation energy of experimentally known Na x P structures and new structure predictions derived from ab initio random structure searching (AIRSS) and structural prototyping (species swapping from known phases of similar compounds). 38In similar systems, such ab initio structure searches have previously determined crystalline products not found in structure databases, 34 as well as identifying possible structural motifs present in disordered phases formed on cycling. 39By the nature of AIRSS and prototyping, there is, however, no information flow between calculations, and thus these techniques are limited to the initial stoichiometries chosen at their outset.
Here, we show that by using a combination of solid-state NMR (ssNMR) and powder X-ray diffraction (PXRD) we can monitor the Na x P phases that are formed on (de)sodiation in black P anodes.Specifically, individual P binding sites within the amorphous Na x P phases that arise during cycling are probed with 2D 31 P phase-adjusted spinning sideband (PASS) ssNMR experiments.In 2D PASS, isotropic shifts are correlated to the corresponding anisotropic chemical shifts, facilitating analysis of the full chemical shift (CS) tensors of individual P sites and ultimately providing information on the geometry of chemical bonding environments and structures of molecular fragments. 37The crystalline structures that occur are characterized by 31 P magic-angle spinning (MAS) ssNMR, variable field 23 Na MAS ssNMR, and XRD.All experimental NMR and XRD data are interpreted through first-principles computational structure prediction.
In this current study, we move beyond the initial constrained AIRSS + prototype search and use a specialized variablecomposition genetic algorithm (GA) code to interpolate to stoichiometries not sampled by the stochastic searches.Combining this theoretical approach with ex situ experimental analyses, we determine how P fragments in amorphous phases change during sodiation vs desodiation and show that the structure of the end member phase, Na 3 P, is a phase previously computationally predicted using the AIRSS technique. 38

■ RESULTS
Genetic Algorithm (GA) Sampling of Theoretical Predictions of Na x P Structures.The GA was initiated from low-lying AIRSS-derived phases to avoid the need to extrapolate to entirely new regions of composition space.Evaluating the fitness of child structures as a function of distance from the convex hull, the GA uncovered several motifs from less ordered phases that could not be sampled with any individual technique alone.These new motifs provided a vast library of possible P binding environments to evaluate as potential local environments and motifs that may arise in the amorphous Na x P structures formed on cycling.These GAderived phases were added to the previously constructed 38 convex hull of thermodynamically stable Na x P structures in Figure 1a.The pristine material, black P, is a layered structure that contains six-membered P rings.The known phosphorus-rich polyphosphide of NaP 7 40 was considered and found to lie on the convex hull.In NaP 7 , one Na coordinates to each [P 7 ] subunit−phosphorus motifs that form infinite tubes arranged in helices about the c axis. 40In previous work, Morris and coworkers found a series of structures containing P cages (Na 3 P 11 and Na 3 P 7 exhibit P 11 and P 7 cages, respectively) (where 0 < x < 0.43 in Na x P) that fall on the convex hull. 38As the sodium content increases, the next on-hull phases are NaP and Na 5 P 4 , which exhibit P helices and four-member P zig-zags, respectively.At the end of sodiation, isolated P ions are found in Na 3 P.
Electrochemical Characterization and Comparison with Calculated Potentials.The theoretical structure models (with energies depicted in Figure 1a) were compared to Na x P phases that arise during cycling of black P in NIBs.Black P electrodes were assembled into Na half-cells for electrochemical measurements.Near-theoretical capacity was reached on the first discharge (2510 mA h g −1 , based on the mass of P alone) and one of the highest reversible capacities observed to date of 2074 ± 80 mA h g −1 without the use of additives is achieved in the first cycle.The improved capacity may be due to the use of 1.0 M NaFSI in 2-MeTHF as the electrolyte, which differs from the electrolytes used previously (1.0 M NaClO 4 in EC/DMC or 1.0 M NaPF 6 in PC).When compared to PC, methyl acetate, or THF, 2-MeTHF shows improved stability against Li, 41 and similar stability is likely present toward Na, thus improving the electrochemical performance of the NIBs.The experimental P reduction profile was compared to the average voltages for on-hull Na x P structures predicted by Morris and co-workers 38 (Figure 1b).When a constant current is applied (C/100), the voltage drops from the open circuit voltage to 1.10 V (Figure 1b).This initial voltage drop suggests that P 11 and P 7 subunits, such as those found in NaP 7 and the computed structures of Na 3 P 11 and Na 3 P 7 , may not form electrochemically during sodiation in a NIB.Their absence may be due to the unfavorable energetics involved with breaking the sixmembered P rings in black P and re-forming P−P bonds to form the P 11 and/or P 7 cages in these structures.Instead a sloping region is seen between 1.10 and 0.63 V, corresponding to 31−420 mA h g −1 and approximate compositions of Na 0.04 P−Na 0.49 P. A steep, but short-lived sloping region is observed from 1.10−0.48V (420−562 mA h g −1 , approximate compositions of Na 0.50 P−Na 0.65 P); which is followed by a more shallow, sloping region from 0.48−0.22V (562−1358 mA h g −1 , approximate compositions of Na 0.65 P−Na 1.57 P), implying gradual sodiation of Na x P (where x ranges from 0.04 to 1.57) from 0.63 to 0.22 V. Following this, a plateau at 0.21 V (ranging from 1358−2340 mA h g −1 , where x corresponds to 1.57−2.70) is observed with a distinct peak in the dQ/dV plot (Figure S2), which is attributed to the formation of c-Na 3 P. 21 Ex Situ 31 P MAS NMR Spectra on Sodiation.Ex situ 31 P MAS NMR spectra were acquired at different stages of (de)sodiation to probe the structural features of various Na x P phases (Figure 2) directly.The 31 P MAS NMR spectrum collected in the first sloping region at 0.90 V (124 mA h g −1 , approximate composition, Na 0.14 P) is nearly identical to that of the pristine material, with the exception that a weak peak at −227 ppm emerges (see Figure S8 for a magnified version of Figure 2b).Assignment of this low frequency resonance is not obvious based on its isotropic 31 P shift alone.However, we note that peaks in this region (ca.− 220 to −250 ppm) persist throughout all stages of (de)sodiation, suggesting this region arises from Na x P species formed in side reactions, minor Na x P intermediates or perhaps Na near P defects (such as the end of a P chain).
The 31 P MAS NMR spectrum for the sample extracted at 0.60 V (approximate composition Na 0.52 P) spans a broad chemical shift range from the resonance of the pristine material (14 ppm) to the small, sharp peak at −224 ppm (Figure 2b).Two pronounced broad peaks with centers-of-mass at −104 and −143 ppm are clearly resolved, these peaks becoming more pronounced as sodiation proceeds (Figure S8).The corresponding ex situ XRD pattern of this material does not show any Bragg reflections that are distinct from the pristine material (Figure S13) indicating that an amorphous or highly disordered phase is present in addition to residual black P. At 0.38 V (780 mA h g −1 ; approximately NaP; x = 1), a broad set of resonances spanning a 31 P chemical shift region from ca. −260 to 40 ppm is again observed, along with the lack of Bragg reflections (Figure S13), suggesting that this phase is also amorphous (a-Na x P).
Following the formation of a-Na x P, a sharp 31 P resonance (fwhm = 10 ppm) emerges at −207 ppm (Figure 2b) that is coincident with the plateau at 0.21 V in the electrochemical profile and is assigned to c-Na 3 P on the basis of previous 21 (and current; see Supporting Information (SI)) ex situ XRD characterization at the corresponding state of charge.The set of broad 31 P resonances decrease in intensity near the end of this plateau (1666 mA h g −1 , approximate composition of Na 1.93 P) consistent with a two-phase reaction of a-Na 3−x P to c-Na 3 P (where x is about 0.6; approximately composition a-Na 2.4 P), and vanish by the end of sodiation at 0.01 V (2510 mA h g −1 , approximate composition of Na 2.90 P).Here x was quantified based on the capacity and integrating the area corresponding to the shifts of black P (14 ppm), c-Na 3 P (−207 ppm), and assuming the integral of the remaining broad resonances corresponds to a-Na 3−x P).Assignment of the Crystal Structure of c-Na 3 P Formed at the End of Sodiation.Previously, the crystallographic structure of Na 3 P formed on discharge has been assigned as Na 3 P-P6 3 /mmc on the basis of both powder 42 and single crystal 43 XRD characterizations of solid state syntheses of Na 3 P.However, our prior DFT calculations predicted a P6 3 cm symmetry phase to be the most stable crystal structure of c-Na 3 P (i.e., located on the convex hull). 38This phase was constructed via species swapping from the Na 3 As structure 44 and can be rationalized as a distorted variant of a previously reported crystal structure of Na 3 P-P6 3 /mmc 42,43 which resides 0.005 eV•f.u.−1 above the DFT hull. 38These two structural models (Figure 3a) are closely related by a group−subgroup transition of the order 6 (Figure S17).In the P6 3 cm structure, the PNa 5 trigonal bipyramids are distorted and tilted, which can be interpreted as frozen vibrational motion of the Na 3 P-P6 3 / mmc crystal structure, consistent with the anisotropic displacement parameters.The simulated XRD pattern of Na 3 P-P6 3 cm shows additional reflections at 31.3°and 33.3°(Cu Kα) that are absent in the otherwise essentially identical pattern of Na 3 P-P6 3 /mmc.Examination of the XRD of Na 3 P samples derived by sodiation of black P in the NIBs revealed similar superstructure reflections and thus lower symmetry (Figure 3c).In addition, we find that Na 3 P synthesized using standard solid-state synthesis techniques in our laboratory also exhibits these superstructure reflections (Figure S16), consistent with the formation of Na 3 P-P6 3 cm.A Rietveld refinement starting from the calculated structure model was performed on both samples.The refined atomic parameters for the sample synthesized via solid state synthesis differ only marginally from those calculated with DFT (Table S2, Figures S14 and S16).The slightly larger deviation observed for Na 3 P formed electrochemically is attributed to the broadening of the Bragg reflections and thus lower signal-to-noise ratio of the data (Figure S15).
Both structural models (P6 3 /mmc and P6 3 cm) feature only one independent P position, with isolated P ions.DFT calculations of the corresponding 31 P NMR shifts appear in a similar spectral region (−180 vs −217 ppm for P6 3 cm vs P6 3 / mmc, respectively) and fall on either side of the experimentally observed value (−207 ppm).Furthermore, both crystal structures show very little anisotropy based on computed CS tensor values (anisotropy, Δ = 38 ppm and asymmetry parameter, η = 0.41 for P6 3 cm; Δ = −39 ppm, η = 0 for P6 3 /mmc), making 31 P NMR an inconclusive diagnostic for structural determination of this species in the NIBs.However, the reduction of symmetry from P6 3 /mmc to P6 3 cm results in four instead of two crystallographically independent Na environmentsall of which exhibit characteristic 23 Na quadrupolar coupling constants and asymmetry parameters (Table 1, C Q and η, respectively).Comparison of the 23 Na MAS NMR of the active materials at the end of sodiation (0.01 V, 2510 mA h g −1 ) to simulations of the 23 Na spectra for both structures reveal significantly better agreement with the Na 3 P-P6 3 cm structure (Figure 3b), supporting the assignment of the formation of Na 3 P-P6 3 cm in NIBs.
Structural Assignment of P Motifs in Amorphous Na x P (x = 1) During Sodiation.With assignment of the crystallographic structure of the final discharge product in hand, the structures of the amorphous intermediates that form during (de)sodiation are now investigated in order to assign the observed 31 P NMR resonances to specific local structures.[N.B. 23 Na MAS NMR spectra were acquired for all the samples studied (Figure S6) and compared with DFT calculations of on-hull structures, but they were less informative than the 31 P spectra (see SI).] NMR parameters for energetically low-lying  compositions corresponding to stoichiometries at or near Na 1 P 1 (NaP, Na 5 P 6 , Na 5 P 4 , Na 7 P 8 , Na 8 P 7 ) were first calculated (Figure S12).The structural diversity computationally accessed in this composition range from known, AIRSS, prototype, and GA-derived structures allows us to sample the 31 P site variation that may be observed in the amorphous Na x P phases more fully.Unfortunately, the 31 P isotropic shifts calculated for these structural models were associated with a narrow spectral range, or the materials themselves did not exhibit a bandgap (which did not allow us to calculate accurate values for δ iso because the Knight shift contribution that describes the influence of free carriers could not be computed), prohibiting assignment based on isotropic shift alone.However, we found that individual 31 P sites within the structural models exhibited distinct differences in chemical shift anistropies (CSAs).Therefore, local 31 P environments in the a-NaP phase were examined in greater detail with 2D 31 P PASS experiments to separate isotropic and anisotropic chemical shifts and thus enable analysis of the full CS tensors of individual P sites.2D PASS is particularly useful here, not only because of ambiguity in the isotropic 31 P dimension, but these spectra are also broadened by strong intramolecular 31 P− 31 P dipolar interactions 45 − potential interferences that are eliminated in the anisotropic dimension in PASS.Deconvolution of the 31 P isotropic projection shows the presence of (at least) three environments with associated 31  The 31 P site at δ iso = −29 ppm agrees well with the CS tensor found for one crystallographically independent P site (P2 of P1−P2) in a metastable phase of NaP + 96 meV above the hull (Figure 4b, Table 2), which contains P arranged in a 3-fold helix that winds through the Na network.(In general, we find that all P atoms that comprise helical fragments are bound to two P atoms and surrounded by four to six Na atoms.)The P1 site of this phase is associated with a chemical shift of −57 ppm, and if present, would likely be obscured by the much more intense δ iso = −78 ppm resonance.Similarly, the 31 P environment at δ iso = −78 ppm agrees with the CS tensor found for the third independent P site (P3, of P1−P7) in NaP + 75 meV above the hull that also contains P helices (Figure 4c, Table 2) but has very different CSA parameters.The 31 P site at δ iso = −143 ppm closely resembles that of the fourth independent P site (P4, of P1−P8) in a second NaP structure + 96 meV above the hull, found via GA (Figure 4d and Table 2), this site corresponding to a terminal unit of a four-member P zig-zag motif, suggesting that this site represents a terminal  2.
unit of a P chain (i.e., site of P−P bond cleavage).Although there are exceptions, Δ is generally larger for sites at the end of P chains or in small clusters, in all the calculated structures examined, again supporting this assignment.It is important to stress that we do not believe that the various crystalline structures are formed, simply that CSA parameters consistent with terminal and helical structures appear to be present in the amorphous phase.We note that (in contrast to the one pulse spectra (Figure 2b) a well-defined peak for δ iso = −104 ppm is not observed in the isotropic projection from 2D PASS.The absence of this and other possible resonances may be due to short spin−spin (T 2 ) relaxation (all sites show T 2 values <1 ms) with respect to the rotor period (100 μs) that significantly lowers the signal over the course of the five consecutive π pulses that are applied during the PASS experiment.Thus, while we can identify individual 31 P sites from different structural models that show similar CS tensors to the experimental PASS data, the remaining sites in these helical units may also not be observed due to fast T 2 relaxation and subsequent dephasing and/or overlap in the isotropic projection.
Similar, broad 31 P shifts near −29, −78, and −143 ppm remain in the 31 P MAS NMR spectra during further sodiation (Figure 2b).Specifically, the resonance near −143 ppm increases in intensity relative to the other broad peaks (Figure S8), suggesting further fragmentation of P helices as more Na is added to the system.Analysis of the chemical shift tensor for the site at −235 ppm at 0.38 V (780 mA h g −1 ; approximately NaP; x = 1, vide inf ra) reveals a small anisotropy value (Δ = 93 ppm, Figure S11) that is similar to the isotropic 31 P shifts and anisotropies associated with isolated P ions (e.g., δ iso = −180 ppm, Δ = 38 ppm for P6 3 cm and δ iso = −217 ppm, Δ = −39 ppm for P6 3 /mmc).
Structural Assignment of P Motifs in a-Na x P (x = 1) during Desodiation.During desodiation, the c-Na 3 P-P6 3 cm 31 P resonance at −207 ppm persists during desodiation until >0.60 V (>1490 mA h g −1 , approximate composition of Na 1.28 P), allowing the assignment of the disappearance of c-Na 3 P to the peak at 0.60 V in the dQ/dV plot (Figure S2).2D 31 P PASS experiments were again performed on the NaP (x = 1) stoichiometry sample (0.80 V, 1730 mA h g −1 ), for comparison to the P environments in a-NaP formed on sodiation.The 31 P peak at δ iso = −143 ppm exhibits lower intensity in the isotropic projection (Figure S9), compared to the MAS spectrum (Figure S8), indicating a change in the T 2 relaxation behavior of this 31 P site on desodiation and/or overlap with a site that exhibits shorter T 2 relaxation.The sideband pattern from 2D PASS for δ iso = −143 ppm could not be reliably simulated with a single CSA pattern (Figure S10); this may be due to overlapping 31 P sites at this frequency and/ or increased structural disorder (either motional or static).These data suggest that although a clear 31 P peak is observed at −143 ppm in the NMR spectra collected at 0.60 and 0.80 V (1490 and 1730 mA h g −1 , Figure 2b) on desodiation, this may be due to different P sites that exhibit fortuitous overlap at this chemical shift and/or that P atoms at the end of a chain may exhibit more variation in coordination environment.
The remaining 31 P sideband patterns from 2D PASS experiments on a-NaP (0.8 V, 1730 mA h g −1 , approximately NaP) during desodiation can again be correlated with P motifs in the structural models (Figure 5 and Table 3), the environments with δ iso = −29 ppm and −78 ppm again showing reasonable agreement with NaP structures containing helices.The CSA parameters, while similar to those observed   3.
on sodiation, are closer to those of structures with longer average Na−P interatomic distances (see SI for an expanded description of the structure models).Increased structural disorder on desodiation may be due to a larger range of clusters that are formed from the isolated P ions in Na 3 P, as opposed to the original P phase, which already contains P−P bonds.At the end of desodiation at 2.5 V, a 31 P shift (2 ppm) close to that of the pristine material (14 ppm) reemerges but is present alongside broad 31 P resonances that span from 2 to −240 ppm (Figure 2b), indicating that the Na ions are not completely removed during desodiation and/or are consumed in deleterious side reactions.Ex situ XRD collected at the end of sodiation shows that black P is not re-formed (Figure S13).Rather, only one reflection at 26°that is also seen in c-black P is observed, suggesting that the material becomes increasingly amorphous during cycling.

■ DISCUSSION
From the experimental and theoretical data presented, we propose a mechanism for sodiation/desodiation in phosphorus anodes for NIBs. 31P NMR spectra acquired during the initial stages of sodiation (Na x P, x = 0.14) are similar to those of the pristine material, consistent with DFT calculations performed by Hembram et al. 47 which showed that the P−P bonds remained intact at low of sodiation (x < 0.25).Only the interlayer distances increase, allowing Na ion intercalation prior to alloying and fragmentation of black P. The P 11 and P 7 cages present in on-hull Na x P phases are not observed here, likely due to the energetic penalty associated with breaking/re-forming P−P bonds required to make P cages.Importantly, our 31 P NMR spectra at early stages of sodiation do not resemble those reported for NaP 7 , which shows four 31 P resonances at −32.9, 20.6, 30.5, and 68.5 ppm that arise from the four independent 31 P sites in this compound. 40Once the composition reaches approximately Na 0.52 P, 31 P NMR shifts consistent with P helices/zig-zags are present, suggesting that at this stage, the initial P−P cleavage has occurred, and P helices are the predominant structural motifs that are subsequently formed.P helices and zig-zags (which are likely truncated helices) may allow some of the original alternating planar/staggered P−P bonds found in pristine black P to persist while also forming an energetically favorable Na x P phase.
As black P is further sodiated, the amorphous Na x P phase that forms at 0.38 V (780 mA h g −1 , approximate composition of NaP) shows 31 P sites consistent with P helices and terminal P units from short P chains (four-member P zig-zags).The calculated voltage profile is in agreement with the measured experimental profile for x > 0.5, again confirming that the structural (helical) motifs found in structures such as NaP and Na 5 P 4 close to the convex hull are likely present.
Crystalline Na 3 P emerges during a two-phase reaction that occurs below 0.22 V, suggesting that isolated P ions are derived from the P helices and shortened P chains (which increase in intensity in samples extracted during the two-phase reaction, Figures 2b and S8).A new crystalline architecture of Na 3 P (P6 3 cm) was discovered from species swapping with Na 3 As that is more thermodynamically stable than the previously determined Na 3 P-P6 3 /mmc structure found in two independent reports. 42,43Theoretical calculations, variable field 23 Na NMR, and XRD are consistent with the formation of this phase in a NIB and in solid-state syntheses performed in our laboratory.Crystalline Na 3 P persists in the batteries during desodiation until a potential >0.6 V (>1490 mA h g −1 ) is reached.In the a-NaP phase accessed during desodiation at 0.80 V (1730 mA h g −1 ), more disordered P helices are found, perhaps due to the increased structural disorder associated with forming helices from isolated P ions in c-Na 3 P-P6 3 cm.Black P does not re-form on Na removal from the system, and this irreversible chemistry may play a role in the poor capacity retention.
The prevalence of P helices in a-NaP during sodiation/ desodiation indicates that this motif is ubiquitous in amorphous/metastable Na x P structures.Four-fold P helices are present in c-NaP, 48 which falls on the convex hull tie-line, suggesting that this structural motif is conserved in the amorphous analogue, where more disorder and thus, lower symmetry is found.Calculations show that the higher symmetry c-NaP is associated with similar isotropic and anisotropic shifts for both independent 31 P sites in the structure (P1, δ iso = −179 ppm, Δ = 288 ppm, and η = 0.40; P2, δ iso = −145 ppm, Δ = 292 ppm, and η = 0.59), neither of which accurately represent the experimental data, suggesting that lower symmetry/more disordered helices are present.DFT calculations by Hembram et al. 47 suggested that P dumbbells are formed at a composition of Na 0.28 P that continue to increase in quantity until individual P ions are observed in Na 3 P.In other Li-39,49−51 and Naion 34,35 Sn, Ge and Si alloying anodes, dumbbell motifs are commonly observed in the amorphous phases formed on cycling.Conversely, our calculations show that the lowest lying structure of Na x P (x = 2) that exhibits dumbbells falls +54 meV above the hull, indicating that dumbbell-containing structures are less stable in the Na−P system.To determine whether dumbbell motifs were also formed during the (de)sodiation of black P, we compared the experimental CS tensors to those found for P−P dumbbells in Na 2 P models with calculated energies ranging from +54 to +199 meV above the hull: while dumbbell motifs are associated with similar asymmetry parameters and sign of anisotropy for some of the 31 P sites found on sodiation, the magnitude of the anisotropy is consistently smaller for the dumbbell motifs than those observed experimentally and the isotropic shift values deviate dramatically (e.g., experimental values for δ iso = −29 ppm, Δ = 179 ppm, and η = 0.72 whereas a dumbbell motif exhibits calculated values of δ iso = −246 ppm, Δ = 171 ppm, and η = 0.68, ultimately indicating better agreement with the P site in helices with δ iso = −21 ppm, Δ = 171 ppm, and η = 0.55, Table Table 3 S1).The lack of experimental evidence combined with the unfavorable formation energies of Na 2 P structures indicates that dumbbell-containing structures do not form in a measurable quantity in these NIBs during cycling.Conversely, the CS tensor for δ iso = −143 ppm is consistent with the terminal P atom in a four-member zig-zag motif, suggesting that this environment may represent the P−P cleavage site of a helix/zig-zag.
■ CONCLUSION Overall, we followed the structural transformations that occur during cycling in high capacity black P anodes for NIBs with solid-state NMR, XRD, and DFT calculations, providing insight into the (de)sodiation pathway.In tracking the various Na x P intermediates that formed on cycling, the discovery of new Na x P structures via the GA was crucial to assigning the P motifs present in a-Na x P intermediates and AIRSS + prototyping structure searches allowed the assignment of the final discharge product.
During the first sodiation, P atoms at the end of a chain are observed as early as 0.60 V (Na 0.52 P), indicating that P−P cleavage begins at this composition.In a-Na x P, P motifs that correspond to both P helices and the terminal unit in P zig-zags are observed, suggesting a sodiation mechanism that is distinct from other alloying anodes including other group 15 elements, e.g.Sb, 35 in which dumbbells are present.Once the potential drops below 0.22 V, the formation of a lower energy c-Na 3 P-P6 3 cm structure, which was predicted to form a priori, appears and persists during desodiation until >0.60 V.By comparing the product formed at the end of sodiation in NIBs with traditional solid-state synthesis, we find that the crystalline architecture is Na 3 P-P6 3 cm in both cases, indicating that both routes access the thermodynamically favorable structure.Slight differences in P environments in a-NaP during sodiation/desodiation may result from different pathways that are accessed to form extended P structures from the pristine black P vs the isolated P atoms in c-Na 3 P-P6 3 cm, respectively.At the end of desodiation, we find that black P is not re-formed, which may contribute to the poor capacity retention in this material.The combined approach of analyzing both experimental and theoretical CSA can be extended to understand other structural transformations on (de)lithiation/(de)sodiation in not only phosphoruscontaining materials, but also systems, such as Si, where distinct structural motifs, such as clusters and chains, play an important role in the battery chemistry of these largely amorphous phases.
Electrode Preparation.Red phosphorus powder (400 mg, 13 mmol) was placed in a ZrO 2 jar in an Ar-filled glovebox and ball-milled using a SPEX 8000 M mixer/mill for 4 h to produce crystalline black phosphorus.To make the electrodes, black P was combined with conductive carbon for a final mass ratio of P:C = 50:50 (where C = 15 graphite:10 carbon nanotubes:10 carbon super P:15 NaCMC).First, black P was mixed with graphite, carbon nanotubes, and carbon super P in a ZrO 2 jar and ball-milled for 1 h under Ar.After ball-milling, NaCMC in water was added to the P:C mixture to form an aqueous slurry which was cast onto etched Cu foil with a 150 μm doctor blade.The resulting film was dried in air overnight, followed by drying at 100 °C under vacuum for an additional 12 h.Supported black P electrodes for assembly in coin cells were made using a 12.7 mm diameter circular punch and used for all electrochemical measurements and ex situ characterization studies with NMR and XRD.Typical loading values of active material were 1−3 mg of black P per electrode.
Electrochemical Measurements.Black P electrodes were assembled into Na half-cells using R2032 coin cells in an Ar-filled glovebox.The electrolyte used for all measurements was 1.0 M NaFSI dissolved in 2-MeTHF.Black P and metallic Na were separated with a glass fiber separator (Whatman, GF/A).All electrochemical measurements were performed using a Biologic MPG2 battery cycler.Cells were galvanostatically cycled in the potential range of 2.5 to 0.01 V at a rate of C/100, unless otherwise noted.C-rate is defined as the inverse of the number of hours it takes to reach a defined theoretical capacityin this case, 2596 mA h g −1 .
Solid-state NMR Spectroscopy.For ex situ NMR measurements, black P electrodes were cycled against Na in half-cells to various states of charge.Batteries were disassembled in an Ar-filled glovebox and active materials were scraped off the Cu foil and washed twice with 2-MeTHF before drying under vacuum for 30 min.The dried powder (approximately 3−8 mg of material) was packed into ZrO 2 rotors for analysis with 31  For these experiments, material collected from a single battery was center-packed using Teflon to fill the void space. 31P and 23 Na shifts were externally referenced to solid ammonium dihydrogen phosphate at 0.8 ppm and solid sodium chloride at 7.2 ppm, respectively.
Solid-state 23 Na MAS NMR measurements were performed on Bruker Avance III 300 and 700 MHz spectrometers (operating at 7.05 and 16.4 T, respectively) equipped with a 1.3 mm HX MAS probehead. 23Na MAS NMR experiments were collected at MAS frequency = 60 kHz using a rotor-synchronized quadrupolar echo sequence and a recycle delay of 5 s (≥5 × T 1 ).Solid-state 31 P MAS NMR measurements were performed on a Bruker Avance III 300 MHz spectrometer (7.05 T) equipped with a 1.3 mm HX probehead.Rotor-synchronized Hahn echo 31 P MAS NMR experiments were collected at MAS frequency = 60 kHz with a recycle delay of 10 s (at least 5 × T 1 , T 1 values ranged from hundreds of ms to approximately 2 s) [N.B.A recycle delay of 60 s was used for the pristine black P due to the longer spin−lattice relaxation time (T 1 = 9.5 s) observed for this sample.]2D 31 P PASS NMR experiments were performed on a Bruker Avance III 500 spectrometer operating at 11.74 T and equipped with a 4 mm HXY MAS probehead.Spinning speeds of 10 kHz and recycle delays of 10 s were used for all PASS experiments.The PASS pulse sequence consists of a π/2 pulse followed by a train of five rotorsynchronized π pulses with interpulse delays that satisfy the PASS equations. 52For each experiment, the indirect dimension was incremented in 16 steps, with either 1200 or 2400 scans collected per step.All 2D data sets were processed by repeating the 2D signal in the indirect dimension eight times to separate spinning sidebands.All spinning sideband patterns were simulated in dmfit 46 to extract the principal components (δ 11 , δ 22 , and δ 33 ) of the chemical shift tensor for comparison to DFT models.All chemical shift tensors reported here use the Haeberlen−Mehring convention 53 to define the isotropic chemical shift (δ iso ), anisotropy (Δ), and asymmetry (η) as follows: In the Haeberlen convention, the principal components are defined as follows: X-ray Diffraction.The air-sensitive Na 3 P powder samples were finely ground in an agate mortar, filled into 0.3 mm diameter glass capillaries, and sealed with two-component glue inside an Ar-filled glovebox.Diffraction data was collected at 296(2) K on a Panalytical Empyrean diffractometer equipped with a Ni filter using Cu Kα radiation (λ = 154.06pm, 154.43 pm).Rietveld refinement was performed using the TOPAS Academic software package (version 4.1). 54Due to inadequate data quality, the displacement parameters of the atoms were fixed to U iso = 0.013 pm 2 .In addition, the z parameter of one atom needed to be fixed owing to strong correlation.After individual refinement of each z parameter while one other was fixed, z(P1) was fixed for the last iteration.Interatomic distances should be treated with care due to the data quality.Crystal structures were visualized using the Diamond software package. 55ata in brief for Na 3 P at 296(2) K: a = 861.10(2)pm, c = 881.88(2)pm, V = 566.30(2)•10 6 pm 3 , P6 3 cm (no.185), Z = 6, ρ calc = 1.76 g cm −3 , μ(Cu Kα) = 7.61 mm −1 , h ≤ 7, k ≤ 7, l ≤ 8, 10°≤ 2θ ≤ 90°, 4705 points, 32 parameters, R wp = 0.037, R exp = 0.023, R p = 0.027, GOF = 1.593,R Bragg = 0.015.Further details of the crystal structure determination are available from the Fachinformationszentrum Karlsruhe, D-76344 Eggenstein-Leopoldshafen, Germany, E-mail: crysdata@fiz-karlsruhe.de, on quoting the depository numbers CSD-433838 (Na 3 P).
Density Functional Theory and Computational Structure Predictions.Following the comprehensive structure prediction study of Mayo et al., 38 a parallel, cross-compositional genetic algorithm (GA) was used to enhance the sampling of metastable Na x P phases.Specifically targeting the region 1 − δ < x < 2 + δ, new low-lying phases were discovered by extrapolating the local structures found using ab initio random structure searching (AIRSS) to larger cells of nontrivial stoichiometries via simple cut-and-splice crossover (70%) and mutation (30%).Aside from physicality of density and bond lengths, two criteria were enforced upon each trial structure; the unit cell must contain fewer than 40 atoms, and it must be unique to other computed structures (as defined via overlap of pair distribution functions).Each trial structure was relaxed with density functional theory (DFT) forces using the planewave pseudopotential code CASTEP (version 16.1) 56 such that formation energies were converged to within 10 meV•atom −1 .This level of accuracy required a planewave cutoff of 400 eV and a Monkhorst−Pack grid with spacing finer than 2π × 0.05 Å −1 to sample the Brillouin zone.The fitness of each structure was evaluated as a logistic function of the distance from the convex hull of previously computed formation energies; fitness proportionate (roulette wheel) selection was used to propagate structures to the next generation, of which 50% were selected from elite structures of earlier generations.Several GA runs were performed with a generation size of 30 relaxed structures, until further generations no longer explored new regions of configuration space.
The 318 unique structures that were found to lie less than 150 meV• atom −1 from the convex hull were re-relaxed using the exchange− correlation functional of Perdew, Burke, and Ernzerhof (PBE) 57 and Vanderbilt ultrasoft pseudopotentials (generated on-the-fly with Na: 2| 1.3|1.3|1.0|16|19|21|20U:30U:21(qc= 8), P: 2|1.8|1.8|1.5|3|5|6|30:31:32LGG).Convergence of formation energies to within 1 meV•atom −1 required truncation of the planewave basis set at 650 eV, and Brillouin zone sampling with a Monkhorst−Pack grid finer than 2π × 0.03 Å −1 . 58he original work of Mayo et al. 38 searched for candidate crystal structures using a combination of the AIRSS and species-swapping structure prediction methods.These techniques have been successfully applied in several previous studies of battery materials, such as in Li− Sn and Li−Sb, 51 Na−Sn, 34 and in Li-FeS 2 . 59As well as the structures arising from our GA searches, we also examined the structures lying within 100 meV•atom −1 of the convex hull in reference 38 with stoichiometries NaP, Na 2 P, Na 5 P 4 , and Na 5 P 6 for our NMR calculations.
NMR Calculations.Chemical shifts and other NMR parameters were calculated using the gauge-including-projector-augmented-wave (GIPAW) algorithm and the Quantum Espresso code, 60 with the GIPAW pseudopotentials Na.pbe-tm-gipaw-dc.UPF and P.pbe-tmnew-gipaw-dc.UPF, 61 and the PBE exchange−correlation functional.A kinetic-energy cutoff of 1360 and 5442 eV was used for the plane-wave expansion of the Kohn−Sham orbitals and the charge density, respectively.
Structures for the NMR calculations were first relaxed using DFT and the CASTEP or Quantum Espresso codes, using a Brillouin zone sampling density of 2π × 0.025 Å −1 .Following this, an NMR calculation was carried out with the aforementioned parameters and pseudopotentials using a stricter sampling density of 2π × 0.008 Å −1 .With this sampling density, we estimate that total NMR chemical shifts are converged to ±1 ppm.We found that this density is necessary as some of our candidate crystal structures display semimetallic or metallic behavior under DFT.
Additional data related to this publication are available at the Cambridge data repository.

Figure 2 .
Figure 2. Ex situ 31 P MAS NMR spectra (MAS = 60 kHz) of black P recorded at various stages of sodiation/desodiation (right) and the corresponding electrochemistry (left).The circles on the electrochemical profile indicate where cycling was arrested and samples were extracted for the NMR experiments.The 31 P chemical shift of the pristine material, black P (dark gray, 14 ppm), and the final discharge product, crystalline Na 3 P (light gray, −207 ppm) are shown in dashed lines. 31P chemical shift regions corresponding to P helical motifs and P near the end of chains are shaded in blue and green, respectively.
P shifts at −29, −78, and −143 ppm that correspond to a-Na x P species (Figure 4a), (with anisotropic projections shown in Figure 4b−d), the width and asymmetry of the CSA patterns differing noticeably for the three resonances.Sideband patterns were simulated in dmfit 46 (Figure 4b−d, gray) and compared to those calculated for 31 P sites in DFT models (Figure 4b−d, purple) with similar compositions (vide supra).

Figure 4 .
Figure 4. (a) Deconvolution of the 31 P isotropic projection from 2D 31 P PASS experiments (black line) performed on a sample of approximate composition NaP, extracted from a NIB at 0.38 V (780 mA h g −1 ) during sodiation.Five regions with centers-of-mass at 14 (black P), −29, −78, −143, and −235 ppm are observed.The fit of the line shape from deconvolution is shown in gray.(b−d) Corresponding experimental anisotropic projections from 2D 31 P PASS (black, bottom) showing sideband patterns for δ iso = −29 (b), −78 (c), and −143 ppm (d) with the respective fits (middle, gray) compared to sideband patterns calculated for individual P sites in DFT models (top, purple) at MAS = 10 kHz.The corresponding P sites (top) and extended P structures (bottom) are displayed to the right of each sideband pattern and are listed in Table2.

Figure 5 .
Figure 5. Experimental anisotropic projections from 2D 31 P PASS (black, bottom) of NIBs at 0.80 V (1730 mA h g −1 ) during desodiation showing sideband patterns for δ iso = −29 (a) and −78 ppm (b) with the respective fits (middle, gray) compared to sideband patterns calculated for individual P sites in DFT models (top, purple) at MAS = 10 kHz.The corresponding P sites (top) and extended P structures (bottom) are displayed to the right of each sideband pattern and are listed in Table3.
P and 23 Na solid-state NMR.Dried powders filled 1.3 mm o.d.(2.5 μL internal volume) ZrO 2 rotors for fast MAS experiments.For 2D 31 P PASS measurements, 3−8 mg of material from a single battery could not fill the entire volume of a 4 mm o.d.ZrO 2 rotor (80 μL internal volume).

Table 1 .
Calculated 23 Na NMR Parameters for Na 3 P Structure Models

Table 2 .
Comparison of the Experimental 31 P CS Tensors Obtained for the "NaP" Sample Extracted at 0.38 V during Sodiation and Calculated CS Tensor Parameters for Individual 31 P Sites in Structural Models with Stoichiometries at or near NaP Structural model was found using AIRSS.b Model was found using GA.c Model does not exhibit a bandgap, and the Knight shift contribution to the isotropic shift was not calculated. a . Comparison of the Experimental 31 P CS Tensors Obtained for Na x P at 0.80 V (1730 mA h g −1 ) during Desodiation and Calculated CS Tensor Parameters for Individual 31 P Sites in Structural Models with Stoichiometries at or near NaP Structural model was found using GA.b Model does not exhibit a bandgap, and the Knight shift contribution to the isotropic shift was not calculated.c Simulation does not accurately represent experimental data due to spectral overlap. a

62 ■
ASSOCIATED CONTENT * S Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/jacs.8b04183.Additional electrochemistry, solid-state NMR spectroscopy, PXRD, and structural modeling data, including Figures S1−S17 and Tables S1 and S2 (PDF) Science for funding under grant number EP/L015552/1.A.J.M. and J.N. acknowledge the Winton Programme for the Physics of Sustainability.J.N. also acknowledges support from the Isaac Newton Fund.L.E.M. thanks Dr. Derrick Kaseman for providing the Matlab script used to process 2D PASS data.We acknowledge Josh Stratford, Dr. Elizabeth Castillo-Marti ́nez, Dr. Michael Gaultois, Dr. Pieter Magusin, and Prof. Michael Ruck (TU Dresden) for helpful discussions.NMR calculations were performed using the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (http://www.csd3.cam.ac.uk/), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council, and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).Structure prediction calculations were performed using the resources of the Center for Functional Nanomaterials, which is a U.S. DOE Office of Science Facility, at Brookhaven National Laboratory under Contract No. DE-SC0012704 and the Thomas Tier 2 facility of the UK national high performance computing service, for which access was obtained via the UKCP consortium and funded by EPSRC grant no.EP/K014560/1.