Imaging electronic states on topological semimetals using scanning tunneling microscopy

The band structure of topological semimetals such as Dirac and Weyl semimetals arises from the discrete number of crossing points of the conduction and valence bands [1–11]. In analogy with relativistic quantum field theory [12–14], the crossing of two non-degenerate or two doubly degenerate bands leads to emergent low-energy, bulk excitations described as two-component Weyl fermions or four-component Dirac fermions, respectively.

In general, the hybridization of electronic bands results in avoided crossings and hinders the realization of such semimetals. The symmetries of the system, therefore, are crucial to restrict the phase space of interactions among the bands and stabilize the semimetal phases. In the presence of time-reversal and space-inversion symmetries, the electronic bands are doubly degenerate at every crystal momentum and a three-dimensional Dirac semimetal phase appears at the quantum critical point between trivial and topological insulators [3, 15–19]. As such realization of Dirac nodes demands a fine-tuning of material parameters, e.g., the strength of spin–orbit interaction, it is difficult to obtain them in real materials. However, the hybridization of bands can be prevented by additional crystallographic symmetries such as discrete rotational symmetries, which lead to protected degeneracies of a three-dimensional Dirac node even in stoichiometric crystalline systems [20–23].

Breaking either the time-reversal or inversion symmetry [3–10] lifts the two-fold degeneracy of the bands and reduces the four-fold crossings to two-fold ones in the bulk band structure. In the vicinity of these Weyl points, the low-energy excitations are equivalent to massless Weyl fermions. Each Weyl node acts as a topological charge since it is a source or sink of the Berry flux depending on its chirality [4, 24–27]. In contrast to Dirac semimetals, the Weyl nodes possess intrinsic topological stability: translationally invariant, small perturbations can only shift the position of Weyl points in the energy–momentum space. However, bringing a pair of opposite Weyl points to the same location may result in a gap opening. These topological charges in the bulk Brillouin zone can lead to protected gapless states on the surface in the form of disjointed Fermi arcs connecting the projected Weyl-nodes with opposite chirality [22, 23, 28, 29].

Dirac and Weyl semimetals have an unusual surface-bulk connectivity that sets them apart from other topological electronic systems. The existence of Weyl nodes leads to exciting physical properties such as the chiral anomaly [30–34], nonlocal transport [35, 36], novel quantum oscillations [37, 38], diverging diamagnetism [39] and unusual optical conductivity [40]. Recent photoemission measurements provided evidence of linearly dispersing bulk bands and the Fermi arc structure of the surface states both on Dirac [41–43] and Weyl semimetals [44–51]. Here we employed scanning tunneling microscopy (STM) to probe the exotic electronic structure of topological semimetals with high energy and spatial resolution, which could shed light on unique aspects of these systems [52–56].

A powerful method to probe the electronic band structure of materials with STM is to use Fourier-transform scanning tunneling spectroscopy (FT-STS) [57–60]. When quasiparticles originating from the surface and/or bulk states of the system scatter from potential barriers caused by lattice imperfections, the incoming and outgoing electronic waves with crystal momenta of k_i and k_f form a quantum interference pattern with q = k_i − k_f  characteristic wavevector. The quasiparticle interference (QPI) results in modulation in the local density of states of the sample, which can be mapped locally by measuring the conductance between the STM tip and the sample surface. The Fourier transform of the obtained conductance map (FT-STS map) can be approximately characterized by wavevectors connecting regions with high density of states in the band structure.

More precisely, if ${\psi }_{n\alpha }({\boldsymbol{k}})$ denotes the Bloch eigenstate of the nth Bloch band with spin component $\alpha ,$ energy of ${\varepsilon }_{n}({\boldsymbol{k}})$ and lifetime ${\rm{\Gamma }},$ we can introduce the spectral function matrix:

$$\begin{eqnarray*}&&{A}_{\alpha \beta }({\boldsymbol{k}},\omega )=-\displaystyle \frac{1}{\pi }\displaystyle \sum _{n}{\rm{Im}}\left(\displaystyle \frac{1}{\omega -{\varepsilon }_{n}({\boldsymbol{k}})+{\rm{i}}{\rm{\Gamma }}}\right){\psi }_{n\alpha }({\boldsymbol{k}}){\psi }_{n\beta }^{+}({\boldsymbol{k}}).\end{eqnarray*}$$

