BEClab

A framework for numerical simulations of trapped BEC behavior.

BEC-specific wrappers

class beclab.System(components, interactions, losses=None, potential=None, linear_terms=None)

Main descriptor of a BEC system with the Hamiltonian

\[\hat{H} = \int d \mathbf{x} \sum_{j=1}^C \sum_{k=1}^C \left\{ \hat{\Psi}_j^{\dagger} \left( -\frac{\hbar^2}{2m_j} \nabla^2 + V_j(\mathbf{x}) \right) \delta_{jk} \hat{\Psi}_k + \frac{U_{jk}}{2} \int d \mathbf{x}^\prime \hat{\Psi}_j^\dagger \hat{\Psi}_k^{\prime\dagger} \hat{\Psi}_j \hat{\Psi}_k^\prime + \hat{\Psi}_j^\dagger \sum_{n=1}^N f_j^{(n)}(\hat{\Psi}_1, \ldots, \hat{\Psi}_C, t) \right\},\]

(where \(f_j\) are linear functions in the \(\hat{\Psi}\) arguments) and the master equation

\[\frac{d \hat{\rho}}{dt} = - \frac{i}{\hbar} \left[ \hat{H}, \hat{\rho} \right] + \sum_{\mathbf{l} \in L} \kappa_{\mathbf{l}} \int d \mathbf{x} \mathcal{L}_{\mathbf{l}} \left[ \hat{\rho} \right].\]
Parameters:
  • components – a list of Component objects.
  • interactions – an array \(U_{jk}\) with shape (len(components), len(components)).
  • losses – a list of tuples (\(\kappa_{\mathbf{l}}\), (\(l_1\), ...)) for every loss source \(\mathbf{l}\). Note that these are not the experimental loss coefficients, which can be calculated as \(\gamma_{j,\mathbf{l}} \equiv 2 l_j \kappa_{\mathbf{l}}\).
  • potential – a beclab.bec.Potential object.
  • linear_terms – a list of modules each containing functions with suffixes 0...``components-1`` (corresponding to the subscript index) that will be passed psi_0,..., psi_C, t, and their return values added to the equation of the corresponding component.
beclab.box_for_tf(system, comp_num, N, pad=1.2)

Returns a tuple with the box sizes that can accommodate the Thomas-Fermi ground state of a BEC with N particles made of component comp_num from the System system.

Grids

class beclab.UniformGrid(shape, box)

Bases: beclab.grid.Grid

Rectangular uniform grid.

Parameters:
  • shape – a tuple with the grid shape.
  • box – a tuple with the simulation box size (same length as shape).
V

The total volume of the simulation box.

class beclab.grid.Grid

The base class for simulation grids.

box

A tuple with simulation box sizes.

dimensions

The number of dimensions (the length of shape or box tuples).

dV

The volume of a grid cell.

dxs

A list of grid steps for each dimension.

shape

A tuple with the grid shape.

size

The total number of grid cells.

xs

A list of dimensions numpy arrays containing grid cell coordinates for each dimension.

Potential

class beclab.bec.Potential

Abstract base class of an external potential.

get_array(grid, components)

Returns an array with the values of the potential for the given grid and components.

Parameters:
  • grid – a Grid object.
  • components – a list of Component objects.
get_module(dtype, grid, components)

Returns the module calculating the value of the potential. The module must contain functions with signatures

INLINE WITHIN_KERNEL ${r_ctype} ${prefix}${comp_num}(
    int idx0, ..., ${r_ctype} t)

where r_ctype is a realy type.

Parameters:
  • grid – a Grid object.
  • components – a list of Component objects.
class beclab.HarmonicPotential(freqs, displacements=None)

Bases: beclab.bec.Potential

A class representing a multidimensional harmonic potential:

\[V_i = \frac{m_i}{2} \sum_{d=1}^D (2 \pi f_d)^2 (x_d + l_{i,d})^2\]
Parameters:
  • freqs – a tuple of trap frequencies (in Hz).
  • displacements – a list of tuples (component, dimension, displacement), specifying the displacement \(l_{i,d}\) for the given component \(i\) and the dimension \(d\).

Cutoffs

