Supercurrent distribution in real-space and anomalous paramagnetic response in a superconducting quasicrystal
Model
Original Hamiltonian
The Hamiltonian is given by $$\hat{\mathcal{H}}=-\sum_{\langle i,j\rangle \sigma} \left\{ t\exp \left( -i\int^{\mathbf{r}_i}_{\mathbf{r}_j} \mathbf{A}(\mathbf{r})\cdot d\mathbf{r}\right) \hat{c}^\dagger_{i\sigma}\hat{c}_{j\sigma}+\text{H.c.}\right\} \\+U\sum_i \hat{n}_{i\uparrow}\hat{n}_{i\downarrow}-\mu\sum_{i\sigma}\hat{n}_{i\sigma}$$
Here, $\hat{c}^\dagger_{i\sigma}$ creates an electron of spin $\sigma$ at site $i$.
then "H.c." is $$t\exp\left(-i\int^{\mathbf{r}_j}_{\mathbf{r}_i} \mathbf{A}(\mathbf{r})\cdot d\mathbf{r}\right) \hat{c}^\dagger_{j\sigma}\hat{c}_{i\sigma}$$
We define a local electron density $n_i=\sum_\sigma\langle \hat{n}_{i\sigma}\rangle$ with $\hat{n}_{i\sigma}=\hat{c}^\dagger_{i\sigma}\hat{c}_{i\sigma}$
The chemical potential $\mu$ is tuned to fix the average electron density $\bar{n}=\sum_i n_i/N$ where $N$ is the system size. (??)
In the case of uniform vector potential $\mathbf{A}$, the Peierls phase can be rewritten as $$-i\int^{\mathbf{r}_i}_{\mathbf{r}_j} \mathbf{A}(\mathbf{r})\cdot d\mathbf{r}=-i\mathbf{A}\cdot\mathbf{r}_{ij},$$ where $\mathbf{r}_{ij}=\mathbf{r}_i-\mathbf{r}_j=a\mathbf{e}_n$ is the bond vector between sites $i$ and $j$.
Mean-Field Approximation
Get effective quadratic (gaussian) Hamiltonian of interaction term.
$$U\sum_i\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}= U\sum_i\langle \hat{n}_{i\uparrow}\rangle\hat{n}_{i\downarrow}+U\sum_i\hat{n}_{i\uparrow}\langle \hat{n}_{i\downarrow}\rangle\\ -\sum_i U\langle \hat{n}_{i\uparrow}\rangle\langle\hat{n}_{i\downarrow}\rangle+\left(\hat{n}_{i\uparrow}-\langle\hat{n}_{i\uparrow}\rangle\right)\left(\hat{n}_{i\downarrow}-\langle\hat{n}_{i\downarrow}\rangle\right)\\ +(\text{other contact term})$$
Third term is contant term, and fourth term is fluctuation term. By neglecting fluctuation, we do "mean-field approximate".
And $$n_{i\sigma}=\langle \hat{n}_{i\sigma}\rangle$$ is site-dependent density order parameter.
We consider other contact term, "site-dependent superconducting order parameter".
$$-U\sum_i\langle \hat{c}^\dagger_{i\uparrow}\hat{c}^\dagger_{i\downarrow}\rangle \hat{c}_{i\downarrow}\hat{c}_{i\uparrow}-U\sum_i\hat{c}^\dagger_{i\uparrow}\hat{c}^\dagger_{i\downarrow}\langle\hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\rangle\cdots$$
changes $\hat{c}$ makes $-$ sign.
Then site-dependent superconducting order parameter is $$\Delta_i=-U\langle \hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\rangle$$
We denote $$\Delta^*_i\delta_{i,j}\hat{c}_{i\uparrow}\hat{c}_{j\downarrow}=(-U\langle \hat{c}^\dagger_{i\uparrow}\hat{c}^\dagger_{i\downarrow}\rangle\delta_{i,j}\hat{c}_{i\uparrow}\hat{c}_{i\downarrow} $$
(Also, there is another mean-field term, but we don't consider it. $\hat{S}$)
Since Hamiltonian is quadratic:
$$\hat{\mathcal{H}}_{eff}=-\sum_{i,j}\left[-t\delta_{\langle i,j\rangle}e^{-i\mathbf{A}\cdot\mathbf{r}_{ij}}+(Un_{i\downarrow}-\mu)\delta_{i,j}\right]\hat{c}^\dagger_{i\uparrow}\hat{c}_{j\uparrow}\\ + \left[-t\delta_{\langle i,j\rangle}e^{-i\mathbf{A}\cdot\mathbf{r}_{ij}}+(Un_{i\uparrow}-\mu)\delta_{i,j}\right]\hat{c}^\dagger_{i\downarrow}\hat{c}_{j\downarrow} \\ +\Delta_i\delta_{i,j}\hat{c}_{i\uparrow}\hat{c}_{j\downarrow}+\Delta^*_i\delta_{i,j}\hat{c}^\dagger_{i\uparrow}\hat{c}^\dagger_{j\downarrow}$$
Then we can make matrix form $$\hat{\mathcal{H}}=\hat{\Psi}^\dagger \mathcal{H}\hat{\Psi} $$ when $$\mathcal{H}=\begin{bmatrix}K_{i,j\uparrow}&\Delta_i\delta_{i,j}\\ \Delta^*_i\delta_{i,j}& -K^*_{i,j\downarrow}\end{bmatrix}$$ with basis $$\hat{\Psi}^\dagger=\begin{bmatrix}\hat{c}^\dagger_{1\uparrow}&\cdots&\hat{c}^\dagger_{N\uparrow}&\hat{c}_{1\downarrow}&\cdots&\hat{c}_{N\downarrow}\end{bmatrix}$$ and $K_\sigma$ matrix $$K_{\sigma i,j}=-t\delta_{\langle i,j\rangle}e^{-i\mathbf{A}\cdot\mathbf{r}_{ij}}+(Un_{i\bar{\sigma}}-\mu)\delta_{i,j}$$
We denote $$-K^*_{i,j\downarrow}\hat{c}_{i\downarrow}\hat{c}^\dagger_{j\downarrow}=K^*_{i,j\downarrow}\hat{c}^\dagger_{j\downarrow}\hat{c}_{i\downarrow}=K^*_{j,i\downarrow}\hat{c}_{i\downarrow}^\dagger\hat{c}_{j\downarrow}=K^\dagger_{i,j\downarrow}\hat{c}^\dagger_{i\downarrow}\hat{c}_{j\downarrow}=K_{i,j\downarrow}\hat{c}^\dagger_{i\downarrow}\hat{c}_{j\downarrow}$$
Bogolubov Equation
Since Hamiltonian matrix is Hermitian, it can be diagonalized.
$$\hat{\mathcal{H}}_{eff}=E_g+\sum_{n,\alpha}\epsilon_n\hat{\gamma}^\dagger_{n\alpha}\hat{\gamma}_{n\alpha}$$
when $\hat{\gamma}$ construct Clifford algebra $$\left\{ \hat{\gamma}^\dagger_{\epsilon\alpha},\hat{\gamma}_{\epsilon'\beta}\right\}=\delta_{\epsilon\epsilon'}\delta_{\alpha\beta},\quad \left\{\hat{\gamma}_{\epsilon\alpha},\hat{\gamma}_{\epsilon'\beta}\right\}=0$$
And unitary trasnformation between $\hat{c}$ and $\hat{\gamma}$ is defined as $$\begin{bmatrix}\hat{c}_{1\uparrow}\\ \vdots\\ \hat{c}_{N\uparrow}\\ \hat{c}^\dagger_{1\downarrow}\\ \vdots\\ \hat{c}^\dagger_{N\downarrow} \end{bmatrix}=\begin{bmatrix}u_{i\epsilon}& -v^*_{i\epsilon}\\ v_{i\epsilon}& u^*_{i\epsilon}\end{bmatrix}\begin{bmatrix}\hat{\gamma}_{1\uparrow}\\ \vdots\\ \hat{\gamma}_{N\uparrow}\\ \hat{\gamma}^\dagger_{1\downarrow}\\ \vdots\\ \hat{\gamma}^\dagger_{N\downarrow} \end{bmatrix}$$
We can calculate commutation relation $$\begin{cases}[\hat{\mathcal{H}}_{eff},\gamma_{n\alpha}]=-\epsilon_n\gamma_{n\alpha}\\ [\hat{\mathcal{H}}_{eff},\gamma^\dagger_{n\alpha}]=\epsilon_n\gamma^\dagger_{n\alpha}\end{cases}$$
We can calculate commutation relation from past Hamiltonian form, $$\begin{cases}[\hat{c}_{i\uparrow},\hat{\mathcal{H}}_{eff}]=K_{i,j\uparrow}\hat{c}_{j\uparrow}+\Delta_i\hat{c}^\dagger_{i\downarrow}\\ [\hat{c}_{i\downarrow},\hat{\mathcal{H}}_{eff}]=K_{i,j\downarrow}\hat{c}_{j\downarrow}+\Delta_i\hat{c}^\dagger_{i\uparrow}\end{cases}$$
Combind that, we get $$\begin{cases}\epsilon u(r)=K_\uparrow u(r)+\Delta(r)v(r)\\ \epsilon v(r)=-K^*_\downarrow v(r)+\Delta^*(r)u(r)\end{cases}$$
Finally we conclude Bogolubov equation, $$\sum_j\hat{\mathcal{H}}_{i,j}\begin{bmatrix}u_{j\epsilon}\\ v_{j\epsilon}\end{bmatrix}= E_\epsilon \begin{bmatrix}u_{i\epsilon}\\ v_{i\epsilon}\end{bmatrix}$$
Compute Expectation Value
Eigenstates of $\hat{\mathcal{H}}_{eff}$ is $|\phi\rangle$. $$\hat{\mathcal{H}}_{eff}|\phi\rangle =E_\phi |\phi\rangle $$
then expectation value of $\hat{\gamma}^\dagger\hat{\gamma}$ is $$\langle \hat{\gamma}^\dagger_{n\alpha}\hat{\gamma}_{m\beta}\rangle = \frac{\sum_\phi\langle \phi |\hat{\gamma}^\dagger_{n\alpha}\hat{\gamma}_{m\beta}|\phi\rangle e^{-\beta E_\phi}}{\sum_\phi e^{-\beta E_\phi}}=\delta_{nm}\delta_{\alpha\beta} f_n$$ since only same $E_n$ term servive of total $\hat{\mathcal{H}}_{eff}$.
Which means, states that have energy $E_n$ servive and eigenvalue is $1$. $$\hat{\gamma}^\dagger_{n\alpha}\hat{\gamma}_{m\beta}|\phi\rangle = \begin{cases}1|\phi\rangle& E_\phi=E_n\\ 0 & E_\phi\ne E_n\end{cases}$$
And $f(E)=1/(e^{E/T}+1)$ is the Fermi-Dirac distribution function at temperature $T$.
We can obtain $\Delta_i$ $$U\langle \hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\rangle \\ =\sum_\epsilon U \left\langle \left(u^*_{i\epsilon}\hat{\gamma}_{\epsilon\downarrow}-v_{i\epsilon}\hat{\gamma}^\dagger_{\epsilon\uparrow}\right)\left(u^*_{i\epsilon}\hat{\gamma}_{\epsilon\uparrow}+v_{i\epsilon}\hat{\gamma}^\dagger_{\epsilon\downarrow}\right)\right\rangle \\ =-\sum_\epsilon U\left( u^*_{i\epsilon}v_{i\epsilon}\langle \hat{\gamma}_{i\downarrow}\hat{\gamma}^\dagger_{i\downarrow}\rangle-v_{i\epsilon}u^*_{i\epsilon}\langle \hat{\gamma}^\dagger_{\epsilon\uparrow}\hat{\gamma}_{\epsilon\uparrow}\rangle \right)\\ =-U\sum_\epsilon u^*_{i\epsilon} v_{i\epsilon}[1-2f(E_{\epsilon})]$$
We can also obtain $n_{i\uparrow}$ $$\langle \hat{c}^\dagger_{i\uparrow}\hat{c}_{i\uparrow}\rangle = \sum_\epsilon |u_{i\epsilon}|^2 f(E_\epsilon)$$
We can also obtain $n_{i\downarrow}$ $$\langle \hat{c}^\dagger_{i\downarrow}\hat{c}_{i\downarrow}\rangle = \sum_\epsilon |u_{i\epsilon}|^2[1- f(E_\epsilon)]$$
So get wavefunction $u$ and $v$ by solving BdG equation is important!
Local Supercurrent
$$\mathbf{J}_{j\rightarrow i} = -\frac{\partial \langle \hat{\mathcal{H}}(\mathbf{A})\rangle}{\partial \mathbf{A}} \\ =2t\text{Im}\left\{ \exp (-i\mathbf{A}\cdot \mathbf{r}_{ij})\sum_\sigma \langle \hat{c}^\dagger_{i\sigma}\hat{c}_{j\sigma}\rangle \right\} \mathbf{r}_{ij}$$
since $\mathbf{A}$ is only in hopping term.
Cooper Pairing
Pair Potential
$$\Delta_{ij}\hat{c}^\dagger_{i\uparrow}\hat{c}^\dagger_{i\downarrow}$$
S-wave
Cooper pair is spherical symmetry.
$$\Delta=\Delta^T$$
Then it can be write $$\Delta_{ij}=\frac{U}{2}\left(\langle \hat{c}_{i\downarrow}\hat{c}_{j\uparrow}\rangle -\langle \hat{c}_{j\downarrow}\hat{c}_{i\uparrow}\rangle\right)$$
which means singlet, spartially symmetric.
for our case, pair potential is more simple. $$\Delta_{ij}=\delta_{i,j}U\langle \hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\rangle$$
P-wave
not yet
Appendix
Second Quantization
first quantization is making Hilbert space of possible states.
State vector is defined as element of Hilbert space $$|r\rangle \in \mathcal{H}$$
And state is evolute through time, it's Schrodinger equation.
To represent that abstract algebra, we bring complete basis.
Ground State (Vaccum)
When Hamiltonian have only hopping term, ground state is $|0\rangle_c$.
When Hamiltonian have BCS term also, ground state is $|0\rangle_\gamma$.
$$\begin{cases}\hat{c}_{i\sigma}|0\rangle_c=0\\ \hat{\gamma}_\alpha|0\rangle_\gamma=0\end{cases}$$ $$|0\rangle_\gamma= \exp \left(S_{ij}\hat{c}^\dagger_{i\uparrow}\hat{c}^\dagger_{j\downarrow}\right) |0\rangle_c$$ which means $$=|0\rangle + |1\rangle +|2\rangle +\cdots$$ Bose condensation of singlet pairing.
By mathematical induction (to series of exp), $$S_{ij}=(vu^{-1})^*$$
Self-Consistency
$$n_{ij}\equiv \langle \hat{c}^\dagger_i \hat{c}^\dagger_j\rangle=\frac{1}{Z}\text{Tr}\left(e^{-\beta \hat{\mathcal{H}}_{eff}}\hat{c}^\dagger_i\hat{c}_j\right)$$
where $$Z=\text{Tr}\left(e^{-\beta \hat{\mathcal{H}}_{eff}}\right)$$
since $\hat{\mathcal{H}}_{eff}$ contains $n_{ij}$, it called "self-consistency condition".
We can understand this as minimizes the free energy $F=H-TS$.
$$0=\frac{d}{dn_{ij}}F=\frac{d}{dn_{ij}}\left(-\frac{1}{\beta}\ln Z\right) \\ =\frac{1}{Z}\text{Tr}\left( e^{-\beta H}\frac{d}{dn_{ij}}H\right) \\ = \frac{1}{Z}\text{Tr}\left(e^{-\beta H}\left(\sum_{i'j'}V_{iji'j'}(\hat{c}^\dagger_{i'}\hat{c}_{j'}-n_{i'j'})\right)\right)\\ = \sum_{i'j'}V_{iji'j'}\left(\langle \hat{c}^\dagger_{i'}\hat{c}_{j'}\rangle -n_{i'j'}\right)$$
what "solve" means?
I think it is finding good parabola location.
Because Mean-field is approximate to quadratic action.
Broken Translation Symmetry: $n$
Translation operator is defined like this: $$\hat{T}(a)|i\rangle = |i+a\rangle$$
Local number operator break translation symmetry.
$$\hat{T}(a)\hat{n}_i\hat{T}^\dagger(a)=\hat{n}_{i+a}\ne \hat{n}_i$$
So nonzero $n_i$ makes change of $\hat{\mathcal{H}}_{eff}$, then ground state are not symmetry under translation.
Broken Gauge Symmtery: $\Delta$
$$\hat{\tilde{c}}_{i\sigma}=\hat{c}_{i\sigma}e^{-i\mathbf{A}\cdot\mathbf{r}_i},\quad \hat{\tilde{c}}_{i\sigma}^\dagger=\hat{c}_{i\sigma}^\dagger e^{i\mathbf{A}\cdot\mathbf{r}_i}$$
Then wavefunctions $u_\epsilon(\mathbf{r}_i)$, $v_\epsilon(\mathbf{r}_i)$ also change $$\tilde{u}_\epsilon(\mathbf{r}_i) = u_\epsilon(\mathbf{r}_i)e^{i\mathbf{A}\cdot\mathbf{r}_i},\quad \tilde{v}_\epsilon(\mathbf{r}_i) = v_\epsilon(\mathbf{r}_i)e^{-i\mathbf{A}\cdot\mathbf{r}_i}$$
The superconducting order parameter is change $$\tilde{\Delta}_i=\Delta_i e^{2i\mathbf{A}\cdot\mathbf{r}_i}$$
Hamiltonian also change, $$\hat{H} = \sum_{i,j} \begin{bmatrix}\hat{c}_{i\uparrow}^\dagger \quad \hat{c}{i\downarrow}\end{bmatrix} \hat{\mathcal{H}}{i,j} \begin{bmatrix} \hat{c}_{j\uparrow}\\ \hat{c}_{j\downarrow}^\dagger \end{bmatrix},$$
where $$\hat{\mathcal{H}}_{i,j} = \begin{bmatrix} \tilde{K}_{i,j} & \tilde{\Delta}_i \delta_{i,j} \\ \tilde{\Delta}_i^* \delta_{i,j} & -\tilde{K}_{i,j}^* \end{bmatrix}$$
where $\tilde{K}_{\sigma,i,j} = -t\delta_{(i,j)} + (Un_{i\bar{\sigma}} - \mu)\delta_{i,j}$.
$$\langle \hat{c}_{i\uparrow}^\dagger \hat{c}_{j\uparrow} \rangle = \langle \hat{\tilde{c}}_{i\uparrow}^\dagger \hat{\tilde{c}}_{j\uparrow} \rangle e^{iA\cdot r_{ij}}$$
$$\langle \hat{c}_{i\downarrow}^\dagger \hat{c}_{j\downarrow} \rangle = \langle \hat{\tilde{c}}_{i\downarrow}^\dagger \hat{\tilde{c}}_{j\downarrow} \rangle e^{iA\cdot r_{ij}}$$
Local supercurrent becomes $$J_{j\rightarrow i} = 2t \text{Im} \left( \sum_{\sigma} \langle \hat{\tilde{c}}_{i\sigma}^\dagger \hat{\tilde{c}}_{j\sigma} \rangle \right) r_{ij}$$
$$J_{j\rightarrow i}^{\text{para}} = 2t \cos (A \cdot r_{ij}) \text{Im} \left\{ e^{iA\cdot r_{ij}} \sum_{\sigma} \langle \hat{\tilde{c}}_{i\sigma}^\dagger \hat{\tilde{c}}_{j\sigma} \rangle \right\} r_{ij}$$
Reference
Bruus & Flensberg - Introduction to Many-body quantum theory in condensed matter physics
Gennes - Superconductivity of Metals Alloys
FULLTEXT01.pdf (diva-portal.org)