Momentum space Aharonov-Bohm interferometry in Rashba spin-orbit coupled Bose-Einstein condensates

Topological orders of matter have recently gained great attention in solid-state and cold-atom physics [1–14] because of their characterization of quantum phases through a scenario different from conventional Ginzburg-Landau orders [15,16] and potential applications in fault-tolerant quantum computation [17–19]. A large variety of topological phases, including quantum Hall states [20–22], topological insulators [23,24], and anomalous Hall states [25–27], can be characterized by a geometric phase [28], or Berry phase, of the underlying band structure of the system. In this context, direct measurements of Berry phases are essential for exploring new physics of topological states of matter.

The Berry phase describes the adiabatic phase shift of the system wave function along a closed loop in parameter space. One well-known example of Berry phases is the Aharonov-Bohm phase [29] in real space that describes the phase obtained for a charge particle encircling a magnetic flux, which can be measured through a real-space interferometry. As an analogy, the Berry phase of an energy band may be measured using a momentum space interferometry. Although it is challenging to realize electronic interferometry in solid-state systems, an Aharonov-Bohm interferometer in reciprocal (lattice-momentum) space has recently been realized with ultracold atoms in optical lattices, where the Berry phase of a two-dimensional (2D) hexagonal lattice [1,30] was measured.

Ultracold atomic gases exhibit great flexibility and controllability in engineering single-particle energy band dispersion for both continuous space and discrete lattice [31–43]. In particular, a distinct type of single-particle band governed by 2D Rashba spin-orbit coupling (SOC) [44–49] has been experimentally realized in both ultracold Fermi [50,51] and Bose [52] gases. The band structure with such 2D Rashba SOC exhibits a Dirac point, which has been detected directly using spin-resolved radio-frequency spectroscopy [50,51]. However, the topological property of the Dirac point, i.e., the π Berry phase, has not been directly measured in experiments.

In this paper, we propose a scheme that conducts momentum space Aharonov-Bohm interferometry in Rashba spin-orbit–coupled Bose-Einstein condensates (BECs), which measures the Berry phase along a loop enclosing the Dirac point. As shown in fig. 1(a), our interferometer exploits the intrinsic ring structure in the Rashba energy band as the interferometer loop and the controllable detuning $\delta_{x}\sigma_{x}$ of ultracold atoms as a trigger for driving the BEC. A sudden change of $\delta_{x}$ sets up the initial BEC state at the highest energy point of the ring (point A). As the BEC naturally pursues the lowest energy (point B), it splits into two, following different halves of the ring (C and C') and exhibiting the interference. Since the whole ring path encloses the region with dense Berry flux, the interference pattern would reflect the Berry phase as a function of $\delta_{z}$, as in fig. 1(b). At $\delta_{z}=0$, the π Berry phase indicates the presence of a Dirac point.

Zoom In Zoom Out Reset image size

We adopt two complementary methods, Gross-Pitaevskii equation (GPE) simulation and variational analysis, for studying the BEC dynamics and interference. Our simulation shows real-time evolution of the interference pattern for ⁸⁷Rb BECs in typical cold-atom experiments, while the variational analysis provides an informative model capturing the key physical features of the interference, including the trajectories in momentum and real spaces as well as the geometric phase acquired during the evolution. We also point out proper conditions for external trapping potential and interatomic interaction under which the interferometry procedure succeeds and discuss why improper trapping frequency or too strong interaction sabotages the desired dynamics for the interference. Our results ought to provide timely guidance for ongoing experimental study on Rashba spin-orbit coupled quantum gases.

We consider a Bose gas with atomic mass m and two hyperfine spin states $(\psi_{\uparrow} \enskip\psi_{\downarrow})^{T}$ subject to synthetic Rashba coupling in the $x\text{-}y$ plane and a tunable Zeeman field. After integrating out the irrelevant z degree of freedom, we write down the effective Hamiltonian,

$$\begin{eqnarray} &&H=\frac{1}{2m}(\hat{p}_{x}^{2}+\hat{p}_{y}^{2})+V+{H_{\mathrm{R}}}+{H_{\mathrm{Z}}}+{H_{\mathrm{I}}}, \end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray} V &= \frac{m}{2}(\omega_{x}^{2}{{\hat{x}}^{2}}+\omega_{y}^{2}{{\hat{y}}^{2}}), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} {H_{\mathrm{R}}} &= {\lambda_{\mathrm{R}}}({{\hat{p}}_{y}}{\sigma_{x}}+a{{\hat{p}}_{x}}{\sigma_{y}}), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} {H_{\mathrm{Z}}} &= {\delta_{x}}{\sigma_{x}}+{\delta_{z}}{\sigma_{z}}/2, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} {H_{\mathrm{I}}} &= \left({\begin{array}{c@{\quad}c} {{g_{\uparrow \uparrow }}{{\left\vert {\psi_{\uparrow }}\right\vert }^{2}}+{g_{\uparrow \downarrow }}{{\left\vert {\psi_{\downarrow }}\right\vert }^{2}}} & 0 \\ 0 & {{g_{\uparrow \downarrow }}{{\left\vert {\psi_{\uparrow }}\right\vert}^{2}}+{g_{\downarrow \downarrow }}{{\left\vert {\psi_{\downarrow }}\right\vert }^{2}}} \end{array}}\right)\!.\qquad \end{eqnarray} \tag{ 5 }$$

