# 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)

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)

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)

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). 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
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)

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: wfs_meta – a WavefunctionSetMetadata object. system – the System object the wavefunction corresponds to.
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: wfs_meta – a WavefunctionSetMetadata object. target_Ns – a tuple of target populations for each component.

### 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: wfs_meta – a WavefunctionSetMetadata object. system – the System object the wavefunction corresponds to.
__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

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).