The total and spin spectral density can be obtained as ${\rho }_{i}({\boldsymbol{k}},\omega )=Tr[{\sigma }_{i}A({\boldsymbol{k}},\omega )]$ for $i=0,$ and $i=1,\,2,\,3$ respectively, where ${\sigma }_{i}$ are the Pauli matrices. In the case that all scattering events are allowed, the FT-STS map can be approximated by the joint-density of states:

$$\begin{eqnarray*}&&{\rm{J}}{\rm{D}}{\rm{O}}{\rm{S}}({\boldsymbol{q}},\omega )=\displaystyle \int {\rm{d}}{\boldsymbol{k}}\,{\rho }_{0}({\boldsymbol{k}},\omega ){\rho }_{0}({\boldsymbol{k}}+{\boldsymbol{q}},\omega ).\end{eqnarray*}$$

Taking into account the effect of the spin texture in the scattering process, we can describe the resulting FT-STS map with the spin-selective scattering probability:

$$\begin{eqnarray*}&&{\rm{S}}{\rm{S}}{\rm{P}}({\boldsymbol{q}},\omega )=\displaystyle \frac{1}{2}\displaystyle \sum _{i=0}^{3}\displaystyle \int {\rm{d}}{\boldsymbol{k}}\,{\rho }_{i}({\boldsymbol{k}},\omega ){\rho }_{i}({\boldsymbol{k}}+{\boldsymbol{q}},\omega ).\end{eqnarray*}$$

In the pure state limit, this expression reduces to a more commonly used form:

$$\begin{eqnarray*}&&{\rm{S}}{\rm{S}}{\rm{P}}({\boldsymbol{q}},\omega )=\displaystyle \int {\rm{d}}{\boldsymbol{k}}\,{\rho }_{0}({\boldsymbol{k}},\omega )T({\boldsymbol{k}},{\boldsymbol{q}}){\rho }_{0}({\boldsymbol{k}}+{\boldsymbol{q}},\omega ),\end{eqnarray*}$$

where $T({\boldsymbol{k}},{\boldsymbol{q}})={| \langle S({\boldsymbol{k}})| S({\boldsymbol{k}}+{\boldsymbol{q}})\rangle | }^{2}$ and $S({\boldsymbol{k}})$ is the spin of the band at momentum ${\boldsymbol{k}}$.

When calculating the integrals above, the integration boundary in the momentum space has to be chosen carefully. Here, we follow a recent theoretical work [61], and we introduce a form factor which envelopes the spectral density ${\tilde{\rho }}_{i}({\boldsymbol{k}},\omega )={\rho }_{i}({\boldsymbol{k}},\omega )\times \exp (-{{\boldsymbol{k}}}^{2}/{\xi }^{2})$ and calculate the integrals over multiple, higher order Brillouin zones without a sharp cutoff at the boundary of the first Brillouin zone.

To enhance the signal-to-noise ratio of the FT-STS maps, generally, a series of data analysis steps are applied to the experimentally obtained real space conductance maps. In this paragraph, we briefly summarize the most common techniques used for this purpose. Firstly, we note that although an impurity is a source of the QPI signal, the conductance value measured on top of it is dominated by tunneling into its bound states rather than by the signal originating from the QPI. Therefore, it is beneficial to replace the real space conductance values obtained at the position of the impurities with the mean conductance value of the map at the given energy to suppress this effect. Drift due to piezo motion or thermal expansion during the long time scale (typically several days) of acquisition of a conductance map also decreases the sharpness of the QPI features. To compensate this artifact, we use a previously developed drift correction method on the real space images [62]. Also, the noise of the scanning piezo can be reduced by normalizing the measured conductance value with the recorded tunneling current at each point of the map. Further improvement of the QPI map can be achieved by suppressing the edge effects arising from the finite size of the conductance map, which can be done by applying certain window functions (e.g. tukey window function). Finally, the Fourier transform of the maps are symmetrized based on the symmetry of the underlying lattice. We note that the reported features are apparent without any of these data analysis steps.

