gpe.bec2

Two component BEC in 1D.

Here we demonstrate the dynamics a two component BEC such as 87Rb.

Simplifications:

• g_aa = g_bb = g_ab = g

• m_a = m_b = m

• Assume that only total particle number is conserved.

Attributes

u

Classes

StateDFTBase	Underlying implementation of features needed for IStateDFT.
_StateBase	Two component state.
GPEMixin	Mixing providing the 2-component GPE equation of state.
StateGPEBase	Two component state with GPE equation of state.
HOMixin	Helper mixin class for backwards compatibility.
StateScaleBase	This state implements the scaling from Castin and Dum.
StateScaleHO	This state implements the scaling from Castin and Dum for

Module Contents

u

class StateDFTBase(**kw)

Bases: gpe.bec.StateDFTBase

Underlying implementation of features needed for IStateDFT.

axpy1(x, a)

axpy where a is a scalar.

axpy2(x, a)

Version of axpy where a has 2 components.

Uses self._axpy for each component.

axpy(x, a=1)

Version of axpy where a has 2 components.

Uses self._axpy for each component.

_axpy_numpy(y, x, a=1)

Default numpy version of axpy.

apply_V(Va, Vb, Vab, exp=False)

fill(value)

braket_GPU(y)

Return the dot product of self.conjugate() with y.

Parameters:

y (State) – The states for which the braket will be computed.

get_Vext_mu_GPU()

Return (Va, Vb, Vab) with the chemical potential subtracted if initializing.

The chemical potential should be subtracted if initializing or minimizing to get the initial state. Minimization might be done with imaginary time evolution, which should be done with negative times.

Thus, we check for self.initializing, self.t < 0, and make sure that self.mu is valid.

get_Vext_GPU()

Return (Va, Vb, Vab), the external potentials.

class _StateBase(basis=None, Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), symmetric_grid=False, twist=None, v_x=0, kwz2=None, t=0, hbar=u.hbar, ms=(u.m, u.m), mus=None, x_TF=None, cooling_phase=1.0, constraint='N', PGPE=False, mu_Bs=(0 * u.mu_B, 0 * u.mu_B), k_r=0.0, rotating_phases=False, basis_args={}, t_final=None)

Bases: StateDFTBase

Two component state.

This state is designed for a spin-orbit (SO) coupled BEC in a harmonic confining potential with the SO field along the x-axis, and a possible magnetic field gradient also along the x-axis.

Parameters:

• basis (IBasis, None) – If provided, then this is used to specify the basis. If not, then a periodic basis is constructed from Nxyz and Lxyz.

• constraint ('N' | 'Nab') – If ‘N’, then constrain the total density allowing particles to convert from one species to the other. If ‘Nab’ Then independently constrain Na and Nb.

• twist ((float, float)) – Twisted boundary conditions for species a and b. This is the twist angle.

• x_TF (float) – Thomas Fermi “radius” for setting the initial state. The initial state will be set so that the density will fall to zero at this point x=x_TF (with y and z in the middle of the trap). If None, then we default to 80% of the length.

• PGPE (bool) – If True, then the code implements the PGPE, imposing a cutoff at kc

t0 = 63.50779925891489

basis = None

data

t = 0

hbar = 1.0

ms

mus = None

x_TF = None

cooling_phase = 1.0

constraint = 'N'

PGPE = False

k_r = 0.0

twist

v_x = 0

kwz2 = None

mu_Bs

property rotating_phases

Get rotating_phase.

_rotating_phases = False

init()

Initialize the state.

This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor __init__() so that the user can change them later and the reinitialize the state. We also call this function from the pre_evolve_hook() so that it is called before any evolution takes place. For this reason, we should not modify the state here.

property x

Flat x abscissa as a numpy array.

set_initial_state(mus=None, x_TF=None)

Set the state using the Thomas Fermi (TF) approximation from either mus or x_TF (pick only one or the other).

Parameters:

• mus ((float, float)) – Fixed chemical potentials.

• x_TF (float) – Position defining the Thomas Fermi “radius”. (The external potential is evaluated at this position and this is used to set mu.)

get_mus_from_Vs_TF(Vs_TF)

Return the corrected chemical potential from Vs_TF.

In some cases, the chemical potential may differ from the value of the external potentials at V(x_TF) due to kinetic energy shifts (in the SOC case for example) or due to the energy of radial excitations (see the tube codes). This function adds the appropriate correction.

Parameters:

Vs_TF ((float, float)) – External potentials at the Thomas Fermi “radius”. (The external potential is evaluated at this position and this is used to get mus.)

get_Vs_TF_from_mus(mus)

Return Vs_TF from the chemical potentials mus.

Parameters:

mus ((float, float)) – Physical chemical potentials (i.e. what you would pass to the minimizer).

get_Vs_TF(x_TF, V_ext=None)

Return the Thomas Fermi chemical potential at x_TF.

Parameters:

x_TF (float) – Position defining the Thomas Fermi “radius”. (The external potential is evaluated at this position and this is used to get mu.)

_get_kx2(kx=None, Lx=None)

Return the effective kx**2 and twist_phase_x for the kinetic energy matrix. This will be passed as the kx2 and twist_phase_x arguments to the basis.lagrangian function.

get_xyz_GPU()

property xyz

get_metric_GPU()

property metric

get_psi_GPU()

Return the physical wavefunction (applying any twist).

set_psi(psi)

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

get_density_GPU()

get_V_GPU()

Return the complete potentials (Va, Vb, Vab) - internal and external.

property dim

Dimension of the state.

property shape

Shape of the state Nxyz (not including components).

property bcast

Return a set of indices suitable for broadcasting masses etc.

property E_max

Return the maximum kinetic energy in the basis.

