Core filling and snaking instability of dark solitons in spin-imbalanced superfluid Fermi gases

TL;DR

This paper investigates how unpaired spins in spin-imbalanced superfluid Fermi gases affect dark solitons, revealing core filling effects and suppression of snake instability through time-dependent Bogoliubov de Gennes simulations.

Contribution

It demonstrates the impact of spin imbalance on dark soliton stability and shape, highlighting the role of unpaired particles in core filling and instability suppression.

Findings

Unpaired spins fill the soliton cores, broadening them.

Spin imbalance suppresses the transverse snake instability.

Potential experimental observation in cold atom systems.

Abstract

We use the time-dependent Bogoliubov de Gennes equations to study dark solitons in three-dimensional spin-imbalanced superfluid Fermi gases. We explore how the shape and dynamics of dark solitons are altered by the presence of excess unpaired spins which fill their low-density core. The unpaired particles broaden the solitons and suppress the transverse snake instability. We discuss ways of observing these phenomena in cold atom experiments.

Equations

$$\begin{split}\hat{H}=&\int{d^{3}\vec{x}\Big{[}\sum_{\sigma=\uparrow,\downarrow}\Psi^{\dagger}_{\sigma}(\vec{x})(-\frac{\hbar^{2}\nabla^{2}}{2m}}-\mu_{\sigma})\Psi_{\sigma}(\vec{x})\\ &-g\Psi^{\dagger}_{\uparrow}(\vec{x})\Psi^{\dagger}_{\downarrow}(\vec{x})\Psi_{\downarrow}(\vec{x})\Psi_{\uparrow}(\vec{x})\Big{]}\end{split}$$

$$1/g=-m/(4\pi {\hslash }^2{a}_s)+\frac{1}{V}\sum _{k}1/(2{ϵ}_k)$$

$$\begin{split}&\hat{H}_{\rm{BdG}}=\\ &\int{\!\!d^{3}\vec{x}}\left(\begin{array}[]{ccc}\Psi_{\uparrow}(\vec{x})\\ \Psi^{\dagger}_{\downarrow}(\vec{x})\end{array}\right)^{\dagger}\!\!\left(\begin{array}[]{ccc}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu_{\uparrow}&\Delta(\vec{x})\\ \Delta^{*}(\vec{x})&\frac{\hbar^{2}\nabla^{2}}{2m}+\mu_{\downarrow}\end{array}\right)\!\!\left(\begin{array}[]{ccc}\Psi_{\uparrow}(\vec{x})\\ \Psi^{\dagger}_{\downarrow}(\vec{x})\end{array}\right)\end{split}$$

$$\left(\begin{array}[]{ccc}\Psi_{\uparrow}(\vec{x})\\ \Psi^{\dagger}_{\downarrow}(\vec{x})\end{array}\right)=\sum_{n}\left(\begin{array}[]{ccc}u_{n}(\vec{x})&-v^{*}_{n}(\vec{x})\\ v_{n}(\vec{x})&u^{*}_{n}(\vec{x})\end{array}\right)\left(\begin{array}[]{ccc}\gamma_{\uparrow,n}\\ \gamma^{\dagger}_{\downarrow,n}\end{array}\right)$$

$$\left(\begin{array}[]{ccc}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu&\Delta(\vec{x})\\ \Delta^{*}(\vec{x})&\frac{\hbar^{2}\nabla^{2}}{2m}+\mu\end{array}\right)\left(\begin{array}[]{ccc}u_{n}(\vec{x})\\ v_{n}(\vec{x})\end{array}\right)=E_{n}\left(\begin{array}[]{ccc}u_{n}(\vec{x})\\ v_{n}(\vec{x})\end{array}\right)$$

$$\Delta (\vec{x})=g\sum _{E_n>0}{u}_n(\vec{x})v_n^{\ast }(\vec{x})(1-\Theta (-E_{\uparrow ,n})-\Theta (-E_{\downarrow ,n}))$$

