Skip to content

Observables

Generated from the docstrings in the source code.

expectation

coefficients

momentum

current

Probability current and angular momentum, with derivatives evaluated spectrally (exact for band-limited wavefunctions on the periodic grid).

gradient(psi: np.ndarray, grid)

Spectral gradient (∂ψ/∂x, ∂ψ/∂y, ∂ψ/∂z) truncated to the grid dimension.

Source code in quantumlab/observables/current.py 
11
12
13
14
def gradient(psi: np.ndarray, grid):
    """Spectral gradient (∂ψ/∂x, ∂ψ/∂y, ∂ψ/∂z) truncated to the grid dimension."""
    psi_k = fftn(psi)
    return tuple(ifftn(1j * K * psi_k) for K in wavevectors(grid))

probability_current(wf, hbar: float = 1.0, m: float = 1.0)

Probability current density j = (ℏ/m) Im(ψ* ∇ψ), returned as a tuple of arrays (j_x[, j_y[, j_z]]). It satisfies ∂|ψ|²/∂t + ∇·j = 0.

Source code in quantumlab/observables/current.py 
17
18
19
20
21
22
23
def probability_current(wf, hbar: float = 1.0, m: float = 1.0):
    """
    Probability current density j = (ℏ/m) Im(ψ* ∇ψ), returned as a tuple of
    arrays (j_x[, j_y[, j_z]]). It satisfies ∂|ψ|²/∂t + ∇·j = 0.
    """
    psi = wf.psi
    return tuple((hbar / m) * np.imag(np.conj(psi) * d) for d in gradient(psi, wf.grid))

angular_momentum_z_expectation(wf, hbar: float = 1.0) -> float

⟨L_z⟩ = ⟨ψ| −iℏ (x ∂_y − y ∂_x) |ψ⟩ for 2D and 3D wavefunctions.

Source code in quantumlab/observables/current.py 
26
27
28
29
30
31
32
def angular_momentum_z_expectation(wf, hbar: float = 1.0) -> float:
    """⟨L_z⟩ = ⟨ψ| −iℏ (x ∂_y − y ∂_x) |ψ⟩ for 2D and 3D wavefunctions."""
    if not isinstance(wf.grid, (Grid2D, Grid3D)):
        raise TypeError('L_z requires a 2D or 3D wavefunction')
    d_dx, d_dy = gradient(wf.psi, wf.grid)[:2]
    Lz_psi = -1j * hbar * (wf.grid.X * d_dy - wf.grid.Y * d_dx)
    return float(np.real(np.sum(np.conj(wf.psi) * Lz_psi) * volume_element(wf.grid)))

autocorrelation(psi0, psi_t, grid) -> complex

C(t) = ⟨ψ(0)|ψ(t)⟩. Its Fourier transform peaks at the eigenenergies.

Source code in quantumlab/observables/current.py 
35
36
37
def autocorrelation(psi0, psi_t, grid) -> complex:
    """C(t) = ⟨ψ(0)|ψ(t)⟩. Its Fourier transform peaks at the eigenenergies."""
    return complex(np.sum(np.conj(psi0) * psi_t) * volume_element(grid))

vortex_positions(wf, min_density: float = 0.2, smoothing: float = 1.0)

Quantised vortices of a 2D wavefunction from the phase winding around each grid plaquette, ∮∇φ·dl = 2πq.

Parameters:

Name Type Description Default
min_density float

Only plaquettes where the smoothed density exceeds this fraction of its maximum are kept. This removes the "ghost" phase singularities in the dilute region outside the condensate.

 0.2
smoothing float

Gaussian smoothing length (physical units) used for that density mask; it should exceed the core size so the cores themselves are not masked.

 1.0

Returns:

Name Type Description
positions (ndarray, shape(n, 2))

(x, y) of the plaquette centres containing a vortex.
charges (ndarray, shape(n))

Winding numbers q (±1 for singly quantised vortices).
Source code in quantumlab/observables/current.py 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def vortex_positions(wf, min_density: float = 0.2, smoothing: float = 1.0):
    """
    Quantised vortices of a 2D wavefunction from the phase winding around each
    grid plaquette, ∮∇φ·dl = 2πq.

    Parameters
    ----------
    min_density : float
        Only plaquettes where the smoothed density exceeds this fraction of
        its maximum are kept. This removes the "ghost" phase singularities in
        the dilute region outside the condensate.
    smoothing : float
        Gaussian smoothing length (physical units) used for that density mask;
        it should exceed the core size so the cores themselves are not masked.

    Returns
    -------
    positions : ndarray, shape (n, 2)
        (x, y) of the plaquette centres containing a vortex.
    charges : ndarray, shape (n,)
        Winding numbers q (±1 for singly quantised vortices).
    """
    from scipy.ndimage import gaussian_filter
    if not isinstance(wf.grid, Grid2D):
        raise TypeError('vortex_positions requires a 2D wavefunction (use a slice for 3D)')
    g = wf.grid
    phase = np.angle(wf.psi)

    def wrapped(a, b):
        return (b - a + np.pi) % (2.0 * np.pi) - np.pi

    circulation = (wrapped(phase[:-1, :-1], phase[1:, :-1]) + wrapped(phase[1:, :-1], phase[1:, 1:])
                   + wrapped(phase[1:, 1:], phase[:-1, 1:]) + wrapped(phase[:-1, 1:], phase[:-1, :-1]))
    charge = np.rint(circulation / (2.0 * np.pi)).astype(int)
    rho = gaussian_filter(wf.probability_density, sigma=(smoothing / g.dx, smoothing / g.dy), mode='constant')
    keep = (charge != 0) & (rho[:-1, :-1] > min_density * rho.max())
    i, j = np.nonzero(keep)
    positions = np.column_stack([g.x[i] + 0.5 * g.dx, g.y[j] + 0.5 * g.dy])
    return positions, charge[i, j]