class beclab.WavelengthCutoff(ksquared)

Bases: beclab.cutoff.Cutoff

Cutoff based on wave lengths in a plane wave basis (essentially, an energy cutoff). All the modes with \(k^2 > \mathrm{cutoff}\) are projected out.

classmethod for_energy(energy, component)

Returns the cutoff expressed in terms of energy for the given component.

classmethod padded(grid, pad=1)

Returns the cutoff calculated to ensure the ellipse of active modes in a rectangular grid is padded. pad = 1 means that the ellipse is touching the sides of the box, pad = 4 corresponds to the “safe” padding used to avoid aliasing problems during FFT transformations.

class beclab.cutoff.Cutoff

The base class for mode cutoffs.

get_mask(grid)

Returns a numpy array with the cutoff mask for the given grid. Elements of the array are 0 in places of the modes that are projected out, and 1 otherwise.

get_modes_number(grid)

Returns the total number of active modes for the given grid.

Ground states

class beclab.ThomasFermiGroundState(thr, dtype, grid, system, cutoff=None)

Thomas-Fermi ground state generator.

Parameters:
  • thr – a Reikna Thread.
  • dtype – the dtype of the generated wavefunction
  • grid – a Grid object.
  • system – a System object.
  • cutoff – a Cutoff object.
wfs_meta

A WavefunctionSetMetadata object representing the generated states.

__call__(Ns)

Gererate a ground state with populations Ns = [N1, N2, ...]. The length of the list must be equal to the number of components in the system given to the constructor. Returns a WavefunctionSet object.

class beclab.ImaginaryTimeGroundState(thr, dtype, grid, system, stepper_cls=<class 'reiknacontrib.integrator.rk46nl_stepper.RK46NLStepper'>, cutoff=None, verbose=True)

Ground state generator based on imaginary time propagation.

Parameters:
  • thr – a Reikna Thread.
  • dtype – the dtype of the generated wavefunction
  • grid – a Grid object.
  • system – a System object.
  • stepper_cls – one of the reiknacontrib.integrator.Stepper classes.
  • cutoff – a Cutoff object.
  • verbose – whether do display additional information about the integration process.
wfs_meta

A WavefunctionSetMetadata object representing the generated states.

__call__(Ns, E_diff=1e-09, E_conv=1e-09, sample_time=0.001, samplers=None, return_info=False)

Gererate a ground state with given populations. The propagation in imaginary time will continue until the difference in energy between two successive samples is larger than the given threshold.

Parameters:
  • Ns – a list of target populations for each component. The length of the list must be equal to the number of components in the system given to the constructor.
  • E_diff – a relative energy difference threshold.
  • E_conv – a convergence threshold for the propagation. Convergence is estimated as the relative difference between the values of energy after propagation with normal and double step for sample_time (that is, the same as in reiknacontrib.integrator.Integrator.adaptive_step).
  • sample_time – time between successive sampling of energy.
  • samplers – additional samplers (reiknacontrib.integrator.Sampler objects) to invoke during propagation.
  • return_info – whether to return additional information about the propagation (see the return section below).
Returns:

if return_info == False, returns a WavefunctionSet object. Otherwise returns a tuple (wfs, result, info), where wfs is a WavefunctionSet object, and the last two are the result data structure and reiknacontrib.integrator.IntegrationInfo object, same as the ones returned by reiknacontrib.integrator.Integrator.adaptive_step.

Integration

class beclab.Integrator(wfs_meta, system, stepper_cls=<class 'reiknacontrib.integrator.rk46nl_stepper.RK46NLStepper'>, cutoff=None, profile=False)

BEC integration class.

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • system – a System object.
  • stepper_cls – one of the reiknacontrib.integrator.Stepper classes. Passed to reiknacontrib.integrator.Integrator.fixed_step or reiknacontrib.integrator.Integrator.adaptive_step.
  • cutoff – a Cutoff object.
  • profile – whether to synchronize with GPU before sampling. Passed to reiknacontrib.integrator.Integrator.fixed_step or reiknacontrib.integrator.Integrator.adaptive_step.
adaptive_step(wfs, *args, **kwds)