In addition to the FT-STS technique, investigating the energy structure of Landau levels formed in high magnetic fields can provide another experimental approach to characterize the electronic band structure. Previously, Landau level scanning tunneling spectroscopy has been utilized to characterize the energy momentum dispersion in two-dimensional systems through the emerging quantized energy levels [63–66]. In three-dimensional systems, however, the presence of magnetic field leads to the formation of Landau bands dispersing along the momenta parallel to the applied field. In an STM experiment, the van-Hove singularities of these Landau bands cause peaks in the measured projected bulk density of states. The ${\varepsilon }_{N}$ Landau level energies for a linearly dispersing band is expected to have ${\varepsilon }_{N}={\varepsilon }_{DP}\pm {v}_{F}\sqrt{2eB(N+\gamma )/\hslash }$ scaling behavior, where ${\varepsilon }_{DP}$ is the Dirac energy, ${v}_{F}$ is the Fermi velocity, e is the electron charge, $B\,$ is the magnetic field, $N$ is the orbital index, $\gamma$ is the phase offset of quantum oscillations, $\hslash$ is the Planck's constant divided by 2π. Hence, the square-root dependence of the singularities in the energy spectrum signals the linear dispersion relation with the slope corresponding to the Dirac velocity.

The first theoretically proposed [22, 23] and experimentally identified [41–43] three-dimensional Dirac semimetals are Na₃Bi and Cd₃As₂. Both materials are characterized by strong spin–orbit interaction and by an inverted band structure. First principle calculations [22, 23] showed that the band inversion occurs in the vicinity of the Γ point between the Cd-5s and As-4p states in Cd₃As₂, while in Na₃Bi, it occurs between Na-3s and Bi-6p_x,y orbitals. The pair of Dirac points in both materials is located along the k_z direction in the first Brillouin zone and stabilized by C₄ and C₃ rotational symmetry of the crystal structure (figures 1(a) and 2(a)) in case of Cd₃As₂ and Na₃Bi, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 1(b) and 2(b) show atomically flat topographic images of the natural (112) cleaved surface of Cd₃As₂ [52] and the (100) surface of the first batch of Na₃Bi(1) [53]. The (112) surface of Cd₃As₂ displays a pseudo-hexagonal structure corresponding to the As layer, whereas the Na₃Bi(1) surface displays stripe-like topographic modulations along the c axis of the crystal. The FT-STS measurements carried out on these surfaces display long wavelength ripples in the real space density of states (figures 1(d)–(g) and 2(d)–(g)). The Fourier transform of the spectroscopic maps (figures 1(h)–(k) and 2(h)–(k)) can be characterized by a rather isotropic set of q scattering wavevectors in both materials. By determining the boundary of the QPI signal in the FT-STS maps as a function of energy, the scattering vectors exhibit an approximately linear energy-evolution even extending to energies far from the Dirac point (figures 1(c) and 2(c)). At low-energies, the diverging quasiparticle scattering wavelength and sample inhomogeneity become comparable, and therefore prevents resolving of the band dispersion close to the Dirac point by FT-STS.

When high magnetic field is applied perpendicular to the (112) face of Cd₃As₂ and the (001) surface of the second batch of Na₃Bi(2) [52, 53], the measured spectra exhibit a series of pronounced peaks in the tunneling spectra as illustrated in figures 3(a) and (c), respectively. In the case of Cd₃As₂, each LL splits into two peaks in a non-uniform way indicating that the spin degeneracy is lifted by the combination of orbital and Zeeman splitting. By mapping the energy of the peaks as a function of magnetic field, we find that they correspond to electron-like Landau levels. The position of the Landau level peaks exhibits an ${\varepsilon }_{N}\propto \sqrt{(N+\gamma )B}$ scaling, providing evidence of Dirac-like dispersion of the bands (figures 3(b) and (d)). For Cd₃As₂, the Dirac energy is located slightly below −200 meV $({v}_{F}\approx 6.2\,\text{eV}{\mathring{\rm{A}}}),$ while in case of Na₃Bi(2) it is located closer to the Fermi energy at −20 meV $({v}_{F}\approx 5.1\,\text{eV}{\mathring{\rm{A}}}).$ We also note that model calculation for Cd₃As₂ showed that the band inversion is expected to be as small as 20 meV [52], which prevents us to grab the signatures of Lifshitz transition in this material.

Zoom In Zoom Out Reset image size

Since studying the physics of Weyl fermions requires either broken inversion or time-reversal symmetry, application of magnetic field, in principle, is one of the possible routes to reach the Weyl semimetal phase from the Dirac phase. In case of Cd₃As₂, the normal vector of the natural cleavage plane is tilted away from the rotational symmetry axis, hence, applying magnetic field perpendicular to the plane eliminates the Weyl nodes [52]. In case of Na₃Bi, the high spatial variability of the measured spectra and surface disorder prevented access to the physics of Weyl fermions [53].