This is useful when using evolvers as convergence should be obtained when the time-step is roughly:

dt = 0.1 * state.hbar / state.E_max

See t_scale.

property t_scale

Return the smallest time-scale for the problem.

Evolvers - especially the ABM evolvers - should use a dt=0.1*t_scale or so. If much smaller dt values are required for convergence, then it usually indicates that your lattice spacing is too larger. Likewise, if you can get away with much larger dt values, then your lattice spacing might be unnecessarily small.

Required by Simulation.

get_Ns_GPU()

get_N_GPU()

integrate_GPU(a)

Integrate over both components, but do not sum.

compute_dy_dt(dy, subtract_mu=True)

Return dy_dt storing the results in dy.

Parameters:

subtract_mu (bool) –

If True, then subtract the chemical potential such that dy_dt is orthogonal to the original state y. This will minimize the evolution of the overall phase during real-time evolution (which can reduce numerical errors) and will ensure that evolution under imaginary or complex time will preserve particle number.

This should not be set if computing physical energy of the state, however, which is why it is a parameter.

linear = False

apply_exp_K(dt)

Apply \(e^{-i K dt}\) in place.

apply_exp_V(dt, state)

Apply \(e^{-i V dt}\) in place using state for any nonlinear dependence in V. (Linear problems should ignore state.)

normalize(s=None)

Normalize the state, return the scale factors and number (s, N).

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

Apply the laplacian multiplied by factor to the state.

Parameters:

• factor (array-like) – The result will be multiplied by this factor.

• exp (bool) – If True then exp(factor*laplacian)(y) will be computed instead.

get_energy(**kw)

get_Hy(subtract_mu=False)

get_energy_density(a=True, b=True, ab=True)

Return the energy density.

Parameters:

• a (bool) – Include the energies associated with this species. This includes the kinetic energy and the self-interaction, but no inter-species interactions.

• b (bool) – Include the energies associated with this species. This includes the kinetic energy and the self-interaction, but no inter-species interactions.

• ab (bool) – Include only the interaction energy.

get_mu()

Compute the chemical potential for convenience only.

plot(log=False)

plot_k()

Plot the momentum distribution.

evolve_to(t_end, dt_t_scale=0.2, Evolver=evolvers.EvolverABM, evolve_steps=200, callback=None)

Evolve state to self.t = t_end.

Parameters:

• t_end (float) – Evolves the state for t_end time units. Remember to convert times by self.t_unit.

• dt_t_scale (float) – Scales dt.

• Evolver (IEvolver) – Pick the relevent evolver, either evolvers.EvolverABM or evolvers.EvolverSplit.

• evolve_steps (float) – Number of evolution steps between callbacks if defined.

• callback (function) – Any function that takes the state as an argument, like plotting or visualization.

class GPEMixin

Mixing providing the 2-component GPE equation of state.

init()

get_Vint_GPU(state=None)

Return (Va, Vb), the “internal” mean-field potential.

get_Eint(state=None)

Return the “internal” mean-field energy-densities (Ea, Eb, Eab).

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.

get_ns_TF(Vs_TF, V_ext=None, gs=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.

• Note (Ignores gab.)

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

Helper the total TF density for a single component.

Parameters:

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

class StateGPEBase(gs=None, **kw)

Bases: GPEMixin, _StateBase

Two component state with GPE equation of state.

gs = (1.0, 1.0, 1.0)

class HOMixin(ws=None, Omega=None, delta=None, **kw)

Bases: gpe.utils.GPUHelper

Helper mixin class for backwards compatibility.

This provides a default get_Vext_GPU().

_Vext = None

init()

get_Vext_GPU(state=None)

Return (V_a, V_b, V_ab), the external potentials.

class StateScaleBase(basis=None, Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), symmetric_grid=False, twist=None, v_x=0, kwz2=None, t=0, hbar=u.hbar, ms=(u.m, u.m), mus=None, x_TF=None, cooling_phase=1.0, constraint='N', PGPE=False, mu_Bs=(0 * u.mu_B, 0 * u.mu_B), k_r=0.0, rotating_phases=False, basis_args={}, t_final=None)

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.

init()

Initialize the state.

This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor __init__() so that the user can change them later and the reinitialize the state. We also call this function from the pre_evolve_hook() so that it is called before any evolution takes place. For this reason, we should not modify the state here.

_get_scale_factor_and_phase_GPU()

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_psi_GPU()

Return the physical wavefunction (applying any twist).

set_psi(psi)

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

get_xyz_GPU()

get_metric_GPU()

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

Apply the laplacian multiplied by factor to the state.

Parameters:

• factor (array-like) – The result will be multiplied by this factor.

• exp (bool) – If True then exp(factor*laplacian)(y) will be computed instead.

get_V_GPU()

Return the complete potential V - internal and external.

This version includes the correction from the scaled coordinates.

class StateScaleHO(basis=None, Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), symmetric_grid=False, twist=None, v_x=0, kwz2=None, t=0, hbar=u.hbar, ms=(u.m, u.m), mus=None, x_TF=None, cooling_phase=1.0, constraint='N', PGPE=False, mu_Bs=(0 * u.mu_B, 0 * u.mu_B), k_r=0.0, rotating_phases=False, basis_args={}, t_final=None)

Bases: StateScaleBase

This state implements the scaling from Castin and Dum for harmonic oscillators. This provides the required get_lambdas() function but now requires the user specify the time-dependence of the trapping frequencies in get_ws().

get_ws(t=None)

Return the trapping frequencies at time t.

init()

Initialize the state.

This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor __init__() so that the user can change them later and the reinitialize the state. We also call this function from the pre_evolve_hook() so that it is called before any evolution takes place. For this reason, we should not modify the state here.

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.

_rhs(q, t)

RHS for lambda(t) ODE.