creating fermionic creation and annihilation operators

doc/changes/2166.feature

@@ -0,0 +1,4 @@
+added fermionic annihilation and creation operators.
+
+Closely followed the protocol proposed in the following `source
+<https://github.com/qutip/qutip/issues/863>`_.

qutip/core/operators.py

@@ -5,11 +5,11 @@
 
 __all__ = ['jmat', 'spin_Jx', 'spin_Jy', 'spin_Jz', 'spin_Jm', 'spin_Jp',
            'spin_J_set', 'sigmap', 'sigmam', 'sigmax', 'sigmay', 'sigmaz',
-           'destroy', 'create', 'qeye', 'qeye_like', 'identity', 'position',
-           'momentum', 'num', 'squeeze', 'squeezing', 'swap', 'displace',
-           'commutator', 'qutrit_ops', 'qdiags', 'phase', 'qzero',
-           'qzero_like', 'enr_destroy', 'enr_identity', 'charge', 'tunneling',
-           'qft']
+           'destroy', 'create', 'fdestroy', 'fcreate', 'qeye', 'qeye_like',
+           'identity', 'position', 'momentum', 'num', 'squeeze', 'squeezing',
+           'swap', 'displace', 'commutator', 'qutrit_ops', 'qdiags', 'phase',
+           'qzero', 'qzero_like', 'enr_destroy', 'enr_identity', 'charge',
+           'tunneling', 'qft']
 
 import numbers
 
@@ -135,11 +135,11 @@ def jmat(j, which=None, *, dtype=None):
         return Qobj(_jplus(j, dtype=dtype).adjoint(), dims=dims, type='oper',
                     isherm=False, isunitary=False, copy=False)
     if which == 'x':
-        A =  _jplus(j, dtype=dtype)
+        A = _jplus(j, dtype=dtype)
         return Qobj(_data.add(A, A.adjoint()), dims=dims, type='oper',
                     isherm=True, isunitary=False, copy=False) * 0.5
     if which == 'y':
-        A =  _data.mul(_jplus(j, dtype=dtype), -0.5j)
+        A = _data.mul(_jplus(j, dtype=dtype), -0.5j)
         return Qobj(_data.add(A, A.adjoint()), dims=dims, type='oper',
                     isherm=True, isunitary=False, copy=False)
     if which == 'z':
@@ -471,6 +471,138 @@ def create(N, offset=0, *, dtype=None):
     return qdiags(data, -1, dtype=dtype)
 
 