The first experimentally confirmed, inversion-symmetry-breaking Weyl semimetal compounds were found in the transition metal monopnictide materials class [27, 28, 44–51]. We here focus on TaAs, which possesses a body-centered tetragonal lattice with four adjacent Ta and As sublayers forming a unit cell (figure 4(a)). The subsequent sublattices are shifted by a fraction of the unit cell, hence the inversion symmetry is naturally broken in the crystal. Previous theoretical calculations [27, 28] showed that the bulk band structure has 12 pairs of Weyl points located close to the mirror planes. Four pairs of the Weyl nodes (W1) sit on the k_z = 0 plane, while the remaining 8 pairs (W2) are located away from it. The projection of the Weyl nodes on the (001) surface, however, results only in 8 pairs of the Weyl points in the surface BZ (figure 4(b)) because the pairs of the W2 nodes have the same in-plane momentum.

Zoom In Zoom Out Reset image size

Similar to our measurements on the Dirac semimetals [52, 53], we utilized quantum interference technique to explore the physics of Weyl fermions in TaAs [54]. Figure 4(c) shows the (001) cleaved, As-terminated surface of single crystal TaAs with a few different crystallographic defects, which give rise to strong modulations in the surface density of states due to scattering of quasiparticles (figure 4(d)). The Fourier transform of this conductance map reveals that the interference pattern can be characterized by a rich variety of scattering wavevectors arising from the rather complex surface band structure (figure 4(e)). We calculated the SSP maps (figure 4(f)) based on our density-functional theory (DFT) results, which unlike in the case of the QPI measured in Dirac semimetals, is undoubtedly insufficient to describe the measured QPI: the SSP displays a large set of wave vectors, which is absent in the measured data. Understanding the discrepancy between the observed and predicted QPI signal, however, reveals an important new aspect of the physics of electrons in a Weyl semimetal: the momentum dependent delocalization of the surface Fermi arcs into the bulk.

Since Fermi arcs connect projected bulk Weyl nodes with opposite chirality, the connectivity of the surface states with the bulk states is a unique aspect of Weyl semimetals. A remarkable consequence of the surface-bulk connectivity of topological semimetals has been possibly illustrated in recent magnetotransport measurements [37, 38]. In contrast to conventional materials, where the electron cyclotron orbit consists of either surface or bulk trajectories on the Fermi surfaces perpendicular to the magnetic field, topological semimetals feature unusual, closed magnetic orbits consisting of the combination of the two. When an electron travels on a disjoint Fermi arc and reaches the end point of the arc (the projection of the Weyl node), it loses its surface characteristic and traverses through the bulk until it reaches the other surface. This unusual quantum path involving electrons sinking around the Weyl nodes leads to anomalous quantum oscillations in transport measurements. In our STM measurements [54], we revealed how this 'sinking effect' influences the QPI pattern on the surface of a Weyl semimetal.

In order to quantitatively understand the surface-bulk connectivity of TaAs, we first study the atomic characteristics of both the bulk and surface states as a function of their momenta using our DFT calculation. Focusing on the bulk properties, we calculate the contribution of the Ta-5d and As-4p orbitals on the bulk bands. As an example, figure 5(a) shows a linecut in the bulk band structure through the W1 and W2 Weyl points, where the red and green colors indicate the corresponding Ta and As weights, respectively. At momenta far from the Weyl points, the As orbitals dominate the wavefunctions, while the Weyl nodes themselves are clearly monopolized by the Ta orbitals (at W1 and W2 the ratio of the contribution of Ta/As orbitals is about 6:1 and 14:1, respectively).

Zoom In Zoom Out Reset image size

In the case of the surface states, we can determine the momentum-dependent atomic character of the Fermi arcs by projecting the DFT calculated spectral density on the top Ta and As layers separately and subtracting the two projections from each other (figure 5(b)). We observe that the bowtie-shaped features are more pronounced on the As layer (green), while the spoon-like features are more visible on the Ta layer (red). Thus, we can see that bowties are generally associated with As orbitals and spoons with Ta orbitals. This behavior is also expected from the momentum-space location of the states: spoons are located closer to the projected Weyl cones, which are dominated by the Ta contribution, and therefore they have a stronger connection to the bulk. Additionally, when we follow the intensity of the bowtie and spoon features as a function of layer position, it is apparent that the contribution of the bowties decreases more rapidly, signaling a weaker connectivity to the bulk as compared to the spoon features (figure 5(c)). A similar trend can be observed by looking at the total spectral weight on the As and Ta layers (figure 5(d)).

