GitHub

Module Hamiltonians

This module contains definitions of Hamiltonians, in particular specific physical models of interest. These are organised by means of an interface around the abstract type AbstractHamiltonian, in the spirit of the AbstractArray interface as discussed in the Julia Documentation.

The Hamiltonians can be used for projector quantum Monte Carlo with ProjectorMonteCarloProblem or for exact diagonalization with ExactDiagonalizationProblem, see Exact Diagonalization.

Rimu.Hamiltonians — Module

The module Rimu.Hamiltonians defines types and functions for working with Hamiltonians.

Exported concrete Hamiltonian types

Real space Hubbard models

Momentum space Hubbard models

Harmonic oscillator models

Other

Wrappers

Observables

Interface for working with Hamiltonians

source

Here is a list of fully implemented model Hamiltonians. There are several variants of the Hubbard model in real and momentum space, as well as some other models.

Real space Hubbard models

Rimu.Hamiltonians.HubbardReal1D — Type
HubbardReal1D(address; u=1.0, t=1.0)

Implements a one-dimensional Bose Hubbard chain in real space.

\[\hat{H} = - \sum_i \left(t a_i^† a_{i+1} + t^* a_{i+1}^† a_i \right) + \frac{u}{2}\sum_i n_i (n_i-1)\]

Arguments

  • address: the starting address, defines number of particles and sites.
  • u: the interaction parameter.
  • t: the hopping strength.

See also

source
Rimu.Hamiltonians.HubbardReal1DEP — Type
HubbardReal1DEP(address; u=1.0, t=1.0, v_ho=1.0)

Implements a one-dimensional Bose Hubbard chain in real space with external potential.

\[\hat{H} = - \sum_i \left(t a_i^† a_{i+1} + t^* a_{i+1}^† a_i \right) + \sum_i ϵ_i n_i + \frac{u}{2}\sum_i n_i (n_i-1)\]

Arguments

  • address: the starting address, defines number of particles and sites.
  • u: the interaction parameter.
  • t: the hopping strength.
  • v_ho: strength of the external harmonic oscillator potential $ϵ_i = v_{ho} i^2$.

The first index is i=0 and the maximum of the potential occurs in the centre of the lattice.

See also

source
Rimu.Hamiltonians.HubbardRealSpace — Type
HubbardRealSpace(address; geometry=PeriodicBoundaries(M,), t=ones(C), u=ones(C, C), v=zeros(C, D))

Hubbard model in real space. Supports single or multi-component Fock state addresses (with C components) and various (rectangular) lattice geometries in D dimensions.

\[ \hat{H} = -\sum_{\langle i,j\rangle,σ} t_σ a^†_{iσ} a_{jσ} + \frac{1}{2}\sum_{i,σ} u_{σσ} n_{iσ} (n_{iσ} - 1) + \sum_{i,σ≠τ}u_{στ} n_{iσ} n_{iτ}\]

If v is nonzero then this calculates $\hat{H} + \hat{V}$ by adding the harmonic trapping potential

\[ \hat{V} = \sum_{i,σ,d} v_{σd} x_{di}^2 n_{iσ}\]

where $x_{di}$ is the distance of site $i$ from the centre of the trap along dimension $d$.

Address types

  • BoseFS: Single-component Bose-Hubbard model.
  • FermiFS: Single-component Fermi-Hubbard model.
  • CompositeFS: For multi-component models.

Note that a single component of fermions cannot interact with itself. A warning is produced if addressis incompatible with the interaction parameters u.

Geometries

Implemented CubicGrids for keyword geometry

Default is geometry=PeriodicBoundaries(M,), i.e. a one-dimensional lattice with the number of sites M inferred from the number of modes in address.

Other parameters

  • t: the hopping strengths. Must be a vector of length C. The i-th element of the vector corresponds to the hopping strength of the i-th component.
  • u: the on-site interaction parameters. Must be a symmetric matrix. u[i, j] corresponds to the interaction between the i-th and j-th component. u[i, i] corresponds to the interaction of a component with itself. Note that u[i,i] must be zero for fermionic components.
  • v: the trap potential strengths. Must be a matrix of size C × D. v[i,j] is the strength of the trap for component i in the jth dimension.
source
Rimu.Hamiltonians.ExtendedHubbardReal1D — Type
ExtendedHubbardReal1D(address; u=1.0, v=1.0, t=1.0, boundary_condition=:periodic, power=nothing)

Implements the extended Hubbard model on a one-dimensional chain in real space. This Hamiltonian can be either real or complex, depending on the choice of boundary_condition.

\[\hat{H} = - \sum_i \left(t a_i^† a_{i+1} + t^* a_{i+1}^† a_i \right) + \frac{u}{2}\sum_i n_i (n_i-1) + v \sum_{i,j>i} f_{j-i} n_i n_j\]

Arguments

  • address: the starting address.
  • u: on-site interaction parameter
  • v: the next-neighbor interaction
  • t: the hopping strength
  • boundary_condition The following values are supported:
    • :periodic: usual period boundary condition realising a ring geometry.
    • :hard_wall: hopping over the boundary is not allowed.
    • :twisted: like :periodic but hopping over the boundary incurs an additional factor of -1.
    • θ <: Number: like :periodic and :twisted but hopping over the boundary incurs a factor $\exp(iθ)$ for a hop to the right and $\exp(−iθ)$ for a hop to the left. With this choice the Hamiltonian will have a complex eltype whereas otherwise the eltype is determined by the type of the parameters t, u, and v.
  • power: the interaction type. The following values are supported:
    • nothing: nearest neighbour interaction (default), i.e. $f_{j-i} = δ_{j-i,1}$.
    • p<:Number: inverse distance interaction, i.e. $f_{j-i} = (j-i)^{-p}$.

See also HubbardRealSpace.

source

Momentum space Hubbard models

Rimu.Hamiltonians.HubbardMom1D — Type
HubbardMom1D(address; u=1.0, t=1.0, dispersion=hubbard_dispersion)