$$\begin{split}&n_{\sigma}=\langle\Psi^{\dagger}_{\sigma}(\vec{x})\Psi_{\sigma}(\vec{x})\rangle\\ &=\sum_{E_{n}>0}\big{[}|u_{n}(\vec{x})|^{2}\Theta(-E_{\sigma,n})+|v_{n}(\vec{x})|^{2}(1-\Theta(-E_{-\sigma,n}))\big{]}\end{split}$$

$$n_I=\frac{n_{\uparrow }(x=0)-n_{\downarrow }(x=0)}{n_o}$$

$$\left(\begin{array}[]{ccc}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu&\Delta(\vec{x},t)\\ \Delta^{*}(\vec{x},t)&\frac{\hbar^{2}\nabla^{2}}{2m}+\mu\end{array}\right)\left(\begin{array}[]{ccc}u_{n}(\vec{x},t)\\ v_{n}(\vec{x},t)\end{array}\right)=i\hbar\frac{\partial}{\partial t}\left(\begin{array}[]{ccc}u_{n}(\vec{x},t)\\ v_{n}(\vec{x},t)\end{array}\right)$$

$$\Delta (\vec{x},t)=g\sum _{E_n>0}{u}_n(\vec{x},t)v_n^{\ast }(\vec{x},t)(1-\Theta (-E_{\uparrow ,n})-\Theta (-E_{\downarrow ,n}))$$

$$\Delta (\vec{x},0)=\Delta (\vec{x})+ϵ\mathrm{sech}(x)\sin \frac{2\pi }{L_y}(y-L_y/4)$$

$$\widetilde{H}\left(\begin{array}[]{ccc}u_{n,k_{z}}(x,y,t)\\ v_{n,k_{z}}(x,y,t)\end{array}\right)=i\hbar\frac{\partial}{\partial t}\left(\begin{array}[]{ccc}u_{n,k_{z}}(x,y,t)\\ v_{n,k_{z}}(x,y,t)\end{array}\right)$$

$$\widetilde{H}=\left(\begin{array}[]{ccc}-\frac{\hbar^{2}}{2m}(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}-k_{z}^{2})-\mu&\Delta(x,y,t)\\ \Delta^{*}(x,y,t)&\frac{\hbar^{2}}{2m}(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}-k_{z}^{2})+\mu\end{array}\right)$$

$$\exp [-i\tilde{H}\delta t]\approx$$

$$\exp [-i{H}_{\Delta }\delta t/2]\exp [-i{H}_{\mathrm{KE}}\delta t]\exp [-i{H}_{\Delta }\delta t/2]+O(\delta t^3)$$

$$H_{\mathrm{KE}}=$$

$$\displaystyle\left(\begin{array}[]{ccc}-\frac{\hbar^{2}}{2m}(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}-k_{z}^{2})-\mu&0\\ 0&\frac{\hbar^{2}}{2m}(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}-k_{z}^{2})+\mu\end{array}\right)$$

$$H_{\rm{\Delta}}=\left(\begin{array}[]{ccc}0&\Delta(x,y,t)\\ \Delta^{*}(x,y,t)&0\end{array}\right)$$

$$\left(\begin{array}[]{ccc}u_{n,k_{z}}(t+\delta t)\\ v_{n,k_{z}}(t+\delta t)\end{array}\right)=\exp[-iH_{\Delta}\delta t/2]\mathcal{F}^{-1}\exp[-iH_{\rm{KE}}\delta t]\mathcal{F}\exp[-iH_{\Delta}\delta t/2]\left(\begin{array}[]{ccc}u_{n,k_{z}}(t)\\ v_{n,k_{z}}(t)\end{array}\right)$$

Full text

I Introduction