To summarize our atomic orbital characterization of the Weyl states, we see that the spoon features associated with Ta have stronger bulk character, while the bowtie surface states residing more strongly on As have a more surface-like character. As the QPI observed in an STM measurement is naturally dominated by the electronic states which remain on the surface, in order to understand our data we need to focus on states which have a smaller probability to connect to the bulk. Therefore, we introduce a weighted Fermi surface [54] which consists of states projected only on the As sites instead of using the projection onto the entire unit cell. As figure 5(e) shows, this projection causes the spoon features to become more faint in the weighted Fermi surface but leaves the bowtie features mostly unaffected.

There are three sets of scattering wavevectors in the first Brillouin zone: between bowties and bowties (P vectors), between bowties and spoons (Q vectors) and finally between spoons and spoons (R vectors) with their intensity rapidly decreasing from P to R owing to the intensity difference between bowties and spoons in the weighted Fermi surface. The possible P, Q and R scattering vectors are shown separately in figures 6(a)–(c), and the corresponding QPI signal in the calculated SSP is highlighted by red and blue rectangles in figures 6(d)–(f). In the case of the P and Q vectors, the same QPI vectors can be clearly identified in the measured FT-STS maps (figures 6(g) and (h)), illustrating an excellent match between theory and measurements. Although some of the S vectors can be found in the SSP calculation (figure 6(f)), their weak intensity or their overlap with other wave vectors prevents us from detecting them in our measurement (figure 6(i)).

Zoom In Zoom Out Reset image size

The strong connection between the atomic orbital character and the Fermi surface features suggests that the different lattice sites carry different parts of the QPI signal. To investigate this idea, we decompose the measured conductance maps according to Ta and As sublattices. By taking the real space conductance map and considering the measured signal only on top of Ta or As sites (associated with local hills and valleys in the topographic image as shown in figure 6(l)), two maps can be constructed: the Ta-site and the As-site selected conductance map. The Fourier transform of each map remarkably shows that they can be independently described by distinct set of scattering vectors corresponding to the P and Q vectors (figures 6(j)–(k)). While the As-site selected conductance map can be understood as scatterings between bowties, the QPI measured on the Ta site can be described as scatterings between bowties and spoons.

The revealed momentum-dependent surface-bulk connectivity is insensitive to the topological nature of the surface states: both topological and trivial states are delocalized into the bulk at momenta close to the projection of the Weyl points. The presented DFT calculation, however, can be used to access additional information about the topological classification of the surface states on TaAs [54]. The calculated band structure encircling W1 and W2 (on a trajectory indicated on figure 7(a)) is shown on figures 7(b) and (c), respectively. In case of the W1 node, the gray shaded regions corresponds to the bulk, whereas yellow and blue lines correspond to the topologically protected surface states connecting the conduction and valence band. The energy–momentum diagram reveals that the surface states give rise to four crossings at the Fermi energy and consequently all four surface states around W1 have topologically nontrivial origin. The number of Fermi arcs naively contradicts the charge of W1, which is two. However, by taking into account the 'chirality' of the states (the sign of their velocities when they cross the Fermi energy), we can see that the three chiral arcs (crossings 2, 3, 4) and the one anti-chiral (crossing 1) arc result in a net number of two crossings in agreement with the total topological charge of W1. We also note that one of the topological arcs is possibly hybridized with a trivial surface state, therefore, the distinction between topological or trivial states is rather vague in this case. The classification of surface states around W2 is clearer. As figure 7(c) shows, the black line corresponds to a trivial, while the green line is a topological surface state. Therefore, the topologically protected Fermi arc state can be associated with the state directly connecting the two W2 points, while the inner bowtie has trivial character.

Zoom In Zoom Out Reset image size

Topological semimetals represent a novel quantum phase of matter, where touching of bulk electronic bands leads to non-trivial band topology. Our STM measurements [52–54] with high energy and spatial resolution provided a crucial insight into the physics of these systems by mapping the electronic behavior in high magnetic field and visualizing the scattering of both the bulk and surface carriers establishing a topological connectivity between those states.