Implements a one-dimensional Bose Hubbard chain in momentum space.

\[\hat{H} = \sum_{k} ϵ_k n_k + \frac{u}{M}\sum_{kpqr} a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k}\]

Arguments

  • address: the starting address, defines number of particles and sites.
  • u: the interaction parameter.
  • t: the hopping strength.
  • dispersion: defines $ϵ_k =$dispersion(t, k)

See also

source
Rimu.Hamiltonians.HubbardMom1DEP — Type
HubbardMom1DEP(address; u=1.0, t=1.0, v_ho=1.0, dispersion=hubbard_dispersion)

Implements a one-dimensional Bose Hubbard chain in momentum space with harmonic external potential.

\[Ĥ = \sum_{k} ϵ_k n_k + \frac{u}{M}\sum_{kpqr} a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k} + V̂_\mathrm{ho} ,\]

where

\[\begin{aligned} V̂_\mathrm{ho} & = \frac{1}{M} \sum_{p,q} \mathrm{DFT}[V_{ext}]_{p-q} \, a^†_{p} a_q ,\\ V_\mathrm{ext}(x) &= v_\mathrm{ho} \,x^2 , \end{aligned}\]

is an external harmonic potential in momentum space, $\mathrm{DFT}[…]_k$ is a discrete Fourier transform performed by fft()[k%M + 1], and M == num_modes(address).

Arguments

  • address: the starting address, defines number of particles and sites.
  • u: the interaction parameter.
  • t: the hopping strength.
  • dispersion: defines $ϵ_k =$dispersion(t, k)
  • v_ho: strength of the external harmonic oscillator potential $v_\mathrm{ho}$.

See also HubbardMom1D, HubbardReal1DEP, Transcorrelated1D, Hamiltonians.

source
Rimu.Hamiltonians.ExtendedHubbardMom1D — Type
ExtendedHubbardMom1D(
    address; 
    u=1.0, t=1.0, v=1.0, dispersion=hubbard_dispersion, boundary_condition = 0.0
)

Implements a one-dimensional extended Hubbard chain, also known as the $t - V$ model, in momentum space.

\[\hat{H} = \sum_{k} ϵ_k n_k + \frac{1}{2M} \sum_{kpqr} (u + 2v \cos(q-p)) a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k}\]

Arguments

  • address: the starting address, defines number of particles and sites.

  • u: the interaction parameter.

  • t: the hopping strength.

  • boundary_condition: θ <: Number: hopping over the boundary incurs a factor $\exp(iθ)$ for a hop to the right and $\exp(−iθ)$ for a hop to the left.

  • dispersion: defines $ϵ_k =$dispersion(t, k + θ)

See also

source

Harmonic oscillator models

Rimu.Hamiltonians.HOCartesianContactInteractions — Type
HOCartesianContactInteractions(addr; S, η, g = 1.0, interaction_only = false, block_by_level = true)

Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions

\[\hat{H} = \sum_{i} \epsilon_\mathbf{i} n_\mathbf{i} + \frac{g}{2}\sum_\mathbf{ijkl} V_\mathbf{ijkl} a^†_\mathbf{i} a^†_\mathbf{j} a_\mathbf{k} a_\mathbf{l}.\]

For a $D$-dimensional harmonic oscillator indices $\mathbf{i}, \mathbf{j}, \ldots$ are $D$-tuples. The energy scale is defined by the first dimension i.e. $\hbar \omega_x$ so that single particle energies are

\[ \frac{\epsilon_\mathbf{i}}{\hbar \omega_x} = (i_x + 1/2) + \eta_y (i_y+1/2) + \ldots.\]

The factors $\eta_y, \ldots$ allow for anisotropic trapping geometries and are assumed to be greater than 1 so that $\omega_x$ is the smallest trapping frequency.

By default the offdiagonal elements due to the interactions are consistent with first-order degenerate perturbation theory

\[ V_{\mathbf{ijkl}} = \delta_{\epsilon_\mathbf{i} + \epsilon_\mathbf{j}} ^{\epsilon_\mathbf{k} + \epsilon_\mathbf{l}} \prod_{d \in x, y,\ldots} \mathcal{I}(i_d,j_d,k_d,l_d),\]

where the $\delta$ function indicates that the total noninteracting energy is conserved meaning all states with the same noninteracting energy are connected by this interaction and the Hamiltonian blocks according to noninteracting energy levels. Setting block_by_level = false will disable this restriction and allow coupling between basis states of any noninteracting energy level, leading to many more offdiagonals and fewer but larger blocks (the blocks are still distinguished by parity of basis states). Alternatively, see HOCartesianEnergyConservedPerDim for a model with the stronger restriction that conserves energy separately per spatial dimension. The integral $\mathcal{I}(a,b,c,d)$ is of four one dimensional harmonic oscillator basis functions, implemented in four_oscillator_integral_general.

Arguments

  • addr: the starting address, defines number of particles and total number of modes.
  • S: Tuple of the number of levels in each dimension, including the groundstate. The allowed couplings between states is defined by the aspect ratio of S .- 1. Defaults to a 1D spectrum with number of levels matching modes of addr. Will be sorted to make the first dimension the largest.
  • η: Define a custom aspect ratio for the trapping potential strengths, instead of deriving from S .- 1. This will only affect the single particle energy scale and not the interactions. The values are always scaled relative to the first dimension, which sets the energy scale of the system, $\hbar\omega_x$.
  • g: the (isotropic) bare interaction parameter. The value of g is assumed to be in trap units.
  • interaction_only: if set to true then the noninteracting single-particle terms are ignored. Useful if only energy shifts due to interactions are required.
  • block_by_level: if set to false will allow the interactions to couple all states without comparing their noninteracting energy.
Warning

num_offdiagonals is a bad estimate for this Hamiltonian. Take care when building a matrix or using QMC methods. Use get_all_blocks first then pass option col_hint = block_size to BasisSetRep to safely build the matrix.

