The Gas class#

The Gas class represents the quantum gas to be simulated.

It is initialized specifying the atomic species, the number of atoms, and details on the region of space to be studied. After setting a first wave function of the gas via the psi or psik attributes, the ground_state() and propagate() methods can be used to evolve the wave function in time. ground_state() performs imaginary time evolution to find the ground state of the system, while propagate() performs real time evolution to study its dynamics. Both methods use the split-step Fourier method to solve the Gross-Pitaevskii equation.

The atomic species is specified by the element parameter, which is a string ('87Rb' by default). If element is a species supported by GPE, its parameters are automatically loaded. The only supported species at the moment are '87Rb' and '39K', which are defined in elements_dict. This is a dictionary where each element is associated with its mass and the angular frequency of its \(d_2\) line, e.g.:

import numpy as np
import scipy.constants as spconsts

# dictionary stored in torchgpe.utils.elements

elements_dict = {
    "87Rb": {
        "m": spconsts.physical_constants["atomic mass constant"][0] * 87,
        "omega d2": 2 * np.pi * 384.2304844685e12,
    },
    "39K": {
        "m": spconsts.physical_constants["atomic mass constant"][0] * 39,
        "omega d2": 2 * np.pi * 391.01617003e12,
    },
}

Additional atomic species can be specified by the user by adding new entries to the dictionary. Similarly, additional properties of the existing species can be added by modifying the corresponding entries.

Adimensionalization#

In order to minimize the numerical errors due to the scale of the system, the wave function is evolved in dimensionless units. The length unit is provided by the adimensionalization_length parameter (hereafter referred to as \(l\)); other adimensionalization parameters are inferred from it:

adimensionalization length: \(l\) (available as adim_length)
Adimensionalization pulse: \(\omega_l = \frac{\hbar}{m l^2}\) (available as adim_pulse)
Adimensionalization time: \(\tau_l = \frac{1}{\omega_l}\)
Adimensionalization energy: \(E_l = \hbar \omega_l\)

Therefore, to transform from adimensionalised units to SI units, one should multiply by the corresponding adimensionalization parameter. For example:

x_SI = x_adim * gas.adim_length        # lengths
w_SI = w_adim * gas.adim_pulse         # frequencies and pulses
t_SI = t_adim / gas.adim_pulse         # times
E_SI = E_adim * hbar * gas.adim_pulse  # energies

Computational grid#

The real space wave function is defined and evolved on a grid specified by the N_grid, grid_size_x and grid_size_y parameters. The grid is centered in the origin, and the grid spacing is given by grid_size/((N_grid-1)*adimensionalization_length). For example, to study a wave function in a region of space \(10\mu m\) wide, one should initialize the Gas class with the argument grid_size=10e-6. Assuming \(l=1\mu m\) and N_grid=101, the resulting grid will span values from \(-5\) to \(5\) in adimensional units, with a grid spacing of \(\Delta = 0.1\) (corresponding to \(100 \, nm\) in SI units).

Similarly, the wave function in momentum space is defined on a grid that is related to the real space one by the application of a Fourier transform. In this case, before applying the transform, we pad the wave function to reduce the presence of boundary effects. Therefore, the wave function in momentum space has 2*N_grid points per dimension, ranging values from \(-\frac{\pi}{\Delta}\) to \(\frac{\pi}{\Delta}\) in adimensional units. The adimensional unit of momentum is given by \(\frac{1}{l}\) and hence the spacing between two sampled momenta is \(\tilde{\Delta} = \pi/(N_{grid} \cdot \Delta)\).

Warning

The grid will have 2*N_grid points per side only if N_grid is even. Otherwise, the grid will have 2*N_grid-1 points per side.

Initialisation of the wave function#

The wave function can be initialized by assigning a value to either the psi or psik attributes. The Gas class is written in such a way that the user can always read the wave function in both real and momentum space, no matter which one of the two is the last updated. In case the user tries to read the wave function in the space that is not the last updated, the wave function is automatically computed via a Fourier transform.