Here V is an external trapping potential, $H_{\mathrm{R}}$ describes the Rashba coupling of strength $\lambda_{\mathrm{R}}$ and anisotropy factor a, $H_{\mathrm{Z}}$ represents the Zeeman field $\delta_{z}~(\delta_{x})$ in the longitudinal (transverse) direction (we set as zero the y component without loss of generality), and $H_{\mathrm{I}}$ results from the spin-dependent mean-field interaction. We use $g_{\uparrow \uparrow}=g_{\downarrow \downarrow}=0.9554g_{\uparrow \downarrow }\equiv g$ in the following simulation for ⁸⁷Rb atoms. (Note that this represents a regime of dominant density-density interaction and negligible spin-spin interaction. One can thus assume that such a condition remains if the system is studied in a dressed spin basis.) The recoil energy $E_{R}= \frac{\hbar^{2}k_{R}^{2}}{2m}=2\pi\hbar \times 1.8\ \text{kHz}$ is chosen as the energy unit.

Given a sufficiently weak trapping potential, the Hamiltonian has two single-particle energy bands in the $p_{x}\text{-}p_{y}$ momentum space, with the lower one being $E_{-}=\frac{\bm{p}^{2}} {2m}-|{{{\vec{d}(\bm{p})}}}|$, where $\bm{p}=(p_{x},p_{y})$ and $\vec{d}=({\lambda_{\mathrm{R}}p_{y}}+ {\delta_{x}},a{\lambda_{\mathrm{R}}}{p_{x}},{\delta_{z}}/2)$. For a = 1 and $\delta_{x}=0$, the set of minima of $E_{-}$ forms a horizontal ring $|\bm{p}|=\sqrt{4\lambda_{\mathrm{R}}^{4}-\delta_{z}^{2}}/ (2\lambda_{\mathrm{R}})$ in the $p_{x}\text{-}p_{y}$-plane. If $\delta_{x}\neq 0$, the ring structure remains but is inclined such that it has only one maximum and one minimum at two intercepts with the p_y-axis, respectively (as points A and B in fig. 1(a)). Our interferometry is performed along this Rashba ring path as we will show below.

The Berry phase γ in the region enclosed by this ring loop can be computed as

$$\begin{equation} \gamma =\oint_{L_{c}}\bm{A}_\text{B}\cdot \text{d}\bm{p}=\int_{S}\mathcal{F}~\text{d}^{2} \bm{p}, \end{equation} \tag{ 6 }$$

where $\bm{A}_\text{B}=-i\langle \zeta (\bm{p})|\nabla_{\bm{p}}|\zeta (\bm{p})\rangle$ is Berry connection for eigenstate $\zeta (\bm{p})$, $\mathcal{F}=\nabla \times \bm{A}_\text{B}$ is Berry curvature, and L_c and S denote the loop and the enclosed region, respectively. For our Hamiltonian, the Berry curvature is related to the unit vector ${\hat{d}}=\vec{d}/|\vec{d}|$ as

$$\begin{equation} \mathcal{F}=\frac{1}{2}\epsilon_{ij}{\hat{d}}\cdot \left( \partial_{i}{\hat{d}}\times \partial_{j}{\hat{d}}\right), \end{equation} \tag{ 7 }$$

where $\epsilon_{ij}$ is the antisymmetric permutation (Levi-Civita) symbol, and the Berry phase is equal to half of the solid angle swept by $\hat{d}$ in the loop integral. The Berry phase depends on both $\delta_{z}$ and $\delta_{x}$ but is insensitive to the latter. For $\delta_{z}=0$, the Berry curvature is a delta function centered at the Dirac point, which gives a Berry phase $\gamma =\pi$. When $\delta_{z}$ increases, the ${\hat{d}}$ vector sweeps a less solid angle, and the Berry phase monotonically decreases as shown in fig. 1(b). Therefore, one can continuously change the Berry phase by tuning $\delta_{z}$, in analogy to changing the magnetic flux in a conventional Aharonov-Bohm interferometer.

In this section, we discuss how to employ our model as an Aharonov-Bohm interferometer for detecting the Berry phase and present the GPE simulation results. In cold-atom systems, the Rashba spin-orbit coupling is generated by a set of Raman lasers that couple different spin and momentum states, and the Zeeman energy shift is determined by the relative detunings of the lasers [50–52]. All parameters are highly tunable in the current experimental setup. We consider a trapped BEC initially prepared at the ground state of fixed $\lambda_{\mathrm{R}}$, a, $\delta_{z}$, and $\delta_{x}=\delta_{0}<0$, which is a momentum space Gaussian wave packet (due to the trap) centered at the minimum of the inclined ring, i.e., $(0,-p_{0})$. Then the detuning $\delta_{x}$ is suddenly changed to $\delta_{x}=\delta_{1}>0$ (with other parameters unchanged). The ring structure is hence inversely inclined such that the current location of the condensate becomes the energy maximum (with a slight deviation ${\sim}O(\delta_{x}\delta_{z}^{2})$, which is negligible in our case), illustrated as point A in fig. 1(a). As a result, the BEC wave packet splits into two parts, which follow separate paths (C and C', respectively) and move toward the new energy minimum at the opposite end of the ring (point B). When two wave packets meet and superpose at point B, the dynamic phase cancels out due to the symmetry between two paths, while the geometric phase (Berry phase), as a function of $\lambda_{\mathrm{R}}$, a, $\delta_{z}$, and $\delta_{x}=\delta_{1}$, can be revealed by the density contrast of the interference pattern.

We remark that a proper external trap is essential for driving the motion of the condensate in the momentum space. If there is no trap, any p state is a stationary state, and the BEC does not move since $\text{d}\hat{p}/\text{d}t=i[H,\hat{p}]=i[V,\hat{p}]= 0$ (given negligible interaction). However, if the trapping potential is too strong (comparable to the Rashba coupling strength), it may also spoil the desired ground state as well as the interference dynamics [53,54]. In addition, if the interaction is too strong, the initial BEC wave packet spontaneously selects one path rather than splitting into two parts. This is because the superposition of two momentum wave packets leads to a real-space density wave that costs too much interaction energy. Furthermore, the Zeeman field $\delta_{z}$ cannot be too large. Otherwise the energy-band tip will be flattened, and consequently, the condensate will not follow the shallow Rashba ring groove.

Given the above constraints, our GPE simulation shows that typical experimental parameters are indeed suited for realizing the interference as shown in fig. 2. In panels (a1)–(a4), we plot density distributions $\rho_{\downarrow }(\bm{p})$ in momentum space at different time frames, with the starting point on the ring path labeled by A and the pursued energy minimum by B. Initially, the condensate locates at A ((a1)), which is lifted from the ground state by the sudden change of $\delta_{x}$. Then the BEC splits into two parts, which separately follow the ring loop (as C and C' in fig. 1(a)) toward B ((a2)). After the two condenates encounter each other at point B, a clear ring-shape interference pattern forms in the momentum space ((a3), (a4)). This ring structure indeed reflects BEC's natural motion along the ring groove toward the lower energy in the Rashba band. Due to this petal pattern in momentum space, the condensate also exhibits exotic circular distributions in real-space ((b1)–(b4)). Our variational analysis of key dynamical variables below will show that the trajectories of splitting BEC's center of mass explain the pattern exhibited by the GPE simulation in both momentum and real space.

Zoom In Zoom Out Reset image size

Next, we turn to study the interference pattern at various Zeeman fields $\delta_{z}$. The GPE simulation results are presented in fig. 3. In panels (a1) and (a2), for $\delta_{z}=0$, the density at point B is constantly zero during the evolution, forming a robust dark fringe that indicates a π phase shift between the splitting BECs upon encountering at point B. Due to the aforementioned dynamic phase canceling, this π phase shift is contributed purely by the accumulated geometric phase around the loop, thus confirming the presence of a Dirac point. Panels (c1) and (c2) show a finite density at B for $\delta_{z}=3\ \text{kHz}$. This indicates that the Berry phase is no longer equal to π as we expect from the smaller solid angle swept by the $\hat{d}$ vector. In figs. 3(b) and (d), we plot the corresponding real-space density distributions. They also exhibit a roughly ring-shaped interference pattern and can be understood by considering the group velocity and phase dynamics of the wave packets, as we will show later.

Zoom In Zoom Out Reset image size

In order to quantitatively relate the interference pattern with the geometric phase, we define a (local) relative brightness as

$$\begin{equation} R_I = \frac{\int_{s_{0}}\text{d}\bm{p}~\rho (\bm{p})}{\int_{s_{n}}\text{d}\bm{p}~\rho (\bm{p})}, \end{equation} \tag{ 8 }$$

where $\rho =\rho_{\uparrow }+\rho_{\downarrow }$ is the total density distribution in the momentum space, s₀ is a proper dark-fringe region around point B, and $s_{n}\supset s_{0}$ includes the adjacent bright-fringe regions (such that $0\leq R_I\leq 1$). While $R_I=0$ for $\delta_{z}=0$, it becomes nonzero as the interference occurs for any $\delta_{z}>0$. Note that since the relative brightness may slightly oscillates with time, in our simulation, we record R_I right after the center of masses of the left and right parts of BECs both move across the p_x-axis. In fig. 4, we plot R_I vs. $\delta_{z}$ obtained from the GPE simulation (blue dots), which fit a monotonically increasing curve (dashed curve). Considering the monotonically decreasing relation between Berry phase and $\delta_{z}$, as in fig. 1(b), this relative brightness R_I may act as a good experimental measurement for directly determining the geometric phase of the energy band. As we will show later, a variational analysis for the interference dynamics also yields a comparable $R_I\text{-}\delta$ relation (green dots).

Zoom In Zoom Out Reset image size

Finally, we study the effects of strong interaction and anisotropic Rashba coupling in the experiment. As mentioned above, a reasonably small interaction strength is favored for our scheme. If the interaction is too strong, splitting of the BEC wave function in momentum space induces real-space density modulation that increases the interaction energy. Consequently, the condensate spontaneously selects one route rather than equally splitting. This is shown by the GPE simulation results in fig. 5(a1) and (a2), in which all parameters are the same as those in fig. 3(c), (d) except for a large $ng=0.2E_{R}$. We see that the axial symmetry is broken in both momentum and real spaces, and the interference fringes become more obscure. In experiments, the Rashba coupling can be anisotropically tuned, i.e., $a\neq 1$ in eq. (3). As a result, the ring structure of the energy band becomes elliptical, so does the interference loop. In the right four panels of fig. 5, we show the GPE simulation results for $a=0.8$, with the other parameters as in fig. 3(a), (b). We still see a dark interference fringe with zero density at point B, indicating the robust π Berry phase against the deformation of the energy band and interference loop, as expected for the property of the Dirac point.

Zoom In Zoom Out Reset image size

In this section, we reveal salient physics features of the dynamical interference process with a simple structured variational wave function. Since we have observed in the GPE simulation that the BEC wave function intends to split into two under weak interatomic interactions, it is natural to consider a superposition that can describe the splitting condensate as

$$\begin{equation} \Psi_{\mathrm{var}}=\phi_{L}e^{i\frac{\varphi }{2}}\cos \alpha +\phi_{R}e^{-i\frac{\varphi }{2}}\sin \alpha, \end{equation} \tag{ 9 }$$

which consists of two Gaussian wave packets $\phi_{j=R,L}$ in the region of p_x > 0 and p_x < 0, respectively. Each Gaussian wave packet takes a general form [55,56] as

$$\begin{eqnarray} \phi_{j}&= \zeta (\mathbf{p}_{j})\prod_{\eta}\biggl[\left(\frac{2}{\pi R_{\eta}^{2}}\right)^{\frac{1} {4}}e^{-(\frac{1}{R_{\eta}^{2}}-\frac{i}{2}\xi_{\eta })(r_{\eta }-A_{j,\eta })^{2}} \nonumber \\ &\quad \times\, e^{ip_{j,\eta }(r_{\eta }-A_{j,\eta })}\biggr], \end{eqnarray} \tag{ 10 }$$

where $\eta =x,y$ stands for the spatial coordinates, $A_{\eta}$ is the center-of-mass position (in real space), $R_{\eta }$ is the width of the wave packet, and $\xi_{\eta }$ is introduced as the conjugate variable for $R_{\eta }$, which is essential for the completeness of this variational method [56]. The axial symmetry of the dynamics allows us to assume that two wave packets have the same $R_{\eta }$ and $\xi_{\eta }$ (independent of $j=L,R$), which have been confirmed by our GPE simulation for reasonably weak interaction. In a semiclassical picture, the system Lagrangian $\mathscr{L}=\int \text{d}^{2}\bm{r}\Psi^{\dagger }(i\frac{\partial }{\partial_{t}}-H)\Psi$ yields the equations of motions (see details in Supplementary Material Supplementarymaterial.pdf ) as

$$\begin{eqnarray} \frac{\text{d}}{\text{d}t}A_{j,\eta } &= \frac{\partial }{\partial p_{j,\eta }}E_{j,-},\qquad\frac{\text{d}}{\text{d}t}p_{j,\eta }=-\omega_{\eta }^{2}A_{j,\eta }, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \frac{\text{d}}{\text{d}t}\xi_{\eta } &= \frac{4}{R_{\eta }^{4}}-\omega_{\eta}^{2}-\xi_{\eta }^{2},\qquad \frac{\text{d}}{\text{d}t} R_{\eta }=R_{\eta }\xi_{\eta }, \end{eqnarray} \tag{ 12 }$$

with $\alpha =\alpha_{0}$ being time-independent and

$$\begin{eqnarray} \frac{\text{d}\varphi }{\text{d}t} &= \frac{1}{2}\left(\frac{\partial E_{L,-}}{\partial\alpha }\cot \alpha -\frac{\partial E_{R,-}}{\partial \alpha }\tan \alpha \right) \nonumber \\ &\quad +\,\sum_{j,\eta }\epsilon_{j}\biggl(p_{j,\eta }\frac{\text{d}A_{j,\eta }}{\text{d}t}-\frac{1}{2}\omega_{\eta }^{2}A_{j,\eta }^{2}-E_{j,-} \nonumber \\ &\quad -\,i\langle \zeta (\mathbf{p})|\partial_{\mathbf{p_{\eta }}}\zeta (\mathbf{p})\rangle \frac{\text{d}p_{\eta }}{\text{d}t}\biggl|_{{p_{\eta }}={p_{j,\eta }}}\biggr), \end{eqnarray} \tag{ 13 }$$

where $E_{j,-}=E_{-}(\bm{p_j})$ and $\epsilon_{j=L,R}=\pm 1$.

We choose the initial condition (at t = 0) of eq. (11) as $A_{x}=A_{y}=0$ for both wave packets (same starting point) and small $p_{L,x} = -p_{R,x}$ for the initial velocity under slight axisymmetric perturbation in momentum space. Such symmetry is actually preserved by the equations of motion. The initial condition for eq. (12) is obtained from the minimization of the system energy functional. The trajectories generated by those equations are presented in fig. 6, in good agreement with the GPE results in fig. 3. Note that we also assume the equal splitting of the condensate, or $\alpha = \pi/4$, given a sufficiently weak interaction.

Zoom In Zoom Out Reset image size

We turn to discussing the dynamics of phase φ of the variational wave function. There are three different contributions to the time derivative of φ —the first comes from energy terms $E_{j,-}$, the second is related to the dynamic parameters like center-of-mass position and momentum, and the last is just the Berry connection $\bm{A}_\text{B}$, determined together by the spinor wave function and time derivative of the momentum. Note that with the axial symmetry to the y (and p_y) axis, all the terms related to energy and other dynamical parameters, such as $A_j^2$ and p_j, vanish, leaving only the last term on the RHS of eq. (13), which becomes

$$\begin{eqnarray} &&\frac{\text{d}}{\text{d}t} \varphi = - i \sum_{j,\eta} \epsilon_j \left( \langle \zeta(\mathbf{p}) | \partial_ {\mathbf{p_ \eta}} \zeta(\mathbf{p}) \rangle \frac{\text{d}p_\eta}{\text{d}t} \biggl|_{{p_\eta}={p_{j,\eta}}} \right). \end{eqnarray} \tag{ 14 }$$

This is exactly the Berry phase defined on the ring loop since

$$\begin{equation} \oint_{L_c} \bm{A}_\text{B} \cdot \text{d}\bm{p} = - \int_{L_{cL}} \bm{A}_\text{B}\cdot \text{d}\bm{p} + \int_{L_{cR}} \bm{A}_\text{B} \cdot \text{d}\bm{p}, \end{equation} \tag{ 15 }$$

where the ring-shaped loop L_c is divided into two parts L_cj, as for each part of the splitting condensate. The integral direction of left-hand side (L_cL) is clockwise, hence carrying a minus sign. The numerical solutions of φ are illustrated in fig. 6(a1) and (b1) as the curve color. In panel (a1), the two wave packets acquire an opposite geometric phase as they encircle the loop. When they encounter each other at point B, the accumulated phases are $\pm\frac{\pi}{2}$, respectively, making a π phase difference, which results in a dark fringe. However, with an extra Zeeman field $\delta_z$, such phase difference never reaches π at point B ((a2)). The relative brightness R_I defined in eq. (8) is also evaluated by the variation analysis, and the results (green dots) are compared with those from GPE (blue dots) in fig. 4. They indeed show the same monotonic trend.

We remark that the variational method well captures the physical features of the interference dynamics as well as the geometric phase with much fewer variables than the GPE simulation. Solving the semiclassical equations of motion is also computationally efficient compared with the GPE simulation.

We have proposed and investigated a realistic approach for conducting momentum-space Aharonov-Bohm interferometry in Rashba spin-orbit–coupled Bose gases and shown that the interference pattern measures the Berry phase of the Rashba energy band. Our approach utilizes the ring structure of the Rashba spectrum as the interferometry loop and the ultracold atoms tunability for triggering the motion of BEC wave packets along the loop. With the real-time GPE simulation for realistic ⁸⁷Rb gases, we have found that the relative brightness of the interference fringes directly indicates the Berry phase as a monotonic function of Zeeman detuning. In particular, the π Berry phase of a Dirac point (without the detuning) is exhibited by a robust dark fringe at the end of interferometry loop. Additionally, we have modeled the interference dynamics with a variational wave function of splitting wave packets and derived semiclassical equations of motion for the most relevant dynamical factors. The variational results have confirmed the trajectories in both momentum and real spaces as well as the local geometric phase acquired by the condensate along the momentum trajectory. The complementary variational analysis and GPE simulation well agree with each other.

Our study would provide guidance for ongoing experimental effort measuring the Berry phase in ultracold atoms with synthetic Rashba spin-orbit coupling. The simulated density pattern in momentum and real space (figs. 2, 3, and 5) can be directly compared with the time-of-flight and direct-imaging measurements, respectively, in the experiment of ⁸⁷Rb gases. Our analysis can be extended to meet different experimental conditions as long as two requirements are fulfilled: 1) the BEC wave packets can be driven to move along a closed path in momentum space and 2) the accumulated dynamic phases along the two halves of the path cancel out (so the interference at the end of the path only reflects the geometric phase). A perfectly flat Rashba ring is not necessary. For example, if several local minima are present in the path (such as in the experimental system in ref. [50]), a sufficiently large detuning may help incline the loop and lift the minima, and a proper choice of parameters for the symmetric dynamic phase is also needed. In a wide scope, the interferometry approach may find applications on various nontrivial energy bands as well as high-spin [57–63] systems.

This work is supported by AFOSR (FA9550-16-1-0387), NSF (PHY-1505496), and ARO (W911NF-17-1-0128).