The Cd₃As₂ crystals were grown from a Cd-rich melt with the stoichiometry Cd₈As₃ in an evacuated quartz ampoule. The sample was heated to 800 °C at 3 °C min⁻¹, then slowly cooled to 400 °C at 0.1 °C min⁻¹. After reaching 400 °C, the sample was subsequently centrifuged to remove excess molten Cd.

In case of the stoichiometric-melt Na₃Bi(1) crystals (measurements on the (100) plane), a 3:1 atomic ratio of Na (99.5%) and Bi (99.999%) was used. The samples were placed in one-end-sealed tubes of C-coated stainless steel (304 l). The open ends of the tubes after starting material insertion were crimped and welded under a continuous Ar-flow. The sealed ampoules were heated to 900 °C at 60 °C h⁻¹ for 6 h, cooled to 750 °C at 10 °C h⁻¹, and further cooled to 300 °C at 1.5 °C h⁻¹, after which they were left to anneal at 300 °C for 24 h and finally furnace cooled to room temperature. For the Na-flux grown crystals (Na₃Bi(2) measurements on the (001) plane), a Na rich composition of Na plus Bi in a 90:10 atomic percent ratio was chosen. The starting materials were placed in C-coated steel tubes and steel wool was used as a plug. The sealed ampoules were heated to 850 °C at a rate of 60 °C h⁻¹ for 12 h to ensure homogeneous mixing of the elements, and then cooled to 750 °C at 10 °C h⁻¹. For the entire process of nucleation and growth, the temperature was lowered slowly to 325 °C, at a rate of 1.5 °C h⁻¹. Due to the extremely air sensitive nature of the Na₃Bi crystals, they were harvested, stored and prepared for the STM measurements in an Ar-filled glove box.

The TaAs single crystals were grown out of Ga flux with the ratio of Ta:As:Ga = 1:1:98. To increase the solubility of Ta in Ga, we used 80% of Ta in the thin foil form and 20% of Ta in the powder form as the starting material. We packed all starting materials in an alumina crucible, together with a catch crucible filled with quartz wool, then sealed them inside a quartz tube under 1/3 ATM Argon gas. The ampoule was heated up to 1100 °C, held for 10 h, cooled down to 1000 °C using 3 °C h⁻¹ and held for 4 h, then cooled down to 700 °C at 1 °C h⁻¹. We then took the ampoule out and spin off the Ga flux at ∼80 °C. Ga flux was removed from the single crystals using concentrated hydrochloride acid.

Prior to the STM measurements the samples were cleaved in ultra-high vacuum and immediately transported into the microscope head. The measurements were carried out in home-built STMs at the following temperatures: 400 mK for the point spectra and 2 K for the QPI on Cd₃As₂, 3 K for Na₃Bi and 40 K for TaAs. The reported QPI maps were acquired on a 208 × 208 grid for Cd₃As₂ and Na₃Bi at lock-in frequency of 798.7 Hz, while on a 256 × 256 grid for TaAs at lock-in frequency of approximately 701 Hz.

Work at Princeton was supported by Army Research Office-Multidisciplinary University Research Initiative (ARO-MURI) program on topological insulators W911NF-12-1-0461, Gordon and Betty Moore Foundation as part of the Emergent Phenomena in Quantum Systems (EPiQS) initiative (GBMF4530), by NSF Materials Resarch Science and Engineering Centers (MRSEC) programs through the Princeton Center for Complex Materials DMR-1420541, NSF-DMR-1104612, NSF-DMR-0819860, NSF CAREER DMR-0952428, DOE DE-FG-02-05ER46200, Packard Foundation, and Keck Foundation. This project was also made possible through use of the facilities at Princeton Nanoscale Microscopy Laboratory supported by grants through ARO-W911NF-11-1-0379, ARO-W911NF-1-0262, ONR-N00014-14-1-0330, ONR-N00014-13-10661, ONR- N00014-11-1-0635, US Department of Energy-Basic Energy Sciences (DOE-BES), Defense Advanced Research Projects Agency—US Space and Naval Warfare Systems Command (DARPA-SPWAR) Meso program N6601-11- 1-4110, LPS and ARO-W911NF-1-0606, and Eric and Wendy Schmidt Transformative Technology Fund at Princeton. Work at University of California-Los Angeles was supported by the DOE-BES (DE-SC0011978). The research at Princeton was also supported by the Department of Energy DE-SC0016239, NSF EAGER Award NOA - AWD1004957, Simons Investigator Award and Schmidt Fund for Innovative Research.