Ultracold atoms have become the best platform for studying collective nonlinear phenomena such as dark solitons. Dark solitons are persistent nonlinear collective excitations in which the density is reduced in a plane. They have been studied in a number of physical settings including atomic Bose-Einstein Condensates (BECs) Frantzeskakis (2010) and superfluid Fermi gases of spin-1/2 atoms Antezza et al. (2007); Liao and Brand (2011); Scott et al. (2011); Yefsah et al. (2013). They are ubiquitous in quenches Lamporesi et al. (2013); Donadello et al. (2014) and can be engineered through phase imprinting protocols Burger et al. (1999); Denschlag et al. (2000); Anderson et al. (2001); Becker et al. (2008); Yefsah et al. (2013); Ku et al. (2014, 2016). Previous experimental Tikhonenko et al. (1996); Anderson et al. (2001); Dutton et al. (2001); Donadello et al. (2014); Ku et al. (2014, 2016) and theoretical work Zakharov and Rubenchik (1974); Jones et al. (1986); Kuznetsov and Turitsyn (1988); Mamaev et al. (1996); Muryshev et al. (1999); Feder et al. (2000); Musslimani and Yang (2001); Brand and Reinhardt (2002); Kamchatnov and Pitaevskii (2008); Brazhnyi and Pérez-García (2011); Cetoli et al. (2013); Mateo and Brand (2014); Bulgac et al. (2014); Reichl and Mueller (2013); Lombardi et al. (2016a) has discovered that these dark solitons are dynamically unstable to a “snaking” instability transverse to the plane of the soliton in both BECs and Fermi gases. In this paper, we theoretically study the dynamics of dark solitons in superfluid Fermi gases in which there is an imbalance between the number of up and down spins in the system. We find that the snaking instability is suppressed by the presence of excess spins which reside within the density depleted plane–or core–of the soliton.

Previous work Cetoli et al. (2013) has studied the snaking instability in spin-balanced Fermi gases using similar approaches as our paper. More recently, the authors of Refs.  Lombardi et al. (2016b, a) applied an effective field theoretic approach Klimin et al. (2015) to studying core filling and snaking instabilities of dark solitons in imbalanced Fermi gases. In this paper, we take a more microscopic approach and model the Fermi gas using the Bogoliubov de Gennes (BdG) equations. Our work also extends recent simulations of the stability of one-dimensional soliton trains Dutta and Mueller (2016) which suggest that excess spin can stabilize dynamical instabilities of dark solitons.

The BdG theory captures the phenomenology of the BEC-BCS crossover Zwerger (2011): at strong attractive interactions the fermions form tightly bound bosonic pairs which condense into a BEC, while at weak interactions the fermions form Cooper pairs which form a neutral analog to a BCS superconductor. Here we study the unitary gas which lies between these two limits. We caution that the BdG theory is a mean-field theory which only approximately models strong correlation physics in the unitary gas. However, the BdG equations has been successfully utilized in previous studies of dark soliton profiles and dynamics in the unitary gas Antezza et al. (2007); Liao and Brand (2011); Scott et al. (2011); Spuntarelli et al. (2011); Scott et al. (2012) and appears to be semi-quantitative.

In Sec. II we discuss the BdG model and find stationary dark solitons in the presence of imbalance. In Sec. III we apply time-dependent BdG equations to simulate the snaking instability. In Sec. IV we discuss how our results might be observed in cold atom experiments.

II Stationary Dark Solitons

II.1 Model

We consider the following Hamiltonian which describes spin-imbalanced spin-1/2 fermions with short-range attractive interactions

[TABLE]

Here $\mu_{\sigma}$ is the chemical potential for spin component $\sigma$ and $g$ is the bare interaction strength, which is related to the s-wave scattering length $a_{s}$ by

[TABLE]

where $V$ is the volume of the system and $\epsilon_{k}=\frac{\hbar^{2}k^{2}}{2m}$. The sum $\sum_{k}1/(2\epsilon_{k})$ is taken over the discrete momenta determined by the grid spacing of our numerics and comes with a momentum cutoff of $k_{c_{\nu}}=N_{\rm{grid}}\pi/L_{\nu}$ (with high energy cutoff $E_{c}=\hbar^{2}k_{c}^{2}/2m$) where $L_{\nu}$ is the length of the system along the direction $\nu$ and $N_{\rm{grid}}$ is the number of grid points along that direction. In this paper we focus our attention on the unitary limit $a_{s}\rightarrow\infty$.

