Bosonic systems, like Fermionic systems, are expressed using the bosonic creation and annihilation operators $b^\dagger_k$ and $b_k$ respectively. Unlike fermions, however, which satisfy the Pauli exclusion principle and thus are distinguished by the canonical fermionic anticommutation relations, the bosonic ladder operators instead satisfy a set of commutation relations:
$$ \begin{align} & [b_i^\dagger, b_j^\dagger] = 0, ~~~ [b_i, b_j] = 0, ~~~ [b_i, b^\dagger_j] = \delta_{ij} \end{align} $$Any weighted sums of products of these operators are represented with the BosonOperator data structure in OpenFermion. Similarly to when we introduced the FermionOperator, the following are examples of valid BosonOperators:
The BosonOperator class is contained in ops/_boson_operators.py. The BosonOperator is derived from the SymbolicOperator, the same class that derives the FermionOperator. As such, the details of the class implementation are identical - as in the fermion case, the class is implemented as hash table (python dictionary). The keys of the dictionary encode the strings of ladder operators and values of the dictionary store the coefficients - the strings are subsequently encoded as a tuple of 2-tuples which we refer to as the "terms tuple".
Each ladder operator is represented by a 2-tuple. The first element of the 2-tuple is an int indicating the quantum mode on which the ladder operator acts. The second element of the 2-tuple is Boole: 1 represents raising and 0 represents lowering. For instance, $b^\dagger_8$ is represented in a 2-tuple as $(8, 1)$. Note that indices start at 0 and the identity operator is an empty list.
$$ \begin{align} I & \mapsto () \nonumber \\ b_1 & \mapsto ((1, 0),) \nonumber \\ b^\dagger_3 & \mapsto ((3, 1),) \nonumber \\ b^\dagger_3 b_1 & \mapsto ((3, 1), (1, 0)) \nonumber \\ b^\dagger_3 b^\dagger_4 b_1 b_9 & \mapsto ((3, 1), (4, 1), (1, 0), (9, 0)) \nonumber \end{align} $$Alternatively, the BosonOperator supports the string-based syntax introduced in the FermionOperator; in this case, the terms are separated by spaces, with the integer corresponding to the quantum mode the operator acts on, and '^' indicating the Hermitian conjugate:
Let's initialize our first term! We do it two different ways below.
In [1]:
from openfermion.ops import BosonOperator
my_term = BosonOperator(((3, 1), (5, 0), (4, 1), (1, 0)))
print(my_term)
my_term = BosonOperator('3^ 5 4^ 1')
print(my_term)
Note the printed order differs from the code, since bosonic operators of different indices commute past each other.
The preferred way to specify the coefficient in openfermion is to provide an optional coefficient argument. If not provided, the coefficient defaults to 1. In the code below, the first method is preferred. The multiplication in the second method actually creates a copy of the term, which introduces some additional cost. All inplace operands (such as +=) modify classes whereas binary operands such as + create copies.
The additive and multiplicative identities can also be created:
BosonOperator(()) and BosonOperator('') initialises the identity (BosonOperator.identity()).BosonOperator() and BosonOperator() initialises the zero operator (BosonOperator.zero()).
In [2]:
good_way_to_initialize = BosonOperator('3^ 1', -1.7)
print(good_way_to_initialize)
bad_way_to_initialize = -1.7 * BosonOperator('3^ 1')
print(bad_way_to_initialize)
identity = BosonOperator('')
print(identity == BosonOperator.identity())
print(identity)
zero_operator = BosonOperator()
print(zero_operator == BosonOperator.zero())
print(zero_operator)
Note that BosonOperator has only one attribute: .terms. This attribute is the dictionary which stores the term tuples.
In [3]:
my_operator = BosonOperator('4^ 1^ 3 9', 1. + 2.j)
print(my_operator)
print(my_operator.terms)
There are various functions and methods that act on the BosonOperator; these include the ability to normal order, double check if the operator is Hermitian, and calculate the Hermitian conjugate.
In [4]:
from openfermion.utils import hermitian_conjugated, is_hermitian, normal_ordered
normal_ordered_boson applies the bosonic commutation relations to write the operator using only normal-ordered terms; that is, that all creation operators are to the left of annihilation operators:
In [5]:
H = BosonOperator('0 0^', 1. + 2.j)
H.is_normal_ordered()
Out[5]:
In [6]:
normal_ordered(BosonOperator('0 0^', 1. + 2.j))
Out[6]:
We can also use a boson operator method to check if the operator converves the particle number - that is, for each qumode, the number of annihilation operators equals the number of creation operators.
In [7]:
H.is_boson_preserving()
Out[7]:
In [8]:
H = BosonOperator('0 0^ 1^ ', 1. + 2.j)
H.is_boson_preserving()
Out[8]:
The Hermitian conjugated function returns the Hermitian conjugate of the operator, and its hermiticity can be checked using is_hermitian:
In [9]:
is_hermitian(H)
Out[9]:
In [10]:
hermitian_conjugated(H)
Out[10]:
In [11]:
H = BosonOperator('0 1^', 1/2.)
H += BosonOperator('1 0^', 1/2.)
print(is_hermitian(H))
print(hermitian_conjugated(H))
Using the bosonic ladder operators, it is common to define the canonical position and momentum operators $\hat{q}$ and $\hat{p}$:
$$ \hat{q}_i = \sqrt{\frac{\hbar}{2}}(\hat{b}_i+\hat{b}^\dagger_i), ~~~ \hat{p}_i = -i\sqrt{\frac{\hbar}{2}}(\hat{b}_i-\hat{b}^\dagger_i)$$These operators are Hermitian, and are referred to as the phase space quadrature operators. They satisfy the canonical commutation relation
$$ [\hat{q}_i, \hat{p}_j] = \delta_{ij}i\hbar$$where the value of $\hbar$ depends on convention, often taking values $\hbar=0.5$, $1$, or $2$.
In OpenFermion, the quadrature operators are represented by the QuadOperator class, and stored as a dictionary of tuples (as keys) and coefficients (as values). For example, the multi-mode quadrature operator $q_0 p_1 q_3$ is represented internally as ((0, 'q'), (1, 'p'), (3, 'q')). Alternatively, QuadOperators also support string input - using string input, the same operator is described by 'q0 p1 q3'.
In [12]:
from openfermion.ops import QuadOperator
H = QuadOperator('q0 p1 q3')
print(H)
print(H.terms)
H2 = QuadOperator('q3 p4', 3.17)
H2 -= 77. * H
print('')
print(H2)
Note that quadrature operators of different indices commute; as such, like the BosonOperator, by default we sort quadrature operators such that the operators acting on the lowest numbered mode appear to the left.
Like the BosonOperator, there are various functions and methods that act on the QuadOperator; these include the ability to normal order, double check if the operator is Hermitian, and calculate the Hermitian conjugate.
In [13]:
from openfermion.utils import hermitian_conjugated, is_hermitian
normal_ordered_quad is an arbitrary convention chosen in OpenFermion that allows us to compare two quadrature operators that might be equivalent, but written in different forms. It is simply defined as a quadrature operator that has all of the position operators $\hat{q}$ to the left of the momentum operators $\hat{q}$. All quadrature operators can be placed in this 'normal form' by making use of the canonical commutation relation.
In [14]:
H = QuadOperator('p0 q0', 1. + 2.j)
H.is_normal_ordered()
Out[14]:
In [15]:
normal_ordered(H)
Out[15]:
By default, we assume the value $\hbar=1$ in the canonical commutation relation, but this can be modified by passing the hbar keyword argument to the function:
In [16]:
normal_ordered(H, hbar=2)
Out[16]:
We can also use a quad operator method to check if the operator is Gaussian - that is, all terms in the quad operator are of quadratic order or lower:
In [17]:
H = QuadOperator('p0 q0', 1. + 2.j)
H.is_gaussian()
Out[17]:
In [18]:
H = QuadOperator('p0 q0 q1', 1. + 2.j)
H.is_gaussian()
Out[18]:
The Hermitian conjugated function returns the Hermitian conjugate of the operator, and its hermiticity can be checked using is_hermitian:
In [19]:
H = QuadOperator('p0 q1 p1', 1-2j)
hermitian_conjugated(H)
Out[19]:
In [20]:
H = QuadOperator('p0 q0', 1/2.)
H += QuadOperator('q0 p0', -1/2.)
print(is_hermitian(H))
print(hermitian_conjugated(H))
In [21]:
H = QuadOperator('p0 q0', 1/2.)
H += QuadOperator('q0 p0', 1/2.)
print(is_hermitian(H))
print(hermitian_conjugated(H))
In [22]:
hermitian_conjugated(H)
Out[22]:
Converting between bosonic ladder operators and quadrature operators is simple - we just apply the definition of the $\hat{q}$ and $\hat{p}$ operators in terms of $\hat{b}$ and $\hat{b}^\dagger$. Two functions are provided to do this automatically; get_quad_operator and get_boson_operator:
In [23]:
from openfermion.transforms import get_boson_operator, get_quad_operator
In [24]:
H = QuadOperator('p0 q0', 1/2.)
H += QuadOperator('q0 p0', 1/2.)
H
Out[24]:
In [25]:
get_boson_operator(H)
Out[25]:
Note that, since these conversions are dependent on the value of $\hbar$ chosen, both accept a hbar keyword argument. As before, if not specified, the default value of $\hbar$ is hbar=1.
In [26]:
H = BosonOperator('0 0^')
normal_ordered(get_quad_operator(H, hbar=0.5), hbar=0.5)
Out[26]:
We also provide support for the Weyl quantization - this maps a polynomial function of the form
$$f(q_0,\dots,q_{N-1},p_0\dots,p_{N-1})=q_0^{m_0}\cdots q_{N-1}^{m_{N-1}} p_0^{m_0}\cdots p_{N-1}^{m_{N-1}}$$on the phase space to the corresponding combination of quadrature operators $\hat{q}$ and $\hat{p}$. To do so, we make use of the McCoy formula,
$$q^m p^n \rightarrow \frac{1}{2^n} \sum_{r=0}^{n} \binom{n}{r} q^r p^m q^{n-r}.$$
In [27]:
from openfermion.transforms import weyl_polynomial_quantization, symmetric_ordering
For weyl_polynomial_quantization, the polynomial function in the phase space is provided in the form of a string, where 'q' or 'p' is the phase space quadrature variable, the integer directly following is the mode it is with respect to, and '^2' is the polynomial power. If the power is not provided, it is assumed to be '^1'.
In [28]:
weyl_polynomial_quantization('q0 p0')
Out[28]:
In [29]:
weyl_polynomial_quantization('q0^2 p0^3 q1^3')
Out[29]:
McCoy's formula is also used to provide a function that returns the symmetric ordering of a BosonOperator or QuadOperator, $S(\hat{O})$. Note that $S(\hat{O})\neq \hat{O}$:
In [30]:
symmetric_ordering(QuadOperator('q0 p0'))
Out[30]:
Consider the symmetric ordering of the square of the bosonic number operator, $\hat{n} = \hat{b}^\dagger \hat{b}$:
In [31]:
from openfermion.utils import number_operator
n2 = number_operator(1, parity=1) * number_operator(1, parity=1)
In [32]:
n2
Out[32]:
In [33]:
Sn2 = symmetric_ordering(n2)
Sn2
Out[33]:
We can use normal_ordered_boson to simplify this result:
In [34]:
Sn2 = normal_ordered(Sn2)
Sn2
Out[34]:
Therefore $S(\hat{n}) = \hat{b}^\dagger \hat{b}^\dagger \hat{b}\hat{b} + 2\hat{b}^\dagger \hat{b} + 0.5$. This is equivalent to $\hat{n}^2+\hat{n}+0.5$:
In [35]:
Sn2 == normal_ordered(n2 + number_operator(1, parity=1) + 0.5*BosonOperator.identity())
Out[35]:
In addition to the bosonic operators discussed above, we also provide Bosonic Hamiltonians that describe specific models. The Bose-Hubbard Hamiltonian over a discrete lattice or grid described by nodes $V=\{0,1,\dots,N-1\}$ is described by:
$$H = - t \sum_{\langle i, j \rangle} b_i^\dagger b_{j + 1} + \frac{U}{2} \sum_{k=1}^{N-1} b_k^\dagger b_k (b_k^\dagger b_k - 1) - \mu \sum_{k=1}^N b_k^\dagger b_k + V \sum_{\langle i, j \rangle} b_i^\dagger b_i b_j^\dagger b_j.$$where
The Bose-Hubbard Hamiltonian function provided in OpenFermion models a Bose-Hubbard model on a two-dimensional grid, with dimensions given by [x_dimension, y_dimension]. It has the form
bose_hubbard(x_dimension, y_dimension, tunneling, interaction,
chemical_potential=0., dipole=0., periodic=True)
where
x_dimension (int): The width of the grid.y_dimension (int): The height of the grid.tunneling (float): The tunneling amplitude $t$.interaction (float): The attractive local interaction $U$.chemical_potential (float, optional): The chemical potential $\mu$ at each site. Default value is 0.periodic (bool, optional): If True, add periodic boundary conditions. Default is True.dipole (float): The attractive dipole interaction strength $V$.Below is an example of a Bose-Hubbard Hamiltonian constructed in OpenFermion.
In [36]:
from openfermion.hamiltonians import bose_hubbard, fermi_hubbard
bose_hubbard(2, 2, 1, 1)
Out[36]:
Like the fermionic operators, OpenFermion contains the capability to represent bosonic operators as a sparse matrix (sparse.csc_matrix). However, as the fermionic operators can be represented as finite matrices, this is not the case of bosonic systems, as they inhabit a infinite-dimensional Fock space. Instead, a integer truncation value $N$ need to be provided - the returned sparse operator will be of size $N^{M}\times N^{M}$, where $M$ is the number of modes in the system, and acts on the truncated Fock basis $\{\left|{0}\right\rangle, \left|{1}\right\rangle, \dots, \left|{N-1}\right\rangle\}$.
In [37]:
from openfermion.utils import boson_operator_sparse
The function boson_operator_sparse acts on both BosonOperators and QuadOperators:
In [38]:
H = boson_operator_sparse(BosonOperator('0^ 0'), 5)
In [39]:
H.toarray()
Out[39]:
In [40]:
H = boson_operator_sparse(QuadOperator('q0'), 5, hbar=1)
H.toarray()
Out[40]: