Ion friction at small values of the Coulomb logarithm

TL;DR

This study investigates ion collision rates in strongly-coupled plasmas with small Coulomb logarithms using experiments on ultracold neutral plasmas and simulations, challenging traditional kinetic theories.

Contribution

It provides the first experimental measurements and simulations of ion collisions in strongly-coupled plasmas at small Coulomb logarithms, using dual-species ultracold neutral plasmas.

Findings

Strong collisional coupling observed between ion species.

Good agreement between experiments and simulations.

Potential for testing kinetic theory expressions for Coulomb logarithms.

Abstract

Transport properties of high-energy-density plasmas are influenced by the ion collision rate. Traditionally, this rate involves the Coulomb logarithm, $\ln\Lambda$. Typical values of $\ln\Lambda$ are $\approx 10~\mbox{to}~20$ in kinetic theories where transport properties are dominated by weak-scattering events caused by long-range forces. The validity of these theories breaks down for strongly-coupled plasmas, when $\ln\Lambda$ is of order one. We present measurements and simulations of collision data in strongly-coupled plasmas when $\ln\Lambda$ is small. Experiments are carried out in the first dual-species ultracold neutral plasma (UNP), using Ca$^+$ and Yb$^+$ ions. We find strong collisional coupling between the different ion species in the bulk of the plasma. We simulate the plasma using a two-species fluid code that includes Coulomb logarithms derived from either a screened Coulomb potential or a the potential of mean force. We find generally good agreement between the experimental measurements and the simulations. With some improvements, the mixed Ca$^+$ and Yb$^+$ dual-species UNP will be a promising platform for testing theoretical expressions for $\ln\Lambda$ and collision cross-sections from kinetic theories through measurements of energy relaxation, stopping power, two-stream instabilities, and the evolution of sculpted distribution functions in an idealized environment in which the initial temperatures, densities, and charge states are accurately known.

Figures

Equations

$$\frac{\partial n_s}{\partial t}+u_s\frac{\partial n_s}{\partial r}$$

$$\frac{\partial T_s}{\partial t}+u_s\frac{\partial T_s}{\partial r}$$

$$\frac{\partial u_s}{\partial t}+u_s\frac{\partial u_s}{\partial r}$$

$$F_{s{s}^{\prime }}=\frac{{\mathbf{R}}^{s{s}^{\prime }}}{n_s}=-\frac{16}{3}\frac{\sqrt{\pi }{e}^4{n}_{s^{\prime }}}{(4\pi {ϵ}_0{)}^2{m}_{s{s}^{\prime }}{\bar{v}}_{s{s}^{\prime }}^3}\Xi (\Delta \overline{V})({\mathbf{u}}_{\mathbf{s}}-{\mathbf{u}}_{{\mathbf{s}}^{\prime }})$$

$$\Xi (\Delta \overline{V})=\frac{3}{16}\frac{1}{\Delta {\overline{V}}^3}\frac{1}{2}{\int }_0^{\infty }d\xi {\xi }^2\frac{{\sigma }_{s{s}^{\prime }}^{(1)}(\xi )}{{\sigma }_0}\mathcal{X},$$

