Here we show the code to calculate transition matrix elements of a $d^4$ system like the one explored in chapter 5 for the ruthenate Ca$_2$RuO$_4$ (CRO). We calculate the possible transitions of the type $d^4d^4\rightarrow d^3d^5$, applying some knowledge regarding the allowed ground and excited states of the CRO system. The transition matrix elements are then used to calculate the amplitude of certain transitions corresponding to the low energy optical excitations found in CRO.

With this in mind, we first explore a small example to test the logic of the computation. Then, we proceed to calculate the transition matrix elements of interest.

For a two sites $(i,j)$ and two levels $\alpha, \beta$ system with one electron per site, the initial states are $d_i^1d_j^1$ configurations where one electron sits on the $\alpha$ level on each site ($4$ in total). The final states are $d_i^0d_j^2$ configurations ($5$ in total).

The representation of the states used is as follows: each level/orbital is composed of two elements of an array, one for each spin (up, down). This is to say that each element of the array corresponds to a creation operator $c^\dagger_{\alpha\sigma}$ when viewed from a second quantization notation. Thus, the first two elements correspond to the $+$ and $-$ spin $\alpha$ level and the following two to the $+$ and $-$ spin $\beta$ level. To build the representation of the states of the combined two sites, we concatenate two single site arrays. This leads to an array where the first 4 elements correspond to the $i$ site and the second 4 ones correspond to the $j$ one.

For example, the state with spin up on both sites at the $\alpha$ level, equivalent to $c^\dagger_{i\alpha\uparrow}c^\dagger_{j\alpha\uparrow}$, is `[1 0 0 0 1 0 0 0]`

. As a second example, we take the configuration where there is one spin up in each level at site $j$ `[0 0 0 0 1 0 1 0]`

. For the remainder of this introduction, we will label those two configurations as `initial`

and `final`

, respectively.

```
In [1]:
```initial = [1, 0, 0, 0, 1, 0, 0, 0]
final = [0, 0, 0, 0, 1, 0, 1, 0]

```
In [2]:
```# Importing necessary extensions
import numpy as np
import itertools
import functools
import operator
# The use of type annotations requires Python 3.6 or newer
from typing import List

The question is, whether there is a transition matrix element between the aforementioned `initial`

and `final`

states. We can easily anser that with a yes, since the receiving level $\beta,+$ is empty and no spin flip is involved when moving from $\alpha,+$. Thus, the question now is how to compute it in a systematic way.

We can start by taking the `XOR`

(exclusive or) operation between the constructed representations of the states. This means that we check where the changes between the states in question are located. Then, we check the positions (index) where we get a `1`

, and if we find that both are odd or both are even, we can say that the transition is allowed. Whereas, if one is in odd and the other in even positions it is not allowed as it would imply a spin flip in the transition. This is equivalent as to write $<f|c^\dagger_{i\alpha'\sigma}c^\dagger_{j\beta'\sigma}|i>$.

Now we can go step by step codifying this rules. First checking for the `XOR`

operation and then asking where the two states differ.

```
In [3]:
```# looking for the positions/levels with different occupation
changes = np.logical_xor(initial, final)
# obtaining the indexes of those positions
np.nonzero(changes)[0].tolist() # .tolist() is called to have a nicer representation

```
Out[3]:
```

```
In [4]:
```np.unique(np.remainder(np.nonzero(changes),2)).size == 1

```
Out[4]:
```

```
In [5]:
```def is_allowed(initial: List[int], final: List[int]) -> bool:
"""Given an initial and final states as represented in a binary list, returns if it is allowed."""
return np.unique(np.remainder(np.nonzero(np.logical_xor(initial,final)),2)).size == 1

`[0 0 0 0 0 1 1 0]`

.

```
In [6]:
```is_allowed(initial, final)

```
Out[6]:
```

```
In [7]:
```is_allowed(initial, [0, 0, 0, 0, 0, 1, 1, 0])

```
Out[7]:
```

Let us first explore the $d^4$ system. In a low spin $d^4$ system, we have only the t$_{2g}$ orbitals (xy yz xz) active leading to a 6 elements representation for a site. Two neighboring states involved in a transition are concatenateed into asingle array consisting of 12 elements.

For this, we create the function to generate the list representation of states given an amount of electrons and levels.

```
In [8]:
```def generate_states(electrons: int, levels: int) -> List[List[int]]:
"""Generates the list representation of a given number of electrons and levels (degeneracy not considered)."""
seed = [1 if position < electrons else 0 for position in range(levels)]
generated_states = list(set(itertools.permutations(seed))) # with set we filter out repeated states
generated_states.sort(reverse=True)
return generated_states

```
In [9]:
```states_d3 = generate_states(3,6)
states_d4 = generate_states(4,6)
states_d5 = generate_states(5,6)

We can consider first the $d^4$ states and take a look at them.

```
In [10]:
``````
states_d4
```

```
Out[10]:
```

It is quite a list of generated states. But from this whoel list, not all states are relevant for the problem at hand. This means that we can reduce the amount of states beforehand by applying the physical constrains we have.

From all the $d^4$ states, we consider only those with a full $d_{xy}$ orbital and those which distribute the other two electrons in different orbitals as possible initial states for the Ca2RuO4 system. In our representation, this means only the states that have a 1 in the first two entries and no double occupancy in $zx$ or $yz$ orbitals.

```
In [11]:
```possible_states_d4 = [list(state) for state in states_d4 # select states that fulfill
if state[0]==1 and state[1]==1 # 1) dxy orbital double occupancy
and state[2] is not state[3] # 2) dzx/dyz orbital single occupancy
]
possible_states_d4

```
Out[11]:
```

We obtain 4 different $d^4$ states that fullfill the conditions previously indicated. From the previous list, the first and last elements correspond to states with $S_z=\pm1$ whereas the ones in the middle correspond to the two superimposed states for the $S=0$ state, namely, a magnon. These four states, could have been easily written down by hand, but the power of this approach is evident when generating and selecting the possible states of the $d^3$ configuration.

For the $d^3$ states, we want at least those which keep one electron in the $d_{xy}$ orbital since we know that the others states are not reachable with one movement as required by optical spectroscopy.

```
In [12]:
```possible_states_d3 = [list(state) for state in states_d3 if state[0]==1 or state[1]==1]
possible_states_d3

```
Out[12]:
```

```
In [13]:
```possible_states_d5 = [list(state) for state in states_d5 if state[0]==1 and state[1]==1]
possible_states_d5

```
Out[13]:
```

We could generate all $d^3d^5$ combinations and check how many of them there are.

```
In [14]:
```def combine_states(first: List[List[int]], second: List[List[int]]) -> List[List[int]]:
"""Takes two lists of list representations of states and returns the list representation of a two-site state."""
# Producing all the possible final states. This has to be read from bottom to top.
final_states = [
functools.reduce(operator.add, combination) # 3) the single site representations are combined into one single two-site representation
for combination # 2) we iterate over all the combinations produced
in itertools.product(first, second) # 1) make the product of the given first and second states lists
]
final_states.sort(reverse=True)
return final_states
print("The number of combined states is: ", len(combine_states(possible_states_d3,possible_states_d5)))

```
```

```
In [15]:
```def transition(initial: List[int], final: List[List[int]], debug = False) -> None:
"""This function takes the list representation of an initial state and a list of
final d3 states. Then, it computes if the transition from the initial state to the
compounded d3d5 final state is possible. The d5 states are implicitly used in
the function from those already generated and filtered.
"""
final_states = combine_states(final, list(possible_states_d5))
print("From initial state {}".format(initial))
# We iterate over all final states and test whether the transition
# from the given initial state is allowed
for state in final_states:
allowed = is_allowed(initial, state)
if allowed:
print(" final state {} is allowed.".format(state))
else:
if debug:
print(" final state {} not allowed".format(state))

With this, we can now explore the transitions between the different initial states and final states ($^4A_2$, $^2E$, and $^2T_1$ multiplets for the $d^3$ sector). Concerning the $d^4$ states, as explained in chapter 5, there is the possibility to be in the $S_z=\pm1$ or $S_z=0$. We will cover each one of them in the following.

First, we will deal with the $^4A_2$ multiplet.

Starting with the pure $S_z=\pm1$ initial states, meaning $d_{\uparrow}^4d_{\uparrow}^4$ (FM) and $d_{\uparrow}^4d_{\downarrow}^4$ (AFM), we have the following representations:

```
In [16]:
```FM = [1,1,1,0,1,0,1,1,1,0,1,0]
AFM_up = [1,1,1,0,1,0,1,1,0,1,0,1]
AFM_down = [1,1,0,1,0,1,1,1,1,0,1,0]

Then, the representation for a $|^4A_2,3/2>$ multiplet, without factors, is given by

```
In [20]:
```A2_32 = [[1,0,1,0,1,0]]
A2_neg_32 = [[0,1,0,1,0,1]]

```
In [22]:
```transition(FM, A2_32)

```
```

Now we proceed with the $|^4A_2,1/2>$ multiplet, which has the following representation:

```
In [33]:
```A2_12 = [[1,0,0,1,0,1],[0,1,1,0,0,1],[0,1,0,1,1,0]] # 4A2 Sz=1/2 d3 multiplet
transition(FM, A2_12)

```
```

Doing the same for the $|^4A_2,-1/2>$ multiplet

```
In [34]:
```A2_neg_12 = [[0,1,1,0,1,0],[1,0,0,1,1,0],[1,0,1,0,0,1]] # 4A2 Sz=-1/2 d3 multiplet
transition(FM, A2_neg_12)

```
```

```
In [37]:
```transition(AFM_up, A2_32)

```
```

We see that the AFM initial state has no transition matrix element for the $|^4A_2,3/2>$ state.

```
In [38]:
```transition(AFM_up, A2_neg_12)

```
```

```
In [39]:
```S0_1 = [1, 1, 1, 0, 0, 1]
S0_2 = [1, 1, 0, 1, 1, 0]
d_zero_down = [1, 1, 0, 1, 0, 1]
d_zero_up = [1, 1, 1, 0, 1, 0]

```
In [40]:
```transition(S0_1 + d_zero_up, A2_32)
transition(S0_2 + d_zero_up, A2_32)

```
```

And for $d^4_\downarrow$ we get

```
In [41]:
```transition(S0_1 + d_zero_down, A2_32)
transition(S0_2 + d_zero_down, A2_32)

```
```

```
In [42]:
```transition(S0_1 + d_zero_up, A2_12)
transition(S0_2 + d_zero_up, A2_12)

```
```

```
In [43]:
```transition(S0_1 + d_zero_up, A2_neg_12)
transition(S0_2 + d_zero_up, A2_neg_12)

```
```

And for completeness we repeat with $d^4_\downarrow$

```
In [44]:
```transition(S0_1 + d_zero_down, A2_12)
transition(S0_2 + d_zero_down, A2_12)
transition(S0_1 + d_zero_down, A2_neg_12)
transition(S0_2 + d_zero_down, A2_neg_12)

```
```

```
In [ ]:
```transition(S0_1 + S0_1, A2_32)
transition(S0_1 + S0_2, A2_32)
transition(S0_2 + S0_1, A2_32)
transition(S0_2 + S0_2, A2_32)

No transitions from the $d^4_0d^4_0$ state to $|^4A_2,3/2>$.

And now repeating the same strategy for the $|^4A_2,1/2>$ multiplet

```
In [ ]:
```transition(S0_1 + S0_1, A2_12)
transition(S0_1 + S0_2, A2_12)
transition(S0_2 + S0_1, A2_12)
transition(S0_2 + S0_2, A2_12)

And $|^4A_2,-1/2>$ multiplet