At zero temperature, up and down spin atoms combine into Cooper pairs and condense to form a superfluid. We rewrite $\hat{H}$ in terms of the Cooper pair field $\Delta(\vec{x})=g\langle\Psi_{\uparrow}(\vec{x})\Psi_{\downarrow}(\vec{x})\rangle$ and neglect quadratic fluctuations. This gives the following mean-field BdG Hamiltonian De Gennes (1966)

[TABLE]

$H_{\rm{BdG}}$ is diagonalized through a Bogoliubov transformation

[TABLE]

where $\gamma^{\dagger}_{\sigma,n}$ is the creation operator for a Bogoliobov excitation of energy $E_{\sigma,n}=E_{n}\pm h$ where $E_{n}$ are the positive eigenvalues of the equation

[TABLE]

and where $h$ and $\mu$ are given by $h=\frac{1}{2}(\mu_{\uparrow}-\mu_{\downarrow})$ and $\mu=\frac{1}{2}(\mu_{\uparrow}+\mu_{\downarrow})$. At zero temperature, $\Delta(\vec{x})$ is expressed in terms of $u$’s and $v$’s as

[TABLE]

where $\Theta(x)$ is the unit step function. The density $n_{\sigma}(\vec{x})$ of fermions with spin $\sigma$ is given by

[TABLE]

II.2 Numerical results

We numerically solve the coupled equations (5) and (6) using an iterative procedure. We first choose an ansatz pair field $\Delta(\vec{x})=\Delta\tanh(x/\xi)$ corresponding to a planar dark soliton fixed at $x=0$. $\xi$ parametrizes the width of the soliton core and is generally chosen in the ansatz to be $\xi\approx k_{F}$, where $k_{F}\equiv(3\pi^{2}n_{o})^{1/3}$ is the Fermi wavevector and $n_{o}$ is the density far from the core of the soliton. We then solve Eq. (5) and calculate a new $\Delta$ from Eq.(6). This process is repeated until $\Delta$ converges to a stationary solution. In all the calculations presented in this paper we check that $\Delta$ converges to same stationary solution after small changes to the initial ansatz.

For simplicity we consider a system in a rectangular box geometry with dimensions $L_{x}\times L_{\perp}\times L_{\perp}$. We impose periodic boundary conditions in the $y$ and $z$ (perpendicular) directions, and in the $x$-direction we impose the conditions: $u_{n}(x+L_{x})=u_{n}(x)$ and $v_{n}(x+L_{x})=-v_{n}(x)$. This boundary condition in the $x$-direction ensures that $\Delta(x+L_{x})=-\Delta(x)$, which is consistent with the profile of a single dark soliton in a finite size box. Because of the homogeneity of the stationary soliton in the perpendicular directions, the solution to Eq. (5) can be expressed in the form $u_{n}(\vec{x})=u_{m,k_{y},k_{z}}(x)\exp(ik_{y}y+ik_{z}z)$, $v_{n}(\vec{x})=v_{m,k_{y},k_{z}}(x)\exp(ik_{y}y+ik_{z}z)$. This effectively reduces the three dimensional problem to a series of one-dimensional problems for each $k_{y}$ and $k_{z}$, and substantially speeds up the calculation.

Fig. 1 shows the total density $n_{\uparrow}+n_{\downarrow}$ (solid blue curve) and density difference $n_{\uparrow}-n_{\downarrow}$ (dashed orange curve) for stationary dark soliton solutions in the presence of spin imbalance. The densities are plotted as functions of $x$ after integrating over the $y$ and $z$ directions. In these calculations the dimensions of the box are set to $L_{x}k_{f}\approx 28$ and $L_{\perp}k_{f}\approx 23$, and we use 60 grid points along the $x$-direction and 50 $k$-space points in both perpendicular directions. These numerical parameters correspond to an energy cutoff $E_{c}\approx 100E_{F}$ in the sum in Eq. (2). Our results are unchanged by using more grid points.

We characterize the spin imbalance using the relative spin imbalance $n_{I}$:

[TABLE]

