The Structure of Polarized Vortices in the Unitary Fermi Gas

Fermionic superfluids do not generally support polarization, and the nature of the ground state of a slightly polarized unitary Fermi gas remains an open question. However, vortices naturally support polarization since the pairing gap vanishes in the core of superfluid vortices. The structure of a polarized vortex is not well understood and may have some interesting properties. To study the microscopic structure of a vortex, we use a density functional theory called the asymmetric superfluid local density approximation (ASLDA) to simulate how vortexes interact, evolve and how energy transfers between paired and unpaired particles. In this talk, I will discuss the structure and properties of polarized vortices using the ASLDA, and how these related to polarized phases through the Thomas-Fermi (TF) approximation.

Some Note from Dr. Forbes

We start with the issue of regularization. In our first version we use an extremely simple approach: we simply use all states in the box, and adjust the fixed coupling constant $g_{c}$ so as to reproduce a particular value of the gap $Δ$ . As a baseline, we use the following known 2D solution.

n_{+} = \frac{k_{F}^{2}}{2 π}, ϵ_{F} = \frac{ℏ^{2} k_{F}^{2}}{2 m}, E_{F G} = \frac{k_{F}^{4}}{8 m π} = \frac{1}{2} n_{+} ϵ_{F}, \tilde{C} = 0, \frac{μ}{ϵ_{F}} = \frac{1}{2}, \frac{Δ}{ϵ_{F}} = \sqrt{2}

The $T = 0$ equations in $d$ dimensions can be expressed in terms of the gap $Δ$ and effective chemical potential $μ_{eff} = μ_{0} - g_{c} n_{+} / 2$ where $g_{c} < 0$ represents an attractive interaction:

Δ = g_{c} ν_{c}, ϵ_{+} = \frac{ℏ^{2} k^{2}}{2 m} - μ_{eff}, E = \frac{E}{V} = \frac{ℏ^{2} τ_{c}}{2 m} + g_{c} (n_{a} n_{b} + ν_{c}^{†} ν_{c}) = \frac{ℏ^{2} τ_{c}}{2 m} + Δ^{†} ν_{c} - \frac{(μ_{eff} - μ_{0}) n_{+}}{2},

through the following integrals:

ν_{c} = - \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{Δ}{2 \sqrt{ϵ_{+}^{2} + Δ^{2}}} \sim - \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{m Δ}{ℏ^{2} k^{2}},

τ_{c} = \int_{c} \frac{d^{d} k}{(2 π)^{d}} k^{2} [1 - \frac{ϵ_{+}}{\sqrt{ϵ_{+}^{2} + Δ^{2}}}] \sim \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{2 m^{2} Δ^{2}}{ℏ^{4} k^{2}},

n_{+} = \int_{c} \frac{d^{d} k}{(2 π)^{d}} [1 - \frac{ϵ_{+}}{\sqrt{ϵ_{+}^{2} + Δ^{2}}}] \sim \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{2 m^{2} Δ^{2}}{ℏ^{4} k^{4}} .

In 1D, these converge, but in 2D and higher, the quantities $τ_{c}$ and $ν_{c}$ diverge as $\ln (k_{c})$ in 2D and as $k_{c}$ in 3D. The density $n_{+}$ converges in both 2D and 3D, although slower in 3D. Note that the divergence in $τ_{c}$ and $ν_{c}$ cancel in the combination

$\begin{aligned} (1) & κ & = \int_{c} \frac{d^{d} k}{(2 π)^{d}} (\frac{ℏ^{2} τ_{c}}{2 m} + Δ^{†} ν_{c}) & \sim \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{2 m^{2} Δ^{2} μ_{eff}}{ℏ^{4} k^{4}}, \end{aligned}$

which enters the energy density $E$ .

To regularize the theory, we must define some finite physical quantities that we can hold fixed as we take the cutoff to infinity so as to defined the parameters of the theory in a way that the physics becomes independent of the cutoff in the large cutoff limit. In 3D we express this as follows (see Bulgac:2011 Eq. (84) but we consider only unit effective mass $α_{+} = 1$ ):

$\begin{aligned} (2) & \tilde{C} (n_{a}, n_{b}) & = \overset{\frac{1}{g_{c}}}{\overset{⏞}{\frac{ν_{c}}{Δ}}} + \overset{Λ_{c}}{\overset{⏞}{\frac{1}{2} \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{1}{ϵ_{+} + i 0^{+}}}} = \frac{1}{g_{c}} + Λ_{c} & \sim \int_{c} \frac{d^{d} k}{(2 π)^{d}} \frac{2 m^{3} Δ^{2}}{ℏ^{6} k^{6}} . \end{aligned}$

where the integral is take in the principal value sense. This is motivated by the familiar Eagles-Leggett model, where this quantity is related to the two-body scattering length $a_{s}$ :

\tilde{C} (n_{a}, n_{b}) = \frac{m}{4 π ℏ^{2} a_{s}} .

Strictly speaking, the denominator in the integral should have $ϵ_{+} = ℏ^{2} k^{2} / 2 m$ without the chemical potential because the calculation of the scattering length should take place in the vacuum. To improve convergence, however, we include the chemical potential in $ϵ_{+} = ℏ^{2} k^{2} / 2 m - μ$ as advocated in BY:2002fk. This does not change the 3D result in the limit $k_{c} \to \infty$ but shifts the finite portion to improve convergence. Without this shift, the integral would behave as $\tilde{C} \sim - \int_{c} d^{d} k 2 m^{2} μ_{eff} / (2 π)^{d} ℏ^{4} k^{4}$ .

The principal value integrals can now be performed analytically:

Λ_{c}^{1 D} = \frac{m}{ℏ^{2}} \frac{1}{2 π k_{0}} \ln \frac{k_{c} - k_{0}}{k_{c} + k_{0}} \to - \frac{m}{ℏ^{2}} \frac{1}{π k_{c}} \to 0,

Λ_{c}^{2 D} = \frac{m}{ℏ^{2}} \frac{1}{4 π} \ln (\frac{k_{c}^{2}}{k_{0}^{2}} - 1) \to \frac{m}{ℏ^{2}} \frac{1}{2 π} \ln \frac{k_{c}}{k_{0}},

Λ_{c}^{3 D} = \frac{m}{ℏ^{2}} \frac{k_{c}}{2 π^{2}} (1 - \frac{k_{0}}{2 k_{c}} \ln \frac{k_{c} + k_{0}}{k_{c} - k_{0}}) \to \frac{m}{ℏ^{2}} \frac{k_{c}}{2 π^{2}},

where

\frac{ℏ^{2} k_{0}^{2}}{2 m} - μ = 0, \frac{ℏ^{2} k_{c}^{2}}{2 m} - μ = E_{c} .

The latter equation allows us to replace the momentum cutoff $k_{c}$ in terms of an energy cutoff $E_{c}$ which is of use in inhomogeneous systems.

Comments:

As we take the cutoff $k_{c} \to \infty$ we notice the following:

In 1D, $Λ_{c} \to 0$ :
- This indicates that $\tilde{C} = - ν_{c} / Δ = - 1 / g_{c}$ is finite, and the theory is well defined in this limit with finite anomalous density $ν$ , and finite cutoff $g = v_{0}$ . In this case, no cutoff is needed, but one might still like to use the cutoff procedure to improve the convergence. In particular, when working on a lattice, there is a maximum momentum that can be represented $k_{max}$ which will limit the convergence roughly according to the behavior discussed above. For example, $ν \sim 1 / k_{max}$ will converge slowly, requiring large lattices to approach the physical limit. Using the cutoff procedure and Thomas-Fermi completion discussed below will speed this convergence.
- Since $g$ is finite, one must also include the Hartree corrections to the chemical potentials:
  $μ_{a, b}^{eff} = μ_{a, b} + g n_{b, a} .$
  This complicates things a bit because now the effective potentials depend on the density. One can still work with fixed $μ$ , but this means that the self-consistent solution is the solution in a funky background potential which has been chosen to absorb these corrections.
In 2D and 3D, $Λ_{c} \to \infty$ divergence. This means that if we fix $\tilde{C}$ to be finite, then $g_{c} \to 0$ . In this case, the Haretree corrections vanish and the effective chemical potentials do not depend on the state:
$μ^{eff} = μ = μ_{0} - V (x) .$
This allows us to easily work in the grand-canonical ensemble.
In 3D, $Λ_{c}$ is independent of $k_{0}$ in the limit of large cutoff. This means that the including the pole, which improves the convergence, does not affect the usual identification of the scattering length.
In 2D, however, there remains a dependence on $k_{0}$ in the large cutoff limit. This means that our procedure for specifying the functional through $\tilde{C}$ depends on $k_{0}$ . This is not ideal for a density function because it means that the physical interpretation of the results depend on the external potential - i.e. the DFT is not independent of the potential as it should be. The previous point that Hartree terms vanish helps a bit since $k_{0}$ then depends only on the external trap and bare chemical potentials. (Otherwise, the DFT would be state-dependent which would be even worse!)