Start integration with adaptive step for the beclab.WavefunctionSet wfs. Other parameters and return values are the same as in reiknacontrib.integrator.Integrator.adaptive_step (except that data is replaced by wfs).

fixed_step(wfs, *args, **kwds)

Start integration with fixed step for the beclab.WavefunctionSet wfs. Other parameters and return values are the same as in reiknacontrib.integrator.Integrator.fixed_step (except that data is replaced by wfs).

Wavefunctions

beclab.REPR_CLASSICAL = 'classical'

str(object=’‘) -> str str(bytes_or_buffer[, encoding[, errors]]) -> str

Create a new string object from the given object. If encoding or errors is specified, then the object must expose a data buffer that will be decoded using the given encoding and error handler. Otherwise, returns the result of object.__str__() (if defined) or repr(object). encoding defaults to sys.getdefaultencoding(). errors defaults to ‘strict’.

beclab.REPR_WIGNER = 'Wigner'

str(object=’‘) -> str str(bytes_or_buffer[, encoding[, errors]]) -> str

Create a new string object from the given object. If encoding or errors is specified, then the object must expose a data buffer that will be decoded using the given encoding and error handler. Otherwise, returns the result of object.__str__() (if defined) or repr(object). encoding defaults to sys.getdefaultencoding(). errors defaults to ‘strict’.

beclab.REPR_POSITIVE_P = 'positive-P'

str(object=’‘) -> str str(bytes_or_buffer[, encoding[, errors]]) -> str

Create a new string object from the given object. If encoding or errors is specified, then the object must expose a data buffer that will be decoded using the given encoding and error handler. Otherwise, returns the result of object.__str__() (if defined) or repr(object). encoding defaults to sys.getdefaultencoding(). errors defaults to ‘strict’.

class beclab.wavefunction.WavefunctionSetMetadata(thread, dtype, grid, components=1, trajectories=1, representation='classical', cutoff=None)

Metadata for a wavefunction container object.

components

The number of components in this wavefunction.

cutoff

A Cutoff object.

data

A Reikna Type object with the container array metadata.

dtype

A numpy.dtype object of the container data type.

grid

A Grid object.

representation

One of REPR_CLASSICAL, REPR_WIGNER, REPR_POSITIVE_P.

shape

A tuple with the full shape of the data container.

thread

A Reikna Thread object the wavefunction is connected to.

trajectories

The number of independent wavefunctions in this container.

class beclab.WavefunctionSet(thread, dtype, grid, components=1, trajectories=1, representation='classical', cutoff=None)

Bases: beclab.wavefunction.WavefunctionSetMetadata

Wavefunction container.

data

A Reikna Array object.

Beam splitter

class beclab.BeamSplitter(wfs_meta, comp1_num=0, comp2_num=1, starting_phase=0, f_detuning=0, f_rabi=0, seed=None)

Beam splitter for a multi-component BEC. The effect of an oscillator is approximated by the application of the matrix

\[\begin{split}\begin{pmatrix} \Psi^\prime_1 \\ \Psi^\prime_2 \end{pmatrix} = \begin{pmatrix} \cos \frac{\theta}{2} & -i e^{-i \phi} \sin \frac{\theta}{2} \\ -i e^{i \phi} \sin \frac{\theta}{2} & \cos \frac{\theta}{2} \end{pmatrix} \begin{pmatrix} \Psi_1 \\ \Psi_2 \end{pmatrix},\end{split}\]

where \(\phi = \delta t + \alpha\).

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • comp1_num – the number of the first affected component.
  • comp2_num – the number of the second affected component.
  • starting_phase – initial phase \(\alpha\) of the oscillator (in radians).
  • f_detuning – detuning frequency \(\delta\) (in Hz).
  • f_rabi – Rabi frequency \(\Omega\) of the oscillator. Connected to the time of the pulse as \(\theta = \Omega t_{\mathrm{pulse}}\).
  • seed – RNG seed for calls with theta_noise > 0.

Samplers

class beclab.samplers.DensityIntegralSampler(wfs_meta, axes=None, beam_splitter=None, theta=0, no_values=False)

Bases: reiknacontrib.integrator.Sampler