where as before $n_{o}$ is the total density far from the core. Fig. 1 shows a range of imbalances from $n_{I}=0$ to $n_{I}=0.33$. As the imbalance is increased, the soliton core (visually represented in the figure as the dip in the total density at $x=0$) fills with excess up spins and widens. This is consistent with previous calculations using different methods Lombardi et al. (2016b, a) and is expected given simple energetic considerations: the most energetically favorable place to store excess unpaired spins is at the core of the soliton where $\Delta=0$ and hence no Cooper pairs need to be broken. On a more microscopic level, the soliton supports a band of midgap Andreev states which are bound to the core of the soliton Antezza et al. (2007) which are filled by excess spins after tuning the chemical potential bias $h$ away from 0.

III Snaking Instability

In this section we discuss time dependent simulations of the snaking instability of dark solitons in the presence of spin imbalance. We find that the instability proceeds slower or, for sufficiently high imbalance, is completely suppressed by the presence of excess spins in the core of the soliton.

We numerically solve the following time-dependent BdG equations:

[TABLE]

where

[TABLE]

The initial set of $u_{n}(\vec{x},0)$ and $v_{n}(\vec{x},0)$ and the $E_{n}$’s in Eq. (10) are stationary solutions of Eqs. (5) and (6). For simplicity we again consider a system at unitarity ($a_{s}=\infty$) in a rectangular box geometry of dimensions $L_{x}\times L_{\perp}\times L_{\perp}$, and we use the same boundary conditions as in Sec. II.2. We assume homogeneity along the $z$-direction and express the $u_{n}$’s and $v_{n}$’s in the form $u_{n}(\vec{x},t)=u_{m,k_{z}}(x,y,t)\exp(ik_{z}z)$ and $v_{n}(\vec{x},t)=v_{m,k_{z}}(x,y,t)\exp(ik_{z}z)$. We use approximately 1000 grid points in the $x-y$ plane and 25 $k_{z}$ points which corresponds to an energy cutoff $E_{c}\approx 40E_{F}$.

In all the simulations described here, we first perturb the stationary $\Delta(x)$ by adding a small term which seeds a snaking instability along the y-direction:

[TABLE]

where $\epsilon\approx 0.1\Delta_{o}$ and $\Delta_{o}$ is the value of $\Delta$ far from the soliton core. We then discretize time and evolve the set of $u_{n}$’s and $v_{n}$’s forward by one time step using a split step method with $\Delta$ calculated from the current time-step. After finding the new $u_{n}$’s and $v_{n}$’s, we calculate $\Delta(t)$ at the next time-step using Eq. (10). Details of the split step method are discussed in the appendix.

Figure 2 shows the dynamics of a dark soliton for three different relative imbalances $n_{I}=0,0.18,0.4$ (columns (a), (b), (c), respectively). The figures show graphs of $|\Delta(x,y,z=0)|/\Delta_{o}$ at different times. In these graphs we have $L_{x}k_{F}\approx 28$ and $L_{\perp}k_{F}\approx 18$. At zero imbalance there is clearly a snaking instability whose rate is consistent with similar calculations in other work Cetoli et al. (2013). The plane of the soliton buckles and eventually breaks leaving behind two vortex cores. However at $n_{I}=0.18$, the instability occurs at a slower rate and finally at $n_{I}=0.4$ the instability is completely suppressed. We have run these simulations up to times of $t=150/E_{F}$, finding no sign of a snaking instability for large imbalances.

Figure 3 shows time scales for the snaking instability at different imbalances and different transverse dimensions $L_{\perp}$. Each time scale $\tau$ was calculated by first extracting the position of the core $x_{\textrm{core}}(t)$ along the $y=0$ line, and fitting it to a exponential function $x_{\textrm{core}}(t)\sim\exp(t/\tau)$. At zero imbalance the time scale gets larger as $L_{\perp}$ increases. This trend is somewhat counter-intuitive, but can be understood by noting that the unstable mode’s wavelength grows as $L_{\perp}$ increases. The instability connects with the Goldstone mode, and hence its frequency $\omega=i\hbar/\tau$ vanishes as $L_{\perp}\to\infty$. Similar results were seen in Ref. Cetoli et al. (2013). For larger $L_{\perp}$ (beyond those shown in this figure) additional decay modes appear. For sufficiently small $L_{\perp}$ ($\frac{\hbar^{2}}{mL_{\perp}^{2}}\sim\mu$) the rate will again decrease as the system becomes quasi-one dimensional. For small enough $L_{\perp}$ the soliton is stable, even without imbalance.