The results of this simple approach must be interpreted instead as the results of a DFT fixing the following combination:
$\tilde{C} - \frac{m}{2 π ℏ^{2}} \ln k_{0}$
which becomes independent of $k_{c}$ in the limit of large cutoff. Note that this diverges at the classical turning points where $k_{0} = 0$ . I am not sure if this poses any issues, and perhaps is expected and related to some of the divergent behaviour in 2D where results depend on the boundary conditions quite strongly. (One should include a dimensionless factor in the log to render $k_{0} / \tilde{k}$ dimensionless, but this will just be absorbed into the constant value.) In any case, this simple regularization procedure should be fine for qualitative and exploratory work keeping in mind that the results might depart strongly from “physical” 2D results near the turning points.

Example: 2D

As an example, consider homogeneous states in 2D. Suppose one knows $Δ$ and $μ_{+}$ , how should $\tilde{C}$ be determined? In principle, one can simply compute:

\tilde{C} (n_{a}, n_{b}) = lim_{k_{c} \to \infty} (\frac{ν_{c}}{Δ} + \frac{1}{2} \int_{0}^{k_{c}} \frac{d k}{2 π} \frac{k}{ϵ_{+} + i 0^{+}}) = \frac{1}{2} \int_{0}^{\infty} \frac{d k}{2 π} (- \frac{k}{\sqrt{ϵ_{+}^{2} + Δ^{2}}} + \frac{k}{ϵ_{+} + i 0^{+}}) .

However, numerically computing this integral can be challenging because of the pole. We can proceed as we would were we performing the computation in a discrete basis.

Thomas Fermi Completion

To facilitate convergence in the case when the solutions vary slowly, we can include states with $E > E_{c}$ by assuming that on length-scales smaller than $1 / k_{c}$ the solution is homogeneous and performing the integrals. This can be done exactly (numerically completing the integrals) or as a series approximation in $1 / k_{c}$ .

As an example, consider the computation of the gap equation. We will illustrate this with symmetric homogeneous matter in 2D. As discussed above, this presents some difficulties due to the divergent behaviour of $Λ_{c}$ . In a typical calculation, this would proceed as follows:

First compute the anomalous density $ν_{c}$ by diagonalizing the Hamiltonian and include only states with $E < E_{c}$ = $ℏ^{2} k_{c}^{2} / 2 m$ . For 2D homogeneous matter this would be
$ν_{c} = - \int_{0}^{k_{c}} \frac{k d k}{2 π} \frac{Δ}{2 \sqrt{ϵ_{+}^{2} + Δ^{2}}} .$
We would then use the regulated equation for $\tilde{C}$ to compute the coupling constant:
$g_{c} = \frac{1}{\tilde{C} - Λ_{c}} .$
Finally, we would use the gap equation to find the self-consistency condition $Δ = g_{c} ν_{c}$ .

Putting these together, we have:

\tilde{C} (n_{a}, n_{b}) = - \frac{ν_{c}}{Δ} + Λ_{c} + O (k_{c}^{- 4}) = - \int_{0}^{k_{c}} \frac{k d k}{2 π} \frac{Δ}{2 \sqrt{ϵ_{+}^{2} + Δ^{2}}} + \frac{1}{2} \int_{0}^{k_{c}} \frac{k d k}{2 π} \frac{1}{ϵ_{+} + i 0^{+}} + O (k_{c}^{- 4}) .

This is missing the contributions from the states above $k_{c}$ which would be corrected by including the full integrals:

\tilde{C} (n_{a}, n_{b}) = \frac{ν_{c}}{Δ} + Λ_{c} + \frac{1}{2} \int_{k_{c}}^{\infty} \frac{k d k}{2 π} (\frac{1}{ϵ_{+} + i 0^{+}} - \frac{1}{\sqrt{ϵ_{+}^{2} + Δ^{2}}}) .

We call this latter correction the “Thomas-Fermi Completion”. The idea is to assume that on length scales shorter than $1 / k_{c}$ , the system is approximately homogeneous and so one can use the homogeneous equations to complete the contributions to various pieces from the missing states. One either computes the integrals above $k_{c}$ numerically, or performs a series expansion in $k_{c}^{- 1}$ to compute the corrections at each point in space. (These corrections will in general depend on position because the quantities in the integrand such as $Δ (x)$ will depend depend on position.)