Collects the density integrated over given axes.

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • axes – indices of axes to integrate over (integrates over all axes if not given).
  • beam_splitter – a BeamSplitter object. If given, it will be applied to the wavefunction before measuring populations.
  • theta – a rotation angle to pass to the beam splitter.
  • no_values – if True, no per-trajectory values will be preserved.
class beclab.samplers.DensitySliceSampler(wfs_meta, fixed_axes={}, beam_splitter=None, theta=0, no_values=False)

Bases: reiknacontrib.integrator.Sampler

Collects a density slice with given fixed indices.

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • fixed_axes – a dictionary {axis: value} of fixed indices for the slice.
  • beam_splitter – a BeamSplitter object. If given, it will be applied to the wavefunction before measuring populations.
  • theta – a rotation angle to pass to the beam splitter.
  • no_values – if True, no per-trajectory values will be preserved.
class beclab.samplers.EnergySampler(wfs_meta, system)

Bases: reiknacontrib.integrator.Sampler

Collects the total energy of a BEC (see beclab.meters.EnergyMeter for details).

Parameters:
class beclab.samplers.InteractionSampler(wfs_meta)

Bases: reiknacontrib.integrator.Sampler

Collects the integral interaction \(I = \int \Psi_1^* \Psi_2 d\mathbf{x}\) (both mean and per-component).

Parameters:wfs_meta – a WavefunctionSetMetadata object.
class beclab.samplers.PopulationSampler(wfs_meta, beam_splitter=None, theta=0)

Bases: reiknacontrib.integrator.Sampler

Collects the component populations (both mean and per-trajectory).

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • beam_splitter – a BeamSplitter object. If given, it will be applied to the wavefunction before measuring populations.
  • theta – a rotation angle to pass to the beam splitter.
class beclab.samplers.PsiSampler

Bases: reiknacontrib.integrator.Sampler

Collects the wavefunction for each trajectory.

class beclab.samplers.StoppingEnergySampler(wfs, system, limit=1e-06)

Bases: reiknacontrib.integrator.Sampler

Same as EnergySampler, but raises reiknacontrib.integrator.StopIntegration when the relative difference of the newly collected sample and the previous one is less than limit.

class beclab.samplers.VisibilitySampler(wfs_meta)

Bases: reiknacontrib.integrator.Sampler

Collects the visibility \(V = 2 \int \Psi_1^* \Psi_2 d\mathbf{x} / (N_1 + N_2)\) (both mean and per-component).

Parameters:wfs_meta – a WavefunctionSetMetadata object.

Filters

class beclab.filters.NormalizationFilter(wfs_meta, target_Ns)

Bases: reiknacontrib.integrator.Filter

Renormalizes the wavefunction to the target population.

Parameters:

Meters

class beclab.meters.DensityIntegralMeter(wfs_meta, axes=None)

Measures the integral of per-component density over chosen axes.

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • axes – indices of axes to integrate over (integrates over all axes if not given).
__call__(wfs_data)

Returns a numpy array with the shape (trajectories, components, size) with the projected density.

__weakref__

list of weak references to the object (if defined)

class beclab.meters.DensitySliceMeter(wfs_meta, fixed_axes={})

Measures the projection of per-component density on a chosen axis (in other words, integrates the density over all the other axes).

Parameters:
  • wfs_meta – a WavefunctionSetMetadata object.
  • fixed_axes – a dictionary {axis: value} of fixed indices for the slice.
__call__(wfs_data)

Returns a numpy array with the shape (trajectories, components, size) with the projected density.

__weakref__

list of weak references to the object (if defined)

class beclab.meters.EnergyMeter(wfs_meta, system)

Measures the energy of a BEC:

\[E = \sum_{j=1}^C \int \Psi_j^* \left( - \frac{\nabla^2}{2 m} + V_j + \sum_{k=1}^C \frac{g_{jk}}{2} \vert \Psi_k \vert^2 \right) \Psi_j d\mathbf{x}.\]
Parameters:
__weakref__

list of weak references to the object (if defined)

Constants

class beclab.constants.Component(atom_name, atom_mass, f, mf)

Represents a BEC component.

atom_name

A string with the type of the atom.