source
Rimu.Hamiltonians.HOCartesianEnergyConservedPerDim — Type
HOCartesianEnergyConservedPerDim(addr; S, η, g = 1.0, interaction_only = false)

Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions

\[\hat{H} = \sum_{i} ϵ_i n_i + \frac{g}{2}\sum_{ijkl} V_{ijkl} a^†_i a^†_j a_k a_l,\]

with the additional restriction that the interactions only couple states with the same energy in each dimension separately. See HOCartesianContactInteractions for a model that conserves total energy.

For a $D$-dimensional harmonic oscillator indices $\mathbf{i}, \mathbf{j}, \ldots$ are $D$-tuples. The energy scale is defined by the first dimension i.e. $\hbar \omega_x$ so that single particle energies are

\[ \frac{\epsilon_\mathbf{i}}{\hbar \omega_x} = (i_x + 1/2) + \eta_y (i_y+1/2) + \ldots.\]

The factors $\eta_y, \ldots$ allow for anisotropic trapping geometries and are assumed to be greater than 1 so that $\omega_x$ is the smallest trapping frequency.

Matrix elements $V_{\mathbf{ijkl}}$ are for a contact interaction calculated in this basis using first-order degenerate perturbation theory.

\[ V_{\mathbf{ijkl}} = \prod_{d \in x, y,\ldots} \mathcal{I}(i_d,j_d,k_d,l_d) \delta_{i_d + j_d}^{k_d + l_d},\]

where the $\delta$-function indicates that the noninteracting energy is conserved along each dimension. The integral $\mathcal{I}(a,b,c,d)$ is of four one dimensional harmonic oscillator basis functions, see four_oscillator_integral_general, with the additional restriction that energy is conserved in each dimension.

Arguments

  • addr: the starting address, defines number of particles and total number of modes.
  • S: Tuple of the number of levels in each dimension, including the groundstate. Defaults to a 1D spectrum with number of levels matching modes of addr. Will be sorted to make the first dimension the largest.
  • η: Define a custom aspect ratio for the trapping potential strengths, instead of deriving from S .- 1. The values are always scaled relative to the first dimension, which sets the energy scale of the system, $\hbar\omega_x$.
  • g: the (isotropic) interparticle interaction parameter. The value of g is assumed to be in trap units.
  • interaction_only: if set to true then the noninteracting single-particle terms are ignored. Useful if only energy shifts due to interactions are required.
source
Rimu.Hamiltonians.HOCartesianCentralImpurity — Type
HOCartesianCentralImpurity(addr; kwargs...)

Hamiltonian of non-interacting particles in an arbitrary harmonic trap with a delta-function potential at the centre, with strength g,

\[\hat{H}_\mathrm{rel} = \sum_\mathbf{i} ϵ_\mathbf{i} n_\mathbf{i} + g\sum_\mathbf{ij} V_\mathbf{ij} a^†_\mathbf{i} a_\mathbf{j}.\]

For a $D$-dimensional harmonic oscillator indices $\mathbf{i}, \mathbf{j}, \ldots$ are $D$-tuples. The energy scale is defined by the first dimension i.e. $\hbar \omega_x$ so that single particle energies are

\[ \frac{\epsilon_\mathbf{i}}{\hbar \omega_x} = (i_x + 1/2) + \eta_y (i_y+1/2) + \ldots.\]

The factors $\eta_y, \ldots$ allow for anisotropic trapping geometries and are assumed to be greater than 1 so that $\omega_x$ is the smallest trapping frequency.

Matrix elements $V_{\mathbf{ij}}$ are for a delta function potential calculated in this basis

\[ V_{\mathbf{ij}} = \prod_{d \in x, y,\ldots} \psi_{i_d}(0) \psi_{j_d}(0).\]

Only even parity states feel this impurity, so all $i_d$ are even. Note that the matrix representation of this Hamiltonian for a single particle is completely dense in the even-parity subspace.

Arguments

  • addr: the starting address, defines number of particles and total number of modes.
  • max_nx = num_modes(addr) - 1: the maximum harmonic oscillator index number in the $x$-dimension. Must be even. Index number for the harmonic oscillator groundstate is 0.
  • ηs = (): a tuple of aspect ratios for the remaining dimensions (η_y, ...). Should be empty for a 1D trap or contain values greater than 1.0. The maximum index in other dimensions will be the largest even number less than M/η_y.
  • S = nothing: Instead of max_nx, manually set the number of levels in each dimension, including the groundstate. Must be a Tuple of Ints.
  • g = 1.0: the strength of the delta impurity in ($x$-dimension) trap units.
  • impurity_only=false: if set to true then the trap energy terms are ignored. Useful if only energy shifts due to the impurity are required.
Warning
Due to use of `SpecialFunctions` with large arguments the matrix representation of
this Hamiltonian may not be strictly symmetric, but is approximately symmetric within
machine precision.

See also HOCartesianContactInteractions andHOCartesianEnergyConservedPerDim.

source

Other model Hamiltonians

Rimu.Hamiltonians.MatrixHamiltonian — Type
MatrixHamiltonian(
    mat::AbstractMatrix{T};
    starting_address::Int = findmin(real.(diag(mat)))[2]
) <: AbstractHamiltonian{T}

Wrap an abstract matrix mat as an AbstractHamiltonian object. Works with stochastic methods of ProjectorMonteCarloProblem() and DVec. Optionally, a valid index can be provided as the starting_address.

Specialised methods are implemented for sparse matrices of type AbstractSparseMatrixCSC. One based indexing is required for the matrix mat.

source
Rimu.Hamiltonians.Transcorrelated1D — Type
Transcorrelated1D(address; t=1.0, v=1.0, v_ho=0.0, cutoff=1, three_body_term=true)

Implements a transcorrelated Hamiltonian for contact interactions in one dimensional momentum space from Jeszenski et al. (2018). Currently limited to two component fermionic addresses.