+def fdestroy(n_sites, site, dtype=None):
+    """
+    Fermionic destruction operator.
+    We use the Jordan-Wigner transformation,
+    making use of the Jordan-Wigner ZZ..Z strings,
+    to construct this as follows (in Latex):
+    a_j = sigma_z^{otimes j}
+          otimes (frac{sigma_x + i sigma_y}{2})
+          otimes I^{otimes N-j-1}
+
+    Parameters
+    ----------
+    n_sites : int
+        Number of sites in Fock space.
+
+    site : int (default 0)
+        The site in Fock space to add a fermion to.
+        Corresponds to j in the above JW transform.
+
+    Returns
+    -------
+    oper : qobj
+        Qobj for destruction operator.
+
+    Examples
+    --------
+    >>> fdestroy(2) # doctest: +SKIP
+    Quantum object: dims=[[2 2], [2 2]], shape=(4, 4), \
+    type='oper', isherm=False
+    Qobj data =
+    [[0. 0. 1. 0.]
+    [0. 0. 0. 1.]
+    [0. 0. 0. 0.]
+    [0. 0. 0. 0.]]
+    """
+    return _f_op(n_sites, site, 'destruction', dtype=dtype)
+
+
+def fcreate(n_sites, site, dtype=None):
+    """
+    Fermionic creation operator.
+    We use the Jordan-Wigner transformation,
+    making use of the Jordan-Wigner ZZ..Z strings,
+    to construct this as follows (in Latex):
+    a_j = sigma_z^{otimes j}
+          otimes (frac{sigma_x - i sigma_y}{2})
+          otimes I^{otimes N-j-1}
+
+
+    Parameters
+    ----------
+    n_sites : int
+        Number of sites in Fock space.
+
+    site : int
+        The site in Fock space to add a fermion to.
+        Corresponds to j in the above JW transform.
+
+    Returns
+    -------
+    oper : qobj
+        Qobj for raising operator.
+
+    Examples
+    --------
+    >>> fcreate(2) # doctest: +SKIP
+    Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), \
+    type = oper, isherm = False
+    Qobj data =
+    [[0. 0. 0. 0.]
+    [0. 0. 0. 0.]
+    [1. 0. 0. 0.]
+    [0. 1. 0. 0.]]
+    """
+    return _f_op(n_sites, site, 'creation', dtype=dtype)
+
+
+def _f_op(n_sites, site, action, dtype=None):
+    """ Makes fermionic creation and destruction operators.
+    We use the Jordan-Wigner transformation,
+    making use of the Jordan-Wigner ZZ..Z strings,
+    to construct this as follows (in Latex):
+    a_j = sigma_z^{otimes j}
+          otimes (frac{sigma_x pm i sigma_y}{2})
+          otimes I^{otimes N-j-1}
+
+    Parameters
+    ----------
+    action : str
+        The type of operator to build.
+        Can only be 'creation' or 'destruction'
+
+    n_sites : int
+        Number of sites in Fock space.
+
+    site : int
+        The site in Fock space to create/destroy a fermion on.
+        Corresponds to j in the above JW transform.
+
+    Returns
+    -------
+    oper : qobj
+        Qobj for destruction operator.
+    """
+    # get `tensor` and sigma z objects
+    from .tensor import tensor
+    s_z = 2 * jmat(0.5, 'z', dtype=dtype)
+
+    # sanity check
+    if site >= n_sites:
+        raise ValueError(f'The specified site {site} is not in \
+                         the range of {n_sites} sites.')
+
+    # figure out which operator to build
+    if action.lower() == 'creation':
+        operator = create(2, dtype=dtype)
+    elif action.lower() == 'destruction':
+        operator = destroy(2, dtype=dtype)
+    else:
+        raise TypeError("Unknown operator '%s'. `action` must be \
+                        either 'creation' or 'destruction.'" % action)
+
+    # build operator
+    if site == 0:
+        return tensor([operator, *([identity(2, dtype=dtype)]*(n_sites-1))])
+    elif n_sites-site-1 > 0:
+        return tensor([*([s_z]*site), operator,
+                       *([identity(2, dtype=dtype)]*(n_sites-site-1))])
+    else:
+        return tensor([*([s_z]*site), operator])
+
+
 def _implicit_tensor_dimensions(dimensions):
     """
     Total flattened size and operator dimensions for operator creation routines

qutip/tests/core/test_operators.py

@@ -251,6 +251,8 @@ def _id_func(val):
     (qutip.spin_Jp, (1,)),
     (qutip.destroy, (5,)),
     (qutip.create, (5,)),
+    (qutip.fdestroy, (5,0)),
+    (qutip.fcreate, (5,0)),
     (qutip.qzero, (5,)),
     (qutip.qeye, (5,)),
     (qutip.position, (5,)),
@@ -342,3 +344,24 @@ def test_qzero_like(dims, superrep, dtype):
     opevo = qutip.QobjEvo(op)
     new = qutip.qzero_like(op)
     assert new == expected
+
+
+@pytest.mark.parametrize('n_sites', [2, 3, 4, 5])
+def test_fcreate_fdestroy(n_sites):
+    identity = qutip.tensor([*([qutip.identity(2)]*(n_sites))])
+    zero_tensor = 0*identity
+    for site_0 in range(n_sites):
+        c_0 = qutip.fcreate(n_sites, site_0)
+        d_0 = qutip.fdestroy(n_sites, site_0)
+        for site_1 in range(n_sites):
+            c_1 = qutip.fcreate(n_sites, site_1)
+            d_1 = qutip.fdestroy(n_sites, site_1)
+            assert qutip.commutator(c_0, c_1, 'anti') == zero_tensor
+            assert qutip.commutator(d_0, d_1, 'anti') == zero_tensor
+            if site_0 == site_1:
+                assert qutip.commutator(c_0, d_1, 'anti') == identity
+                assert qutip.commutator(c_1, d_0, 'anti') == identity
+            else:
+                assert qutip.commutator(c_0, d_1, 'anti') == zero_tensor
+                assert qutip.commutator(c_1, d_0, 'anti') == zero_tensor
+    assert qutip.commutator(identity, c_0) == zero_tensor