# Load packages


```python
import numpy as np
from sympy import *
from sympsi import *

from sympy.physics.secondquant import F
from sympy.physics.quantum.dagger import Dagger
from sympy.physics.quantum import TensorProduct
from sympy.physics.quantum import Operator
from sympy.matrices import Matrix, banded

# This is required to get around a bug in sympy
Matrix.adjoint = lambda self: self.T.applyfunc(Dagger) 

```


```python
# We start with creating the full pairing Hamiltonian for N sites
N = 4

a = symbols(f"a_:{N}")
a = [Operator(ai) for ai in a]
```


```python
# We construct the vector containing the annihilation and creation operators for fermions on each site
# This follows the explanation at https://topocondmat.org/w1_topointro/1D.html
c = Matrix(a + [Dagger(ai) for ai in a])
c
```




$\displaystyle \left[\begin{matrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{0}^{\dagger}\\a_{1}^{\dagger}\\a_{2}^{\dagger}\\a_{3}^{\dagger}\end{matrix}\right]$




```python
# Check whether taking the hermitin conjugate works properly
Dagger(c)
```




$\displaystyle \left[\begin{matrix}a_{0}^{\dagger} & a_{1}^{\dagger} & a_{2}^{\dagger} & a_{3}^{\dagger} & a_{0} & a_{1} & a_{2} & a_{3}\end{matrix}\right]$




```python
# We create the Bogoliubov Hamiltonian matrix
epsilon = symbols(f"epsilon_:{N}")
delta = symbols('Delta')
E = banded({0: epsilon})
C = banded({-1: (N-1)*[-delta], 1: (N-1)*[delta]})

G = Matrix([[E,C], [-C,-E]])
G
```




$\displaystyle \left[\begin{matrix}\epsilon_{0} & 0 & 0 & 0 & 0 & \Delta & 0 & 0\\0 & \epsilon_{1} & 0 & 0 & - \Delta & 0 & \Delta & 0\\0 & 0 & \epsilon_{2} & 0 & 0 & - \Delta & 0 & \Delta\\0 & 0 & 0 & \epsilon_{3} & 0 & 0 & - \Delta & 0\\0 & - \Delta & 0 & 0 & - \epsilon_{0} & 0 & 0 & 0\\\Delta & 0 & - \Delta & 0 & 0 & - \epsilon_{1} & 0 & 0\\0 & \Delta & 0 & - \Delta & 0 & 0 & - \epsilon_{2} & 0\\0 & 0 & \Delta & 0 & 0 & 0 & 0 & - \epsilon_{3}\end{matrix}\right]$




```python
# Next we show that the Bogoliubov matrix is indeed like the original matrix
H = Dagger(c)*G*c
h1 = H[0]
h1
```




$\displaystyle \left(- \Delta a_{1}^{\dagger} - \epsilon_{0} a_{0}\right) a_{0}^{\dagger} + \left(\Delta a_{2}^{\dagger} - \epsilon_{3} a_{3}\right) a_{3}^{\dagger} + \left(\Delta a_{1} + \epsilon_{0} a_{0}^{\dagger}\right) a_{0} + \left(- \Delta a_{2} + \epsilon_{3} a_{3}^{\dagger}\right) a_{3} + \left(\Delta a_{0}^{\dagger} - \Delta a_{2}^{\dagger} - \epsilon_{1} a_{1}\right) a_{1}^{\dagger} + \left(\Delta a_{1}^{\dagger} - \Delta a_{3}^{\dagger} - \epsilon_{2} a_{2}\right) a_{2}^{\dagger} + \left(- \Delta a_{0} + \Delta a_{2} + \epsilon_{1} a_{1}^{\dagger}\right) a_{1} + \left(- \Delta a_{1} + \Delta a_{3} + \epsilon_{2} a_{2}^{\dagger}\right) a_{2}$




```python
# We will now continue showing that the this big matrix can actually be written as the sum of a smaller matrix of only two sites
# To do this we actually first create the system for two sites with an undefined Bogoliubov matrix
Nr = 2

ar = symbols(f"a_:{Nr}")
ar = [Operator(ai) for ai in ar]

cr = Matrix(ar + [Dagger(ai) for ai in ar])

xr = np.array(symbols(f"x_:{2*Nr}:{2*Nr}")).reshape((2*Nr,2*Nr))
Gr = Matrix(xr)

# Execute the formulism to get the Hamiltonian in reduced space
Hr = Dagger(cr)*Gr*cr
h2 = Hr[0]
h2
```




$\displaystyle \left(x_{00} a_{0}^{\dagger} + x_{10} a_{1}^{\dagger} + x_{20} a_{0} + x_{30} a_{1}\right) a_{0} + \left(x_{01} a_{0}^{\dagger} + x_{11} a_{1}^{\dagger} + x_{21} a_{0} + x_{31} a_{1}\right) a_{1} + \left(x_{02} a_{0}^{\dagger} + x_{12} a_{1}^{\dagger} + x_{22} a_{0} + x_{32} a_{1}\right) a_{0}^{\dagger} + \left(x_{03} a_{0}^{\dagger} + x_{13} a_{1}^{\dagger} + x_{23} a_{0} + x_{33} a_{1}\right) a_{1}^{\dagger}$




```python
# We will now match the coefficients
# First step: create a library of coefficients of the full hamiltonian 
lib = {e.args[-1]: e/e.args[-1] for e in h1.args}
lib
```




    {Dagger(a_3): Delta*Dagger(a_2) - epsilon_3*a_3,
     a_0: Delta*a_1 + epsilon_0*Dagger(a_0),
     a_3: -Delta*a_2 + epsilon_3*Dagger(a_3),
     Dagger(a_0): -Delta*Dagger(a_1) - epsilon_0*a_0,
     Dagger(a_1): Delta*Dagger(a_0) - Delta*Dagger(a_2) - epsilon_1*a_1,
     Dagger(a_2): Delta*Dagger(a_1) - Delta*Dagger(a_3) - epsilon_2*a_2,
     a_1: -Delta*a_0 + Delta*a_2 + epsilon_1*Dagger(a_1),
     a_2: -Delta*a_1 + Delta*a_3 + epsilon_2*Dagger(a_2)}




```python
# Now, for each term in the reduced Hamiltonian, we will match the element of the terms in the reduced and full Hamiltonian
Gr1 = Gr.copy()
for i in range(2*Nr):
    # get the operator of the term
    term, ci = h2.args[i].args

    # create a lookup for operators and elements for this particular term 
    store = {e.args[-1]: e/e.args[-1] for e in ncollect(lib[ci], cr).args}
    
    # for each operator in the term, match with the full Hamiltonian and substitute
    for term in ncollect(term, cr).args:
        xii, ki = term.args
        ki = store.get(ki, 0)
        Gr1 = Gr1.subs(xii, ki)
Gr1
```




$\displaystyle \left[\begin{matrix}\epsilon_{0} & 0 & 0 & \Delta\\0 & \epsilon_{1} & - \Delta & 0\\0 & - \Delta & - \epsilon_{0} & 0\\\Delta & 0 & 0 & - \epsilon_{1}\end{matrix}\right]$



OK, great! We have a reduced matrix for the first two sites. Next we need to check whether the matrix of the next two sites is structually identical.


```python
# We will now continue showing that the this big matrix can actually be written as the sum of a smaller matrix of only two sites
# To do this we actually first create the system for two sites with an undefined Bogoliubov matrix
ar = symbols(f"a_:{N}")
ar = ar[2:]
ar = [Operator(ai) for ai in ar]

cr = Matrix(ar + [Dagger(ai) for ai in ar])

xr = np.array(symbols(f"x_:{2*Nr}:{2*Nr}")).reshape((2*Nr,2*Nr))
Gr = Matrix(xr)

# Execute the formulism to get the Hamiltonian in reduced space
Hr = Dagger(cr)*Gr*cr
h2 = Hr[0]

# Now, for each term in the reduced Hamiltonian, we will match the element of the terms in the reduced and full Hamiltonian
for i in range(2*Nr):
    # get the operator of the term
    term, ci = h2.args[i].args

    # create a lookup for operators and elements for this particular term 
    store = {e.args[-1]: e/e.args[-1] for e in ncollect(lib[ci], cr).args}
    
    # for each operator in the term, match with the full Hamiltonian and substitute
    for term in ncollect(term, cr).args:
        xii, ki = term.args
        ki = store.get(ki, 0)
        Gr = Gr.subs(xii, ki)
Gr
```




$\displaystyle \left[\begin{matrix}\epsilon_{2} & 0 & 0 & \Delta\\0 & \epsilon_{3} & - \Delta & 0\\0 & - \Delta & - \epsilon_{2} & 0\\\Delta & 0 & 0 & - \epsilon_{3}\end{matrix}\right]$



It is!
We can now proceed and check whether the sum of these two is actually the same as the full Hamiltonian


```python
# We will now continue showing that the this big matrix can actually be written as the sum of a smaller matrix of only two sites
# To do this we actually first create the system for two sites with an undefined Bogoliubov matrix
result = Add()
for i in range(0, N-1, 1):
    ar = a[i:i+2]
    #print(ar)
    cr = Matrix(ar + [Dagger(ai) for ai in ar])

    xr = np.array(symbols(f"x_:{4}:{4}")).reshape((4,4))
    Gr = Matrix(xr)

    # Execute the formulism to get the Hamiltonian in reduced space
    Hr = Dagger(cr)*Gr*cr
    h2 = Hr[0]

    # Now, for each term in the reduced Hamiltonian, we will match the element of the terms in the reduced and full Hamiltonian
    for i in range(2*Nr):
        # get the operator of the term
        term, ci = h2.args[i].args

        # create a lookup for operators and elements for this particular term 
        store = {e.args[-1]: e/e.args[-1] for e in ncollect(lib[ci], cr).args}

        # for each operator in the term, match with the full Hamiltonian and substitute
        for term in ncollect(term, cr).args:
            xii, ki = term.args
            ki = store.get(ki, 0)
            Gr = Gr.subs(xii, ki)
            
    result += (Dagger(cr)*Gr*cr)[0]
result
```




$\displaystyle \left(\Delta a_{0}^{\dagger} - \epsilon_{1} a_{1}\right) a_{1}^{\dagger} + \left(- \Delta a_{1}^{\dagger} - \epsilon_{0} a_{0}\right) a_{0}^{\dagger} + \left(\Delta a_{1}^{\dagger} - \epsilon_{2} a_{2}\right) a_{2}^{\dagger} + \left(- \Delta a_{2}^{\dagger} - \epsilon_{1} a_{1}\right) a_{1}^{\dagger} + \left(\Delta a_{2}^{\dagger} - \epsilon_{3} a_{3}\right) a_{3}^{\dagger} + \left(- \Delta a_{3}^{\dagger} - \epsilon_{2} a_{2}\right) a_{2}^{\dagger} + \left(- \Delta a_{0} + \epsilon_{1} a_{1}^{\dagger}\right) a_{1} + \left(- \Delta a_{1} + \epsilon_{2} a_{2}^{\dagger}\right) a_{2} + \left(\Delta a_{1} + \epsilon_{0} a_{0}^{\dagger}\right) a_{0} + \left(- \Delta a_{2} + \epsilon_{3} a_{3}^{\dagger}\right) a_{3} + \left(\Delta a_{2} + \epsilon_{1} a_{1}^{\dagger}\right) a_{1} + \left(\Delta a_{3} + \epsilon_{2} a_{2}^{\dagger}\right) a_{2}$




```python
H = Dagger(c)*G*c
h1 = H[0]
h1
```




$\displaystyle \left(- \Delta a_{1}^{\dagger} - \epsilon_{0} a_{0}\right) a_{0}^{\dagger} + \left(\Delta a_{2}^{\dagger} - \epsilon_{3} a_{3}\right) a_{3}^{\dagger} + \left(\Delta a_{1} + \epsilon_{0} a_{0}^{\dagger}\right) a_{0} + \left(- \Delta a_{2} + \epsilon_{3} a_{3}^{\dagger}\right) a_{3} + \left(\Delta a_{0}^{\dagger} - \Delta a_{2}^{\dagger} - \epsilon_{1} a_{1}\right) a_{1}^{\dagger} + \left(\Delta a_{1}^{\dagger} - \Delta a_{3}^{\dagger} - \epsilon_{2} a_{2}\right) a_{2}^{\dagger} + \left(- \Delta a_{0} + \Delta a_{2} + \epsilon_{1} a_{1}^{\dagger}\right) a_{1} + \left(- \Delta a_{1} + \Delta a_{3} + \epsilon_{2} a_{2}^{\dagger}\right) a_{2}$




```python
h1.expand()
```




$\displaystyle \Delta a_{0}^{\dagger} a_{1}^{\dagger} - \Delta a_{1}^{\dagger} a_{0}^{\dagger} + \Delta a_{1}^{\dagger} a_{2}^{\dagger} - \Delta a_{2}^{\dagger} a_{1}^{\dagger} + \Delta a_{2}^{\dagger} a_{3}^{\dagger} - \Delta a_{3}^{\dagger} a_{2}^{\dagger} - \Delta a_{0} a_{1} + \Delta a_{1} a_{0} - \Delta a_{1} a_{2} + \Delta a_{2} a_{1} - \Delta a_{2} a_{3} + \Delta a_{3} a_{2} + \epsilon_{0} a_{0}^{\dagger} a_{0} - \epsilon_{0} a_{0} a_{0}^{\dagger} + \epsilon_{1} a_{1}^{\dagger} a_{1} - \epsilon_{1} a_{1} a_{1}^{\dagger} + \epsilon_{2} a_{2}^{\dagger} a_{2} - \epsilon_{2} a_{2} a_{2}^{\dagger} + \epsilon_{3} a_{3}^{\dagger} a_{3} - \epsilon_{3} a_{3} a_{3}^{\dagger}$




```python
result.expand()
```




$\displaystyle \Delta a_{0}^{\dagger} a_{1}^{\dagger} - \Delta a_{1}^{\dagger} a_{0}^{\dagger} + \Delta a_{1}^{\dagger} a_{2}^{\dagger} - \Delta a_{2}^{\dagger} a_{1}^{\dagger} + \Delta a_{2}^{\dagger} a_{3}^{\dagger} - \Delta a_{3}^{\dagger} a_{2}^{\dagger} - \Delta a_{0} a_{1} + \Delta a_{1} a_{0} - \Delta a_{1} a_{2} + \Delta a_{2} a_{1} - \Delta a_{2} a_{3} + \Delta a_{3} a_{2} + \epsilon_{0} a_{0}^{\dagger} a_{0} - \epsilon_{0} a_{0} a_{0}^{\dagger} + 2 \epsilon_{1} a_{1}^{\dagger} a_{1} - 2 \epsilon_{1} a_{1} a_{1}^{\dagger} + 2 \epsilon_{2} a_{2}^{\dagger} a_{2} - 2 \epsilon_{2} a_{2} a_{2}^{\dagger} + \epsilon_{3} a_{3}^{\dagger} a_{3} - \epsilon_{3} a_{3} a_{3}^{\dagger}$



The two are nearly identical. The center energy terms are bing counted twice...


```python
# When we remove the double counted energies teh two results are equal...
out = result.copy()
for i in [1, 2]:
    out += epsilon[i]*a[i]*Dagger(a[i]) - epsilon[i]*Dagger(a[i])*a[i]
    
h1.expand() == out.expand()
```




    True



# Programmetic construction of the BdG matrix
Considering the construction of the Bogoliubov de Gennes Matrix is procedural we should be able to construct it using some programmatic procedures. 

The Hamiltonian that we will take here is

$\Delta \sum_{k=0}^{N} \big(c_{k+1}^{\dagger}c_{k}^{\dagger} + \mathrm{h.c.}\big)
- \mu \sum_{k=0}^{N} \big( c_{k}^{\dagger}c_{k} + \mathrm{h.c.} \big) 
- \tau \sum_{k=0}^{N} \big(c_{k+1}^{\dagger}c_{k} + \mathrm{h.c.} \big)$ 

There is a conveniant way to construct these summations in Sympy using the Sum functions. This shows a nice mathematical representation of the sum until you actually execute the summation. Unfortunately the Indexed object are always commutation despite setting those to non-commutative. These makes it useles for operators. 


```python
# Checking the commutative behavior of Indexed objects
Indexed.is_commutative = False
c = IndexedBase('c', commutative=False)
# This is not what we expect -> these commute
c[0]*c[1] + (c[1]*c[0])
```




$\displaystyle {c}_{0} {c}_{1} + {c}_{1} {c}_{0}$




```python
# This is what we are supposed to get
c = symbols('c:2', commutative=False)
c[0]*c[1] + c[1]*c[0]
```




$\displaystyle c_{0} c_{1} + c_{1} c_{0}$



Because of this we need to create a Sum function ourself. We will use this to construct the Hamiltonian


```python
# Simple summation function
def Sum(fnc, n, N):
    result = Add()
    for i in range(n+N):
        result += fnc(i)
    return result
        
# Create the non-commuting operators
n = 2
a = symbols(f"a_:{n}")
a = [Operator(ai) for ai in a]

# Construct the Hamiltonian
H0 = -Sum(lambda i: Dagger(a[i])*a[i] , 0, n) + \
     -tau* Sum(lambda i: Dagger(a[i+1])*a[i] + Dagger(a[i])*a[i+1] , 0, n-1) + \
     delta*Sum(lambda i: Dagger(a[i]*a[i+1]) + a[i+1]*a[i] , 0, n-1)

# NOTE: we already limit the index of summations that include i+1 terms to n-1
# Otherwise we would need include an extra operator term and later remove terms included that extra operator with
# H = drop_terms_containing(H0.expand(), [a[-1], Dagger(a[-1])])
# a = a[:-1]
H0
```




$\displaystyle \Delta \left(a_{1}^{\dagger} a_{0}^{\dagger} + a_{1} a_{0}\right) - \tau \left(a_{0}^{\dagger} a_{1} + a_{1}^{\dagger} a_{0}\right) - a_{0}^{\dagger} a_{0} - a_{1}^{\dagger} a_{1}$



So far we just created the original Hamiltonian. To construct a BdG Hamiltonian we need to use some commutation identities and rewrite the Hamiltonian. 

The $\mu$-term

$a_{k}^{\dagger}a_{k} + a_{k+1}^{\dagger}a_{k} = 1 \to a_{k}^{\dagger}a_{k} = \frac{1}{2}a_{k}^{\dagger}a_{k} + \frac{1}{2}(1-a_{k}a_{k}^{\dagger})$

The $\Delta$-term

$a_{k+1}^{\dagger}a_{k}^{\dagger} + a_{k}^{\dagger}a_{k+1}^{\dagger} = 0 \to a_{k+1}^{\dagger}a_{k}^{\dagger} = \frac{1}{2}a_{k+1}^{\dagger}a_{k}^{\dagger} - \frac{1}{2}a_{k}^{\dagger}a_{k+1}^{\dagger}$

$a_{k+1} a_{k} + a_{k} a_{k+1} = 0 \to a_{k+1}a_{k} = \frac{1}{2}a_{k+1}a_{k}- \frac{1}{2}a_{k}a_{k+1}$

The $\tau$-term

$a_{k+1}^{\dagger}a_{k} + a_{k}^{\dagger}a_{k+1} = 0 \to a_{k+1}^{\dagger}a_{k} = \frac{1}{2}a_{k+1}^{\dagger}a_{k} - \frac{1}{2}a_{k}^{\dagger}a_{k+1}$

$a_{k}^{\dagger}a_{k+1} + a_{k+1}^{\dagger}a_{k} = 0 \to a_{k}^{\dagger}a_{k+1} = \frac{1}{2}a_{k}^{\dagger}a_{k+1} - \frac{1}{2}a_{k+1}^{\dagger}a_{k}$

After making these replacement we will be able to construct the BdG Hamiltonian


```python
# To replace them we need to trick because Sympy's subs method we do the replacement sequentially which means, 
# parts that are added duning early replacements will we replaced later on during later replacements.
# The trick is to do the replcement is two stages using dummies in between.
submap = [[Dagger(a[i])*a[i], 1/2*(Dagger(a[i])*a[i] + 1 - a[i]*Dagger(a[i]))] for i in range(2)] + \
         [[Dagger(a[i]*a[j]), 1/2*(Dagger(a[j]*a[i]) - Dagger(a[i]*a[j]))] for i in range(2) for j in range(2)] + \
         [[a[i]*a[j], 1/2*(a[i]*a[j] - a[j]*a[i])] for i in range(2) for j in range(2)] + \
         [[Dagger(a[i+1])*a[i], 1/2*(Dagger(a[i+1])*a[i] - a[i]*Dagger(a[i+1]))] for i in range(2-1)] + \
         [[Dagger(a[i])*a[i+1], 1/2*(Dagger(a[i])*a[i+1] - a[i+1]*Dagger(a[i]))] for i in range(2-1)]

H = H0.expand()
for item in submap:
    key = Dummy()
    old, new = item
    
    H = H.subs(old, key)
    item[0] = key
    
H = H.subs(submap)
H.expand()
```




$\displaystyle 0.5 \Delta a_{0}^{\dagger} a_{1}^{\dagger} - 0.5 \Delta a_{1}^{\dagger} a_{0}^{\dagger} - 0.5 \Delta a_{0} a_{1} + 0.5 \Delta a_{1} a_{0} - 0.5 \tau a_{0}^{\dagger} a_{1} - 0.5 \tau a_{1}^{\dagger} a_{0} + 0.5 \tau a_{0} a_{1}^{\dagger} + 0.5 \tau a_{1} a_{0}^{\dagger} - 1.0 - 0.5 a_{0}^{\dagger} a_{0} - 0.5 a_{1}^{\dagger} a_{1} + 0.5 a_{0} a_{0}^{\dagger} + 0.5 a_{1} a_{1}^{\dagger}$




```python
# We construct the operator vectors!
c = Matrix(a + [Dagger(ai) for ai in a])
Dagger(c)
```




$\displaystyle \left[\begin{matrix}a_{0}^{\dagger} & a_{1}^{\dagger} & a_{0} & a_{1}\end{matrix}\right]$




```python
# Before we can create the 
expr = ncollect(H.expand())
expr
```




$\displaystyle -1.0 + a_{0}^{\dagger} \left(0.5 \Delta a_{1}^{\dagger} - 0.5 \tau a_{1} - 0.5 a_{0}\right) + a_{1}^{\dagger} \left(- 0.5 \Delta a_{0}^{\dagger} - 0.5 \tau a_{0} - 0.5 a_{1}\right) + a_{0} \left(- 0.5 \Delta a_{1} + 0.5 \tau a_{1}^{\dagger} + 0.5 a_{0}^{\dagger}\right) + a_{1} \left(0.5 \Delta a_{0} + 0.5 \tau a_{0}^{\dagger} + 0.5 a_{1}^{\dagger}\right)$




```python
# This routine actually creates the BdG matrix Hamiltonian. 
# First we express the Hamiltonian in terms of the leading (!!) operators.
expr = ncollect(H.expand())

# Next we collect the terms belonging to each of the elements in the operator vector Dagger(c). 
# Those select the rows in the matrix
G = Matrix.zeros(len(c))
expr = ncollect(H.expand())
for i, l_ops in enumerate(Dagger(c)):
    h = get_coefficient(expr, l_ops)
    
    # For each term, we again match the operators but now to the operator vector c (!!)
    # This selects the columns. 
    for j, r_ops in enumerate(c):
        coeff = get_coefficient(h, r_ops)
        G[i,j] = coeff 

# This gives us the BdG Hamiltonian in matrix form
G
```




$\displaystyle \left[\begin{matrix}-0.5 & - 0.5 \tau & 0 & 0.5 \Delta\\- 0.5 \tau & -0.5 & - 0.5 \Delta & 0\\0 & - 0.5 \Delta & 0.5 & 0.5 \tau\\0.5 \Delta & 0 & 0.5 \tau & 0.5\end{matrix}\right]$



The original Hamiltonain in expanded sum-form can be obtained again (up to constant terms) using

$\mathrm{C^{\dagger} g C}$

where $\mathrm{C^{\dagger}} = [a_0^{\dagger}, a_1^{\dagger}, ... a_N^{\dagger}, a_0, a_1, ..., a_N]$.


```python
# Note the constant term '1' is missing
ncollect((Dagger(c)*G*c)[0].expand())
```




$\displaystyle a_{0}^{\dagger} \left(0.5 \Delta a_{1}^{\dagger} - 0.5 \tau a_{1} - 0.5 a_{0}\right) + a_{1}^{\dagger} \left(- 0.5 \Delta a_{0}^{\dagger} - 0.5 \tau a_{0} - 0.5 a_{1}\right) + a_{0} \left(- 0.5 \Delta a_{1} + 0.5 \tau a_{1}^{\dagger} + 0.5 a_{0}^{\dagger}\right) + a_{1} \left(0.5 \Delta a_{0} + 0.5 \tau a_{0}^{\dagger} + 0.5 a_{1}^{\dagger}\right)$



# alternative expression using Pauli matrices
Finally following https://topocondmat.org/w1_topointro/1D.html it should also be possible to express the BdG Hamiltonian using Pauli-matrixes. Ofcourse, this would help you actually constuct the Hamiltonian, it's just a nice way to write it. 


```python
# We define the Pauli-matrices
sigma0 = Matrix([[1,0],[0,1]])
sigmax = Matrix([[0,1],[1,0]])
sigmay = Matrix([[0,-1j],[1j,0]])
sigmaz = Matrix([[1,0],[0,-1]])
```

The approach makes use of a 'selecting' vector $|n\rangle$ of size n consisting of a single non-zero element which is 1 i.e.

$|n\rangle = [..., 1,...]^{T}$ e.g. $|2\rangle = [0, 0, 1, ....]$

Note, we index from 0 here.


```python
# Conveniance function to create those special 'selecting vectors'
def nvec(N, k = 0):
    v = Matrix(N*[0])
    v[k] = 1
    return v
nvec(4, 1)
```




$\displaystyle \left[\begin{matrix}0\\1\\0\\0\end{matrix}\right]$



At the website they give the folowing identity:

$\mathrm{C^{\dagger}} \sigma_z |n\rangle \langle n| \mathrm{C} = 2 a_n^{\dagger} a_n - 1$

This seems weird because $\sigma_z$ is of size [2,2] but $|n\rangle \langle n|$ will be of size [n,n]. Moreover $\mathrm{C} is on size [2*n,1]$. To resolve this we need to take into account the Tensor product of $\sigma_z$ and $|n\rangle \langle n|$ which will be of size [2n, 2n]

$\mathrm{C^{\dagger}} \Big( \sigma_z \otimes |n\rangle \langle n| \Big) \mathrm{C} = 2 a_n^{\dagger} a_n - 1$


```python
# We check the identity
# We just reuse the c-vector of before
n0 = nvec(int(len(c)/2), 0)
eq = Dagger(c) * TensorProduct(mu*sigmaz, n0*Dagger(n0)) * c
eq[0]
```




$\displaystyle \mu a_{0}^{\dagger} a_{0} - \mu a_{0} a_{0}^{\dagger}$



Taking into account the commutation relation $a_0^{\dagger}a_0 + a_0 a_0^{\dagger} = 1$, this identity is indeed correct.

Next we can check the formation of the BdG Hamiltonian i.e.


```python
n1 = nvec(2, 1)
TensorProduct(tau*sigmaz + 1j*delta*sigmay, n0*Dagger(n1))
```




$\displaystyle \left[\begin{matrix}0 & \tau & 0 & 1.0 \Delta\\0 & 0 & 0 & 0\\0 & - 1.0 \Delta & 0 & - \tau\\0 & 0 & 0 & 0\end{matrix}\right]$




```python
Dagger(TensorProduct(tau*sigmaz + 1j*delta*sigmay, n0*Dagger(n1)))
```




$\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\tau^{\dagger} & 0 & - 1.0 \Delta^{\dagger} & 0\\0 & 0 & 0 & 0\\1.0 \Delta^{\dagger} & 0 & - \tau^{\dagger} & 0\end{matrix}\right]$