\[\begin{aligned} \tilde{H} &= t \sum_{kσ}k^2 n_{k,σ} \\ &\quad + \sum_{pqkσσ'} T_{pqk} a^†_{p-k,σ} a^†_{q+k,σ'} a_{q,σ'} a_{p,σ} \\ &\quad + \sum_{pqskk'σσ'} Q_{kk'}a^†_{p-k,σ} a^†_{q+k,σ} a^†_{s+k-k',σ'} a_{s,σ'} a_{q,σ} a_{p,σ} \\ &\quad + V̂_\mathrm{ho} \end{aligned}\]

where

\[\begin{aligned} \tilde{u}(k) &= \begin{cases} -\frac{2}{k^2} &\mathrm{if\ } |k| ≥ k_c\\ 0 & \mathrm{otherwise} \end{cases} \\ T_{pqk} &= \frac{v}{M} + \frac{2v}{M}\left[k^2\tilde{u}(k) - (p - q)k\tilde{u}(k)\right] + \frac{2v^2}{t}W(k)\\ W(k) &= \frac{1}{M^2}\sum_{q} (k - q)q\, \tilde{u}(q)\,\tilde{u}(k - q) \\ Q_{kl} &= -\frac{v^2}{t M^2}k \tilde{u}(k)\,l\tilde{u}(l), \end{aligned}\]

Arguments

  • address: The starting address, defines number of particles and sites.
  • v: The interaction parameter.
  • t: The kinetic energy prefactor.
  • v_ho: Strength of the external harmonic oscillator potential $V̂_\mathrm{ho}$. See HubbardMom1DEP.
  • cutoff controls $k_c$ in equations above. Note: skipping generating off-diagonal elements below the cutoff is not implemented - zero-valued elements are returned instead.
  • three_body_term: If set to false, generating three body excitations is skipped. Note: when disabling three body terms, cutoff should be set to a higher value for best results.

See also

source
Rimu.Hamiltonians.FroehlichPolaron — Type
FroehlichPolaron(address::OccupationNumberFS{M}; kwargs...) <: AbstractHamiltonian

The Froehlich polaron Hamiltonian for a 1D lattice with M momentum modes is given by

\[H = (p̂_f - p)^2/m + ωN̂ - v Σₖ(âₖ^† + â₋ₖ)\]

where $p$ is the total momentum, $p̂_f = Σ_k k âₖ^† âₖ$ is the momentum operator for the bosons, and $k$ part of the momentum lattice with separation $2π/l$. $N̂$ is the number operator for the bosons.

Keyword Arguments

  • p=0.0: the total momentum $p$.
  • v=1.0: the coupling strength $v$.
  • mass=1.0: the particle mass $m$.
  • omega=1.0: the oscillation frequency of the phonons $ω$.
  • l=1.0: the box size in real space $l$. Provides scale parameter of the momentum lattice.
  • momentum_cutoff=nothing: the maximum boson momentum allowed for an address.
  • mode_cutoff: the maximum number of bosons in each momentum mode. Defaults to the maximum value supported by the address type OccupationNumberFS.

Examples

julia> fs = OccupationNumberFS(0,0,0)
OccupationNumberFS{3, UInt8}(0, 0, 0)

julia> ham = FroehlichPolaron(fs; v=0.5)
FroehlichPolaron(fs"|0 0 0⟩{8}"; v=0.5, mass=1.0, omega=1.0, l=1.0, p=0.0, mode_cutoff=255)

julia> dimension(ham)
16777216

julia> dimension(FroehlichPolaron(fs; v=0.5, mode_cutoff=5))
216

See also OccupationNumberFS, dimension, AbstractHamiltonian.

source

Convenience functions

Rimu.Hamiltonians.momentum — Function
momentum(ham::AbstractHamiltonian)

Momentum as a linear operator in Fock space. Pass a Hamiltonian ham in order to convey information about the Fock basis. Returns an AbstractHamiltonian that represents the momentum operator.

Note: momentum is currently only defined on HubbardMom1D.

Example

julia> add = BoseFS((1, 0, 2, 1, 2, 1, 1, 3));


julia> ham = HubbardMom1D(add; u = 2.0, t = 1.0);


julia> mom = momentum(ham);


julia> diagonal_element(mom, add) # calculate the momentum of a single configuration
-1.5707963267948966

julia> v = DVec(add => 10; capacity=1000);


julia> rayleigh_quotient(mom, v) # momentum expectation value for state vector `v`
-1.5707963267948966

Part of the AbstractHamiltonian interface.

source

Hamiltonian wrappers

The following Hamiltonians are constructed from an existing Hamiltonian instance and change its behaviour:

Rimu.Hamiltonians.GutzwillerSampling — Type
GutzwillerSampling(H::AbstractHamiltonian; g)

Wrapper over any AbstractHamiltonian that implements Gutzwiller sampling. In this importance sampling scheme the Hamiltonian is modified as follows

\[\tilde{H}_{ij} = H_{ij} e^{-g(H_{ii} - H_{jj})} .\]

This way off-diagonal spawns to higher-energy configurations are discouraged and spawns to lower-energy configurations encouraged for positive g while keeping the spectrum of the Hamiltonian intact.

Constructor

  • GutzwillerSampling(::AbstractHamiltonian, g)
  • GutzwillerSampling(::AbstractHamiltonian; g)

After construction, we can access the underlying Hamiltonian with G.hamiltonian and the g parameter with G.g.

Example

julia> H = HubbardMom1D(BoseFS(1,1,1); u=6.0, t=1.0)
HubbardMom1D(fs"|1 1 1⟩"; u=6.0, t=1.0)

julia> G = GutzwillerSampling(H, g=0.3)
GutzwillerSampling(HubbardMom1D(fs"|1 1 1⟩"; u=6.0, t=1.0); g=0.3)

julia> Matrix(H; sort=true)
4×4 Matrix{Float64}:
 9.0      0.0       4.89898  0.0
 0.0      0.0       4.89898  0.0
 4.89898  4.89898  12.0      4.89898
 0.0      0.0       4.89898  9.0

julia> Matrix(G; sort=true)
4×4 Matrix{Float64}:
 9.0      0.0        12.0495  0.0
 0.0      0.0       179.294   0.0
 1.99178  0.133858   12.0     1.99178
 0.0      0.0        12.0495  9.0

julia> eigen(Matrix(H)).values
4-element Vector{Float64}:
 -2.3661456273236645
  4.9594958589580465
  8.999999999999996
 18.406649768365643

julia> eigen(Matrix(G)).values
4-element Vector{Float64}:
 -2.366145627323686
  4.959495858958046
  8.999999999999998
 18.40664976836564

Observables

See AllOverlaps for calculation of observables with a transformed Hamiltonian.

source
Rimu.Hamiltonians.GuidingVectorSampling — Type
GuidingVectorSampling

Wrapper over any AbstractHamiltonian that implements guided vector a.k.a. guided wave function sampling. In this importance sampling scheme the Hamiltonian is modified as follows.

\[\tilde{H}_{ij} = v_i H_{ij} v_j^{-1}\]

and where v is the guiding vector. v_i and v_j are also thresholded to avoid dividing by zero (see below).

Constructors

  • GuidingVectorSampling(::AbstractHamiltonian, vector, eps)
  • GuidingVectorSampling(::AbstractHamiltonian; vector, eps)

eps is a thresholding parameter used to avoid dividing by zero; all values below eps are set to eps. It is recommended that eps is in the same value range as the guiding vector. The default value is set to eps=norm(v, Inf) * 1e-2

After construction, we can access the underlying hamiltonian with G.hamiltonian, the eps parameter with G.eps, and the guiding vector with G.vector.

Example

julia> H = HubbardMom1D(BoseFS(1,1,1); u=6.0, t=1.0);

julia> v = DVec(starting_address(H) => 10);

julia> G = GuidingVectorSampling(H, v, 0.1);

julia> Matrix(H; sort=true)
4×4 Matrix{Float64}:
 9.0      0.0       4.89898  0.0
 0.0      0.0       4.89898  0.0
 4.89898  4.89898  12.0      4.89898
 0.0      0.0       4.89898  9.0

julia> Matrix(G; sort=true)
4×4 Matrix{Float64}:
   9.0      0.0     0.0489898    0.0
   0.0      0.0     0.0489898    0.0
 489.898  489.898  12.0        489.898
   0.0      0.0     0.0489898    9.0

julia> eigen(Matrix(H)).values
4-element Vector{Float64}:
 -2.3661456273236645
  4.9594958589580465
  8.999999999999996
 18.406649768365643

julia> eigen(Matrix(G)).values
4-element Vector{Float64}:
 -2.366145627323689
  4.9594958589580465
  8.999999999999998
 18.406649768365643

Observables

See AllOverlaps for calculation of observables with a transformed Hamiltonian.

source
Rimu.Hamiltonians.ParitySymmetry — Type
ParitySymmetry(ham::AbstractHamiltonian{T}; even=true) <: AbstractHamiltonian{T}

Impose even or odd parity on all states and the Hamiltonian ham as controlled by the keyword argument even. Parity symmetry of the Hamiltonian is assumed. For some Hamiltonians, ParitySymmetry reduces the size of the Hilbert space by half.

ParitySymmetry performs a unitary transformation, leaving the eigenvalues unchanged and preserving the LOStructure. This is achieved by changing the basis set to states with defined parity. Effectively, a non-even address $|α⟩$ is replaced by $\frac{1}{√2}(|α⟩ ± |ᾱ⟩)$ for even and odd parity, respectively, where ᾱ == reverse(α).

Notes

  • This modifier currently only works on starting_addresss with an odd number of modes.
  • For odd parity, the starting_address of the underlying Hamiltonian cannot be symmetric.
  • If parity is not a symmetry of the Hamiltonian ham then the result is undefined.
  • ParitySymmetry works by modifying the offdiagonals iterator.
julia> ham = HubbardReal1D(BoseFS(0,2,1))
HubbardReal1D(fs"|0 2 1⟩"; u=1.0, t=1.0)

julia> size(Matrix(ham))
(10, 10)

julia> size(Matrix(ParitySymmetry(ham)))
(6, 6)

julia> size(Matrix(ParitySymmetry(ham; odd=true)))
(4, 4)

julia> eigvals(Matrix(ham))[1] ≈ eigvals(Matrix(ParitySymmetry(ham)))[1]
true

See also TimeReversalSymmetry.

source
Rimu.Hamiltonians.TimeReversalSymmetry — Type
TimeReversalSymmetry(ham::AbstractHamiltonian{T}; even=true) <: AbstractHamiltonian{T}

Impose even or odd time reversal on all states and the Hamiltonian ham as controlled by the keyword argument even. If time reversal is a symmetry of the Hamiltonian it will block (reducing Hilbert space dimension) preserving the eigenvalues and LOStructure.

Notes

  • This modifier only works two component starting_addresses.
  • For odd time reversal symmetry, the starting_address of the underlying Hamiltonian must not be symmetric.
  • If time reversal is not a symmetry of the Hamiltonian ham then the result is undefined.
  • TimeReversalSymmetry works by modifying the offdiagonals iterator.
julia> ham = HubbardMom1D(FermiFS2C((1,0,1),(0,1,1)));

julia> size(Matrix(ham))
(3, 3)

julia> size(Matrix(TimeReversalSymmetry(ham)))
(2, 2)

julia> size(Matrix(TimeReversalSymmetry(ham, even=false)))
(1, 1)

julia> eigvals(Matrix(TimeReversalSymmetry(ham)))[1] ≈ eigvals(Matrix(ham))[1]
true

See also ParitySymmetry.

source
Rimu.Hamiltonians.Stoquastic — Type
Stoquastic(ham <: AbstractHamiltonian) <: AbstractHamiltonian

A wrapper for an AbstractHamiltonian that replaces all off-diagonal matrix elements v by -abs(v), thus making the new Hamiltonian stoquastic.

A stoquastic Hamiltonian does not have a Monte Carlo sign problem. For a hermitian ham the smallest eigenvalue of Stoquastic(ham) is ≤ the smallest eigenvalue of ham.

source
Rimu.Hamiltonians.TransformUndoer — Type
TransformUndoer(transform::AbstractHamiltonian, op::AbstractObservable) <: AbstractHamiltonian

Create a new operator for the purpose of calculating overlaps of transformed vectors, which are defined by some transformation transform. The new operator should represent the effect of undoing the transformation before calculating overlaps, including with an optional operator op.

Not exported; transformations should define all necessary methods and properties, see AbstractHamiltonian. An ArgumentError is thrown if used with an unsupported transformation.

Example

A similarity transform $\hat{G} = f \hat{H} f^{-1}$ has eigenvector $d = f \cdot c$ where $c$ is an eigenvector of $\hat{H}$. Then the overlap $c' \cdot c = d' \cdot f^{-2} \cdot d$ can be computed by defining all necessary methods for TransformUndoer(G) to represent the operator $f^{-2}$ and calculating dot(d, TransformUndoer(G), d).

Observables in the transformed basis can be computed by defining TransformUndoer(G, A) to represent $f^{-1} A f^{-1}$.

Supported transformations

source
Rimu.Hamiltonians.HamiltonianSum — Type
HamiltonianSum(A::AbstractHamiltonian, B::AbstractHamiltonian; weight=0.5) <: AbstractHamiltonian
add(A::AbstractHamiltonian, B::AbstractHamiltonian, [a=1, b=1]; weight=0.5) -> HamiltonianSum
+(A::AbstractHamiltonian, B::AbstractHamiltonian)

The sum of two AbstractHamiltonians, $A + B$. The two Hamiltonians must act on the same address space. The keyword argument weight affects random spawning with random_offdiagonal and determines the probability of random spawns from A, with 1 - weight the probability of spawning from B.

If coefficients a and b are given, the Hamiltonians are scaled with ScaledHamiltonian, to represent $aA + bB$.

See also ScaledHamiltonian, HamiltonianProduct, AbstractHamiltonian.

source

Observables

Rimu.jl offers two other supertypes for operators that are less restrictive than AbstractHamiltonian. AbstractObservable and AbstractOperators both can represent a physical observable. Their expectation values can be sampled during a ProjectorMonteCarloProblem simulation by passing them into a suitable ReplicaStrategy, e.g. AllOverlaps. Some observables are also AbstractHamiltonians. The full type hierarchy is

AbstractHamiltonian{T} <: AbstractOperator{T} <: AbstractObservable{T}
Rimu.Hamiltonians.ParticleNumberOperator — Type
ParticleNumberOperator() <: AbstractOperator{Float64}

The number operator in Fock space. This operator is diagonal in the Fock basis and returns the number of particles in the Fock state. It works with any address type that is a subtype of AbstractFockAddress.

julia> p = ExactDiagonalizationProblem(FroehlichPolaron(fs"|0 0⟩{}"; mode_cutoff=5, v=3));

julia> gs = solve(p).vectors[1]; # normalised ground state vector

julia> dot(gs, ParticleNumberOperator(), gs) # particle number expectation value
2.8823297252925917

See also AbstractHamiltonian.

source
Rimu.Hamiltonians.G2RealCorrelator — Type
G2RealCorrelator(d::Int) <: AbstractOperator{Float64}

Two-body operator for density-density correlation between sites separated by d with 0 ≤ d < M.

\[ \hat{G}^{(2)}(d) = \frac{1}{M} \sum_i^M \hat{n}_i (\hat{n}_{i+d} - \delta_{0d}).\]

Assumes a one-dimensional lattice with periodic boundary conditions where

\[ \hat{G}^{(2)}(-M/2 \leq d < 0) = \hat{G}^{(2)}(|d|),\]

\[ \hat{G}^{(2)}(M/2 < d < M) = \hat{G}^{(2)}(M - d),\]

and normalisation

\[ \sum_{d=0}^{M-1} \langle \hat{G}^{(2)}(d) \rangle = \frac{N (N-1)}{M}.\]

For multicomponent basis, calculates correlations between all particles equally, equivalent to stacking all components into a single Fock state.

Arguments

  • d::Integer: distance between sites.

See also

source
Rimu.Hamiltonians.G2RealSpace — Type
G2RealSpace(geometry::CubicGrid, σ=1, τ=1; sum_components=false) <: AbstractOperator{SArray}

Two-body operator for density-density correlation for all Displacements $d⃗$ in the specified geometry.

\[ \hat{G}^{(2)}_{σ,τ}(d⃗) = \frac{1}{M} ∑_{i⃗} n̂_{σ,i⃗} (n̂_{τ,i⃗+d⃗} - δ_{0⃗,d⃗}δ_{σ,τ}).\]

For multicomponent addresses, σ and τ control the components involved. Alternatively, sum_components can be set to true, which treats all particles as belonging to the same component.

Examples

julia> geom = CubicGrid(2, 2);

julia> g2 = G2RealSpace(geom)
G2RealSpace(CubicGrid((2, 2), (true, true)), 1,1)

julia> diagonal_element(g2, BoseFS(2,0,1,1))
2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.5  1.0
 0.5  1.0

julia> g2_cross = G2RealSpace(geom, 1, 2)
G2RealSpace(CubicGrid((2, 2), (true, true)), 1,2)

julia> g2_sum = G2RealSpace(geom, sum_components=true)
G2RealSpace(CubicGrid((2, 2), (true, true)); sum_components=true)

julia> diagonal_element(g2, fs"|⇅⋅↓↑⟩")
2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 0.0  0.5

julia> diagonal_element(g2_cross, fs"|⇅⋅↓↑⟩")
2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.25  0.25
 0.25  0.25

julia> diagonal_element(g2_sum, fs"|⇅⋅↓↑⟩")
2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.5  1.0
 0.5  1.0

See also

source
Rimu.Hamiltonians.G2MomCorrelator — Type
G2MomCorrelator(d::Int) <: AbstractOperator{ComplexF64}

Two-body correlation operator representing the density-density correlation at distance d. It returns a Complex value.

Correlation within a single component:

\[\hat{G}^{(2)}(d) = \frac{1}{M}\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} a^†_{p} a_q a_r δ_{s+p,q+r}\]

The diagonal element, where (p-q)=0, is

\[\frac{1}{M}\sum_{k,p=1}^M a^†_{k} b^†_{p} b_p a_k .\]

Arguments

  • d::Integer: the distance between two particles.

See also

source
Rimu.Hamiltonians.SuperfluidCorrelator — Type
SuperfluidCorrelator(d::Int) <: AbstractOperator{Float64}

Operator for extracting superfluid correlation between sites separated by a distance d with 0 ≤ d < M:

\[ \hat{C}_{\text{SF}}(d) = \frac{1}{M} \sum_{i}^{M} a_{i}^{\dagger} a_{i + d}\]

Assumes a one-dimensional lattice with $M$ sites and periodic boundary conditions. $M$ is also the number of modes in the Fock state address.

Usage

Superfluid correlations can be extracted from a Monte Carlo calculation by wrapping SuperfluidCorrelator with AllOverlaps and passing into ProjectorMonteCarloProblem with the replica keyword argument. For an example with a similar use of G2RealCorrelator see G2 Correlator Example.

See also HubbardReal1D, G2RealCorrelator, AbstractOperator, and AllOverlaps.

source
Rimu.Hamiltonians.StringCorrelator — Type
StringCorrelator(d::Int; address=nothing, type=nothing) <: AbstractOperator{T}

Operator for extracting string correlation between lattice sites on a one-dimensional Hubbard lattice separated by a distance d with 0 ≤ d < M

\[ Ĉ_{\text{string}}(d) = \frac{1}{M} \sum_{j}^{M} δ n̂_j (e^{i π \sum_{j ≤ k < j + d} δ n̂_k}) δ n̂_{j+d}\]

Here, $δ n̂_j = n̂_j - n̄$ is the boson number deviation from the mean filling number and $n̄ = N/M$ is the mean filling number of lattice sites with $N$ particles and $M$ lattice sites (or modes).

Assumes a one-dimensional lattice with periodic boundary conditions. For usage see SuperfluidCorrelator and AllOverlaps.

The default element type T is ComplexF64. This can be overridden with the type keyword argument. If an address is provided, then T is calculated from the address type. It is set to ComplexF64 for non-integer filling numbers, and to Float64 for integer filling numbers or if d==0.

See also HubbardReal1D, G2RealCorrelator, SuperfluidCorrelator, AbstractOperator, and AllOverlaps.

source
Rimu.Hamiltonians.ReducedDensityMatrix — Type
ReducedDensityMatrix{T=Float64}(p) <: AbstractObservable{Hermitian{T, Matrix{T}}}

A matrix-valued operator that can be used to calculate the p-particle reduced density matrix. The matrix elements are defined as:

\[\hat{ρ}^{(p)}_{j_1,...,j_1,k_1,...,k_p} = \prod_{i=1}^{p} â^†_{j_i} \prod_{l=p}^{1} â_{k_{l}}\]

The integer indices j_i and k_i represent single particle modes. For efficiency they are chosen to be distinct and ordered:

\[j_1 < j_2 < \ldots < j_{p} \quad \land \quad k_1 < k_2 < \ldots < k_{p}\]

ReducedDensityMatrix can be used to construct the single-particle reduced density matrix (with p == 1) for fermionic and bosonic Fock spaces with address types <: SingleComponentFockAddress. For higher order reduced density matrices with p > 1 only fermionic Fock addresses (FermiFS) are supported due to the ordering of indices.

ReducedDensityMatrix can be used with dot or AllOverlaps to calculate the whole matrix in one go.

Examples

julia> dvec_b = PDVec(BoseFS(1,1) => 0.5, BoseFS(2,0) => 0.5)
2-element PDVec: style = IsDeterministic{Float64}()
  fs"|2 0⟩" => 0.5
  fs"|1 1⟩" => 0.5

julia> Op1 = ReducedDensityMatrix(1)
ReducedDensityMatrix{Float64}(1)

julia> dot(dvec_b, Op1, dvec_b)
2×2 Hermitian{Float64, Matrix{Float64}}:
 0.75      0.353553
 0.353553  0.25

julia> Op2 = ReducedDensityMatrix{Float32}(2)
ReducedDensityMatrix{Float32}(2)

julia> dvec_f = PDVec(FermiFS(1,1,0,0) => 0.5, FermiFS(0,1,1,0) => 0.5)
2-element PDVec: style = IsDeterministic{Float64}()
  fs"|↑↑⋅⋅⟩" => 0.5
  fs"|⋅↑↑⋅⟩" => 0.5

julia> dot(dvec_f, Op2, dvec_f)
6×6 Hermitian{Float32, Matrix{Float32}}:
 0.25  0.0  0.25  0.0  0.0  0.0
 0.0   0.0  0.0   0.0  0.0  0.0
 0.25  0.0  0.25  0.0  0.0  0.0
 0.0   0.0  0.0   0.0  0.0  0.0
 0.0   0.0  0.0   0.0  0.0  0.0
 0.0   0.0  0.0   0.0  0.0  0.0

See also single_particle_density, SingleParticleDensity, SingleParticleExcitation, TwoParticleExcitation.

source
Rimu.Hamiltonians.Momentum — Type
Momentum(component=0; fold=true) <: AbstractHamiltonian

The momentum operator $P̂$.

The component argument controls which component of the address is taken into consideration. A value of 0 sums the contributions of all components. If fold is true, the momentum is folded into the Brillouin zone.

julia> address = BoseFS((1, 0, 2, 1, 2, 1, 1, 3))
BoseFS{11,8}(1, 0, 2, 1, 2, 1, 1, 3)

julia> v = DVec(address => 10);

julia> rayleigh_quotient(Momentum(), DVec(address => 1))
-2.0

julia> rayleigh_quotient(Momentum(fold=false), DVec(address => 1))
14.0
source
Rimu.Hamiltonians.AxialAngularMomentumHO — Type
AxialAngularMomentumHO(S; z_dim = 3, addr = BoseFS(prod(S))) <: AbstractHamiltonian

Angular momentum operator for application to Cartesian harmonic oscillator basis, see HOCartesianContactInteractions or HOCartesianEnergyConservedPerDim. Represents the projection of angular momentum onto z-axis:

\[\hat{L}_z = i \hbar \sum_{j=1}^N \left( b_x b_y^\dag - b_y b_x^\dag \right),\]

where $b_x^\dag$ and $b_x$ are raising and lowering (ladder) operators for a harmonic oscillator in the $x$ dimension, and simlarly for $y$.

This is implemented for an $N$ particle Fock space with creation and annihilation operators as

\[\frac{1}{\hbar} \hat{L}_z = i \sum_{n_x=1}^{M_x} \sum_{n_y=1}^{M_y} \left( a_{n_x-1,n_y+1}^\dag - a_{n_x+1,n_y-1}^\dag \right) a_{n_x, n_y}.\]

in units of $\hbar$.

Argument S is a tuple defining the range of Cartesian modes in each dimension and their mapping to Fock space modes in a SingleComponentFockAddress. If S indicates a 3D system the z dimension can be changed by setting z_dim; S should be be isotropic in the remaining x-y plane, i.e. must have S[x_dim] == S[y_dim]. The starting address addr only needs to satisfy num_modes(addr) == prod(S).

Example

Calculate the overlap of two Fock addresses interpreted as harmonic oscillator states in a 2D Cartesian basis

julia> S = (2,2)
(2, 2)

julia> Lz = AxialAngularMomentumHO(S)
AxialAngularMomentumHO((2, 2); z_dim = 3, addr = BoseFS{0,4}(0, 0, 0, 0))

julia> v = DVec(BoseFS(prod(S), 2 => 1) => 1.0)
DVec{BoseFS{1, 4, BitString{4, 1, UInt8}},Float64} with 1 entry, style = IsDeterministic{Float64}()
  fs"|0 1 0 0⟩" => 1.0

julia> w = DVec(BoseFS(prod(S), 3 => 1) => 1.0)
DVec{BoseFS{1, 4, BitString{4, 1, UInt8}},Float64} with 1 entry, style = IsDeterministic{Float64}()
  fs"|0 0 1 0⟩" => 1.0

julia> dot(w, Lz, v)
0.0 + 1.0im
source

Geometry

Lattices in higher dimensions are defined here and can be passed with the keyword argument geometry to HubbardRealSpace and G2RealSpace.

Rimu.Hamiltonians.CubicGrid — Type
CubicGrid(dims::NTuple{D,Int}, fold::NTuple{D,Bool})

Represents a D-dimensional grid. Used to define a cubic lattice and boundary conditions for some AbstractHamiltonians, e.g. with the keyword argument geometry when constructing a HubbardRealSpace. The type instance can be used to convert between cartesian vector indices (tuples or SVectors) and linear indices (integers). When indexed with vectors, it folds them back into the grid if the out-of-bounds dimension is periodic and 0 otherwise (see example below).

  • dims controls the size of the grid in each dimension.
  • fold controls whether the boundaries in each dimension are periodic (or folded in the case of momentum space).
julia> geo = CubicGrid((2,3), (true,false))
CubicGrid{2}((2, 3), (true, false))

julia> geo[1]
(1, 1)

julia> geo[2]
(2, 1)

julia> geo[3]
(1, 2)

julia> geo[(1,2)]
3

julia> geo[(3,2)] # 3 is folded back into 1
3

julia> geo[(3,3)]
5

julia> geo[(3,4)] # returns 0 if out of bounds
0

See also PeriodicBoundaries, HardwallBoundaries and LadderBoundaries for special-case constructors. See also HubbardRealSpace and G2RealSpace.

source
Rimu.Hamiltonians.Directions — Type
Directions(D) <: AbstractVector{SVector{D,Int}}
Directions(geometry::CubicGrid) <: AbstractVector{SVector{D,Int}}

Iterate over axis-aligned direction vectors in D dimensions.

julia> Directions(3)
6-element Directions{3}:
 [1, 0, 0]
 [0, 1, 0]
 [0, 0, 1]
 [-1, 0, 0]
 [0, -1, 0]
 [0, 0, -1]

See also CubicGrid.

source
Rimu.Hamiltonians.Displacements — Type
Displacements(geometry::CubicGrid) <: AbstractVector{SVector{D,Int}}

Return all valid offset vectors in a CubicGrid. If center=true the (0,0) displacement is placed at the centre of the array.

julia> geometry = CubicGrid((3,4));

julia> reshape(Displacements(geometry), (3,4))
3×4 reshape(::Displacements{2, CubicGrid{2, (3, 4), (true, true)}}, 3, 4) with eltype StaticArraysCore.SVector{2, Int64}:
 [0, 0]  [0, 1]  [0, 2]  [0, 3]
 [1, 0]  [1, 1]  [1, 2]  [1, 3]
 [2, 0]  [2, 1]  [2, 2]  [2, 3]

julia> reshape(Displacements(geometry; center=true), (3,4))
3×4 reshape(::Displacements{2, CubicGrid{2, (3, 4), (true, true)}}, 3, 4) with eltype StaticArraysCore.SVector{2, Int64}:
 [-1, -1]  [-1, 0]  [-1, 1]  [-1, 2]
 [0, -1]   [0, 0]   [0, 1]   [0, 2]
 [1, -1]   [1, 0]   [1, 1]   [1, 2]
source
Rimu.Hamiltonians.LadderBoundaries — Function
LadderBoundaries(dims...) -> CubicGrid
LadderBoundaries(dims) -> CubicGrid

Return a CubicGrid where the first dimension is dimensions non-periodic and the rest are periodic. Equivalent to CubicGrid(dims, (true, false, ...)).

source

Index