```
In [ ]:
```transition(S0_1 + S0_1, A2_neg_12)
transition(S0_1 + S0_2, A2_neg_12)
transition(S0_2 + S0_1, A2_neg_12)
transition(S0_2 + S0_2, A2_neg_12)

There are allowed transitions.

The prefactors for the S=0 and the $|^4A_2,1/2>$ states are $\frac{1}{\sqrt{2}}$ and $\frac{1}{\sqrt{3}}$. Multiplying each S=0 state, we get a prefactor of $1/2$ and checking the transitions we get two jumps to $d_{yz}$ and two to $d_{xz}$. [NOT TRUE! ket to recycle symbols later This leads to a transition matrix element of $\frac{1}{2\sqrt{3}}(2t_{xy/yz}+2t_{xy/xz})$ which due to the symmetry we get $\frac{2}{\sqrt{3}}t_{xy/yz}$.]

First we encode the $|^2E,a>$ multiplet

```
In [ ]:
```Ea = [[0,1,1,0,1,0],[1,0,0,1,1,0],[1,0,1,0,0,1]]
transition(AFM_down, Ea)
transition(AFM_up, Ea)
transition(FM_down, Ea)
transition(FM_up, Ea)

For the $|^2E,a>$ multiplet, only transitions from the AFM ground state are possible. Collecting the prefactors we get that the transition matrix element in $-\sqrt{2/3}t_{xy,xz}$ and $-\sqrt{2/3}t_{xy,yz}$ as could be easily checked by hand.

Then, for the $|^2E,b>$ multiplet

```
In [ ]:
```Eb = [[1,0,1,0,0,1],[1,0,0,1,1,0]]
transition(AFM_down, Eb)
transition(AFM_up, Eb)
transition(FM_down, Eb)
transition(FM_up, Eb)

From the $S=\pm1$ initial states, no transitions possible to Eb.

We follow with the situation when considering the $S=0$. In this case, each initial state is decomposed in two parts leading to 4 terms.

```
In [ ]:
```transition(S0_1 + S0_1, Ea)
transition(S0_1 + S0_2, Ea)
transition(S0_2 + S0_1, Ea)
transition(S0_2 + S0_2, Ea)

Each one of the combinations is allowed, thus considering the prefactors of the $S_0$ and $|^2E,a>$ we obtain $\sqrt{\frac{2}{3}}t_{xy,xz}$ and $\sqrt{\frac{2}{3}}t_{xy,yz}$.

Doing the same for $|^2E,b>$

```
In [ ]:
```transition(S0_1 + S0_1, Eb)
transition(S0_1 + S0_2, Eb)
transition(S0_2 + S0_1, Eb)
transition(S0_2 + S0_2, Eb)

Adding all the contributions of the allowed terms we obtain, that due to the - sign in the $|^2E,b>$ multiplet, the contribution is 0.

We sill have to cover the ground state of the kind $d_0^4d_\uparrow^4$. As done previously, we again will split the $d_0^4$ in the two parts.

```
In [ ]:
```S0_1 = [1, 1, 1, 0, 0, 1]
S0_2 = [1, 1, 0, 1, 1, 0]

```
In [ ]:
```transition(S0_1 + d_zero_up, Ea)
transition(S0_2 + d_zero_up, Ea)

Here, both parts of the $S_z=0$ state contribute. Checking the prefactors for $S_z=0$ ($1/\sqrt{2}$) and $|^2E, Ea>$ ($1/\sqrt{6}$) we get a matrix element $\sqrt{\frac{2}{3}}t_{xy/xz}$.

Following for transitions into the $|^2E, Eb>$

```
In [ ]:
```transition(S0_1 + [1, 1, 0, 1, 0, 1], Eb)
transition(S0_2 + [1, 1, 0, 1, 0, 1], Eb)

we see that there is no allowed transition.

```
In [ ]:
```T1_p_xy = [[1,0,1,1,0,0],[1,0,0,0,1,1]]
transition(AFM_down, T1_p_xy)
transition(AFM_up, T1_p_xy)
transition(FM_down, T1_p_xy)
transition(FM_up, T1_p_xy)

