---
# System prepended metadata

title: Excitation Diffusion in doped NaYF$_4$
tags: [PhD projects, Diffusion, 'but this is not soft matter :(']

---


$$
\newcommand{\ket}[1]{\left|{#1}\right\rangle}
\newcommand{\bra}[1]{\left\langle{#1}\right|}
\newcommand{\proj}[1]{\left|{#1}\middle\rangle\middle\langle{#1}\right|}
\newcommand{\mel}[3]{\left\langle{#1}\middle|{#2}\middle|{#3}\right\rangle}
\newcommand{\quot}[2]{{{#1}\big/{ #2}}}
\newcommand{\barrow}[]{\rightarrowtail\!\!\!\!\!\rightarrow}
\newcommand{\floor}[1]{\lfloor{#1}\rfloor}
\newcommand{\lb}[2]{\left[{#1},{#2}\right]}
\newcommand{\alb}[2]{\left\{{#1},{#2}\right\}}
\newcommand{\avg}[1]{\left\langle{#1}\right\rangle}
\newcommand{\tr}[1]{\mathrm{Tr}\left\{{#1}\right\}}
$$

# Excitation Diffusion in doped NaYF$_4$

###### tags: `PhD projects` `Diffusion` `but this is not soft matter :(`

---

### Owners (the only one with the permission to edit the main text, others can comment)

Alptug, Laura, Freddy(coming soon?) 

---
## To do
- [x] Think
- [ ] Calculate :)
- [ ] Write the introduction of "Total System" section

---
This document is Alptug's theoretical thoughts on Rene's project

## Short description
We have a sample full of $X\%$ doped $\text{NaYF}_4$ microcrystals. We shine laser on the sample, excite the doped sites, and measure how the excitation diffuses.

### What do we know?
* Only the doped sites can be excited
* hopping probability scales with $1/r^{6}$ (derivation?)
* Hopping between microcrystals is possible (How probable?)

### What is Rene up to?
* He carries out the experiments (I am still amazed how they can measure individual excitations)
* he does modelling:
    * inter-microcrystal hopping (current)
    * ideal gas model on the lattice
    * coupled N-level systems, write the DE and fit (I do not like this approach, too many parameters, once a wise man said  "With four parameters I can fit an elephant, and with five I can make him wiggle his trunk") 

## Theory
In my opinion Hubbard model fits perfectly for this system where the lattice is formed only by doping atoms of a 100% doped $\text{NaYF}_4$. (I ignore the inter-microcrystal interaction for now, I believe we can add that interaction in a mean-fieldesque fashion in the future).We have three options to consider:
1. Doping atoms are 2-level systems, i.e., they are either at ground state $\left|g\right>$ or excited state $\left|e\right>$.
2. Doping atoms are N-level systems, i.e., they can occupy any state $\left|n\right>$ where $n\in\mathbb{Z}/N\mathbb{Z}$.
3. Doping atoms have (countable) infinitely many levels, i.e., they can occupy any state $\left|n\right>$ where $n\in\mathbb{Z}$.

### Option 1
If we consider option (1), we can model the excitations as spinless fermions where the Hamiltonian is given by:
$$
H = -\sum_{i\neq j}t_{ij}c^\dagger_i c^{\vphantom{\dagger}}_j ,
$$
where $t_{ij}$ is the hopping matrix, and $c_i^{(\dagger)}$ is the annihilation(creation) operator for lattice site $i$. The only non-trivial anti-commutator is given by
$$
\left\{ c^{\vphantom{\dagger}}_i,c^\dagger_j \right\} = \delta_{ij}.
$$
This model reduces to tight binding model if we restrict the sum into nearest neighbours (NN) and take $t_{ij}$ to be constant.

### Option 2
In this option, we have to assume N is odd. Then, we can model the excitations as spin-$(N/2-1)$ fermions where the Hamiltonian is given by:
$$ 
H=
-\sum_{i\neq j,\sigma}t_{ij}c^\dagger_{i,\sigma} c^{\vphantom{\dagger}}_{j,\sigma}  + \frac{U}{2} \sum_{i,\sigma \neq \sigma^\prime}c^\dagger_{i,\sigma} c^{\vphantom{\dagger}}_{i,\sigma}c^\dagger_{i,\sigma^\prime} c^{\vphantom{\dagger}}_{i,\sigma^\prime},
$$
where $t_{ij}$ is the hopping matrix, $U$ is the on-site interaction energy, and $c_{i,\sigma}^{(\dagger)}$ is the annihilation(creation) operator for lattice site $i$ and spin $\sigma$. The only non-trivial anti-commutator is given by
$$
\left\{ c^{\vphantom{\dagger}}_{i,\sigma\vphantom{\sigma^\prime}},c^\dagger_{j,\sigma^\prime} \right\} = \delta_{ij}\delta_{\sigma\sigma^\prime}.
$$

### Option 3
In this option, we model the excitations as (non-)interacting bosons. The Hamiltonian is given by:

$$
H= - \sum_{i\neq j}t_{ij} b^\dagger_ib^{\vphantom{\dagger}}_j + \frac{U}{2}\sum_i b_i^\dagger b_i^\vphantom{\dagger} \left(b_i^\dagger b_i^\vphantom{\dagger} - 1\right),
$$
where $t_{ij}$ is the hopping matrix, $U$ is the on-site interaction energy, and $b_{i}^{(\dagger)}$ is the bosonic annihilation(creation) operator for lattice site $i$. The only non-trivial anti-commutator is given by
$$
\left[ b^{\vphantom{\dagger}}_i,b^\dagger_j \right] = \delta_{ij}.
$$

## What about doping ratio?
In the above mentioned models, we worked in a lattice of doping atoms. We model 100% doped lattice as the perfect crystal and we can lower the percentage of doping by introducing vacancy impurities since the excitations cannot live in non-doped sites. Furthermore, we can deal with these "impurities" using disorder averaged Green's functions.


# Toy Model (Markov Chain)
Let $\mathbf{T}$ be the transition matrix that advances the system by $\Delta t$ amount of time. Assuming time is discritized in $\Delta t$ steps, we define the Green's funtion as
$$
G(i,j,t-t_0) = \left(\mathbf{T}^{(t-t_0)/\Delta t}\right)_{ij}.
$$
Given a particle located in the lattice site $j$ at time $t_0$. The Green's function gives the probability of finding that particle at lattice site &i$ at time $t$. Therefore, the mean square dispalecement of this particle is given by:
$$
\mathbb{E}\left[\left| \mathbf{r}(t) - \mathbf{r}_j\right|^2\right]= \sum_{i=1}^N G(i,j,t-t_0)\left|\mathbf r_i - \mathbf r_j \right|^2
$$
Then, we can find then ensemble average by introducing an initial state vector $\mathbf X$ whose $i^{th}$ entry is $1$ if there is a particle in lattice site $i$, and $0$ otherwise. 
$$
\begin{aligned}
\langle \mathbf{r}(t) - \mathbf{r}_\mathbf{0}\rangle ={}& \frac{1}{\sum_iX_i}\sum_{j=1}^N \mathbb{E}\left[\left| \mathbf{r}(t) - \mathbf{r}_j\right|^2\right] X_j\\
={}& \frac{1}{\sum_iX_i}\sum_{j=1}^N \sum_{i=1}^N G(i,j,t)\left|\mathbf r_i - \mathbf r_j \right|^2 X_j\\
={}& 2Dt
\end{aligned}
$$
where we set $t_0=0$ WLOG. We can find the diffusion coefficient by differentiating the equation wrt time. Therefore,
$$
\begin{aligned}
D
={}& \frac{1}{2\sum_iX_i}\sum_{j=1}^N \sum_{i=1}^N \partial_tG(i,j,t)\left|\mathbf r_i - \mathbf r_j \right|^2 X_j\\
\approx{}& \frac{1}{2\Delta t\sum_iX_i}\sum_{j=1}^N \sum_{i=1}^N \left(G(i,j,t+\Delta t)-G(i,j,t)\right)\left|\mathbf r_i - \mathbf r_j \right|^2 X_j.
\end{aligned}
$$
Before moving on, it is beneficial to investigate the finite difference approximation of the Green's function:
\begin{aligned}
\mathbf G(t+\Delta t)-\mathbf G(t) 
={}& \left(\mathbf T - \mathbb{I}\right) \mathbf{T}^{t/\Delta t}\\
={}&\left(\frac{1}{2^{\Delta t / \tau_{\mathrm{hl}}}}\mathbb{I} + \left(1-\frac{1}{2^{\Delta t/ \tau_{\mathrm{hl}}}}\right)\mathbf{H}-\mathbb{I}\right)\left(\frac{1}{2^{\Delta t/ \tau_{\mathrm{hl}}}}\mathbb{I} + \left(1-\frac{1}{2^{\Delta t/ \tau_{\mathrm{hl}}}}\right)\mathbf{H}\right)^{t/\Delta t}\\
={}&\frac{1}{2^{ t/ \tau_{\mathrm{hl}}}}\left(1-\frac{1}{2^{\Delta t/ \tau_{\mathrm{hl}}}}\right)
\left( \mathbf{H}-\mathbb{I}\right)\left(\mathbb{I} + \left(2^{\Delta t/ \tau_{\mathrm{hl}}}-1\right)\mathbf{H}\right)^{t/\Delta t}\\
\approx{} & 2^{-t/ \tau_{\mathrm{hl}}}\ln(2) \frac{\Delta t}{ \tau_{\mathrm{hl}}} \left( \mathbf{H}-\mathbb{I}\right)
\left(\mathbb{I} +\ln(2)\frac{\Delta t}{ \tau_{\mathrm{hl}}} \mathbf{H}\right)^{t/\Delta t}\\
\approx{} & 2^{-t/ \tau_{\mathrm{hl}}}\ln(2) \frac{\Delta t}{ \tau_{\mathrm{hl}}} \left( \mathbf{H}-\mathbb{I}\right)
\exp\left(\frac{\ln(2)}{\tau_{\mathrm{hl}}}\mathbf H t\right)\\
 ={} &\Delta t \frac{\ln(2)}{ \tau_{\mathrm{hl}}} \left( \mathbf{H}-\mathbb{I}\right)
\exp\left(\frac{\ln(2)}{\tau_{\mathrm{hl}}}\left( \mathbf{H}-\mathbb{I}\right) t\right)\\
\end{aligned}
where $\mathbf H$ is the hopping matrix whose elements are given by:
$$
H_{ij} = \left(1-\delta_{ij}\right)\frac{\mathcal{N}_i}{\left|\mathbf r_i - \mathbf r_j \right|^6},
$$
and $\mathcal{N}_i$ are chosen such that $\sum_jH_{ij}=1$. Then,

$$
\begin{aligned}
D
\approx{}& \frac{\ln(2)}{2\tau_{\mathrm{hl}}\sum_iX_i} \sum_{i\neq j} H_{ij} \left|\mathbf r_i - \mathbf r_j \right|^{2} X_j.
\end{aligned}
$$


![](https://i.imgur.com/jZcxDq7.png)

# Microscopic Theory
Dear reader, remember the three theory options I gave a section ago... Yeah... I lied. The real theory will be much more complicated (I hope I am wrong). So

>Come with me
And you'll be
In a world of
Pure imagination
Take a look
And you'll see
Into your imagination

Let's start with the facts:

## Crystal Structure
![](https://i.imgur.com/ndHpxj0.png)
Hexagonal lattice with 2 point basis $\mathbf{e}_1=a(1,0,0)^T$, $\mathbf{e}_2=a(-\sqrt{3}/2, 1/2, 0)^T$, $\mathbf{e}_3=b(0,0,1)^T$, where $a=5.9148\times 10^{-10}\mathrm{m}$ and $b=3.4960\times 10^{-10}\mathrm{m}$. The two points in the lattice are 1a and 1f whose positions are given by $\mathbf{0}$ and $2/3\, \mathbf{e}_1 + 1/3\, \mathbf{e}_2 + 1/2\, \mathbf{e}_3$ respectively.

Furthermore, for doping concentration $r_d \in [0,1]$, the filling probability is $r_d$ for site 1a and $r_d/2$ for site 1f. Notice that this means the lattice is never fully filled.

In our last discussion, we claimed that we can treat the unfilled sites as defects and introduce a disorder potential. This approach is perturbative and it will fail even in 100% doped case. The vast number of non-doped sites makes our system non-perturbative. Therefore, we need a non-perturbative theory to deal with the system. To construct that theory, we first need to think about single site states.

## Local Hilbert Space

Let's focus on a single site. We have 3 options to consider: The first option is the site is not doped, i.e., it is empty; the second option is the site hosts an Erbium ion, i.e., it is a 2-level system; and finally the third option is the site host a Ytterbium ion, i.e., it is a 3-level system. Luckily, we can capture all three flavors with only one *generalized* local Hilbert space. If we consider empty option as a 1-level system, the corresponding local Hilbert space for each option is given by:


| System | $\mathcal H_s$ |
| -------- | -------- | 
| Undoped     | $\mathbb C$     |
| $\mathrm{Yb}^{3+}$     | $\mathbb C^2$     |
| $\mathrm{Er}^{3+}$      | $\mathbb C^3$     |

We can describe all the spaces by the generalized space $\mathcal H \cong \mathbb C \oplus \mathbb C \oplus \mathbb C \cong \mathbb C^3$. To get individual spaces, we can project $\mathcal H$ onto any subspace $\mathcal H_s$. Therefore, we have a 3-level system to deal with. Individual levels can be described by:


| State | Notation |Energy | Undoped | $\mathrm{Yb}^{3+}$ | $\mathrm{Er}^{3+}$ |
| -------- | -------- | -------- | -------- | -------- | -------- |
| Ground     | $\ket{-1}$      |0| $+$     | $+$     | $+$     |
| $1^\mathrm{st}$ excited     | $\ket{0}$      |$\epsilon$| $-$     | $+$     | $+$     |
| $2^\mathrm{nd}$ excited     | $\ket{1}$      |$2\epsilon$| $-$     | $-$     | $+$     |

where the last three columns describe if the state is available to individual subsystems. Following above, we can define the projection operators canonically:
\begin{aligned}
P_e ={}& \proj{-1},\\
P_\mathrm{Yb} ={}& \proj{-1} + \proj{0}, \\
P_\mathrm{Er} ={}&1\\
\end{aligned}
where $P_S$ is the projection operator of subsystem $S$ which is Hermitian, i.e., $P_S^\dagger = P_S$. Therefore any operator, $O$, can be projected into a subspace using the formula:
$$
O_S = P_S^\dagger O P_S = P_S O P_S.
$$

Furthermore, our choice of notation is not a coincidence. 
>We'll begin
With a spin
Traveling in
The world of my creation
What we'll see
Will defy
Explanation

The Hilbert space $\mathcal H$ is isomorphic to a spin-1 particle's Hilbert space. Then, we can construct the local Hamiltonian of the system by making use of spin operators: 

$$
H = \epsilon \left( \frac{1}{\hbar}S_z + 1 \right)
$$
where
\begin{aligned}
S_z\ket{n} ={}& \hbar n\ket{n}\\
S_+\ket{n} ={}& \hbar \sqrt{2-n(n+1)}\ket{n+1}\\
S_-\ket{n} ={}& \hbar \sqrt{2-n(n-1)}\ket{n-1}\\
S_+S_-\ket{n} ={}& \hbar^2 \left(2-n(n-1)\right)\ket{n} = 2\hbar^2\left(1-P_e\right)\ket{n}\\
S_-S_+\ket{n} ={}& \hbar^2 \left(2-n(n+1)\right)\ket{n} = 2\hbar ^2P_\mathrm{Yb}\ket{n}\\
\end{aligned}
## The Total System
:::warning
TODO: Write an intro here
:::
$$
H = 
-\frac{1}{2\hbar^2}\sum_{i\neq j}\sum_{\sigma,\sigma^\prime}t_{i,\sigma;j,\sigma^\prime}S_+^{(i,\sigma)} S_-^{(j,\sigma^\prime)}
+\epsilon\sum_i\sum_{\sigma}\left( \frac{1}{\hbar}S_z^{(i,\sigma)} + 1 \right)
$$
where the Latin index denotes the index of the unit cell and Greek index denotes the site (1a, 1f) in the unit cell, and $t_{i,\sigma;j,\sigma^\prime} = t^*_{j,\sigma^\prime ;i,\sigma}$. 

Furthermore, for states $\ket{\psi}$ and $\ket{\phi}$, the transition rate, $\Gamma_{\psi\rightarrow\phi}$ is proportional to the squared norm of the Hamiltonian matrix element, i.e.,
$$
\Gamma_{\psi\rightarrow\phi} \propto \left|\mel{\phi}{H}{\psi} \right|^2.
$$
Combining this with our knowledge of hopping rate scale, we can simplify the Hamiltonian:
$$
H = 
-\frac{t}{2\hbar^2}\sum_{i\neq j}\sum_{\sigma,\sigma^\prime}
\frac{
S_+^{(i,\sigma)} S_-^{(j,\sigma^\prime)}
}{
\left|\mathbf r_{i,\sigma} - \mathbf r_{j,\sigma^\prime} \right|^3
}
+\epsilon\sum_i\sum_{\sigma}\left( \frac{1}{\hbar}S_z^{(i,\sigma)} + 1 \right)
$$

::: info
Note that the kinetic term (the first term) knows nothing about the identity (undoped/Er/Yb) of the sites. 
:::

### The Lattice

We can encode all the lattice info into a projection operator. To do this, we take the operators $P_e,\, P_
\mathrm{Yb},\, P_\mathrm{Er}$ from the previous section, and promote them to local projectors, $P_e^{(i,\sigma)},\, P_
\mathrm{Yb}^{(i,\sigma)},\, P_\mathrm{Er}^{(i,\sigma)}$ which only act on unit cell $i$ and site $\sigma$.

Now, we would like to construct a mathematical structure to describe our doped lattice. The macroscopic properties that we can measure about the doped lattice are the number of unit cells, $N$, Er concerntration $r_\mathrm{Yb} = \frac{2N_\mathrm{Yb}}{3N}$, and Yb concerntration $r_\mathrm{Er} = \frac{2N_\mathrm{Er}}{3N}$ where $N_\mathrm{Yb}$ and $N_\mathrm{Er}$ are number of Er and Yb ions respectively.

We can describe the microscopic properties of the doped lattice by using two permutation functions, $\alpha_a,\alpha_f \colon \quot{\mathbb Z}{N\mathbb Z} \barrow \quot{\mathbb Z}{N\mathbb Z}$. The way to do it is to construct a projection operator that projects the generalized Hilbert space $\mathcal{H}^{2N}$ into the Hilbert space of one particular realization of a doped crystal.


$$
P_a(\alpha_a) =
\left( 
\prod_{i=1}^{{r_\mathrm{Yb}N}} P_\mathrm{Yb}^{(\alpha_a(i),a)} 
\right)
\left( 
\prod_{i=1+{r_\mathrm{Yb}N}}^{{(r_\mathrm{Er}+r_\mathrm{Yb})N}} P_\mathrm{Er}^{(\alpha_a(i),a)} 
\right)
\left( 
\prod_{i=1+{(r_\mathrm{Er}+r_\mathrm{Yb})N}}^{N} P_\mathrm{e}^{(\alpha_a(i),a)} 
\right),
$$

$$
P_f(\alpha_f) =
\left( 
\prod_{i=1}^{{r_\mathrm{Yb}N}} P_\mathrm{Yb}^{(\alpha_f(i),f)} 
\right)
\left( 
\prod_{i=1+{r_\mathrm{Yb}N}}^{{(r_\mathrm{Er}+r_\mathrm{Yb})N}} P_\mathrm{Er}^{(\alpha_f(i),f)} 
\right)
\left( 
\prod_{i=1+{(r_\mathrm{Er}+r_\mathrm{Yb})N}}^{N} P_\mathrm{e}^{(\alpha_f(i),f)} 
\right),
$$

$$
P_\alpha = P_a(\alpha_a) P_f(\alpha_f),
$$

Therefore a doped lattice can be described by the tuple $(N,r_\mathrm{Yb},r_\mathrm{Er},\alpha_a,\alpha_f)$ where the first three elements are macroscopic and last two are microscopic quantities. In addition, the following identities hold
\begin{aligned}
\tr{P_\alpha} ={}& 2N_\mathrm{Er} + N_\mathrm{Yb} + 2N \\
\tr{P_1 P_\alpha} ={}& N_\mathrm{Er} \\
\tr{P_0 P_\alpha} ={}& N_\mathrm{Er}+N_\mathrm{Yb}  \\
\tr{P_{-1} P_\alpha} ={}& 2N
\end{aligned}
where
\begin{aligned}
P_1  ={}& \prod_{i=1}^{N}\prod_{\sigma\in\{a,f\}}\left(P^{(i,\sigma)}_\mathrm{Er}-P^{(i,\sigma)}_\mathrm{Yb}\right),\\
P_0  ={}& \prod_{i=1}^{N}\prod_{\sigma\in\{a,f\}}\left(P^{(i,\sigma)}_\mathrm{Yb}-P^{(i,\sigma)}_e\right),\\
P_{-1} ={}& \prod_{i=1}^{N}\prod_{\sigma\in\{a,f\}}P^{(i,\sigma)}_e.
\end{aligned}
The projected Hamiltonian is
$$
H_\alpha = P^\dagger_\alpha H P_\alpha = P_\alpha H P_\alpha,
$$

from this definition $\lb{H_\alpha}{P_\alpha}=0$ immediatelly follows (use $P^2_\alpha = P_\alpha$) which ensures the constraint we put on the Hilbert space is a consistent one, i.e., the constraint does not evolve in time.
### Lattice Average
Due to the linearity of the EoMs, we can take a lattice average of the Hamiltonian to define an effective Hamiltonian:
$$
H_{eff} = \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N} P_\alpha H P_\alpha,
$$
where $S_N$ denotes the permutation group of order $N!$. This sum is particularly difficult to handle analytically, in contrast it is straigtforward with computer aid. 

To gain physical insight, we propose an an expansion of the projectors such that $P_\alpha = \avg{P} + \delta P_\alpha$. We can justify the expansion around a common _effective_ projector by stating the family of projectors, $P_\alpha$, share the same macroscopic properties; $N,r_\mathrm{Yb},r_\mathrm{Er}$. Still, this proposal is a bold one. Let's see if it works and captures the macroscopic properties as we claimed:

\begin{aligned}
H_{eff} ={}& \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N} P_\alpha H P_\alpha\\
={}& \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}
\left(\avg{P} + \delta P_\alpha\right) 
H 
\left(\avg{P} + \delta P_\alpha\right)\\
={}& \avg{P}H\avg{P} + \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}\left[
\avg{P}
H 
\delta P_\alpha
+
\delta P_\alpha
H 
\avg{P} 
+
\delta P_\alpha
H 
\delta P_\alpha
\right],\\
\end{aligned}
up to this step, we did not impose anything about the expansion. In other words, $\avg{P}$ can be chosen as an arbitrary operator and the expansion would still be valid. Now, we will break this property, and demand the the part of the expansion that is linear in $\avg{P}$ has to vanish, i.e, it is the operator that contributes the most to the lattice average. Then, our demand yields
\begin{aligned}
0 ={}& \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}\left[
\avg{P}
H 
\delta P_\alpha
+
\delta P_\alpha
H 
\avg{P} 
\right]\\
={}& \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}\left[
\avg{P}
H 
\left(P_\alpha - \avg{P}\right)
+
\left(P_\alpha - \avg{P}\right)
H 
\avg{P} 
\right]\\
={}&
-2\avg{P}H\avg{P}+
\frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}\left[
\avg{P}
H 
P_\alpha 
+
P_\alpha 
H 
\avg{P} 
\right]\\
\Leftrightarrow \avg{P} ={}& \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}P_\alpha.
\end{aligned}
Therefore, $\avg{P}$ really is the average of the projectors, and the effective Hamiltonian becomes
$$
H_{eff}=\avg{P}H\avg{P} + \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}
\delta P_\alpha
H 
\delta P_\alpha .
$$
:::danger
$\avg{P}$ is not a real projector since $\avg{P}^2\neq\avg{P}$.
:::

Furthermore, for operators $P_s$ where $s = 1,0,-1$, we get
$$
\tr{P_s\avg{P}}= \frac{1}{(N!)^2}\sum_{\alpha\in S_N\times S_N}\tr{P_sP_\alpha} = \tr{P_sP_\alpha},
$$
then we conclude,
$$
\tr{P_s\delta P_\alpha}=0.
$$
:::info
Can we find a nice basis for $\delta P_\alpha$ that can quantify the heterogeneity? I hope so :eyes: 
:::

### The nice basis
Let's start with the trivial basis:
\begin{aligned}
\Sigma^{(i,\sigma)}_{-1} ={}&-\frac{1}{2\hbar^3}S_z^{(i,\sigma)}S_-^{(i,\sigma)}S_+^{(i,\sigma)}\prod_{j\neq i}\prod_{\sigma^\prime\neq\sigma} \varnothing^{(j,\sigma^\prime)}\\
\Sigma^{(i,\sigma)}_{0} ={}&
\left(1-\frac{1}{\hbar^2}S_z^{(i,\sigma)}S_z^{(i,\sigma)}\right)
\prod_{j\neq i}\prod_{\sigma^\prime\neq\sigma} \varnothing^{(j,\sigma^\prime)}\\
\Sigma^{(i,\sigma)}_{1} ={}&\frac{1}{2\hbar^3}S_z^{(i,\sigma)}S_+^{(i,\sigma)}S_-^{(i,\sigma)}\prod_{j\neq i}\prod_{\sigma^\prime\neq\sigma} \varnothing^{(j,\sigma^\prime)}
\end{aligned}
which satisfy the orthogonality relation:
$$
\Sigma^{(i,\sigma)}_{a}\Sigma^{(j,\lambda)}_{b} = \delta_{ij}\delta_{\sigma\lambda}\delta_{ab}\Sigma^{(i,\sigma)}_{a},
$$
in addition, the basis operators satisfy $\tr{\Sigma^{(i,\sigma)}_{a}}=1$. therefore we can write any operator, $D$, that is diagonal in $S_z$ basis as
$$
D = \sum_i\sum_\sigma\sum_a d^{(i,\sigma)}_{a}\Sigma^{(i,\sigma)}_{a}
$$
and $d^{(i,\sigma)}_{a}$ is given by:
$$
d^{(i,\sigma)}_{a} = \tr{\Sigma^{(i,\sigma)}_{a} D}.
$$
:::warning
Here I will change my notation a little bit. Generally, I dislike this kind of inconsistencies but this document is a flow of thoughts and I have just realized I have to make this change to introduce my _nice_ basis. for any labelled operator $O^{(i,\sigma)}$, I upgrade my notation $O^{(i,\sigma)}\rightarrow O^{(\mathbf i,\sigma)}$ where $i \rightarrow \mathbf i=i_1\mathbf e_1 + i_2\mathbf e_2 + i_3\mathbf e_3$ became a 3D vector. Notice that this does not change any physics, previously we were indexing unit cells with integers, now we are indexing with vectors, which gives us more insight about where the unit cell is.
:::
:::info
Reminder to myself:
construct a reciprocal lattice kind of structure where the wavelenghts are multiples of lattice constants, and use a fourier type of coefficients to construct another basis from the trivial basis, the "energy" of the "wave" will be a measure of heterogeneity. Then we can expand $\delta P_\alpha$ in this basis and include increasing heterogeneity perturbatively.

also prove that this procedure gives a complete basis, find the orthogonality relation.
:::
## Dynamics

$$
\partial_t n(t,\mathbf r) = \frac{i}{\hbar} \mathrm{Tr}\left\{ \rho
\lb{H_\alpha}{\mathcal N(\mathbf r)} 
\right\}  + \mathrm{Tr}\left\{\rho\mathcal D(\mathcal N(\mathbf r))\right\}
$$
where




$$
\mathcal N(\mathbf r) =\sum_{i,\sigma}\delta\left(\mathbf r - \mathbf r_{(i,\sigma)}\right)\left(\frac{1}{\hbar}S_z^{(i,\sigma)}+1\right)
$$

$$
\mathcal N_{-1}(\mathbf r) =-\frac{1}{2\hbar^3}\sum_{i,\sigma}\delta\left(\mathbf r - \mathbf r_{(i,\sigma)}\right)
S_z^{(i,\sigma)}S_-^{(i,\sigma)}S_+^{(i,\sigma)}
$$

$$
\mathcal N_{0}(\mathbf r) =\sum_{i,\sigma}\delta\left(\mathbf r - \mathbf r_{(i,\sigma)}\right)
\left(1-\frac{1}{\hbar^2}S_z^{(i,\sigma)}S_z^{(i,\sigma)}\right)
$$

$$
\mathcal N_{1}(\mathbf r) =\frac{1}{2\hbar^3}\sum_{i,\sigma}\delta\left(\mathbf r - \mathbf r_{(i,\sigma)}\right)
S_z^{(i,\sigma)}S_+^{(i,\sigma)}S_-^{(i,\sigma)}
$$



$$
\mathcal D(\rho) = \gamma \sum_{i,\sigma}\left( L_{i,\sigma}\rho L_{i,\sigma}^\dagger - \frac{1}{2}\left\{L^\dagger_{i,\sigma}L_{i,\sigma},\rho\right\}  \right)
$$
where
$$
L_{i,\sigma} = \frac{1}{2\hbar^2}S_-^{(i,\sigma)}S_-^{(i,\sigma)}
$$


$$
\mathrm{Tr}\left\{\rho\mathcal D(\mathcal N(\mathbf r))\right\} = -\gamma\tr{\rho\mathcal N_1} = -\gamma n_1(t,\mathbf r)
$$

### Conservative term
\begin{aligned}
\lb{H_{eff}}{\mathcal O} ={}& \lb{\avg{P}H\avg{P}}{\mathcal O}\\
={}& \avg{P}\lb{H\avg{P}}{\mathcal O} + \lb{\avg{P}}{\mathcal O}H\avg{P}\\
={}& \avg{P}H\lb{\avg{P}}{\mathcal O} + \avg{P}\lb{H}{\mathcal O}\avg{P} + \lb{\avg{P}}{\mathcal O}H\avg{P}\\
\end{aligned}


Is $\mathcal N_1$ Hermitian?

\begin{aligned}
\left(S_z S_+ S_-\right)^\dagger ={}& S_-^\dagger S_+^\dagger S_z^\dagger\\
={}& S_+ S_- S_z\\
={}& S_+ \left(S_z S_- - \lb{S_z}{S_-}\right)\\
={}& S_+ S_z S_- + \hbar S_+ S_-\\
={}& \left(S_z S_+ - \lb{S_z}{S_+} \right) S_- + \hbar S_+ S_-\\
={}& S_z S_+S_- - \hbar S_+S_-+ \hbar S_+ S_-\\
={}& S_z S_+S_- 
\end{aligned}

yes it is :smile: 

\begin{aligned}
\lb{P_e^{(j,\alpha)}}{\mathcal N_{-1}} ={}& 0
\end{aligned}
## Useful Identities

\begin{aligned}
L^\dagger_{i,\sigma}L_{i,\sigma} ={}&
\frac{1}{2\hbar^3}S_+^{(i,\sigma)}S_z^{(i,\sigma)}S_-^{(i,\sigma)}
+\frac{1}{2\hbar^2}S_+^{(i,\sigma)}S_-^{(i,\sigma)}\\
={}& \frac{1}{2\hbar^3}S_z^{(i,\sigma)}S_+^{(i,\sigma)}
S_-^{(i,\sigma)}\\
L^\dagger_{i,\sigma}S_z^{(j,\alpha)}L_{i,\sigma} ={}&
\left(S_z^{(j,\alpha)}-2\hbar\delta_{ij}\delta_{\sigma\alpha} \right)
L^\dagger_{i,\sigma}L_{i,\sigma}
\end{aligned}

\begin{aligned}
\lb{P_e^{(j,\alpha)}}{S_z^{(i,\sigma)}S_z^{(i,\sigma)}} ={}& 
S_z^{(i,\sigma)}\lb{P_e^{(j,\alpha)}}{S_z^{(i,\sigma)}} + \lb{P_e^{(j,\alpha)}}{S_z^{(i,\sigma)}}S_z^{(i,\sigma)}\\
={}& 
\alb{S_z^{(i,\sigma)}}{\lb{P_e^{(j,\alpha)}}{S_z^{(i,\sigma)}} }
\end{aligned}

\begin{aligned}
\lb{P_{-1}^{(j,\alpha)}}{S_z^{(i,\sigma)}} 
={}& - \frac{1}{2\hbar^3}\delta_{ij}\delta_{\sigma\alpha}
\lb{S_zS_-S_+}{S_z}\\
={}& - \frac{1}{2\hbar^3}\delta_{ij}\delta_{\sigma\alpha}
\left(
S_zS_-\lb{S_+}{S_z} + \lb{S_z}{S_z} S_-S_+ + S_z\lb{S_-}{S_z} S_+
\right)\\
={}&  \frac{1}{2\hbar^3}\hbar\delta_{ij}\delta_{\sigma\alpha}
\left( 
S_zS_-S_+ - S_zS_-S_+
\right)\\
={}&\hbar\delta_{ij}\delta_{\sigma\alpha}
\left( 
P^{(i,\sigma)}_{1} + P^{(i,\sigma)}_{-1}
\right)\\
={}&\hbar\delta_{ij}\delta_{\sigma\alpha}
S_z^{(i,\sigma)} S_z^{(i,\sigma)} 
\\
\end{aligned}

\begin{aligned}
\lb{P_{0}^{(j,\alpha)}}{S_z^{(i,\sigma)}} 
={}& 0
\end{aligned}

\begin{aligned}
\lb{P_{1}^{(j,\alpha)}}{S_z^{(i,\sigma)}} 
={}& \frac{1}{2\hbar^3}\delta_{ij}\delta_{\sigma\alpha}
\lb{S_zS_+S_-}{S_z}\\
={}& \frac{1}{2\hbar^3}\delta_{ij}\delta_{\sigma\alpha}
\lb{S_z
\left(\lb{S_+}{S_-}+S_-S_+\right)
}{S_z}\\
={}& \frac{1}{2\hbar^3}\delta_{ij}\delta_{\sigma\alpha}
\lb{S_z
\lb{S_+}{S_-}
}{S_z}
+
\frac{1}{2\hbar^3}\delta_{ij}\delta_{\sigma\alpha}
\lb{S_z
S_-S_+
}{S_z}\\
={}& \frac{1}{\hbar^2}\delta_{ij}\delta_{\sigma\alpha}
\lb{S_z
S_z
}{S_z} -
\lb{P_{-1}^{(j,\alpha)}}{S_z^{(i,\sigma)}}\\
={}&  -
\lb{P_{-1}^{(j,\alpha)}}{S_z^{(i,\sigma)}}\\
\end{aligned}

\begin{aligned}
\lb{S_+^{(i,\sigma)} S_-^{(j,\sigma^\prime)}}{S_z^{(k,\alpha)}}
={}& S_+^{(i,\sigma)}\lb{ S_-^{(j,\sigma^\prime)}}{S_z^{(k,\alpha)}} + \lb{S_+^{(i,\sigma)} }{S_z^{(k,\alpha)}}S_-^{(j,\sigma^\prime)}\\
={}& \hbar \delta_{jk}\delta_{\sigma^\prime\alpha} S_+^{(i,\sigma)}S_-^{(j,\sigma^\prime)} - 
\hbar  \delta_{ik}\delta_{\sigma\alpha}S_+^{(i,\sigma)} S_-^{(j,\sigma^\prime)}\\
={}& \hbar
\left( \delta_{jk}\delta_{\sigma^\prime\alpha} - \delta_{ik}\delta_{\sigma\alpha}\right)
 S_+^{(i,\sigma)}S_-^{(j,\sigma^\prime)} \\
\end{aligned}