$$\mathcal{X}=[$$

$$(2\xi \Delta \overline{V}-1)e^{-(\xi -\Delta \overline{V}{)}^2},$$

$${\sigma }_0=\frac{\pi e^4}{(4\pi {ϵ}_0{)}^2{m}_{s{s}^{\prime }}^2{\bar{v}}_{s{s}^{\prime }}^4}.$$

$$Q_{s{s}^{\prime }}=-\frac{16\sqrt{\pi }{n}_s{n}_{s^{\prime }}{e}^4{k}_B}{(4\pi {ϵ}_0{)}^2{m}_s^2{\bar{v}}_{s{s}^{\prime }}^3}\tilde{\Xi }(\Delta \overline{V})(T_s-T_{s^{\prime }})$$

$$-\frac{v_{Ts}^2}{{\bar{v}}_{s{s}^{\prime }}^2}\Delta \mathbf{V}\cdot {\mathbf{R}}^{s{s}^{\prime }}$$

$$\tilde{\Xi }(\Delta \overline{V})=\frac{1}{8\Delta \overline{V}}{\int }_0^{\infty }d\xi {\xi }^4\frac{{\sigma }_{s{s}^{\prime }}^{(1)}(\xi )}{{\sigma }_0}\left[e^{-(\xi -\Delta \overline{V}{)}^2}-e^{-(\xi +\Delta \overline{V}{)}^2}\right]$$

$$n_i^{m+1}$$

$$T_i^{m+1}$$

$$\frac{Q}{(2/3)\nabla \cdot (u_s\hat{r})}(1-e^{-(2/3)\nabla \cdot (u_s\hat{r})\tau })$$

$$u_i^{m+1}$$

$$+\left[-\frac{k_B}{n_s{m}_s}\frac{\partial n_s{T}_s}{\partial r}-\frac{k_B{T}_e}{n_s{m}_s}\frac{\partial n_s}{\partial r}+\frac{F_{s,drag}}{m_s}\right]\tau$$

$$\dot{r}\approx u_i+u_i^{\prime }(r-r_i)$$

$$u_i^{\prime }=\frac{u_{i+1}-u_{i-1}}{2\Delta r}$$

$$\delta r=\frac{u_i}{u_i^{\prime }}(1-e^{-u_i^{\prime }\tau })\approx u_i\tau -\frac{1}{2}{u}_i{u}_i^{\prime }{\tau }^2$$

$$G(r_i-\delta r)$$

$$\frac{\delta r^2}{2\Delta r^2}(G_{i+1}-2G_i+G_{i-1})$$

$$\Xi (\Delta \overline{V})=\frac{3\sqrt{\pi }}{4}\frac{\psi (\Delta {\overline{V}}^2)}{\Delta {\overline{V}}^3}\ln \Lambda ,$$

$$\psi (x)=\mathrm{erf}(\sqrt{x})-\frac{2}{\sqrt{\pi }}\sqrt{x}{e}^{-x}.$$

$$\tilde{\Xi }(\Delta \overline{V})=\frac{\sqrt{\pi }}{2\Delta \overline{V}}\mathrm{erf}(\Delta \overline{V})\ln \Lambda .$$

$$\ln \Lambda =\ln (1+C/g)$$

$$g=\frac{e^2}{4\pi {ϵ}_0}{\left[{\lambda }_D\left(\frac{1}{2}{m}_{s{s}^{\prime }}\right)\left({\bar{v}}_{s{s}^{\prime }}^2+\frac{2}{3}{\left|u_s-u_{s^{\prime }}\right|}^2\right)\right]}^{-1}$$

$$\frac{1}{{\lambda }_D^2}=\frac{1}{{\lambda }_e^2}+\sum _i\frac{1}{{\lambda }_i^2}\left(\frac{1}{1+(u_s-u_{s^{\prime }}{)}^2/v_{th,i}^2+3{\Gamma }_i}\right),$$

$${\Gamma }_i\equiv \frac{e^2}{4\pi {ϵ}_0{a}_{ws}}\frac{1}{k_B{T}_i}=2.3,$$

$$f(v_z)\propto$$

$${\int }_0^{\infty }\frac{n(r)}{u(r)}\left[\mathrm{erf}\left(\frac{u(r)-v_z}{v_{th}}\right)+\mathrm{erf}\left(\frac{u(r)+v_z}{v_{th}}\right)\right]r^{2}dr,$$

Full text

strongly-coupled plasmas, warm-dense plasmas, two-fluid & multi-fluid model, time-resolved light scattering spectroscopy, optical plasma measurements, randon & disordered media, low-temperature plasmas, ion-inertial confinement, atom & ion cooling

\doi

10.1103/PhysRevE.99.053206

I Introduction

Understanding energy transport and relaxation processes is an important aspect of optimizing plasma fusion experiments and determining their equations of state Fletcher et al. (2014); Sperling et al. (2017); Stanton and Murillo (2016). In some systems, the electron and ion temperatures can be different by an order of magnitude or more Lindl et al. (2014). As fusion proceeds, this energy difference is exacerbated when fusion products asymmetrically deposit their energy into the electron system due to the large mass ratio Vorberger and Gericke (2014). The resulting two-temperature problem is a long-standing plasma physics issue in many systems Scott et al. (2018); Dharma-wardana and Murillo (2008); Cooley et al. (2016); Dharuman et al. (2017); Haack et al. (2017); Benedict et al. (2017); Zastrau et al. (2014); Mabey et al. (2017); Shaffer et al. (2017); Christopherson et al. (2018a, b).

Kinetic calculations, based on the Boltzmann equation or one of its simplifications, rely on a statistical assumption about Coulomb collisions between charged particles. They are most accurate when the collisions are frequent and weak, corresponding to large impact parameters and small-angle scattering. Mathematical expressions for two-body processes such as the electron-ion collision rate or the ion-ion collision cross-section are modified by the Coulomb logarithm to account for the many-body, long-range nature of Coulomb collisions. In the Landau-Spitzer treatment, the Coulomb logarithm is conveniently written as $\ln\Lambda=\ln(\lambda_{D}/r_{0})$, where the electron Debye length is $\lambda_{D}=[\epsilon_{0}k_{B}T/(ne^{2})]^{1/2}$ and the classic distance of closest approach is $r_{0}=e^{2}/(4\pi\epsilon_{0}k_{\rm B}T)$. The Coulomb logarithm is a function of density and temperature and it multiplies every cross-section and collision rate in kinetic calculations Huba (2013). Its value ranges from 10 to 20 for weakly coupled plasmas, and represents an average over many long-range binary collisions in the plasma.

Recent theoretical and computational work extend the Landau-Spitzer treatment described above to higher density plasmas. These treatments use the screened Coulomb interactions in which the many-body physics lacking in binary models is accounted for through the choice of an ad-hoc effective screening length Liboff (1959); Muchmore (1984); Paquette et al. (1986); Stanton and Murillo (2016). Recently, the effective potential theory which models many-body correlation effects by treating binary interactions as arising through the potential of mean force rather than the screened Coulomb potential has been proposed in Ref. Baalrud and Daligault (2016). All these models have been used to calculate transport properties for plasmas across coupling regimes Stanton and Murillo (2015, 2016); Dimonte and Daligault (2008); Daligault (2012); Baalrud and Daligault (2013); Moldabekov et al. (2017). These are typically compared to results from molecular dynamics simulations, and, where possible, experimental data Affolter et al. (2016); Strickler et al. (2016); Liu and Goree (2014); Bannasch et al. (2012); Affolter et al. (2018).

Ultracold neutral plasmas (UNPs) span the phase space region where the Coulomb logarithm values are small. These systems have enabled studies of plasma dynamics and evolution in a highly idealized environment Killian et al. (2007); Lyon and Rolston (2017), serving, in a way, as high-energy-density plasma simulators Murillo (2004, 2006). They are generated by resonantly photo-ionizing mK-temperature atoms Killian et al. (1999); Rolston and Roberts (2003); Simien et al. (2004); Cummings et al. (2005a); Murphy et al. (2014) or molecules Morrison et al. (2008); Haenel et al. (2017). The initial plasma density ($10^{7}$ to $10^{12}~{}\mbox{cm}^{-3}$) and electron temperature (5 to 500 K) are selected with small uncertainties. Simulations Robicheaux and Hanson (2002); Pohl et al. (2004), laser spectroscopy Cummings et al. (2005b); Bergeson and Robicheaux (2008); Bannasch et al. (2012), radio-frequency measurementsRoberts et al. (2004); Wilson et al. (2013), charged particle detection and imaging Zhang et al. (2008); Kulin et al. (2000) are all used to characterize these systems. They have deepened our understanding of the time-evolving density Bergeson and Spencer (2003); Twedt and Rolston (2012), electron and ion temperatures Wilson et al. (2013); Lyon et al. (2013); Pohl et al. (2005); Kuzmin and O’Neil (2002); McQuillen et al. (2015); Langin et al. (2016), collision and recombination rates Fletcher et al. (2007); Lyon et al. (2015), expansion, velocity relaxation Bannasch et al. (2012), localization Sous and Grant (2018), and self-diffusion Strickler et al. (2016) in strongly-coupled Coulomb systems. The recent realization of laser-cooling ions in an ultracold neutral plasma open the possibility of extending all of these studies farther into the strongly coupled plasma regime Langin et al. (2019).

In this article we present the first dual-species UNP, using Ca+ and Yb+ ions. We report simulations and measurements of the Ca+ ion velocity distribution in the dual-species UNP. We also present a two-fluid model with two representations of the friction force between the ions. Those simulations reproduce the main features of the measured ion velocity distribution. This system provides a unique platform for future studies of collision physics in strongly-coupled plasmas. In this system it should be possible to study idealized versions of classic plasma problems such as inter-species diffusion Amendt et al. (2010), multi-species plasma expansion Bin et al. (2018), two-stream instabilities, the sensitivity of bump-on-tail evolution to electron screening Aslanyan et al. (2017), shock evolution Marciante and Murillo (2017); Shang et al. (2017), and evaluations of the Coulomb logarithm when the plasma approaches the non-ideal state Baalrud and Daligault (2013); Stanton and Murillo (2016).

II Experimental description

We simultaneously trap $10^{7}$ neutral 40Ca and 174Yb atoms in a magneto-optical trap (MOT) Raab et al. (1987) at a temperature of a few mK Cummings et al. (2005a). The spatially-overlapped MOTs operate on the strong resonance transitions at 423 and 399 nm for Ca and Yb, respectively. Unlike dual-species MOTs with alkali atoms or combinations of alkali and alkaline-earth atoms, the 40Ca and 174Yb atoms occupy the same physical space in the MOT without adversely influencing the number of trapped atoms of either species because there is no ground-state hyperfine structure Busch et al. (2006); McCarron et al. (2011); Extavour et al. (2006); Ridinger et al. (2011). The spatial density profile is Gaussian, $n=n_{0}\exp(-r^{2}/2\sigma^{2})$. In order to minimize spatial inhomogeneities in the neutral atom clouds stemming from imperfect laser beams, the neutral atoms are allowed to expand for $100~{}\mu\mbox{s}$ before formation of the plasma.

The neutral atoms in the MOT are resonantly-ionized using ns-duration laser pulses in a two-step process. The initial electron temperature in the plasma is determined by varying the wavelengths of the 390 nm (Ca) and 395 nm (Yb) laser pulses. In experiments reported here, the electron temperature is $T_{e}=96~{}\mbox{K}$. The ion densities are determined by varying the intensity of these same laser pulses. With our few-mJ pulses, we can ionize all of the Ca atoms and up to 60% of the Yb atoms. The peak density of the Ca+ plasma is $n_{0}=1.8\times 10^{10}~{}\mbox{cm}^{3}$ with an initial rms size of $\sigma_{0}=0.29~{}\mbox{mm}$. The peak density of the Yb+ plasma in the experiments reported here varies from $n_{0}=0.2\times 10^{10}~{}\mbox{cm}^{-3}$ to $n_{0}=1.8\times 10^{10}~{}\mbox{cm}^{-3}$ with an initial rms size of $\sigma_{0}=0.37~{}\mbox{mm}$.

The electron temperature, $T_{e}$, is determined by the excess photon energy above the atomic ionization potential. Because the electron energy drives the plasma expansion rate, we must characterize $T_{e}$ accurately. A partial energy level diagram for Ca and Yb is shown in Fig. 1. We adjust the wavelength of the 390 nm laser ($\lambda^{-1}=23754~{}\mbox{cm}^{-1}$) so that it ionizes Ca atoms out of the $4s4p~{}^{1}P_{1}$ level ($E_{4s4p}^{\rm Ca}=23654~{}\mbox{cm}^{-1}$), imparting $100~{}\mbox{cm}^{-1}$ of kinetic energy to the electrons. However, this same laser, if it is coincident with the 399 nm Yb excitation laser, would ionize Yb atoms out of the $6s6p~{}^{1}P_{1}$ level ($E_{6s6p}^{\rm Yb}=25068~{}\mbox{cm}^{-1}$), imparting $279~{}\mbox{cm}^{-1}$ of kinetic energy to the electrons. We minimize this problem in two ways. First, we use fast optical modulators to turn off all of the MOT laser beams 100 $\mu$s before ionizating the MOT. This ensures that none of the Ca and Yb atoms are in the ${}^{1}P_{1}$ states when the ionizing laser pulses arrive at the MOT. Second, we delay the Ca ionization pulses by 40 ns relative to those for Yb. Because of the larger Yb mass, this delay is short enough that the Yb+ plasma does not expand before the dual-species plasma is formed. We verify that this delay is equal to 5 times the measured pulse width of the laser pulses and that these precautions prevent spurious ionization to a level below our detection sensitivity.

The time-evolving rms width of the Ca+ ion velocity distribution is determined using laser-induced fluorescence measurements at 397 nm Kramida et al. (2018). A linearly-polarized probe laser beam at this wavelength passes through the plasma and is retro-reflected. The single-beam intensity is $I=50~{}\mbox{mW/cm}^{2}\approx I_{\rm sat}$. A strong laser beam ($I=2000~{}\mbox{mW/cm}^{2}$) at 850 nm is used to minimize optical pumping of the Ca+ ions into dark states. The size of these laser beams is large compared to the size of the plasma. Both laser beams illuminate the entire plasma for the duration of the experiment.

When the frequency of the probe laser is detuned by a frequency $\Delta\omega=2\pi\times(f-f_{0})$ from the atomic resonance frequency, $f_{0}$, the fluorescence signal is proportional to the number of ions Doppler-shifted into resonance with the laser beam. By repeating fluorescence measurements for a range of different probe laser frequencies, we are able to map out the $v_{z}$ velocity distribution as a function of time, averaged over the entire plasma.

III Two-fluid simulation

The dual-species UNP environment is highly collisional. The Vlasov equation, which has modeled single-species UNP expansion with high accuracy Laha et al. (2007) is not appropriate. When the low-mass Ca+ ions expand in the presence of the heavier Yb+ plasma, ion friction transfers momentum from calcium to ytterbium, dramatically changing the behavior of the plasma.

To interpret the experimental results, we have built a two-fluid 1-d code in spherical coordinates. A kinetic treatment is almost certainly required, but such a calculation is difficlt and the mean free path in the bulk of the plasma is small enough that a fluid treatment should give some insight. The two species are denoted by the subscript $s$ and it is assumed that the important physical effects are convection, adiabatic expansion, pressure acceleration, acceleration by an ambipolar electric field, and interspecies friction, including Joule heating due to the relative velocity between the two species. Our plasma is not very strongly coupled, so we assume that the monatomic ideal gas law is the equation of state for both species, so that each one has a distribution function approximated by a drifing Maxwellian and so that each species has adiabatic exponent $\gamma=5/3$. For our conditions viscous effects, the thermal force, and ion thermal conduction are small and are not included in the fluid equations. With these assumptions the three equations to be solved for each species are

[TABLE]

In these equations $(n_{s},T_{s},u_{s})$ are, respectively, the density, temperature, and radial fluid velocity for species $s$. The quantity $F_{ss^{\prime}}$ is the interspecies friction force and $Q_{ss^{\prime}}$ is a term representing frictional heating and temperature equilibration between the two species.

To compute the interspecies friction force $F_{ss^{\prime}}$ and the heating term $Q_{ss^{\prime}}$, we follow the treatment of Baalrud and Daligault in Ref. Baalrud and Daligault (2016), including the energy exchange density in Eqs (44) through (51) of that reference. Our friction force $F_{ss^{\prime}}$ is given in terms of their fluid friction force density ${\bf R}^{ss^{\prime}}$ by

[TABLE]

where $m_{ss^{\prime}}$ is the reduced mass $m_{ss^{\prime}}=m_{s}m_{s^{\prime}}/(m_{s}+m_{s^{\prime}})$ and where $\bar{v}_{ss^{\prime}}=(2k_{B}T_{s}/m_{s}+2k_{B}T_{s^{\prime}}/m_{s^{\prime}})^{1/2}$. The quantity $\Delta\bar{V}=|{\bf u}_{s}-{\bf u}_{s^{\prime}}|/\bar{v}_{ss^{\prime}}$, where ${\bf u}_{s}$ and ${\bf u}_{s^{\prime}}$ are the species fluid velocities. The particle velocities ${\bf v}_{s}$ and ${\bf v}_{s^{\prime}}$ of the two species are assumed to be distributed according to two Maxwellians flowing relative to each other with relative velocity $\Delta V={\bf u}_{s}-{\bf u}_{s^{\prime}}$.

The quantity $\Xi(\Delta\overline{V})$ is a generalized Coulomb logarithm and is given by

[TABLE]

where the function ${\cal{X}}$ is

[TABLE]

where $\sigma_{ss^{\prime}}^{(1)}(\xi)$ is the usual first momentum transfer cross-section Baalrud and Daligault (2016), $\xi$ is the ratio of the particle velocity $v_{s}$ to the thermal velocity $v_{Ts}=\sqrt{2k_{B}T_{s}/m_{s}}$, and where

[TABLE]

Once this friction force is computed we use it in Eqs. (2) and (3) of the fluid model.

Baalrud and Daligault compute the energy exchange and frictional heating term $Q_{ss^{\prime}}$ similarly. They find

[TABLE]

[TABLE]

where

[TABLE]

This term may then be used in Eq. (2) of the fluid model.

The code is built on a cell-centered spherical grid with $r_{i}=(i-\textstyle{\frac{1}{2}})\Delta r$, $i=1,2,3,...$, with $r$ the spherical radial coordinate and with $\Delta r$ the constant grid spacing. We solve these equations using the method of characteristics.

In what follows $i$ denotes the spatial position on the radial grid and $m$ indicates time step in equal time increments $\tau$.

[TABLE]

We reach back in time from $r_{i}$ to $r_{i}-\delta r$ in order to find the quantity to be convected forward using the approximate characteristic equation

[TABLE]

where $u_{i}$ is the fluid velocity at radial grid point $i$ and where $u_{i}^{\prime}$ is the centered approximation to the radial derivative of the fluid velocity at grid point $i$,

[TABLE]

Solving Eq. (13) to find the radius from which the density is convected to $r_{i}$ at $t^{m+1}$ yields for $\delta r$ in $r_{i}-\delta r$

[TABLE]

If $G$ represents any of the quantities $(n,T,v)$ evaluated at the retarded position $r_{i}-\delta r$, then

[TABLE]

To handle the non-convective parts of the time advance a simple two step predictor-corrector method is used. In the first step old values of $(n_{s},T_{s},u_{s})$ are used to advance to time level $m+1/2$. In the second step these intermediate values are used to advance $(n_{s},T_{s},v_{s})$ to time level $t^{m+1}$.

III.1 The Coulomb logarithm and momentum transfer

We have studied two treatments of momentum transfer. The first uses the usual Coulomb cross section, modified by a suitable generalization of $\ln\Lambda$ Dimonte and Daligault (2008); Daligault and Dimonte (2009); Baalrud and Daligault (2013). As shown in Ref. Baalrud and Daligault (2016), for this case the friction force generalized Coulomb logarithm $\Xi(\Delta\overline{V})$ is given by

[TABLE]

where

[TABLE]

Similarly, the Baalrud-Daligault effective Coulomb logarithm for energy exchange $\tilde{\Xi}(\Delta\overline{V})$ in the case of Coulomb scattering with a Coulomb logarithm multiplier is given by

[TABLE]

For the case of electron-ion temperature relaxation, molecular dynamics simulations Dimonte and Daligault (2008) indicate that a Coulomb logarithm of the form,

[TABLE]

is appropriate, where $C=0.7$, and where $g=(e^{2}/4\pi\epsilon_{0})[1/(\lambda_{De}k_{B}T_{e})]$ is the so-called plasma parameter. Effective potential theory calculations suggest that this might be appropriate for our dual-species plasma as well Baalrud and Daligault (2014).

Because we are calculating ion-ion momentum transfer, some caution is in order. In the NRL Plasma Formulary, the plasma parameter is $g=r_{\rm min}/r_{\rm max}$. For ion-ion collisions in flowing Maxwellians,

[TABLE]

where the Debye length $\lambda_{D}$ includes both the ion contribution and a correction due to ion flow and strong coupling, as given in Ref. Stanton and Murillo (2016),

[TABLE]

where the summation is over the ion species and where $v_{th,i}=(2k_{B}T_{i}/m_{i})^{1/2}$. Consistent with Ref. McQuillen et al. (2015) and many other UNP studies, we take the ion strong coupling parameter to be,

[TABLE]

where $a_{ws}=[3/(4\pi n_{0})]^{1/3}$ is the Wigner-Seitz radius. Near the center of the plasma, where the density is the highest, the value of the plasma parameter is $g=2.6$.

The second form for the momentum transfer cross section uses the Debye-screened Coulomb potential described by Stanton and Murillo in Section III and Appendix C, Eq. (41), of Ref. Stanton and Murillo (2016). The collision integrals in this reference are expressed as convenient functions of the plasma parameter, $g$, as discussed above. Using the screened Coulomb cross section cited aboveStanton and Murillo (2016), the integral in Eq. (5) was performed numerically and fit to an analytic form for use in the fluid code. In this treatment the energy exchange term is density-weighted, as opposed to velocity-weighted in the Baalrud-Daligault treatment, and the temperature equilibration term was neglected since its effect turned out to be small when comparing the simulation to the experiment.

III.2 The velocity distribution $\bm{f}(v_{z})$

In order to compare the fluid code directly with the experimental measurements, we calculate the $v_{z}$ velocity distribution from the simulation results for $u(r)$, $n(r)$, and $v_{th}(r)$. In doing so we assume that the particles of each species are drifting Maxwellians with the parameters given in the previous sentence and we integrate over all 3 dimensions in space and over $v_{x}$ and $v_{y}$ in velocity space to obtain the following distribution in $v_{z}$:

[TABLE]

where $v_{th}=(2k_{B}T/m)^{1/2}$. We evaluate this integral numerically. Note that normalizing constants have been omitted since the experimental data are not normalized.

IV Comparison of the simulation to the lab data

Before comparing the simulation results to the laboratory measurements we modify the raw data in two ways. First, the simulated velocity distribution is convolved with a Lorentzian line shape. The laboratory data infers the velocity distribution using fluorescence measurements, and those measurements necessarily include a contribution from the natural line width of the atomic transition. The second modification corrects the laboratory measurements for optical pumping. At a given time $t$, the data are multiplied by $\exp(t/\tau)$, where $t$ is the time for which the comparison is made and $\tau$ is estimated to be $30~{}\mu\mbox{s}$.

In Figs. 2 and 3 we show the simulated and measured $f(v_{z})$ distributions at 0.2, 0.5, 1, 2, and 3 $\mu$s after the plasma is formed. The laboratory data are shown in black circles with errorbars indicating the rms noise in the measurements. The red solid line shows the simulated distribution using the effective Coulomb logarithm in Eq. (20) with $C=0.7$. The blue dashed line shows the distribution from the momentum transfer treatment of Ref. Stanton and Murillo (2016). As can be seen in Fig. 2, these two treatments are in good agreement with each other. The uncertainties in the experimental measurements stemming from optical pumping allow both of these simulated distributions to be in agreement with the measurements. We have run the simulation for a range of $C$ values in Eq. (20). As the value of $C$ increases, the velocity distributions at late times become narrower and the relative values of $f(v_{z}=0)$ fall less quickly. When $C=1.1$, the velocity distributions generated using the two simulations appear to be in perfect agreement.

A comparison of the data on a logarithmic scale, shown in Fig. 3, gives a better view of the wings of the distributions. At 0.2 $\mu$s the simulated and measured distributions agree. After 1.0 $\mu$s, the wings of the measured distribution clearly rise above that of the simulations. At later times, that difference becomes less pronounced for the range of velocities that can be measured at present. This difference in the wings is most likely due to kinetic effects not included in the fluid simulation. This is expected, since the mean free path in the edge of the plasma is an appreciable fraction (about 1/3) of the plasma radius. In the wings of the spatial density distribution, the lighter Ca+ ions would be accelerated quickly outwards by the persistent density gradient of the heavier Yb+ ions. Because of the lower density, the friction force would be small. This hypothesis could be tested in the lab using spatial imaging techniques, and that experiment is currently underway.

From the simulation, we can extract information that is not experimentally accessible. During the $3~{}\mu\mbox{s}$ expansion, the central Yb and Ca densities falls from the initial value of $1.8\times 10^{10}~{}\mbox{cm}^{-3}$ to $0.8\times 10^{10}~{}\mbox{cm}^{-3}$ for Yb and $0.6\times 10^{10}~{}\mbox{cm}^{-3}$ for Ca. The electron temperature, which is assumed to be spatially uniform, falls from 96 K to 45 K.

The ion radial velocity, density, and temperature depend on the radial coordinate $r$. In the simulation, all three of these parameters develop features near the outside edge of the density distribution. In Fig. 4(a) we show the (radial) flow velocity of the Ca+ and Yb+ ions as a function of $r$ after $3.1~{}\mu\mbox{s}$. Near the center of the plasma, where the densities are high, the flow velocities match. As the density falls off, the flow velocities diverge.

At the location of the flow-velocity divergence, we also see changes in the ion temperature and density. These quantities are plotted in Fig. 4(b) and (c). The limitations in the simulations suggests these data should not be taken literally. However, these data suggest a cold interior surrounded by a warm shell that is itself surrounded by a cold exterior. We are in the process of setting up an experiment to measure this directly. If these predictions are confirmed, it may be possible to determine $\Gamma$-dependent transport properties in this binary plasma mixture.

V Discussion

We present measurements and simulations of the Ca+ ion velocity distribution in an dual-species UNP of Ca+ and Yb+ ions. The simulation uses two treatments for momentum transfer. One is based on the effective potential approach using the potential of mean force. It uses a Coulomb logarithm extracted from MD simulations. The other uses the momentum transfer treatment of Ref. Stanton and Murillo (2016), based on the interaction of Yukawa-screened charges. These momentum transfer treatments are included in the friction force between two flowing Maxwellian distributions of Ca+ and Yb+ ions at the same temperature.

Each of these treatments result from different assumptions. The effective potential approach using the potential of mean force is most appropriate for near-equilibrium processes. It assumes that the ion pair distribution function can be calculated using thermodynamic considerations. After the initial disorder-induced heating process in the UNP, this assumption is almost certainly valid. In the present work, we have used a Coulomb logarithm derived from MD simulations based on the effective potential approach. Because that work studied electron-ion energy relaxation, that Coulomb logarithm should be appropriate for momentum transfer processes. However, because we are studying momentum transfer between ions of different mass and not between electrons and ions, we have modified the parameters used in the Debye length and in the plasma parameter, as described previously. Future work should use the full theory in order to remove the approximations used in the present study.

The treatment of Ref. Stanton and Murillo (2016) assumes that the ion-ion potential is accurately represented using a Yukawa-screened interaction. For small values of $\Gamma$ (or large values of $g$) when the plasma is not strongly-coupled, this approach works nicely. As the plasma becomes strongly-coupled, the screening length is modified and Yukawa screening is assumed to be largely correct.

It might be useful to compare these approaches with laboratory data in which the coupling parameter $\Gamma$ is larger, where the underlying assumptions in the theoretical approaches could be tested more directly. For example, one can imagine an experiment in which the Ca plasma is generated and allowed to expand. As it expands, the ion-ion coupling parameter increases to values near 5 McQuillen et al. (2015). The co-located Yb plasma could then be generated and the interaction between the cold Ca+ ions and the hot Yb+ ions could be measured. Alternatively, very recent work demonstrated the successful laser-cooling of ions in a strontium UNP Langin et al. (2019). That method could be used to reach $\Gamma=11$ in one species while observing collisions and interactions due to the presence of the other, with the perturbing plasma either at the same or at higher temperatures.

On the other hand, experiments at higher $\Gamma$ values might produce only trivial transport results. Strong coupling should result in small values of the transport coefficients due to the greater collisional locking or caging of the plasma ions. Such predictions could be verified in our dual-species ultracold neutral plasma. Either way, the laser cooling or heating in UNPs shown in Ref. Langin et al. (2019) could provide a convenient way to measure $\Gamma$-dependent transport in a tightly-controlled environment.

The simulations show that the Ca+ flow velocity matches the Yb+ flow velocity in the center of the plasma. One can envision an experiment in which that flow could be interrupted, or in which two somewhat spatially-offset plasmas could approach equilibrium, flowing over each other. Perhaps flow-related instabilities could be observed. The simulations also show ion heating in regions where the spatial density gradient increases. Experiments are underway now to study this effect.

VI Conclusion

In conclusion, we report the first experimental realization of a two-ion species ultracold neutral plasma, comprised of Yb+ and Ca+ ions and electrons. The mm-sized spherical plasma is not confined but expands radially under the influence of the ambipolar field. We measure the spatially-averaged velocity distribution of the Ca+ ions as the plasma expands and observe that the Yb+ ions significantly slow the rate at which the distribution broadens. This results from momentum transfer between the Yb+ and Ca+ ions in the plasma.

We compare these measurements with the output of a fluid-code simulation. In the simulation we use two different expressions for the momentum transfer cross section. One derives from the potential of mean force. The other derives from a coupling-corrected Yukawa interaction. In the fluid equations, momentum transfer manifests most strongly as a friction force between flowing (assumed) Maxwellian distributions of Ca+ and Yb+ ions at the same temperature. The momentum transfer cross section traditionally includes a Coulomb logarithm, that, in our system, has a value less than 1.

We find that both formulations of the momentum transfer cross section, when included in our fluid simulation, produce nearly identical radial velocity distributions. The main features of the simulated distribution match the measured velocity distribution well. Some behavior in the wings of the distribution are noted, perhaps due to kinetic effects outside of the fluid code assumptions. It is possible that differences between these two momentum transfer treatments might appear if spatially-resolved measurements could be made in plasmas of varying levels of strong Coulomb coupling. It is possible that higher fidelity simulations might also reveal differences. Experiments and calculations are currently underway to test these ideas.

VII Acknowledgements

This project was supported by grants from the National Science Foundation (Grant No. PHY-1500376) and the Air Force Office of Scientific Research (AFOSR FA9550-17-1-0302). The authors express appreciation to Scott Baalrud and Jerome Daligault for helpful discussions concerning the effective potential approach.