nonlinear

interaction_energy(wf, nonlinearity: float) -> float

Mean-field interaction energy E_int = (g/2) ∫ |ψ|⁴ dV.

Source code in quantumlab/observables/nonlinear.py 
6
7
8
def interaction_energy(wf, nonlinearity: float) -> float:
    """Mean-field interaction energy E_int = (g/2) ∫ |ψ|⁴ dV."""
    return float(0.5 * nonlinearity * np.sum(wf.probability_density ** 2) * volume_element(wf.grid))

gross_pitaevskii_energy(wf, potential, nonlinearity: float, hbar: float = 1.0, m: float = 1.0) -> float

GP energy functional E[ψ] = ⟨T⟩ + ⟨V⟩ + (g/2)∫|ψ|⁴ (conserved by real-time GPE dynamics).

Source code in quantumlab/observables/nonlinear.py 
11
12
13
14
def gross_pitaevskii_energy(wf, potential, nonlinearity: float, hbar: float = 1.0, m: float = 1.0) -> float:
    """GP energy functional E[ψ] = ⟨T⟩ + ⟨V⟩ + (g/2)∫|ψ|⁴ (conserved by real-time GPE dynamics)."""
    return float(kinetic_energy_expectation(wf, hbar, m) + potential_energy_expectation(wf, potential)
                 + interaction_energy(wf, nonlinearity))

chemical_potential(wf, potential, nonlinearity: float, hbar: float = 1.0, m: float = 1.0) -> float

μ = ⟨T⟩ + ⟨V⟩ + g∫|ψ|⁴, the eigenvalue of the stationary GP equation.

Source code in quantumlab/observables/nonlinear.py 
17
18
19
20
def chemical_potential(wf, potential, nonlinearity: float, hbar: float = 1.0, m: float = 1.0) -> float:
    """μ = ⟨T⟩ + ⟨V⟩ + g∫|ψ|⁴, the eigenvalue of the stationary GP equation."""
    return float(kinetic_energy_expectation(wf, hbar, m) + potential_energy_expectation(wf, potential)
                 + 2.0 * interaction_energy(wf, nonlinearity))

coupled_gp_energy(states, potentials, couplings, hbar: float = 1.0, m: float = 1.0) -> float

Energy functional of a multi-component condensate, E = Σ_i (⟨T⟩_i + ⟨V_i⟩) + ½ Σ_ij g_ij ∫ n_i n_j dV, with n_i = |ψ_i|².

Source code in quantumlab/observables/nonlinear.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
def coupled_gp_energy(states, potentials, couplings, hbar: float = 1.0, m: float = 1.0) -> float:
    """
    Energy functional of a multi-component condensate,
    E = Σ_i (⟨T⟩_i + ⟨V_i⟩) + ½ Σ_ij g_ij ∫ n_i n_j dV, with n_i = |ψ_i|².
    """
    g = np.asarray(couplings, dtype=float)
    dV = volume_element(states[0].grid)
    E = sum(kinetic_energy_expectation(s, hbar, m) + potential_energy_expectation(s, p)
            for s, p in zip(states, potentials))
    n = [s.probability_density for s in states]
    E += 0.5 * sum(g[i, j] * float(np.sum(n[i] * n[j]) * dV)
                   for i in range(len(states)) for j in range(len(states)))
    return float(E)

miscibility_overlap(wf_a, wf_b) -> float

Spatial overlap Λ = ∫ n_a n_b / sqrt(∫ n_a² ∫ n_b²) of two components: 1 for identical density profiles, → 0 for fully separated ones.

Source code in quantumlab/observables/nonlinear.py 
38
39
40
41
42
43
44
def miscibility_overlap(wf_a, wf_b) -> float:
    """
    Spatial overlap Λ = ∫ n_a n_b / sqrt(∫ n_a² ∫ n_b²) of two components:
    1 for identical density profiles, → 0 for fully separated ones.
    """
    na, nb = wf_a.probability_density, wf_b.probability_density
    return float(np.sum(na * nb) / np.sqrt(np.sum(na ** 2) * np.sum(nb ** 2)))