```python
# Lets construct the whole thing

result = Matrix.zeros(2*n, 2*n) #Add()
for ni in range(n):
    n0 = nvec(n, ni)
    result += -0.5*TensorProduct(mu*sigmaz, n0*Dagger(n0))
    
    if ni < n-1:
        n1 = nvec(n, ni+1)
        result += -0.5*TensorProduct(tau*sigmaz - 1j*delta*sigmay, n0*Dagger(n1))
        result += -0.5*TensorProduct(Dagger(tau*sigmaz - 1j*delta*sigmay), n1*Dagger(n0))
result
```




$\displaystyle \left[\begin{matrix}- 0.5 \mu & - 0.5 \tau & 0 & 0.5 \Delta\\- 0.5 \tau^{\dagger} & - 0.5 \mu & - 0.5 \Delta^{\dagger} & 0\\0 & - 0.5 \Delta & 0.5 \mu & 0.5 \tau\\0.5 \Delta^{\dagger} & 0 & 0.5 \tau^{\dagger} & 0.5 \mu\end{matrix}\right]$




```python
# The Hamiltonian in sum-form can again be obtained by multiplying with the operator vector. 
ncollect((Dagger(c)*result*c)[0].expand())
```




$\displaystyle a_{0}^{\dagger} \left(- 0.5 \Delta a_{1}^{\dagger} - 0.5 \mu a_{0} - 0.5 \tau a_{1}\right) + a_{1}^{\dagger} \left(- 0.5 \mu a_{1} + 0.5 \Delta^{\dagger} a_{0}^{\dagger} - 0.5 \tau^{\dagger} a_{0}\right) + a_{0} \left(0.5 \Delta a_{1} + 0.5 \mu a_{0}^{\dagger} + 0.5 \tau a_{1}^{\dagger}\right) + a_{1} \left(0.5 \mu a_{1}^{\dagger} - 0.5 \Delta^{\dagger} a_{0} + 0.5 \tau^{\dagger} a_{0}^{\dagger}\right)$




```python

```

50f9883218aa22ee4748bb305c643727fe8c2650
ipynb
Jupyter Notebook
sympsi/examples/BogoliubovFormalism.ipynb
dmtvanzanten/sympsi
d717d7ceac5ae2b928939f1eda8cbd06f4e18c29
BSD-3-Clause
sympsi/examples/BogoliubovFormalism.ipynb
dmtvanzanten/sympsi
d717d7ceac5ae2b928939f1eda8cbd06f4e18c29
BSD-3-Clause
sympsi/examples/BogoliubovFormalism.ipynb
dmtvanzanten/sympsi
d717d7ceac5ae2b928939f1eda8cbd06f4e18c29
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES