r"""Separable Potentials.
This module contains various three-dimensional separable potentials of the form
.. math::
\op{V} = g V_k V_q
The relevant one-body Schrodinger equation becomes in momentum space (for bound
states):
.. math::
\left(\frac{k}{2m} - E\right)\psi_k
= -gV_k \int\dbar^{3}{\vect{q}}\;V_q\psi_q.
The effective range expansion can be expressed succinctly for separable
potentials:
.. math::
k\cot\delta_k = -4\pi\left(\frac{1}{2 m g V_k^2} +
\sum_{p}\frac{\frac{V_p^2}{V_k^2} - \frac{k^2}{p^2}}
{p^2 - k^2}\right)
= \frac{-1}{a} + \frac{r_{\text{eff}}k^2}{2}
+ P r_{\text{eff}}^{3}k^{4} + \order(k^6).
The first two coefficients are computed by the methods
:meth:`Separable.ainv` and :meth:`Separable.reff`.
.. note:: Everything here is expressed in units where
:math:`\hbar=1`. To restore these constants make the following
substitutions. Also, these units are for the one-body problem: i.e. $m$ is
the mass of a single species. If you wish to use these for the two-body
problem, be sure to use the reduced mass in place of $m$.
.. math::
a, r_{\text{eff}}, k & \rightarrow a, r_{\text{eff}}, k\\
g &\rightarrow \frac{g}{\hbar^2},\\
\Lambda, p & \rightarrow \frac{\Lambda}{\hbar}, \frac{p}{\hbar}.
The dimensions are thus:
[$\hbar$] :
$ET = ML^2/T = 1$ which implies $T \equiv ML^2$.
[$v$] :
$1$
[$g$] :
$EV = ML^5/T^2 = 1/T = L/M$.
[$a$] = [$r_e$] :
$1/L$
"""
from __future__ import division
__all__ = ['VPoschlTeller', 'Separable', 'SharpCutoff', 'SmoothCutoff',
'Yamaguchi', 'EffectiveRange']
from numpy import pi, inf
import numpy as np
import scipy.integrate
sp = scipy
from mmf.objects import StateVars, process_vars
from mmf.objects import Required, Computed
import mmf.math.differentiate
import mmf.math.special
import mmf.math.integrate.integrate_1d.quadrature
from mmf.math.integrate.integrate_1d.quadrature import Decay, PowerDecay
from mmf.interfaces import Interface, implements, Attribute
_ABS_TOL = 1e-12
_REL_TOL = 1e-12
_EPS = np.finfo(float).eps
[docs]class VPoschlTeller(StateVars):
r"""Poschl-Teller potential used by GFMC groups."""
_state_vars = [
('V0', -1.0, "Strength of the potential."),
('r0', Required,\
"""Range of potential. (This is only the effective range
if V0 = -1)""")]
process_vars()
[docs] def V(self, r, m=1.0):
return 4.0*self.V0/m/self.r0**2/np.cosh(2.0*r/self.r0)**2
[docs] def Vp(self, p, m=1.0):
r"""Fourier transform of potential.
Notes
-----
.. math::
\begin{aligned}
\tilde{V}_{p} &= \frac{V_{0}\hbar^2\pi^3 r_{0}}{m}\left(
\frac{\coth(x) - x^{-1}}{\sinh(x)}
\right), &
x &= \frac{\pi p r_{0}}{4}
\end{aligned}
"""
x = pi*p*self.r0/4.0
x2 = x*x
x1 = np.where(x2>0.01,x,1.0)
fact = np.where(abs(x2)>0.01,
(1./np.tanh(x1) - 1./x1)/np.sinh(x1),
1./3. + x2*(-7./90. + x2*31./2520.))
return self.V0*pi**3*self.r0/m*fact
[docs] def Vpq(self, p, q, m=1.0):
r"""Integrated potential over angles that appears in the BCS
s-wave gap equation:
.. math::
\Delta_{q} = \int_{0}^{\infty}\frac{\d{p}}{2\pi^2}
\frac{V_{qp}\Delta_{p}}{2\sqrt{E_{p}^2 + \Delta_{p}^2}}
Notes
-----
.. math::
V_{pq} &= \frac{V_{0}\hbar^2\pi^3 r_{0}}{2 m p q}\left(
\frac{p - q}{\sinh(p - q)} - \frac{p + q}{\sinh(p + q)}
\right)
Examples
--------
>>> V = VPoschlTeller(r0=1)
>>> p = np.array([[0,1,2]])
>>> q = p.T
>>> V.Vpq(p,q)
array([[-6.57973627, -6.21059264, -5.25788171],
[-6.21059264, -5.88551639, -5.03455853],
[-5.25788171, -5.03455853, -4.4270986 ]])
"""
# These factors should be 1 according to text, but that does
# not work... need to figure out!
_a_ = 2.0 # Could be problem with reduced mass
_b_ = 2./pi
const = _a_*self.V0*pi**2*self.r0/2.0/m
p = np.asarray(p)
q = np.asarray(q)
x = _b_*pi*p*self.r0/4.0 + 0*q
y = _b_*pi*q*self.r0/4.0 + 0*p
res = np.where(x == 0,
np.where(y == 0,
2./3.,
2*(1./np.tanh(y) - 1./y)/np.sinh(y)),
np.where(y == 0,
2*(1./np.tanh(x) - 1./x)/np.sinh(x),
np.where(x == y,
(1./x - 2.0/np.sinh(2.0*x))/x,
((x - y)/np.sinh(x - y) -
(x + y)/np.sinh(x + y))/x/y)))
return const*res
[docs]class Separable(StateVars):
r"""Base class for separable potentials, including methods to
compute their scattering properties: :math:`a^{-1}` and
:math:`r_{eff}` are computed by :meth:`potentials.ainv` and
:meth:`potentials.reff`.
Notes
-----
These potentials have the form:
.. math::
V(p,q) = gv(p)v(q)
Our convention for normalization is that v(0) = 1. By default,
the coupling constant is set to -1, so the potential is
attractive. The form-factors :math:`v(p)` should be functions
that decay for large momenta and are roughly constant for small
momenta.
"""
_state_vars = [('g', -1.0, "Coupling constant"),
('m', 1.0, "Mass. Set to reduced mass for two-body"),
('inf', Computed,\
"""Infinite endpoint wrapped with a Decay
subclass to encode the nature of the decay.""")]
process_vars()
######################################################################
# Required methods for subclasses to implement
[docs] def v(self,p):
r"""Return the form factor."""
raise NotImplementedError
# Suggested method for subclasses to implement
def _get_p_c(self):
r"""Return a list of typical/critical momenta for integration.
These points are explicitly included in integrations to ensure
that important features of the potential are not missed."""
return []
def _get_inf(self):
r"""Return a form of `inf`, preferably wrapped with
:class:`Decay` to describe the asymptotic behaviour of
(v_p)^2.
"""
return inf
def _kappa(self):
r"""Return diff(v(p),p^2)/v(0) where p=0."""
def f(p2):
return self.v(np.sqrt(p2))
kappa = mmf.math.differentiate.differentiate(f,
dir=1.0)/self.v(0)
return kappa
[docs] def set_g(self, ainv, reff=None):
r"""Set the coupling and Lambda to fix the scattering length
and range."""
if reff is None:
def _f0(g):
self.g = g
return self.ainv() - ainv
self.g = sp.optimize.fsolve(_f0,self.g)
else:
def _f1(x):
self.g, self.Lambda = x
return [self.ainv() - ainv, self.reff() - reff]
self.g, self.Lambda = sp.optimize.fsolve(_f1,[self.g, self.Lambda])
######################################################################
[docs] def __init__(self, *v, **kw):
r"""Subclasses should call to add p_c to `self.V`"""
self.v.__dict__['p_c'] = self.p_c
self.inf = self._get_inf()
@property
[docs] def p_c(self):
r"""List of typical/critical momenta. The points should be
carefully inspected during integration etc."""
return self._get_p_c()
def _v2p2(self, p2):
r"""Return $v(p)^2$ as a function of $p^2$. This is used by
:meth:`k_cot_delta`."""
if p2 < 0:
p2 += 0j
return self.v(np.sqrt(p2))**2
def _int(self, f, a, b, abs_tol=_ABS_TOL, rel_tol=_REL_TOL):
r"""Integrate f including self.p_c as special points."""
res, err = mmf.math.integrate.integrate_1d.quadrature.integrate(
f, a, b, points=self.p_c,
abs_tol=abs_tol, rel_tol=rel_tol)
return res
[docs] def k_cot_delta(self, p2, abs_tol=_ABS_TOL, rel_tol=_REL_TOL):
r"""Return the phase shift $p\cot\delta_p$ as a function of $p^2$.
.. math::
p\cot\delta_p = -4\pi\left(\frac{1}{2 m g V_p^2} +
\sum_{k}\frac{\frac{V_k^2}{V_p^2} - \frac{p^2}{k^2}}
{k^2 - p^2}\right)
"""
v2p2 = self._v2p2(p2)
def f(k):
k2 = k*k
return k2*(self._v2p2(k2)/v2p2 - p2/k2)/(k2-p2)
mg = self.m*self.g
return -4.0*np.pi*(
0.5/mg/v2p2 + self._int(f, 0, self.inf,
abs_tol, rel_tol)/(2*np.pi**2))
[docs] def ainv(self, abs_tol=_ABS_TOL, rel_tol=_REL_TOL):
r"""Return the inverse s-wave scattering-length.
Notes
-----
This is computed by
.. math::
a^{-1} = \frac{2\pi}{gmV_{0}^2}
+\frac{2}{\pi V_{0}^2}\int_{0}^{\infty}\d{k}\;V_{k}^2
"""
v0 = self.v(0)
def f(k):
return (self.v(k)/v0)**2
mg = self.m*self.g
return 2*pi/mg/v0/v0 + 2/pi*self._int(f, 0, self.inf,
abs_tol, rel_tol)
[docs] def reff(self, abs_tol=_ABS_TOL, rel_tol=_REL_TOL):
r"""Return the s-wave effective range.
Notes
-----
This is computed by:
.. math::
r_{\text{eff}} = \frac{8\pi \kappa}{gmV_{0}^2}
- \frac{8}{\pi}\int_{0}^{\infty}\d{k}
\frac{\frac{V_{k}^{2}}{V_{0}^{2}}\left(1-2k^2\kappa\right)
-1
}{k^2}
where :math:`\kappa = \left.\diff{\ln V_{p}}{p^2}\right|_{p=0}`
"""
v0 = self.v(0)
kappa = self._kappa()
if isinstance(self.inf, Decay):
if 0 == kappa:
inf_ = PowerDecay(inf,
x0=self.inf.x0,
n=2)
else:
inf_ = PowerDecay(inf,
x0=max(self.inf.x0,
1./np.sqrt(abs(kappa))),
n=2)
else:
inf_ = PowerDecay(inf,
x0=max(self.p_c),
n=2)
def f(k):
return ((self.v(k)/v0)**2*(1.0-2*k*k*kappa)-1.0)/k/k
mg = self.g*self.m
return 8*pi*kappa/mg/v0/v0 - 4/pi*self._int(f, 0, inf_,
abs_tol, rel_tol)
def _test(self):
r"""Automatic test for consistency of the various calculations."""
ai, re = self.ainv(), self.reff()
h = 1e-6
ainvs = [-1.0,0.0,1.0]
reffs = [0.1,0.01]
for _ai in ainvs:
for _re in reffs:
self.set_g(ainv=_ai, reff=_re)
assert(np.allclose(self.ainv(), -self.k_cot_delta(0)))
assert(np.allclose(self.ainv(), _ai))
reff = 2*(self.k_cot_delta(h) - self.k_cot_delta(-h))/(2*h)
assert(np.allclose(self.reff(), reff, rtol=1e-7, atol=1e-5))
assert(np.allclose(self.reff(), _re))
self.set_g(ainv=ai, reff=re)
[docs]class SharpCutoff(Separable):
r"""Separable potential with a sharp momentum cutoff.
Notes
-----
.. math::
v(p) &= \theta(p - \Lambda),\\
\frac{1}{a} &= \frac{2\pi}{mg} + \frac{2\Lambda}{\pi},\\
r_{\text{eff}} &= \frac{4}{\pi\Lambda},\\
Pr_{\text{eff}}^{3} &= \frac{2}{3\pi\Lambda^3}.
::
assume(L>0);assume(k>0);assume(p>0);additionally(L>k):
kcd:=(v)->-4*Pi*(1/2/m/g/v(k)**2 + int((v(p)**2/v(k)**2
- k^2/p^2)/(p^2-k^2)*p^2/2/Pi^2,p=0..infinity)):
v:=p->Heaviside(L-p):
ainv=-subs(k=0,kcd(v));
re=2*coeff(series(kcd(v),k,6),k,2);
P*re^3=coeff(series(kcd(v),k,6),k,4);
Examples
--------
>>> V = SharpCutoff(Lambda=10.0)
>>> V.set_g(ainv=1.0, reff=0.1)
>>> V.ainv() # doctest: +ELLIPSIS
1.00000000...
>>> V.reff() # doctest: +ELLIPSIS
0.10000000...
"""
_state_vars = [('Lambda', 1.0, "Momentum cutoff")]
process_vars()
[docs] def v(self,p):
return mmf.math.special.step(self.Lambda - p)
def _v2p2(self,p2):
return mmf.math.special.step(self.Lambda**2 - p2)
def _get_p_c(self):
return [self.Lambda]
def _get_inf(self):
return inf
def _kappa(self):
r"""Return diff(v(p),p^2)/v(0) where p=0."""
return 0.0
[docs] def set_g(self, ainv, reff=None):
r"""Set the coupling to fix the scattering length.
Notes
-----
.. math::
\Lambda &= \frac{4}{r_{\text{eff}}\pi},\\
g &= m^{-1}\frac{2\pi}{a^{-1} - \frac{2\Lambda}{\pi}}.
"""
if reff is not None:
self.Lambda = 4./reff/pi
self.g = 2.0*pi/(ainv - 2.0*self.Lambda/pi)/self.m
[docs]class SmoothCutoff(Separable):
r"""Separable potential with a smooth momentum cutoff.
.. math::
v(p) &= e^{-\tfrac{\pi}{8}\frac{p^2}{\Lambda^2}},\\
\frac{1}{a} &= \frac{2\pi}{mg} + \frac{2\Lambda}{\pi},\\
r_{\text{eff}} &= \frac{1}{\Lambda} - \frac{\pi^2}{mg\Lambda^2},\\
Pr_{\text{eff}}^{3} &= \frac{\pi}{48\Lambda^3}
- \frac{2\pi^3}{8g\Lambda^4}.
::
assume(L>0);assume(k>0);assume(p>0);additionally(L>k):
kcd:=(v)->-4*Pi*(1/2/m/g/v(k)**2 + int((v(p)**2/v(k)**2
- k^2/p^2)/(p^2-k^2)*p^2/2/Pi^2,p=0..infinity)):
v:=p->exp(-Pi/8*p^2/L^2):
eqs:={ainv=-expand(simplify(eval(subs(k=0,kcd(v))))),
re=2*expand(simplify(coeff(series(kcd(v),k,6),k,2)))};
P*re^3=coeff(series(kcd(v),k,6),k,4);
solve(eqs,[g,L]);
Examples
--------
>>> a = 1.0
>>> V = SmoothCutoff(Lambda=10.0)
>>> V.set_g(ainv=1.0, reff=0.1)
>>> V.ainv() # doctest: +ELLIPSIS
1.000000000...
>>> V.reff() # doctest: +ELLIPSIS
0.100000000...
"""
_state_vars = [('Lambda', 1.0, "Momentum cutoff")]
process_vars()
[docs] def v(self, p):
return np.exp(-pi/8*(p/self.Lambda)**2)
def _v2p2(self, p2):
return np.exp(-pi/4*p2/self.Lambda**2)
def _get_p_c(self):
# Add the point where V.v(p) < eps
return [self.Lambda, np.sqrt(abs(2*np.log(_EPS)))*self.Lambda]
def _get_inf(self):
return inf
def _kappa(self):
r"""Return diff(v(p),p^2)/v(0) where p=0."""
return -pi/8/self.Lambda**2
[docs] def set_g(self, ainv, reff=None):
r"""Set the coupling to fix the scattering length.
.. math::
\Lambda &= \frac{1 + \sqrt{1-\frac{\pi}{2}
\frac{r_{\text{eff}}}{a}}}
{r_{\text{eff}}},\\
g &= m^{-1}\frac{2\pi}{a^{-1} - \frac{2\Lambda}{\sqrt{\pi}}}.
"""
if reff:
d = 1.0 - pi*reff*ainv/2.0
if d < 0:
raise ValueError("reff must be less than 2a/pi")
self.Lambda = (1.0 + np.sqrt(d))/reff
self.g = 2*pi/(ainv - 2*self.Lambda/pi)/self.m
[docs]class EffectiveRange(Separable):
r"""Separable Effective Range Interaction, completely
determined by the scattering length and effective range.
.. math::
v(p) &= \frac{1}{\sqrt{\left(\frac{\pi p}{2\Lambda}\right)^2 + 1}},\\
\frac{1}{a} &= \frac{2\pi}{mg} + \frac{2\Lambda}{\pi},\\
r_{\text{eff}} &= -\frac{\pi^3}{mg\Lambda^2},\\
Pr_{\text{eff}}^{3} &= 0
::
assume(L>0);assume(k>0);assume(p>0);additionally(L>k):
kcd:=(v)->-4*Pi*(1/2/m/g/v(k)**2 + int((v(p)**2/v(k)**2
- k^2/p^2)/(p^2-k^2)*p^2/2/Pi^2,p=0..infinity)):
v:=p->1/sqrt((Pi*p/2/L)^2+1):
eqs:={ainv=-expand(simplify(eval(subs(k=0,kcd(v))))),
re=2*expand(simplify(coeff(series(kcd(v),k,6),k,2)))};
P*re^3=coeff(series(kcd(v),k,6),k,4);
simplify(expand(rationalize(subs(allvalues(solve(eqs,[g,L]))[1][1],L)),symbolic));
solve(eqs,[g,L]);
Note that this as a closed effective range expansion stopping at
the second term.
Examples
--------
>>> V = EffectiveRange(Lambda=10.0)
>>> V.set_g(ainv=1.0,reff=0.5)
>>> V.ainv() # doctest: +ELLIPSIS
1.000000...
>>> V.reff() # doctest: +ELLIPSIS
0.500000...
"""
_state_vars = [('Lambda', 1.0, "Momentum cutoff")]
process_vars()
[docs] def v(self, p):
return 1./np.sqrt((pi*p/2/self.Lambda)**2 + 1.0)
def _v2p2(self, p2):
return 1./((pi/2/self.Lambda)**2*p2 + 1.0)
def _get_p_c(self):
# Add the point where V.v(p) < eps
return [2*self.Lambda/pi]
# 2*self.Lambda/pi*np.sqrt(abs(1./_EPS**2 - 1.0))]
def _get_inf(self):
return PowerDecay(inf, x0=2*self.Lambda/pi, n=2)
def _kappa(self):
r"""Return diff(v(p),p^2)/v(0) where p=0."""
return -(pi/self.Lambda)**2/8
[docs] def set_g(self, ainv, reff=None):
r"""Set the coupling and Lambda to fix the scattering length
and range.
.. math::
\Lambda &= \frac{1+\sqrt{1 - \frac{2r_{\text{eff}}}{a}}}
{2r_{0}/\pi},\\
g &= \frac{4\pi}{a^{-1} - \frac{\Lambda}{\sqrt{\pi}}}.
"""
if reff:
r_a = reff*ainv
d = (1.0 - 2.0*r_a)
if d < 0:
raise ValueError("reff must be less than a/2")
self.Lambda = (1.0 + np.sqrt(d))/(2*reff/pi)
self.g = 2*pi/(ainv - 2*self.Lambda/pi)/self.m
[docs]class Yamaguchi(Separable):
r"""Separable Yamaguchi potential.
.. math::
v(p) &= \frac{1}{\left(\frac{\pi p}{4\Lambda}\right)^2 + 1},\\
\frac{1}{a} &= \frac{2\pi}{mg} + \frac{2\Lambda}{\pi},\\
r_{\text{eff}} &= \frac{\pi}{4\Lambda} - \frac{\pi^3}{2mg\Lambda^2}\\
Pr_{\text{eff}}^{3} &= -\frac{\pi^5}{128 mg \Lambda^4}
::
assume(L>0);assume(k>0);assume(p>0);additionally(L>k):
kcd:=(v)->-4*Pi*(1/2/m/g/v(k)**2 + int((v(p)**2/v(k)**2
- k^2/p^2)/(p^2-k^2)*p^2/2/Pi^2,p=0..infinity)):
v:=p->1/(1+(Pi*p/4/L)^2):
eqs:={ainv=-expand(simplify(eval(subs(k=0,kcd(v))))),
re=2*expand(simplify(coeff(series(kcd(v),k,6),k,2)))};
P*re^3=coeff(series(kcd(v),k,6),k,4);
solve(eqs,[g,L]);
simplify(expand(rationalize(allvalues(solve(eqs,[g,L]))[1][1])));
Note: this has a closed effective range expansion with three
terms.
Examples
--------
>>> V = Yamaguchi(Lambda=10.0)
>>> V.set_g(ainv=1.0,reff=0.5)
>>> V.ainv() # doctest: +ELLIPSIS
1.0...
>>> V.reff() # doctest: +ELLIPSIS
0.499999...
"""
_state_vars = [('Lambda', 1.0, "Momentum cutoff")]
process_vars()
[docs] def v(self, p):
return 1./((pi*p/4.0/self.Lambda)**2 + 1.0)
def _v2p2(self, p2):
return 1./((pi/4.0/self.Lambda)**2*p2 + 1.0)**2
def _get_inf(self):
return PowerDecay(inf, x0=self.Lambda, n=4)
def _get_p_c(self):
# Add the point where V.v(p) < eps
return [self.Lambda]
# np.sqrt(abs(1.0/_EPS - 1.0))*self.Lambda]
def _kappa(self):
r"""Return diff(v(p),p^2)/v(0) where p=0."""
return -pi*pi/self.Lambda**2/16
[docs] def set_g(self, ainv, reff=None):
r"""Set the coupling and Lambda to fix the scattering length
and range.
.. math::
\Lambda = \frac{\pi}{8r_{\text{eff}}}\left(
3 + \sqrt{9-16\frac{r_{\text{eff}}}{a}}\right).
"""
if reff:
d = (9.0 - 16.0*reff*ainv)
if d < 0:
raise ValueError("reff must be less than 9*a/16")
self.Lambda = pi*(3 + np.sqrt(d))/8/reff
self.g = 2*pi/(ainv - 2*self.Lambda/pi)/self.m