Once the soliton core is filled with excess spins the rate of the snaking instability becomes slower and is eventually suppressed altogether (up to times of at least $tE_{F}=150$). At smaller $L_{\perp}$ the snaking instability is suppressed for smaller values of spin imbalance.

IV Discussion

In this paper we have studied the dynamics of a dark soliton in an imbalanced Fermionic superfluid. We have found that the snaking instability of the soliton is suppressed by the presence of excess spins which reside at the low density core of the soliton. Although we have focused on the unitary limit, we expect these phenomena to persist in the BEC regime ($a_{s}\gg 0$) and BCS regime ($a_{s}\ll 0$). The behavior in the BEC regime should be similar to that of a two-component BEC: the paired superfluid maps onto one component, while the unpaired fermions in the core can crudely be modeled by a second component. Indeed, calculations Musslimani and Yang (2001); Brazhnyi and Pérez-García (2011) and experiments Anderson et al. (2001) of solitons in two-component BECs observe a suppression of the snaking instability.

We feel that experimentally observing the physics presented in this paper is feasible given existing tools. In a trapped imbalanced Fermi gas at equilibrium, excess spins reside along the edge of the trap Olsen et al. (2015). One naive idea is therefore to phase imprint a soliton onto the system and allow for the excess spins to diffuse from the edge of the trap into the core of the soliton. Unfortunately, the time scales for this process are prohibitively long. Instead we suggest first using a laser to create a potential barrier across the center of the trap and separating the imbalanced superfluid into two disjoint regions. This geometry was produced in Ref. Valtolina et al. (2015). Excess spins will then reside at the center of the trap between the two superfluid halves. Phase imprinting should then result in a soliton whose core is at the location of the excess spins. Varying the shape and dynamics of the applied potential barrier should allow experimentalists to control the relative imbalance $n_{I}$. One can image the soliton in time of flight after a ramp to the BEC limit as done in Refs. Yefsah et al. (2013); Ku et al. (2014, 2016).

V Acknowledgements

This material was based upon work supported by the National Science Foundation Grant No. PHY-1508300 and the ARO-MURI Non-Equilibrium Many-Body Dynamics Grant W9111NF-14-1-0003.

VI Appendix

In this appendix we discuss the details of the split step method used to perform the time-dependent BdG simulations discussed in Sec. III. Assuming homogeneity in the z-direction, the time-dependent BdG equations are given by

[TABLE]

where

[TABLE]

and $\Delta(x,y,t)$ is given by Eq. (10) in the main text. We approximate the operator $\exp[-i\widetilde{H}\delta t]$ which evolves the $u$’s and $v$’s forward by a time step $\delta t$ as:

[TABLE]

where

[TABLE]

and

[TABLE]

Because $H_{\rm{\Delta}}$ is diagonal in real space, applying the operator $\exp(-iH_{\Delta}\delta t/2)$ to \left(\begin{array}[]{ccc}u_{n,k_{z}}(x,y,t)\\ v_{n,k_{z}}(x,y,t)\end{array}\right) can be done with $O(N)$ operations, where $N$ is the number of grid points in the $x-y$ plane. Similarly, because $H_{\rm{KE}}$ is diagonal in $k$-space, the operator $\exp(-iH_{\rm{KE}}\delta t)$ can be applied to the Fourier transform of the $u$’s and $v$’s in $O(N)$ operations. Thus, our split-step algorithm can is represented formally as

[TABLE]

where $\mathcal{F}$ represents the application of a Fast Fourier Transform which takes $O(N\log N)$. Note that in each time step we reevaluate $\Delta(x,y,t)$ using Eq. (10) with the updated $u$’s and $v$’s.