Aside: this procedure also allows us to compute $\tilde{C}$ in two dimensions. Although the integral is finite, numerically computing through the pole is a challenge since one really needs a PV integral. If we follow this procedure, however, we can compute $ν_{c}$ numerically, $Λ_{c}$ analytically, then apply the TF completion.

$n_{+} = \int d k \overset{f_{+}}{\overset{⏞}{(1 - \frac{ϵ_{k}^{+}}{E_{k}} \underset{f_{ν}}{\underset{⏟}{(f (ω_{-}) - f (ω_{+}))}})}}$ $n_{-} = \int d k \overset{f_{-}}{\overset{⏞}{(f (ω_{+}) - f (- ω_{-}))}}$

τ_{+} = \int d k k^{2} f_{+}, τ_{-} = \int d k k^{2} f_{-}, ν = - \int d k \frac{Δ}{2 E_{k}} \overset{f_{ν}}{\overset{⏞}{(f (ω_{-}) - f (ω_{+}))}}

κ = \int d k (\frac{ℏ^{2} k^{2}}{2 m_{+}} f_{+} - \frac{Δ^{2}}{2 E_{k}} f_{ν}), \tilde{C} = \int d k \frac{1}{2} (\frac{1}{ϵ_{+} + i 0^{+}} - \frac{f_{ν}}{E_{k}}) .

n_{+}^{3 D} \approx \int_{E < E_{c}} \frac{d^{d} k}{(2 π)^{d}} [1 - \frac{ϵ_{+}}{\sqrt{ϵ_{+}^{2} + Δ^{2}}}] + \frac{m^{2} Δ^{2}}{k_{c} ℏ^{4} π^{2}} + \frac{4 m^{3} Δ^{2} μ}{3 ℏ^{6} π^{2} k_{c}^{3}} + O (k_{c}^{- 5}),

κ^{3 D} \approx \int_{E < E_{c}} \frac{d^{d} k}{(2 π)^{d}} (\frac{ℏ^{2} τ_{c}}{2 m} + Δ^{†} ν_{c}) + \frac{m^{2} μ Δ^{2}}{ℏ^{4} π^{2} k_{c}} + \frac{m^{3} Δ^{2} (8 μ^{2} - Δ^{2})}{6 ℏ^{6} π^{2} k_{c}^{3}} + O (k_{c}^{- 5}),

3D

Here we check the results for homogeneous matter in various dimensions. In 3D can compare these with the results we know for the BdG equation (Eagles-Leggett model) at unitarity:

n_{+} = n_{a} + n_{b} = \frac{k_{F}^{3}}{3 π^{2}}, ϵ_{F} = \frac{ℏ^{2} k_{F}^{2}}{2 m}, E_{F G} = \frac{k_{F}^{5}}{10 m π^{2}} = \frac{3}{5} n_{+} ϵ_{F}, \tilde{C} = 0,

ξ = \frac{E}{E_{F G}} = \frac{μ}{ϵ_{F}} = 0.59060550703283853378393810185221521748413488992993 \dots

\frac{Δ}{ϵ_{F}} = 0.68640205206984016444108204356564421137062514068346 \dots

$\frac{Δ}{μ} = 1.162200561790012570995259741628790656202543181557689 \dots .$

2D

In 2D, we must compute results but we seem to have:

n_{+} = \frac{k_{F}^{2}}{2 π}, ϵ_{F} = \frac{ℏ^{2} k_{F}^{2}}{2 m}, E_{F G} = \frac{k_{F}^{4}}{8 m π} = \frac{1}{2} n_{+} ϵ_{F}, \tilde{C} = 0, \frac{μ}{ϵ_{F}} = \frac{1}{2}, \frac{Δ}{ϵ_{F}} = \sqrt{2}

1D

In 1D, we have the following:

n_{+} = \frac{k_{F}}{π}, ϵ_{F} = \frac{ℏ^{2} k_{F}^{2}}{2 m}, E_{F G} = \frac{k_{F}^{3}}{3 m π} = \frac{2}{3} n_{+} ϵ_{F},

\tilde{C} = - \frac{m}{ℏ^{2} k_{F}}, \frac{μ}{ϵ_{F}} = 0.28223521359741266, \frac{Δ}{ϵ_{F}} = 0.41172622996179004 .

Some Note from Dr. Forbes

Example: 2D

Thomas Fermi Completion

3D

2D

1D