And for the $|^2T_1,->$

```
In [ ]:
```T1_n_xy = [[0,1,1,1,0,0],[0,1,0,0,1,1]]
transition(AFM_down, T1_n_xy)
transition(AFM_up, T1_n_xy)
transition(FM_down, T1_n_xy)
transition(FM_up, T1_n_xy)

In this case, there is no possible transition.

```
In [ ]:
```T1_p_xz = [[1,1,1,0,0,0],[0,0,1,0,1,1]]
transition(AFM_down, T1_p_xz)
transition(AFM_up, T1_p_xz)
transition(FM_down, T1_p_xz)
transition(FM_up, T1_p_xz)

```
In [ ]:
```T1_n_xz = [[1,1,0,1,0,0],[0,0,0,1,1,1]]
transition(AFM_down, T1_n_xz)
transition(AFM_up, T1_n_xz)
transition(FM_down, T1_n_xz)
transition(FM_up, T1_n_xz)

$t_{yz,xy}$ and $t_{yz,yz}$ transitions are allowed

```
In [ ]:
```T1_p_yz = [[1,1,0,0,1,0],[0,0,1,1,1,0]]
transition(AFM_down, T1_p_yz)
transition(AFM_up, T1_p_yz)
transition(FM_down, T1_p_yz)
transition(FM_up, T1_p_yz)

```
In [ ]:
```T1_n_yz = [[1,1,0,0,0,1],[0,0,1,1,0,1]]
transition(AFM_down, T1_n_yz)
transition(AFM_up, T1_n_yz)
transition(FM_down, T1_n_yz)
transition(FM_up, T1_n_yz)

$t_{xz,xz}$ and $t_{xz,yz}$ transitions are allowed

We corroborate the straightforward observation that for a ferromagnetic initial state, the excitation to a $^2T_1$ multiplet is forbiden when considering the $\textit{xy}$ orbital order.

Thus, considering the prefactor of the $^2T_1$ multiplets, the hopping amplitudes are $t_{yz,xy}/\sqrt{2}$ and $t_{yz,yz}/\sqrt{2}$ $t_{xz,xz}/\sqrt{2}$, and $t_{xz,yz}/\sqrt{2}$.

Now the challenge of addressing this multiplet when considering the $S=0$ component in the ground state.

```
In [ ]:
```S0_1 = [1, 1, 1, 0, 0, 1]
S0_2 = [1, 1, 0, 1, 1, 0]
T1_p_xz = [[1,1,1,0,0,0],[0,0,1,0,1,1]]

```
In [ ]:
```transition(S0_1 + d_zero_up, T1_p_xz)
transition(S0_2 + d_zero_up, T1_p_xz)
transition(S0_1 + d_zero_up, T1_n_xz)
transition(S0_2 + d_zero_up, T1_n_xz)

```
In [ ]:
```transition(S0_1 + d_zero_down, T1_p_xz)
transition(S0_2 + d_zero_down, T1_p_xz)
transition(S0_1 + d_zero_down, T1_n_xz)
transition(S0_2 + d_zero_down, T1_n_xz)

Thus, for initial states of the type $|0,{\uparrow,\downarrow}>$ considering the singly occupied $\textit{xz}$ multiplet, we obtain transitions involving $t_{yz,xz}$ and $t_{yz,yz}$ for $|0,{\uparrow}>$ and $t_{xz,xz}$ and $t_{xz,yz}$ $|0,{\downarrow}>$.

And at last $d^4_0d^4_0$

```
In [ ]:
```transition(S0_1 + S0_1, T1_p_xz)
transition(S0_1 + S0_2, T1_p_xz)
transition(S0_2 + S0_1, T1_p_xz)
transition(S0_2 + S0_2, T1_p_xz)
print("------------------------")
transition(S0_1 + S0_1, T1_n_xz)
transition(S0_1 + S0_2, T1_n_xz)
transition(S0_2 + S0_1, T1_n_xz)
transition(S0_2 + S0_2, T1_n_xz)

```
In [ ]:
```