atom_mass

Atom mass (in a.m.u).

f

Total spin of the atom.

mf

Total spin projection.

m

Atom mass (in kg).

class beclab.constants.FeshbachResonance(value, background_scattering, width, decay_width)

Represents a Feshbach resonance.

value

Resonance magnetic field (in G).

background_scattering

Background scattering length (in Bohr radii).

width

Resonance width (in G).

decay width

Resonance decay width (in G).

beclab.constants.HBAR = 1.054571628e-34

Planck constant (in J*s).

beclab.constants.M_PROTON = 1.6588053425287357e-27

Proton mass (in kg).

beclab.constants.R_BOHR = 5.2917720859e-11

Bohr radius (in m).

beclab.constants.effective_area(m, f1, f2)

Returns the effective area (in squared meters) for a pseudo-1D trap with transverse frequencies f1 and f2 (in Hz) and atoms with mass m (in kg).

beclab.constants.feshbach_interaction(m, B, resonance)

Returns a tuple of the scattering length (in Bohr radii) and the corresponding two-body loss rate \(\kappa\) (in \(\mathrm{m}^3/\mathrm{s}\)) for the given Feshbach resonance.

Parameters:
  • m – atom mass (kg).
  • B – magnetic field (G).
  • resonance – a FeshbachResonance object.
beclab.constants.loss_rate(*comps, **kwds)

Returns the loss rate (in SI units) for a list of one or more Component objects in the magnetic field B (in G). The loss rate returned is the theoretical one (\(\kappa\)).

beclab.constants.magical_field_Rb87_1m1_2p1 = 3.228

The “magical” magnetic field (in G) for \(^{87}\mathrm{Rb}\), \(\vert F=1,\, m_F=-1 \rangle\) and \(^{87}\mathrm{Rb}\), \(\vert F=2,\, m_F=+1 \rangle\) interaction.

beclab.constants.mu_tf_1d(freq, N, m, g)

Returns the TF-approximated chemical potential (in J) for N atoms of mass m and elastic interaction coefficient g (in \(\mathrm{J} \cdot \mathrm{m}\)) in a 1D trap with frequency freq (in Hz).

beclab.constants.mu_tf_3d(freqs, N, m, g)

Returns the TF-approximated chemical potential (in J) for N atoms of mass m and elastic interaction coefficient g (in \(\mathrm{J} \cdot \mathrm{m}^3\)) in a 3D trap with frequencies freqs (in Hz).

beclab.constants.rb87_1_1 = Component(Rb,87,1,1)

\(^{87}\mathrm{Rb}\), \(\vert F=1,\, m_F=+1 \rangle\).

beclab.constants.rb87_1_minus1 = Component(Rb,87,1,-1)

\(^{87}\mathrm{Rb}\), \(\vert F=1,\, m_F=-1 \rangle\).

beclab.constants.rb87_2_1 = Component(Rb,87,2,1)

\(^{87}\mathrm{Rb}\), \(\vert F=2,\, m_F=-1 \rangle\).

beclab.constants.rb87_2_2 = Component(Rb,87,2,2)

\(^{87}\mathrm{Rb}\), \(\vert F=2,\, m_F=+2 \rangle\).

beclab.constants.rb87_2_minus1 = Component(Rb,87,2,-1)

\(^{87}\mathrm{Rb}\), \(\vert F=2,\, m_F=-1 \rangle\).

beclab.constants.resonance_Rb_9_1 = <beclab.constants.FeshbachResonance object>

The Feshbach resonance for Rb87 at 9.1G. See Kaufman et al., PRA 80 050701 (2009).

beclab.constants.scattering_3d(a, m)

Returns the 3D elastic interaction coefficient (in \(\mathrm{J} \cdot \mathrm{m}^3\)) for the scattering length a (in Bohr radii) and atoms with mass m (in kg).

beclab.constants.scattering_length(comp1, comp2, B=None)

Returns the scattering length (in Bohr radii) for two Component objects in the magnetic field B (in G).

beclab.constants.scattering_matrix(comps, B=None)

Returns a numpy array of the elastic interaction coefficients for the given list of Component objects and magnetic field B (in G).

Indices and tables