gpe.bec

Code for simulating a single component BEC.

This is the reference implementation of the GPE for a single component BEC. The current parameters are designed for easily working with the Rb87 experiments in Peter Engels lab at WSU, but the framework is quite general and should be easy to adapt to other experiments.

Note: for performance, the classes here use the AsNumpyMixin class which turns “sunder” methods - those starting and ending with an underscore - into proper public methods using either numpy or cupy arrays.

Classes

Units	Units and physical constants.
GPEMixin	Mixing providing the GPE equation of state.
StateBase	State for the GPE in a homogeneous box.
HOMixin	Helper mixin class for harmonically trapped systems.
StateTwist_x	Minimal modification of the State class that implements twisted boundary
StateScaleBase	This state implements the scaling from Castin and Dum.

Module Contents

class Units

Units and physical constants.

This class is simply a container for physical constants and units. We have chosen values here that are relevant for 87Rb which is studied at WSU in Peter Engels’ lab.

hbar = 1.0

micron = 1.0

mm = 1000.0

cm = 10000.0

nm = 0.001

meter = 1000000.0

u = 1.0

kg = 6.0221408585491615e+26

G = 1.0

mG = 0.001

m_K39 = 38.96370648

m_Li6 = 6.0151228874

a_B = 5.291772105498088e-05

a_KK = 0.010583544210996176

Hz = 1.5746097513521148e-05

kHz = 0.015746097513521146

MHz = 15.746097513521146

s = 63507.79925891489

ms = 63.50779925891489

us = 0.0635077992589149

ns = 6.35077992589149e-05

c = 4720561277.486194

mW

microK = 2.061483743918704

nK = 0.002061483743918704

scattering_lengths

masses

a

mu_B

static get_magnetic_moment_mu_B(species=(1, -1))

Return the magnetic moment for the specified species in units of mu_B.

See also

https

//steck.us/alkalidata/

https

//link.aps.org/doi/10.1103/RevModPhys.49.31

Parameters:

species ((int, int)) – (F, mF) specifying the hypdrfine states of the 87Rb ground state (5^2 S_{1/2}).

Examples

>>> u.get_magnetic_moment_mu_B((1, -1)), u.magnetic_moment[(1, -1)]/u.mu_B
(0.5018..., 0.5018...)
>>> u.get_magnetic_moment_mu_B((1, 0))
-0.0
>>> u.get_magnetic_moment_mu_B((2, -2))
-0.99967...
>>> u.get_magnetic_moment_mu_B((2, 0))
0.0
>>> u.get_magnetic_moment_mu_B((2, 2))
0.99967...

class GPEMixin(g=None, **kw)

Bases: gpe.utils.GPUHelper

Mixing providing the GPE equation of state.

g

init()

get_n_TF(V_TF, V_ext=None, g=None, state=None)

Return the Thomas Fermi density profile n from V_TF.

Parameters:

V_TF (float) – Value of V(x_TF) where the density should vanish in the TF limit.

get_Vint_GPU(state=None)

Return the “internal” mean-field potential.

This version implements the standard GPE where the energy-density has \(gn^2/2\), so we have the derivative gn here.

get_Eint(state=None)

Return the “internal” mean-field energy-density.

This version implements the standard GPE where the energy-density has \(gn^2/2\). The method get_Vint() should return the appropriate derivative of this.

class StateBase(g=None, **kw)

Bases: GPEMixin, _StateBase

State for the GPE in a homogeneous box.

class HOMixin(ws=None, **kw)

Bases: gpe.utils.GPUHelper

Helper mixin class for harmonically trapped systems.

It is common in experiments for the clouds to be trapped in an external harmonic oscillator potential. Often, components of this trapping potential are turned off or relaxed to allow for expansion before imaging. This class provides an implementation of get_Vext_GPU() that implements the trapping potential. The user should call this (with super().get_Vext_GPU()) and then add any additional potentials.

This provides a default get_Vext_GPU().

g

m = 86.909187

ws

memoize = True

get_ws(t=None)

Return the trapping frequencies at time t.

init()

get_Vext_GPU(state=None)

Return the external potential.

Overload this method if you want to change the external potential. If the potential should be time dependent, use self.t which will be updated by the evolvers.

If a chemical potential self.mu is defined, then this is subtracted from Vext. This allows gpe.minimize to find states at a constant chemical potential. Note: for general evolution, it is better not to set the chemical potential as this is automatically set by compute_dy_dt and will then cause get_energy to return the thermodynamic potential instead.

get_mu_HO(N)

Return the chemical potential required to set the particle number in the Thomas Fermi (TF) approximation for a 3D harmonic oscillator.

class StateTwist_x(twist=None, v_x=0, kwz2=None, **kw)

Bases: StateBase

Minimal modification of the State class that implements twisted boundary conditions along the x direction and a boost with velocity v_x.

Parameters:

• twist (float) – Twisted boundary conditions. This is the twist angle from left to right over the length of the box. If there is a twist, then the wavefunction stored in the state is the periodic version with the twist removed. For this reason, we include the functions get_psi() and set_psi() which allow the user to access the physical non-periodic wavefunction.

• v_x (float) – Boost velocity along the x axis.

• kwz2 (float) – Angular velocity of about z-axis in expressed as kwz2 = m*omega_z/hbar.

twist

v_x = 0

kwz2 = None

property Lx

property k_B

Bloch wave-vector. This can be redefined by subclasses if needed.

property kx

Momentum along x direction.

Twisted boundary conditions are implemented as a shift in the momenta by the “Bloch” momentum

property twist_phase_x

Return the current twist phase.

get_psi_GPU()

Return the physical wavefunction (applying any twist).

set_psi(psi)

Set the state from a physical wavefunction (removing any twist).

property x_lab

Return the abscissa in the non-moving (lab) frame for comparison.

property x_v

Return the abscissa in the moving frame. This is just a shortcut for the abscissa since the problem is formulated in this frame.

property kx2

Return the equivalent of \(k_x^2\) which enters the dispersion.

For normal dispersion, we have:

K = (hbar*k)**2/2/m = - K_factor * kx2

Hence, for a modified dispersion E(kx2) we should have:

K = E(kx2) = - K_factor * kx2 kx2 = -E(kx2) / K_factor

class StateScaleBase(g=None, **kw)

Bases: StateBase

This state implements the scaling from Castin and Dum.

To use this class, you must provide get_lambdas() which defines the scaling factors at the current time.

get_lambdas(t=None)

Return (lams, dlams, ddlams): the scale factors and derivatives.

These should be computed at time t which should be self.t if t is None. There should be exactly self.dim scale factors - one for each dimension. If a dimension is not to be scaled, then the scale factor should be kept constant at 1.

_get_scale_factor_and_phase_GPU()

get_psi_GPU()

Return the physical wavefunction. Use this so classes can overload.

set_psi(psi)

Set the state from a physical wavefunction.

get_xyz_GPU()

Return the scaled physical coordinates (XYZ in notes).

property metric

apply_laplacian(factor, exp=False, **_kw)

Apply the laplacian multiplied by factor to the state.

_get_Vcorr_GPU()

Return the correction to the potential induced by coordinate transform.

get_V_GPU()

Return the complete potential V - internal and external.

This version includes the correction from the scaled coordinates.