For example, to initialize the real space wave function to a Gaussian wave packet of size \(1\mu m\), one can use:

from torchgpe.bec2D import Gas
import torch

bec = Gas(N_particles=2e5, grid_size=1e-5)
bec.psi = torch.exp( -(bec.X**2 + bec.Y**2)/(2*(1e-6/bec.adim_length)**2) )

Note that since the grid is expressed in adimensional units, the Gaussian width must be divided by the adimensionalization length. Also, normalizing the wave function is not necessary, as the library takes care of it.

Evolution#

GPE evolves the wave function of the gas integrating the Gross-Pitaevskii equation in time. The equation is solved using the split-step Fourier method, which consists in alternating the application of the kinetic and potential operators in real and momentum space.

Consider an Hamiltonian \(\hat{H} = \hat{T} + \hat{V}\), where \(\hat{T}\) is the kinetic term (diagonal in momentum space) and \(\hat{V}\) is the potential one, assumed diagonal in real space. The time evolution from time \(t\) to time \(t+\Delta t\) is given by the operator

\[U = \exp(-\frac{i}{\hbar}\hat{H}\Delta t).\]

Since \(\hat{T}\) and \(\hat{V}\) don’t commute, the operator \(U\) cannot be factorized. However, defining the dimensionless operators

\[\tilde{T} = -\frac{i}{\hbar}\hat{T}\Delta t \quad \text{and} \quad \tilde{V} = -\frac{i}{\hbar}\hat{V}\Delta t\]

and making use of the Baker-Campbell-Hausdorff formula, one can write

\[\begin{split}U &= \exp(\tilde{T}+\tilde{V}) = \exp(\tilde{T}/2 + \tilde{V} + \tilde{T}/2)\\ &= e^{\tilde{T}/2}e^{\tilde{V} + \tilde{T}/2}\exp(-\frac{1}{2}\left[\tilde{T}/2,\tilde{V}\right] + \mathcal{O}\left(\Delta t^3\right))\\ &= e^{\tilde{T}/2}e^{\tilde{V}}e^{\tilde{T}/2}\exp(-\frac{1}{2}\left[\tilde{T}/2,\tilde{V}\right] -\frac{1}{2}\left[\tilde{V},\tilde{T}/2\right] + \mathcal{O}\left(\Delta t^3\right))\\ &= e^{\tilde{T}/2}e^{\tilde{V}}e^{\tilde{T}/2}e^{\mathcal{O}\left(\Delta t^3\right)}.\end{split}\]

That is, the subsequent application of the kinetic propagator \(e^{\tilde{T}/2}\) followed by the potential one \(e^{\tilde{V}}\) and the kinetic propagator one last time is equivalent to the application of the full operator \(U\) up to third order in \(\Delta t\).

Ground state#

The ground_state() method uses the split-step Fourier method to evolve the wave function in imaginary time to find the ground state of the system. The method takes as input the number of iterations to perform, the time step to use, the list of the potentials acting on the system and, optionally, a set of callbacks to monitor the evolution. Note that imaginary time evolution does not support the use of a time-dependent potential. After checking that the time step is indeed a purely imaginary number, the method integrates the Gross Pitaevskii equation taking care of renormalising the wave function after each step. At the end of the propagation, the psi and psik attributes are updated to the final wave function.

Real time propagation#

The propagate() method uses the split-step Fourier method to evolve the wave function in real time. The method takes as input the time step to use, the final time, the list of the potentials acting on the system and, optionally, a set of callbacks to monitor the evolution. Note that the time step can be slightly adjusted to ensure that the wave function is evaluated at the final time. The method integrates the Gross Pitaevskii equation taking care of renormalising the wave function after each step. At the end of the propagation, the psi and psik attributes are updated to the final wave function.