Kernel: Python3.6

Neutron scattering spectrum for a superexchange and double exchange model

J. K. Glasbrenner

August 3, 2017

Python boilerplate

For our analysis, we need to import the following Python modules.

In [1]:

%matplotlib inline

import IPython
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np
import sympy
from sympy import Matrix, Rational, SparseMatrix
from sympy.physics.quantum import TensorProduct
import tabulate

import copy
import functools
import operator

Exchange model Hamiltonian

We adopt the following ansatz as our Hamiltonian for a two spin system,

\begin{darray}{rcl} \text{H} = 2 J^{SE} \textbf{S}_{1} \cdot \textbf{S}_{2} - t \sum_{n=0}^{\infty} J_{n}^{DE} \left( S \right) \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{n}, \end{darray}

where the first term corresponds to superexchange as defined in Ref. 1 and the second term to double exchange as defined in Ref. 2.

We are modeling the exchange interaction of a pair of manganese atoms, which in the atomic limit — and ignoring spin-orbit coupling — can be idealized as a system of two identical spin-5/2 particles. It is convenient to work with Pauli matrices when taking a computational approach. The Pauli matrices for a spin-5/2 particle 3 are:

\begin{align*} \mathbf{\sigma}_x &= \frac{1}{2} \begin{pmatrix} 0 & \sqrt{5} & 0 & 0 & 0 & 0 \\ \sqrt{5} & 0 & \sqrt{8} & 0 & 0 & 0 \\ 0 & \sqrt{8} & 0 & \sqrt{9} & 0 & 0 \\ 0 & 0 & \sqrt{9} & 0 & \sqrt{8} & 0 \\ 0 & 0 & 0 & \sqrt{8} & 0 & \sqrt{5} \\ 0 & 0 & 0 & 0 & \sqrt{5} & 0 \end{pmatrix}, \\ \mathbf{\sigma}_y &= \frac{1}{2i} \begin{pmatrix} 0 & \sqrt{5} & 0 & 0 & 0 & 0 \\ -\sqrt{5} & 0 & \sqrt{8} & 0 & 0 & 0 \\ 0 & -\sqrt{8} & 0 & \sqrt{9} & 0 & 0 \\ 0 & 0 & -\sqrt{9} & 0 & \sqrt{8} & 0 \\ 0 & 0 & 0 & -\sqrt{8} & 0 & \sqrt{5} \\ 0 & 0 & 0 & 0 & -\sqrt{5} & 0 \end{pmatrix}, \\ \mathbf{\sigma}_z &= \begin{pmatrix} -5/2 & 0 & 0 & 0 & 0 & 0 \\ 0 & -3/2 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1/2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3/2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 5/2 \end{pmatrix} \end{align*}

The total spin operator is a three-vector as follows:

\begin{align*} \mathbf{\sigma} &= \mathbf{\sigma}_x \hat{x} + \mathbf{\sigma}_y \hat{y} + \mathbf{\sigma}_z \hat{z} \\ &= \begin{pmatrix} \mathbf{\sigma}_x \\ \mathbf{\sigma}_y \\ \mathbf{\sigma}_z \end{pmatrix}. \end{align*}

Let \(\textbf{S}_1\) be the spin operator for particle 1 and \(\textbf{S}_1\) be the spin operator for particle 2. These operators are tensor products as follows:

\begin{align*} \textbf{S}_1 &= \left( \mathbf{\sigma}_x \otimes \textbf{I} \right) \hat{x} + \left( \mathbf{\sigma}_y \otimes \textbf{I} \right) \hat{y} + \left( \mathbf{\sigma}_z \otimes \textbf{I} \right) \hat{z}, \\ \textbf{S}_2 &= \left( \textbf{I} \otimes \mathbf{\sigma}_x \right) \hat{x} + \left( \textbf{I} \otimes \mathbf{\sigma}_y \right) \hat{y} + \left( \textbf{I} \otimes \mathbf{\sigma}_z \right) \hat{z}, \end{align*}

where I is the identity matrix. For shorthand, we define:

\begin{align*} \textbf{S}_{1i} &\equiv \mathbf{\sigma}_i \otimes \textbf{I}, \\ \textbf{S}_{2i} &\equiv \textbf{I} \otimes \mathbf{\sigma}_i. \end{align*}

We then write the rest of the relevant two-spin operators, which are

\begin{align*} \textbf{S}_z &= \textbf{S}_{1z} + \textbf{S}_{2z}, \\ \textbf{S} &= \textbf{S}_1 + \textbf{S}_2, \\ \textbf{S}^2 &= \left( \textbf{S}_{1x} + \textbf{S}_{2x} \right)^2 + \left( \textbf{S}_{1y} + \textbf{S}_{2y} \right)^2 + \left( \textbf{S}_{1z} + \textbf{S}_{2z} \right)^2, \\ \textbf{S}_1 \cdot \textbf{S}_2 &= \textbf{S}_{1x} \cdot \textbf{S}_{2x} + \textbf{S}_{1y} \cdot \textbf{S}_{2y} + \textbf{S}_{1z} \cdot \textbf{S}_{2z}. \end{align*}

We implement the above operators in Python as the following SpinOperators class. This class also solves the eigenvalue problem for the two-spin operators \(\textbf{S}_1^2\), \(\textbf{S}_2^2\), \(\textbf{S}^2\), \(\textbf{S}_z\), and \(2 \textbf{S}_1 \cdot \textbf{S}_2\).

In [2]:

class SpinOperators(object):
    """Defines the one- and two-particle spin operators for a pair of
    identical spin-5/2 particles using Pauli matrices and solves for
    the eigenvalues and eigenvectors.
    
    The __init__ method uses Pauli matrices to explicitly construct
    the  one- and two-spin operators for a spin-5/2 system. It then
    finds the eigenvalues and eigenvectors for the |S1^2|, |S2^2|,
    |S^2|, |Sz|, and |2S1*S2| two-spin operators, and stores them in
    the attributes described below.
    
    Attributes
    ----------
    one_spin_operators : `dict`
        The one-spin operators |σx|, |σy|, |σz|, and |σ| for spin-5/2
        particles as Pauli matrices.
    two_spin_operators : `dict`
        The two-spin operators |S1x|, |S1y|, and |S1z| for
        key=`"particle1"`, |S2x|, |S2y|, and |S2z| for
        key=`"particle2"`, and |S|, |S1^2|, |S2^2|, |S^2|, and |2S1*S2|
        for key=`"system"` constructed with Pauli matrices.
    two_spin_eigensystem : `dict`
        The eigenvalues (key=`"eigenvalues"`) and eigenvectors
        (key=`"eigenvectors"`) for the |S|, |S1^2|, |S2^2|, |S^2|, and
        |2S1*S2| two-spin operators.
        
    Methods
    -------
    show_two_spin_eigenvalues
        Display the two-spin eigenvalues as nicely rendered tables.
        
    .. |σx| replace:: :math:`\mathbf{\sigma}_x`
    .. |σy| replace:: :math:`\mathbf{\sigma}_y`
    .. |σz| replace:: :math:`\mathbf{\sigma}_z`
    .. |σ| replace:: :math:`\mathbf{\sigma}`
    .. |Sx| replace:: :math:`\textbf{S}_x`
    .. |Sy| replace:: :math:`\textbf{S}_y`
    .. |Sz| replace:: :math:`\textbf{S}_z`
    .. |S| replace:: :math:`\textbf{S}`
    .. |S1^2| replace:: :math:`\textbf{S}_1^2`
    .. |S2^2| replace:: :math:`\textbf{S}_2^2`
    .. |S^2| replace:: :math:`\textbf{S}^2`
    .. |2S1*S2| replace:: :math:`2 \textbf{S}_1 \cdot \textbf{S}_2`
    """
    
    def __init__(self):
        
        self._one_spin_operators_list = ["Sx", "Sy", "Sz"]
        self._two_spin_operators_list = [
            "S1^2", "S2^2", "S^2", "Sz", "2S1*S2"
        ]

        self.one_spin_operators = (
            self._construct_one_spin_operators(self._one_spin_operators_list)
        )

        self.two_spin_operators = (
            self._construct_two_spin_operators(
                copy.deepcopy(self.one_spin_operators),
                "5/2",
                self._one_spin_operators_list
            )
        )
        
        self.two_spin_eigensystem = {
            "eigenvalues": self._calculate_two_spin_eigenvalues(
                self.two_spin_operators,
                self._two_spin_operators_list
            ),
            "eigenvectors": self._calculate_two_spin_eigenvectors(
                self.two_spin_operators,
                self._two_spin_operators_list
            )
        }
        
    def show_two_spin_eigenvalues(self):
        """Displays the eigenvalues found in attribute `two_spin_eigensystem`
        as nicely rendered tables.
        """

        output_markdown = []

        eigenvalues = copy.deepcopy(
            self.two_spin_eigensystem["eigenvalues"]
        )

        latex_rep_list = [
            "\\\\(\\textbf{S}_1^2\\\\)",
            "\\\\(\\textbf{S}_2^2\\\\)",
            "\\\\(\\textbf{S}^2\\\\)",
            "\\\\(\\textbf{S}_z\\\\)",
            "\\\\(2 \\textbf{S}_1 \\cdot \\textbf{S}_2\\\\)"
        ]

        
        for spin_op, latex_rep in zip(self._two_spin_operators_list, latex_rep_list):

            output_eigenvalue = (
                self._format_eigenvalues_for_output(eigenvalues[spin_op])
            )

            output_table = (
                self._construct_output_table_for_eigenvalues(output_eigenvalue)
            )
            
            output_markdown.append(
                self._format_markdown_string(latex_rep, output_table)
            )
            
        # Join markdown strings and render output
        IPython.display.display(
            IPython.display.Markdown("\n\n".join(output_markdown))
        )

    @staticmethod
    def _operator_dot_product(matrix1, matrix2):
        
        dot_product = matrix1.conjugate().transpose() * matrix2
        
        return dot_product
        
    @staticmethod
    def _construct_one_spin_operators(one_spin_operators_list):
        
        cls = SpinOperators
        
        one_spin_operators = {
            "5/2": {
                "Sx": (
                    Rational(1, 2) *
                    Matrix(
                        [[0, sympy.sqrt(5), 0, 0, 0, 0],
                         [sympy.sqrt(5), 0, sympy.sqrt(8), 0, 0, 0],
                         [0, sympy.sqrt(8), 0, sympy.sqrt(9), 0, 0],
                         [0, 0, sympy.sqrt(9), 0, sympy.sqrt(8), 0],
                         [0, 0, 0, sympy.sqrt(8), 0, sympy.sqrt(5)],
                         [0, 0, 0, 0, sympy.sqrt(5), 0]])
                ),
                "Sy": (
                    Rational(1, 2) * sympy.I ** (-1) *
                    Matrix(
                        [[0, sympy.sqrt(5), 0, 0, 0, 0],
                         [-sympy.sqrt(5), 0, sympy.sqrt(8), 0, 0, 0],
                         [0, -sympy.sqrt(8), 0, sympy.sqrt(9), 0, 0],
                         [0, 0, -sympy.sqrt(9), 0, sympy.sqrt(8), 0],
                         [0, 0, 0, -sympy.sqrt(8), 0, sympy.sqrt(5)],
                         [0, 0, 0, 0, -sympy.sqrt(5), 0]])
                ),
                "Sz": (
                    Rational(1, 2) *
                    Matrix(
                        [[5, 0, 0, 0, 0, 0],
                         [0, 3, 0, 0, 0, 0],
                         [0, 0, 1, 0, 0, 0],
                         [0, 0, 0, -1, 0, 0],
                         [0, 0, 0, 0, -3, 0],
                         [0, 0, 0, 0, 0, -5]])
                )
            }
        }
        
        one_spin_operators["5/2"]["S"] = (
            cls._construct_s_operator(
                copy.deepcopy(one_spin_operators),
                "5/2",
                one_spin_operators_list
            )
        )
        
        return one_spin_operators
    
    @staticmethod
    def _construct_s_operator(spin_operators,
                              dict_key,
                              one_spin_operators_list):

        s_operator = []        
        for Si in one_spin_operators_list:
            s_operator.append(spin_operators[dict_key][Si])
            
        return SparseMatrix(s_operator)

    @staticmethod
    def _construct_two_spin_operators(one_spin_operators,
                                      spin_quantum_number,
                                      one_spin_operators_list):

        cls = SpinOperators
        nrows, ncols = (
            one_spin_operators[spin_quantum_number]["Sz"].shape
        )
        identity = sympy.eye(nrows)
        
        two_spin_operators = (
            cls._construct_two_particle_tensor_products(one_spin_operators,
                                                        identity,
                                                        spin_quantum_number,
                                                        one_spin_operators_list)
        )

        for spin_i in ["spin1", "spin2"]:
            two_spin_operators[spin_i]["S"] = (
                cls._construct_s_operator(
                    copy.deepcopy(two_spin_operators),
                    spin_i,
                    one_spin_operators_list
                )
            )
            
        for two_spin_op in ["Sz", "S"]:
            two_spin_operators["system"][two_spin_op] = (
                two_spin_operators["spin1"][two_spin_op] +
                two_spin_operators["spin2"][two_spin_op]
            )
            
        for (op1,
             op2,
             prefactor,
             two_spin_op) in [["spin1", "spin1", 1, "S1^2"],
                              ["spin2", "spin2", 1, "S2^2"],
                              ["system", "system", 1, "S^2"],
                              ["spin1", "spin2", 2, "2S1*S2"]]:

            two_spin_operators["system"][two_spin_op] = (
                cls._operator_dot_product(
                    prefactor * two_spin_operators[op1]["S"],
                    two_spin_operators[op2]["S"]
                )
            )
            
        return two_spin_operators
        
    @staticmethod
    def _construct_two_particle_tensor_products(one_spin_operators,
                                                identity_matrix,
                                                spin_quantum_number,
                                                one_spin_operators_list):
        two_spin_operators = {
            "spin1": {},
            "spin2": {},
            "system": {}
        }

        for Si in one_spin_operators_list:
            two_spin_operators["spin1"][Si] = (
                SparseMatrix(
                    TensorProduct(one_spin_operators[spin_quantum_number][Si],
                                  identity_matrix)
                )
            )
            two_spin_operators["spin2"][Si] = (
                SparseMatrix(
                    TensorProduct(identity_matrix,
                                  one_spin_operators[spin_quantum_number][Si])
                )
            )
            
        return two_spin_operators
    
    @staticmethod
    def _calculate_two_spin_eigenvalues(two_spin_operators,
                                        two_spin_operators_list):
        
        eigenvalues = {}

        for spin_op in two_spin_operators_list:
            spin_op_eigenvalues = (
                two_spin_operators["system"][spin_op].eigenvals()
            )
            spin_op_eigenvalues = [
                [eigenvalue, degeneracy]
                for eigenvalue, degeneracy in spin_op_eigenvalues.items()
            ]
            spin_op_eigenvalues = sorted(
                spin_op_eigenvalues, key=operator.itemgetter(0)
            )

            eigenvalues[spin_op] = (
                spin_op_eigenvalues
            )
            
        return eigenvalues
    
    @staticmethod
    def _calculate_two_spin_eigenvectors(two_spin_operators,
                                         two_spin_operators_list):
        
        eigenvectors = {}

        for spin_op in two_spin_operators_list:
            eigenvectors[spin_op] = (
                two_spin_operators["system"][spin_op].eigenvects()
            )
            
        return eigenvectors
   
    @staticmethod
    def _format_eigenvalues_for_output(eigenvalues):
        
        output_eigenvalue = [
            [f"{sympy.nsimplify(eigenvalue, tolerance=0.001, rational=True)}",
             f"{sympy.nsimplify(degeneracy, tolerance=0.001, rational=True)}"]
            for eigenvalue, degeneracy in eigenvalues
        ]
        
        return output_eigenvalue

    @staticmethod
    def _construct_output_table_for_eigenvalues(output_eigenvalue):
        
        output_table = (
            tabulate.tabulate(tabular_data=output_eigenvalue,
                              headers=("Eigenvalue", "Degeneracy"),
                              tablefmt='pipe',
                              numalign='center',
                              stralign='center')
        )

        return output_table

    @staticmethod
    def _format_markdown_string(latex_rep, output_table):
        
        output_markdown = (
            f"{latex_rep} eigenvalues\n\n"
            f"{output_table}"
        )
            
        return output_markdown

Our goal is to model the inelastic neutron scattering (INS) spectra in the figure below

The main features we would like to explain using a model and our DFT calculations are

The peak around 13.3 meV for Ba(Zn_0.85Mn_0.15)₂As₂.
The broadened peak around 13.3 meV for Ba_0.8K_0.2(Zn_0.85Mn_0.15)₂As₂.

The outline for the rest of this notebook is as follows. First, we focus on the superexchange term and derive the energy differences for the \(\Delta{}S_z = \pm 1\) transitions, which are what the INS experiment detects. Next, we focus on the double exchange term, show that we can eliminate the infinite sum for a system of two spin-5/2 particles, and obtain the energy for different eigenvalues of \(\textbf{S}z\). Finally, we combine the results together by assuming that the superexchange and double exchange terms are independent such that the energy contributions from each can be simply summed together. We then use our DFT calculations to fit the superexchange \(J^{SE}\) and double exchange \(t\) parameters. We then assume that the superexchange parameter \(J^{SE}\) remains constant under doping and that we should multiply \(J^{SE}\) and \(t\) by 0.548 and 0.529 respectively to take into account magnetic moment softness. This allows us to get an approximate estimate of how double exchange and moment softness modifies the \(\Delta{}S{z} = \pm 1\) spectra. The notebook concludes with an interactive model of the INS spectra, assuming that Gaussian peaks, whose widths and amplitudes are tunable empirical parameters, are centered around the model's predicted \(\Delta{}S_{z} = \pm 1\) spectra.

Superexchange term

It's well-known that the eigenvalues for the superexchange term \(2 J^{SE} \textbf{S}_1 \cdot \textbf{S}_2\) are given by the expression \(J^{SE} \left[ S_z \left( S_z + 1 \right) - S_1 \left( S_1 + 1 \right) - S_2 \left( S_2 + 1 \right) \right]\). Even though the result is known, there is pedagogical worth in carrying out the analysis and reproducing the result, so let's do that.

The eigenvalue problem is already solved in the SpinOperators class from the previous section, so the reader is referred to the code for the details. The eigenvalues for the different two-particle operators are as follows

In [3]:

SpinOperators().show_two_spin_eigenvalues()

\(\textbf{S}_1^2\) eigenvalues

Eigenvalue	Degeneracy
35/4	36

\(\textbf{S}_2^2\) eigenvalues

Eigenvalue	Degeneracy
35/4	36

\(\textbf{S}^2\) eigenvalues

Eigenvalue	Degeneracy
0	1
2	3
6	5
12	7
20	9
30	11

\(\textbf{S}_z\) eigenvalues

Eigenvalue	Degeneracy
-5	1
-4	2
-3	3
-2	4
-1	5
0	6
1	5
2	4
3	3
4	2
5	1

\(2 \textbf{S}_1 \cdot \textbf{S}_2\) eigenvalues

Eigenvalue	Degeneracy
-35/2	1
-31/2	3
-23/2	5
-11/2	7
5/2	9
25/2	11

Using the above tables, you can quickly confirm that the \(2 \textbf{S}_1 \cdot \textbf{S}_2\) eigenvalues match the predictions of the analytic expression.

The two-spin operators \(\textbf{S}^2\) and \(\textbf{S}_z\) commute with \(\textbf{S}_1 \cdot \textbf{S}_2\), and so we know that the eigenvectors for \(\textbf{S}_1 \cdot \textbf{S}_2\) will diagonalize \(\textbf{S}^2\) and \(\textbf{S}_z\). We can use this to associate the ferromagnetic and antiferromagnetic configurations from our DFT calculations with two eigenvectors from this model.

First, we claim that the following eigenvector can be associated with the antiferromagnetic configuration

In [4]:

# Solve eigenvalue problem
calc_spinoperators = SpinOperators()

In [5]:

# S1⋅S2 singlet eigenvector
calc_afm = calc_spinoperators.two_spin_eigensystem["eigenvectors"]["2S1*S2"][0][2][0].normalized()
print(calc_afm.transpose())

Matrix([[0, 0, 0, 0, 0, -sqrt(6)/6, 0, 0, 0, 0, sqrt(6)/6, 0, 0, 0, 0, -sqrt(6)/6, 0, 0, 0, 0, sqrt(6)/6, 0, 0, 0, 0, -sqrt(6)/6, 0, 0, 0, 0, sqrt(6)/6, 0, 0, 0, 0, 0]])

This eigenvector has the following eigenvalue for \(2 \textbf{S}_1 \cdot \textbf{S}_2\)

In [6]:

# Confirm the eigenvalue
calc_afm.transpose().conjugate() * calc_spinoperators.two_spin_operators["system"]["2S1*S2"] * calc_afm

Matrix([[-35/2]])

the following eigenvalue for \(\textbf{S}^2\)

In [7]:

# Shared basis, so what is the eigenvalue of S^2 with this singlet eigenvector?
calc_afm.transpose().conjugate() * calc_spinoperators.two_spin_operators["system"]["S^2"] * calc_afm

Matrix([[0]])

and the following eigenvalue for \(\textbf{S}_z\)

In [8]:

# Shared basis, so what is the eigenvalue of Sz with this singlet eigenvector?
calc_afm.transpose().conjugate() * calc_spinoperators.two_spin_operators["system"]["Sz"] * calc_afm

Matrix([[0]])

This is in line with what we would expect with an antiferromagnetic state. This state is a singlet, it has the lowest energy, and the eigenvalues for \(\textbf{S}^2\) and \(\textbf{S}_z\) are consistent with a state where half the spins point up and the other half point down.

Next, we claim that the following eigenvector can be associated with the ferromagnetic configuration

In [9]:

# One of the S1⋅S2 eigenvectors with a 25/4 eigenvalue
# The 25/4 eigenvalue is 11-fold degenerate
calc_fm = calc_spinoperators.two_spin_eigensystem["eigenvectors"]["2S1*S2"][5][2][0].normalized()
print(calc_fm.transpose())

Matrix([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

This eigenvector has the following eigenvalue for \(2 \textbf{S}_1 \cdot \textbf{S}_2\)

In [10]:

# Confirm the eigenvalue
calc_fm.transpose().conjugate() * calc_spinoperators.two_spin_operators["system"]["2S1*S2"] * calc_fm

Matrix([[25/2]])

the following eigenvalue for \(\textbf{S}^2\)

In [11]:

# Shared basis, so what is the eigenvalue of S^2 with this particular eigenvector?
calc_fm.transpose().conjugate() * calc_spinoperators.two_spin_operators["system"]["S^2"] * calc_fm

Matrix([[30]])

and the following eigenvalue for \(\textbf{S}_z\)

In [12]:

# Shared basis, so what is the eigenvalue of Sz with this particular eigenvector?
calc_fm.transpose().conjugate() * calc_spinoperators.two_spin_operators["system"]["Sz"] * calc_fm

Matrix([[5]])

In [13]:

# Remove calc_* variables from scope
IPython.get_ipython().magic("reset_selective -f calc_*")

This is in line with what we would expect with a ferromagnetic state. It has the highest energy and the eigenvalues for \(\textbf{S}^2\) and \(\textbf{S}_z\) are consistent with a state where all the spins are pointing in the same direction. During subsequent analysis, we will refer to the state with \(2 S_1 S_2 = -35/2 \), \(S^2 = 0\), and \(S_z = 0\) as the antiferromagnetic state and the state with \(2 S_1 S_2 = 25/2 \), \(S^2 = 30\), and \(S_z = 5\) as the ferromagnetic state.

Now that we have the eigenvalues, we want to find the change in energy for all \(\Delta{}S_z = \pm 1\) transitions. The analytic expression for calculating the transition energy is \(J^{SE} \left[ S_z' \left( S_z' + 1 \right) - S_z \left( S_z + 1 \right) \right]\). This results in the following transition energies

S_z → S^'_z	ΔE
-5 → -4	-8J^SE
-4 → -3	-6J^SE
-3 → -2	-4J^SE
-2 → -1	-2J^SE
-1 → 0	0
0 → 1	2J^SE
1 → 2	4J^SE
2 → 3	6J^SE
3 → 4	8J^SE
4 → 5	10J^SE

Double exchange term

We start with the two spin double exchange Hamiltonian ansatz defined in equation (5.213) in Ref. 2,

\begin{darray}{rcl} \text{H}^{DE} = - t \sum_{n=0}^{\infty} J_{n}^{DE} \left( S \right) \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{n}. \end{darray}

In order to show that this is a valid way to model doube exchange, we need to complete two tasks:

The infinite summation can be exactly replaced with a six-term sum for a system of two spin-5/2 particles.
There is a solution to the six-term reformulated model that yields the double exchange energy \(E^{DE} = -t \dfrac{\left\lvert S_z \right\rvert}{S_1 + S_2 + 1}\), where the eigenvalues \(S_1 = S_2 = 5/2\).

For task 1, we follow the procedure on pages 223-226 of Nolting and Ramakanth's book [2]. For convenience, we define the eigenvalue \(S_1 = S_2 \equiv S_{1,2}\). We start by calculating the values of \(\dfrac{1}{\hbar{}^{2}} \left( \textbf{S}1 \cdot \textbf{S}2 \right) = \dfrac{1}{2} S_z \left(S_z + 1\right) - S{1,2} \left( S{1,2} + 1 \right) \) for all values of \(S\) using the following helper functions,

In [14]:

def product_of_two_identical_spins(s, s_t):
    s1_dot_s2 = sympy.Rational(1, 2) * s_t * (s_t + 1) - s * (s + 1)
    
    return s1_dot_s2


def calculate_spin_products(s, list_of_s_t):
    s_t_and_spin_product_list = [
        [str(s_t),
         str(product_of_two_identical_spins(s, sympy.Rational(s_t)))]
        for s_t in list_of_s_t]
    
    return s_t_and_spin_product_list


def write_s_t_and_spin_product_markdown_table(s_t_and_spin_product_list,
                                              table_headers):
    s_t_spin_product_table = (
        tabulate.tabulate(tabular_data=s_t_and_spin_product_list,
                          headers=table_headers,
                          tablefmt='pipe',
                          numalign='center',
                          stralign='center'))
    
    IPython.display.display(IPython.display.Markdown(
        '{0}'.format(s_t_spin_product_table)))

We then run the calculation to get the following table

In [15]:

def run_s1_dot_s2_calculation():
    s = sympy.Rational(5, 2)
    list_of_s_t = range(0, 6, 1)
    s_t_and_spin_product_list = calculate_spin_products(s, list_of_s_t)
    table_headers = ("*S*",
                     "ℏ<sup>-2</sup> "
                     "(***S***<sub>1</sub> ⋅ ***S***<sub>2</sub>)")
    
    write_s_t_and_spin_product_markdown_table(s_t_and_spin_product_list,
                                              table_headers)


# Run calculation and output Markdown compatible table
run_s1_dot_s2_calculation()

S	ℏ^-2 (S₁ ⋅ S₂)
0	-35/4
1	-31/4
2	-23/4
3	-11/4
4	5/4
5	25/4

Following the procedure of Nolting and Ramakanth, we write down a second ansatz

\begin{darray}{rcl} \frac{1}{\hbar^{12}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{6} = \alpha + \frac{\beta}{\hbar^{2}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right) + \frac{\gamma}{\hbar^{4}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{2} + \frac{\delta}{\hbar^{6}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{3} + \frac{\epsilon}{\hbar^{8}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{4} + \frac{\zeta}{\hbar^{10}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{5}. \end{darray}

If a set of coefficients \(\left( \alpha, \beta, \gamma, \delta, \epsilon, \zeta \right)\) exists for a system of two spin-5/2 particles, then all terms with powers higher than \( \left( \textbf{S}{1} \cdot \textbf{S}{2} \right)^5 \) can be re-expressed as products of the lower-order terms, resulting in a six-term Hamiltonian. To show this, we implement the second ansatz and the method of solving for the coefficients in the functions below:

In [16]:

def build_spin_product_powers_expression(number_of_terms):
    hbar = sympy.var("ℏ")
    s1s2 = sympy.var("S1⋅S2")
    coefficients = [sympy.var(coefficient) for coefficient in
                    ["α",
                     "β",
                     "γ",
                     "δ",
                     "ϵ",
                     "ζ"]]
    lhs = []
    rhs = []
    
    lhs.append(
        (1 / hbar ** (2 * number_of_terms)) *
        (s1s2) ** (number_of_terms)
    )
    
    if len(coefficients) >= number_of_terms:
        for n, coefficient in zip(range(0, number_of_terms, 1),
                                  coefficients[:number_of_terms]):
            rhs.append(
                (coefficient / hbar ** (2 * n)) * (s1s2) ** (n)
            )
    else:
        raise
        
    lhs = functools.reduce(sympy.Add, lhs)
    rhs = functools.reduce(sympy.Add, rhs)
    equality = sympy.Eq(lhs, rhs)
    
    return equality


def substitute_s1s2_in_expression(equality, substitution):
    s1s2 = sympy.var("S1⋅S2")
    
    return equality.subs(s1s2, substitution)


def series_of_equations(equality):
    s = sympy.Rational(5, 2)
    list_of_s_t = range(0, 6, 1)
    s_t_and_spin_product_list = [
        [s_t,
         product_of_two_identical_spins(s, sympy.Rational(s_t))]
        for s_t in list_of_s_t]
    
    equation_series = []
    for s_t_pair in s_t_and_spin_product_list:
        equation_series.append(substitute_s1s2_in_expression(equality, s_t_pair[1]))
    
    return equation_series

As a sanity check, let's test out the above functions against the \(S_{1,2} = 1\) case, where we know that the coefficients should be

coeff	value
α	2
β	1
γ	-2

In [17]:

def run_spin_one_coefficients():

    def series_of_equations(equality):
        s = sympy.Rational(1)
        list_of_s_t = range(0, 3, 1)
        s_t_and_spin_product_list = [
            [s_t,
             product_of_two_identical_spins(s, sympy.Rational(s_t))]
            for s_t in list_of_s_t]

        equation_series = []
        for s_t_pair in s_t_and_spin_product_list:
            equation_series.append(substitute_s1s2_in_expression(equality, s_t_pair[1]))

        return equation_series

    
    threeterms = build_spin_product_powers_expression(3)
    equation_system = series_of_equations(threeterms)
    A, b = sympy.linear_eq_to_matrix(
        [equations.rhs - equations.lhs * sympy.var("ℏ") ** 6 for equations in equation_system],
        [sympy.var(coefficient) / sympy.var("ℏ") ** n for coefficient, n in zip(["α", "β", "γ"], [0, 2, 4])])
    threetermssolution = sympy.linsolve(
        (A, b),
        [sympy.var(coefficient) for coefficient in ["α", "β", "γ"]])
    print(threetermssolution)
    
    
run_spin_one_coefficients()

{(2, 1, -2)}

The coefficients match as expected, so now we apply the same procedure to our ansatz of interest:

In [18]:

def run_spin_five_half_coefficients():
    sixterms = build_spin_product_powers_expression(6)
    equation_system = series_of_equations(sixterms)
    A, b = sympy.linear_eq_to_matrix(
        [equations.rhs - equations.lhs * sympy.var("ℏ") ** 12 for equations in equation_system],
        [sympy.var(coefficient) / sympy.var("ℏ") ** n
         for coefficient, n in zip(
             ["α", "β", "γ", "δ", "ϵ", "ζ"],
             [0, 2, 4, 6, 8, 10])
        ]
    )
    sixtermssolution = sympy.linsolve(
        (A, b),
        [sympy.var(coefficient)
         for coefficient in ["α", "β", "γ", "δ", "ϵ", "ζ"]
        ]
    )
    
    print(sixtermssolution)


run_spin_five_half_coefficients()

{(-34313125/4096, 768325/512, 885385/256, 10343/16, -707/16, -35/2)}

The solution exists, which means that we can satisfy the second ansatz as follows

\begin{darray}{rcl} \frac{1}{\hbar^{12}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{6} = -\frac{34313125}{4096} + \frac{768325}{512 \hbar^{2}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right) + \frac{885385}{256 \hbar^{4}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{2} + \frac{10343}{16 \hbar^{6}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{3} - \frac{707}{16 \hbar^{8}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{4} - \frac{35}{2 \hbar^{10}} \left( \textbf{S}_{1} \cdot \textbf{S}_{2} \right)^{5}. \end{darray}

Thus, the series for \(H^{DE}\) with \(S_{1,2} = \dfrac{5}{2}\) terminates after the quintic term.

For task 2, we continue to follow the lead of Nolting and Ramakanth, writing down the following new six-term ansatz

\begin{darray}{rcl} H = -t \left[ J_0(5/2) + J_1(5/2) \left( \textbf{S}_1 \cdot \textbf{S}_2 \right) + J_2(5/2) \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^2 + J_3(5/2) \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^3 + J_4(5/2) \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^4 + J_5(5/2) \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^5 \right]. \end{darray}

In task 2, we specified that the eigenvalues for the double exchange interaction need to be

\begin{darray}{rcl} E = -t \frac{\left\lvert S_z \right\rvert}{2 S_{1,2} + 1}, \end{darray}

which means our Hamiltonian should yield the following eigenvalues when \(S_z\) equals the following values:

\begin{darray}{rcl} E = - \frac{t}{6} \begin{cases} 5 & \text{for } S_z = -5 \\ 4 & \text{for } S_z = -4 \\ 3 & \text{for } S_z = -3 \\ 2 & \text{for } S_z = -2 \\ 1 & \text{for } S_z = -1 \\ 0 & \text{for } S_z = 0 \\ 1 & \text{for } S_z = 1 \\ 2 & \text{for } S_z = 2 \\ 3 & \text{for } S_z = 3 \\ 4 & \text{for } S_z = 4 \\ 5 & \text{for } S_z = 5 \end{cases} \end{darray}

We follow a similar procedure as before. This time, we will solve a series of equations and obtain the constants \(J_0(5/2)\) through \(J_5(5/2)\). The functions below construct and solve this series of equations:

In [19]:

def ansatz2_build_spin_product_powers_expression(number_of_terms, s_t, s1s2_eigenvalue):
    hbar = sympy.var("ℏ")
    s1s2 = sympy.var("S1⋅S2")
    coefficients = [sympy.var(coefficient) for coefficient in
                    ["J0",
                     "J1",
                     "J2",
                     "J3",
                     "J4",
                     "J5"]]
    lhs = [sympy.Rational(s_t, 6)]
    rhs = []
    
    if len(coefficients) >= number_of_terms:
        for n, coefficient in zip(range(0, number_of_terms, 1),
                                  coefficients[:number_of_terms]):
            rhs.append(
                coefficient * (hbar ** 2 * s1s2_eigenvalue) ** (n)
            )
    else:
        raise
        
    lhs = functools.reduce(sympy.Add, lhs)
    rhs = functools.reduce(sympy.Add, rhs)
    equality = sympy.Eq(lhs, rhs)
    
    return equality


def ansatz2_series_of_equations():
    s = sympy.Rational(5, 2)
    list_of_s_t = range(0, 6, 1)
    s_t_and_spin_product_list = [
        [s_t,
         product_of_two_identical_spins(s, sympy.Rational(s_t))]
        for s_t in list_of_s_t]
    
    equation_series = []
    for s_t_pair in s_t_and_spin_product_list:
        equation_series.append(ansatz2_build_spin_product_powers_expression(6, s_t_pair[0], s_t_pair[1]))
    
    return equation_series

As a sanity check, let's again test out the above functions against the \(S_{1,2} = 1\) case, where we know that the constants should be

coeff	value
J₀(1)	5/9
J₁(1)	ℏ^-2/6
J₂(1)	- ℏ^-4/18

In [20]:

def run_spin_one_coefficients_for_ansatz2():

    def ansatz2_build_spin_product_powers_expression(number_of_terms, s_t, s1s2_eigenvalue):
        hbar = sympy.var("ℏ")
        s1s2 = sympy.var("S1⋅S2")
        coefficients = [sympy.var(coefficient) for coefficient in
                        ["J0",
                         "J1",
                         "J2",
                         "J3",
                         "J4",
                         "J5"]]
        lhs = [sympy.Rational(s_t, 3)]
        rhs = []

        if len(coefficients) >= number_of_terms:
            for n, coefficient in zip(range(0, number_of_terms, 1),
                                      coefficients[:number_of_terms]):
                rhs.append(
                    coefficient * (hbar ** 2 * s1s2_eigenvalue) ** (n)
                )
        else:
            raise

        lhs = functools.reduce(sympy.Add, lhs)
        rhs = functools.reduce(sympy.Add, rhs)
        equality = sympy.Eq(lhs, rhs)

        return equality


    def ansatz2_series_of_equations():
        s = sympy.Rational(1)
        list_of_s_t = range(0, 3, 1)
        s_t_and_spin_product_list = [
            [s_t,
             product_of_two_identical_spins(s, sympy.Rational(s_t))]
            for s_t in list_of_s_t]

        equation_series = []
        for s_t_pair in s_t_and_spin_product_list:
            equation_series.append(
                ansatz2_build_spin_product_powers_expression(3, s_t_pair[0], s_t_pair[1]))

        return equation_series


    equation_system = ansatz2_series_of_equations()
    A, b = sympy.linear_eq_to_matrix(
        [equations.rhs - equations.lhs for equations in equation_system],
        [sympy.var(coefficient)
         for coefficient in
         ["J0", "J1", "J2"]]
    )
    threetermssolution = sympy.linsolve(
        (A, b),
        [sympy.var(coefficient)
         for coefficient in
         ["J0", "J1", "J2"]
        ]
    )
    
    print(threetermssolution)


run_spin_one_coefficients_for_ansatz2()

{(5/9, 1/(6*ℏ**2), -1/(18*ℏ**4))}

The coefficients match as expected, so now we apply the same procedure to our ansatz of interest:

In [21]:

def run_spin_five_half_coefficients_for_ansatz2():
    equation_system = ansatz2_series_of_equations()
    A, b = sympy.linear_eq_to_matrix(
        [equations.rhs - equations.lhs for equations in equation_system],
        [sympy.var(coefficient)
         for coefficient in
         ["J0", "J1", "J2", "J3", "J4", "J5"]]
    )
    sixtermssolution = sympy.linsolve(
        (A, b),
        [sympy.var(coefficient)
         for coefficient in
         ["J0", "J1", "J2", "J3", "J4", "J5"]
        ]
    )
    
    print(sixtermssolution)


run_spin_five_half_coefficients_for_ansatz2()

{(1267325/1990656, 601247/(17418240*ℏ**2), -86357/(10886400*ℏ**4), -1633/(2721600*ℏ**6), 43/(272160*ℏ**8), 1/(48600*ℏ**10))}

The coefficients are therefore

coeff	value
J₀(5/2)	1267325/1990656
J₁(5/2)	601247ℏ^-2/17418240
J₂(5/2)	-86357ℏ^-4/10886400
J₃(5/2)	-1633ℏ^-6/2721600
J₄(5/2)	43ℏ^-8/272160
J₅(5/2)	ℏ^-10/48600

and the double exchange Hamiltonian for a pair of spin \(\textbf{S} = \frac{5}{2}\) particles is

\begin{darray}{rcl} H = -t \left[ \frac{1267325}{1990656} + \frac{601247}{17418240 \hbar^2} \left( \textbf{S}_1 \cdot \textbf{S}_2 \right) - \frac{86357}{10886400 \hbar^4} \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^2 - \frac{1633}{2721600 \hbar^6} \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^3 + \frac{43}{272160 \hbar^8} \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^4 + \frac{1}{48600 \hbar^{10}} \left( \textbf{S}_1 \cdot \textbf{S}_2 \right)^5 \right] \end{darray}

We implement the double exchange Hamiltonian into the function below

In [22]:

def double_exchange_model_energy(s_z):
    t = sympy.var("t")
    s1s2 = sympy.var("S1⋅S2")
    model = (
        -t * (
            Rational(1267325, 1990656) +
            Rational(601247, 17418240) * s1s2 -
            Rational(86357, 10886400) * s1s2 ** 2 -
            Rational(1633, 2721600) * s1s2 ** 3 +
            Rational(43, 272160) * s1s2 ** 4 +
            Rational(1, 48600) * s1s2 ** 5
        )
    )
    
    s1s2_energies = {
        -5: Rational(5, 4),
        -4: -Rational(11, 4),
        -3: -Rational(23, 4),
        -2: -Rational(31, 4),
        -1: -Rational(35, 4),
        0: -Rational(35, 4),
        1: -Rational(31, 4),
        2: -Rational(23, 4),
        3: -Rational(11, 4),
        4: Rational(5, 4),
        5: Rational(25, 4)
    }
    
    return model.subs(s1s2, s1s2_energies[s_z])

How double exchange affects the superexchange transition spectra

We now combine the superexchange term and the double exchange term and analyze how this affects the \(\Delta{}S_z = \pm 1\) transition spectra.

First, we will proceed by assuming that we can treat the superexchange and double exchange terms as independent. We've already solved both systems separately in the previous sections, so all we need to do is add together the energies. We construct a table showing how the energy spectra changes with \(S_z\):

S_z	E^SE	E^DE	E^SE + E^DE
-5	5J^SE/2	-5t/6	\( \dfrac{5J^{SE}}{2} - \dfrac{5t}{6} \)
-4	-11J^SE/2	-2t/3	\( -\dfrac{11J^{SE}}{2} - \dfrac{2t}{3} \)
-3	-23J^SE/2	-t/2	\( -\dfrac{23J^{SE}}{2} - \dfrac{t}{2} \)
-2	-31J^SE/2	-t/3	\( -\dfrac{31J^{SE}}{2} - \dfrac{t}{3} \)
-1	-35J^SE/2	-t/6	\( -\dfrac{35J^{SE}}{2} - \dfrac{t}{6} \)
0	-35J^SE/2	0	\( -\dfrac{35J^{SE}}{2} \)
1	-31J^SE/2	-t/6	\( -\dfrac{31J^{SE}}{2} - \dfrac{t}{6} \)
2	-23J^SE/2	-t/3	\( -\dfrac{23J^{SE}}{2} - \dfrac{t}{3} \)
3	-11J^SE/2	-t/2	\( -\dfrac{11J^{SE}}{2} - \dfrac{t}{2} \)
4	5J^SE/2	-2t/3	\( \dfrac{5J^{SE}}{2} - \dfrac{2t}{3} \)
5	25J^SE/2	-5t/6	\( \dfrac{25J^{SE}}{2} - \dfrac{5t}{6} \)

In the superexchange section, we showed that the \(2 S_1 S_2 = -35/2 \), \(S^2 = 0\), and \(S_z = 0\) state can be associated with the antiferromagnetic configuration and the state \(2 S_1 S_2 = 25/2 \), \(S^2 = 30\), and \(S_z = 5\) with the ferromagnetic one. This means that the energy cost for a full spin is:

\begin{darray}{rcl} \Delta E(0 \to 5) = 30J^{SE} - \frac{5t}{6} \end{darray}

We apply this model to (Ba_1-xK_x)(Zn, Mn)₂As₂ and consider the case for a nearest-neighbor pair of Mn atoms for doping levels of x = 0, x = 0.2, and x = 0.4. The VASP calculated spin flip energies are

x	E(AFM) - E(FM) (meV)
0.0	407
0.2	207
0.4	148

For x = 0, there is no double exchange and thus t = 0. We can calculate the superexchange constant \(J^{SE}\) for \(0 \to 5\):

In [23]:

j = sympy.var("J")
t = sympy.var("t")
model1 = j * sympy.Rational(30) - t * sympy.Rational(5, 6)
round(sympy.solveset(sympy.Eq(407, model1.subs(t, 0)), j, domain=sympy.S.Reals).n().args[0], 1)

13.6

This yields a superexchange parameter of \(J^{SE} = 13.6\) meV. Let's assume that \(J^{SE}\) is relatively unchanged by doping. This allows us to solve for \(t\) for the x = 0.2 and x = 0.4 cases.

In [24]:

print(f"x = 0.2: t = "
      f"{round(sympy.solveset(sympy.Eq(207, model1.subs(j, 13.6)), t, domain=sympy.S.Reals).n().args[0], 0):.0f} "
      f"meV\n"
      f"x = 0.4: t = "
      f"{round(sympy.solveset(sympy.Eq(148, model1.subs(j, 13.6)), t, domain=sympy.S.Reals).n().args[0], 0):.0f} "
      f"meV")

x = 0.2: t = 241 meV
x = 0.4: t = 312 meV

Now that we have \(J^{SE}\) and \(t\), we can predict the \(\Delta S_z = \pm 1\) transition energies. The function below implements this calculation:

In [25]:

def calculate_transition_energies(xdoping, tvalue, jvalue):

    s_z = sympy.var("Sz")
    j = sympy.var("J")
    t = sympy.var("t")
    model_energy_vs_s_z = (
        j *
        (s_z * (s_z + 1) -
         2 * Rational(5, 2) * (Rational(5, 2) + 1)) -
        t * np.abs(s_z) * Rational(1, 6)
    )
    
    initial_s_z = list(range(-5, 5))
    final_s_z = list(range(-4, 6))
    transition_list = []

    for initial, final in zip(initial_s_z, final_s_z):
        initial_energy = model_energy_vs_s_z.subs({s_z: initial, t: tvalue, j: jvalue})
        final_energy = model_energy_vs_s_z.subs({s_z: final, t: tvalue, j: jvalue})
        delta_energy = final_energy - initial_energy
        superexchange_only_energy = (
            model_energy_vs_s_z.subs({s_z: final, t: 0, j: jvalue}) -
            model_energy_vs_s_z.subs({s_z: initial, t: 0, j: jvalue})
        )
        transition_list.append([f"{initial}", f"{final}", f"{round(delta_energy, 1):.1f}",
                                f"{superexchange_only_energy}"])
    
    table_headers = ("*S*<sub>z</sub>",
                     "*S*<sup>'</sup><sub>z</sub>",
                     "ΔE<sup>SE+DE</sup> (meV)", "ΔE<sup>SE</sup> (meV)")

    transition_list_table = (
            tabulate.tabulate(tabular_data=transition_list,
                              headers=table_headers,
                              tablefmt='pipe',
                              numalign='center',
                              stralign='center'))
    
    IPython.display.display(IPython.display.Markdown(
        f"For x = {xdoping}, J = {jvalue}, t = {tvalue}\n\n"
        f"{transition_list_table}"
    ))

We run the calculation using parameters from above, and get

In [26]:

calculate_transition_energies(0.0, 0, 13.6)
calculate_transition_energies(0.2, 241, 13.6)
calculate_transition_energies(0.4, 312, 13.6)

For x = 0.0, J = 13.6, t = 0

S_z	S^'_z	ΔE^SE+DE (meV)	ΔE^SE (meV)
-5	-4	-108.8	-108.8
-4	-3	-81.6	-81.6
-3	-2	-54.4	-54.4
-2	-1	-27.2	-27.2
-1	0	0	0
0	1	27.2	27.2
1	2	54.4	54.4
2	3	81.6	81.6
3	4	108.8	108.8
4	5	136	136

For x = 0.2, J = 13.6, t = 241

S_z	S^'_z	ΔE^SE+DE (meV)	ΔE^SE (meV)
-5	-4	-68.6	-108.8
-4	-3	-41.4	-81.6
-3	-2	-14.2	-54.4
-2	-1	13	-27.2
-1	0	40.2	0
0	1	-13	27.2
1	2	14.2	54.4
2	3	41.4	81.6
3	4	68.6	108.8
4	5	95.8	136

For x = 0.4, J = 13.6, t = 312

S_z	S^'_z	ΔE^SE+DE (meV)	ΔE^SE (meV)
-5	-4	-56.8	-108.8
-4	-3	-29.6	-81.6
-3	-2	-2.4	-54.4
-2	-1	24.8	-27.2
-1	0	52	0
0	1	-24.8	27.2
1	2	2.4	54.4
2	3	29.6	81.6
3	4	56.8	108.8
4	5	84	136

Quantitatively, these don't align all that well with the neutron scattering peaks, especially for x = 0.0. In that case, we're off by a factor of 2. However, an interesting outcome of this is that including double exchange lifts the degeneracy between \(S_z = -1\) and \(S_z = 0\). And, since INS picks up both \(\Delta{}S_z = 1\) and \(\Delta{}S_z = -1\), both negative and positive energies will be detected. We already see that, in the x = 0.2 case, two different peaks can be centered nearby each other, leading to overlap. Overlapping peaks can be a potential source of the strong peak broadening we see in the INS data.

We should note that the previous analysis did not take into account magnetic moment softening. The DFT local Mn moments in the undoped case are approximately 3.7 \(\mu_B\) each, which is significantly reduced from the ideal 5 \(\mu_B\) case. We can make a back of the envelope-style estimate for how much this will reduce the eigenvalues of the two-spin operators:

In [27]:

# S^2 eigenvalue
test_s_squared_no_reduced_moment = 5 * (5 + 1)
test_s_squared_reduced_moment = 3.7 * (3.7 + 1)
# S1^2, S2^2 eigenvalues
test_s1_s2_squared_no_reduced_moment = (5/2) * (5/2 + 1)
test_s1_s2_squared_reduced_moment = (3.7/2) * (3.7/2 + 1)
# S1 dot S2
test_2_s1_dot_s2_no_reduced_moment = (
    test_s_squared_no_reduced_moment -
    2 * test_s1_s2_squared_no_reduced_moment
)
test_2_s1_dot_s2_reduced_moment = (
    test_s_squared_reduced_moment -
    2 * test_s1_s2_squared_reduced_moment
)
# Double exchange
test_double_exchange_no_reduced_moment = [
    sz / (2 * (5 / 2) + 1) for sz in range(1, 6)
]
test_double_exchange_reduced_moment = [
    (sz * (3.7 / 5)) / (2 * (3.7) + 1) for sz in range(1, 6)
]

# Reduction in S^2
print(test_s_squared_reduced_moment /
      test_s_squared_no_reduced_moment)

# Reduction in 2(S1 dot S2)
print(test_2_s1_dot_s2_reduced_moment /
      test_2_s1_dot_s2_no_reduced_moment)

# Reduction in double exchange
test_ratio_double_exchange = [
    reduced / unreduced
    for reduced, unreduced in zip(test_double_exchange_reduced_moment,
                                  test_double_exchange_no_reduced_moment)
]
print(test_ratio_double_exchange)

0.5796666666666667
0.5475999999999999
[0.5285714285714286, 0.5285714285714286, 0.5285714285714285, 0.5285714285714286, 0.5285714285714286]

What this means is, if we assume that the basic physics of the two-spin model with spin-5/2 particles is correct, but that hybridization reduces the effective exchange coupling strength, then a direct fit of this model to the DFT energies will yield parameters \(J^{SE}\) and \(t\) that are too large by a factor of nearly two. A quick fix then is to reduce the superexchange parameter \(J\) and the hopping parameter \(t\) using the above ratios:

In [28]:

calculate_transition_energies(0.0, 0, 13.6*0.548)
calculate_transition_energies(0.2, 241*0.529, 13.6*0.548)
calculate_transition_energies(0.4, 312*0.529, 13.6*0.548)

For x = 0.0, J = 7.452800000000001, t = 0

S_z	S^'_z	ΔE^SE+DE (meV)	ΔE^SE (meV)
-5	-4	-59.6	-59.6224
-4	-3	-44.7	-44.7168
-3	-2	-29.8	-29.8112
-2	-1	-14.9	-14.9056
-1	0	0	0
0	1	14.9	14.9056
1	2	29.8	29.8112
2	3	44.7	44.7168
3	4	59.6	59.6224
4	5	74.5	74.528

For x = 0.2, J = 7.452800000000001, t = 127.489

S_z	S^'_z	ΔE^SE+DE (meV)	ΔE^SE (meV)
-5	-4	-38.4	-59.6224
-4	-3	-23.5	-44.7168
-3	-2	-8.6	-29.8112
-2	-1	6.3	-14.9056
-1	0	21.2	0
0	1	-6.3	14.9056
1	2	8.6	29.8112
2	3	23.5	44.7168
3	4	38.4	59.6224
4	5	53.3	74.528

For x = 0.4, J = 7.452800000000001, t = 165.048

S_z	S^'_z	ΔE^SE+DE (meV)	ΔE^SE (meV)
-5	-4	-32.1	-59.6224
-4	-3	-17.2	-44.7168
-3	-2	-2.3	-29.8112
-2	-1	12.6	-14.9056
-1	0	27.5	0
0	1	-12.6	14.9056
1	2	2.3	29.8112
2	3	17.2	44.7168
3	4	32.1	59.6224
4	5	47	74.528

The spectra seem more plausible now for x = 0.0. And the trend for x = 0.2 of overlapping peaks becomes even more plausible with these renormalized energies. It's worth modeling the INS peaks as Gaussians and using the above energies as their respective means. The widths and amplitudes are left as empirical parameters to be tuned. The goal is to show that this explanation of overlapping peaks is a plausible explanation for the INS data. We implement this visualization in the next section.

Interactive model of INS spectra

We implement an interactive model visualization for the INS spectra, using our above results for input. This visualization will plot two main curves, the predicted peaks if only the superexchange term is used, and the predicted peaks if both superexchange and double exchange are used. Some things to note about this visual model:

We take the absolute value of all transition energies, as we assume that both \(\Delta{}S_z = 1\) and \(\Delta{}S_z = -1\) transitions will be detected.
The width of the Gaussian peaks is an empirical parameter that we can tune. Based on the INS data, peak widths seem to be in the general neighborhood of 4 to 6 meV.
The amplitudes of the Gaussian peaks are an empirical parameter that we can tune. Based on the INS data, peak amplitudes seem to be in the general neighborhood of 3 to 4 units.
The default widths assume that the doped spectra have a slightly larger Gaussian width than the undoped spectra.
A major assumption is that the amplitude for transitions between states with larger \(\left\lvert S_z \right\rvert\) eigenvalues is reduced compared to the smaller ones, for example the transition \(0 \to 1\) has a larger amplitude than \(1 \to 2\).
- For this model, we reduce the amplitude of each higher transition — such as going from \(0 \to 1\) to \(1 \to 2\) — by a factor of 2.

The code below implements the interactive dashboard for the INS spectra plot.

In [29]:

def create_gaussian_curve(amplitude, mean, sd, emin, emax):
    x = np.linspace(emin, emax, 5000)
    gaussian = amplitude * np.exp(- ((x - mean)**2)/(2*sd**2))
    return np.column_stack((x, gaussian))


def widget_calculate_transition_energies(tvalue, jvalue):

    s_z = sympy.var("s_z")
    j = sympy.var("J")
    t = sympy.var("t")
    model_energy_vs_s_z = (
        j *
        (s_z * (s_z + 1) -
         2 * Rational(5, 2) * (Rational(5, 2) + 1)) -
        t * np.abs(s_z) * Rational(1, 6)
    )
    
    initial_s_z = list(range(-5, 5))
    final_s_z = list(range(-4, 6))
    transition_list = []

    for initial, final in zip(initial_s_z, final_s_z):
        initial_energy = model_energy_vs_s_z.subs({s_z: initial, t: tvalue, j: jvalue})
        final_energy = model_energy_vs_s_z.subs({s_z: final, t: tvalue, j: jvalue})
        delta_energy = np.abs(np.float(final_energy - initial_energy))
        superexchange_only_energy = (
            model_energy_vs_s_z.subs({s_z: final, t: 0, j: jvalue}) -
            model_energy_vs_s_z.subs({s_z: initial, t: 0, j: jvalue})
        )
        transition_list.append(delta_energy)
    
    return np.array(transition_list)


def update_plot(hopping_parameter,
                superexchange_parameter,
                emax,
                gaussian_sd,
                gaussian_sd_no_hop,
                gaussian_amplitude,
                gaussian_amplitude_no_hop,
                use_envelope,
                env_amp,
                tau,
                vertical_axis_max):
    
    gaussian_mean = widget_calculate_transition_energies(
        hopping_parameter, superexchange_parameter
    )
    gaussian_mean_no_hop = widget_calculate_transition_energies(
        0, superexchange_parameter
    )
    
    gaussian_curves = []
    amplitude_scaling = [
        0.0625,
        0.125,
        0.25,
        0.5,
        1,
        1,
        0.5,
        0.25,
        0.125,
        0.0625
    ]
    summed_curve = np.zeros(5000)
    summed_curve_no_hop = np.zeros(5000)

    for mean_value, mean_value_no_hop, amp_scale in zip(
        gaussian_mean,
        gaussian_mean_no_hop,
        amplitude_scaling
    ):
    
        curve_for_mean = create_gaussian_curve(
            amp_scale * gaussian_amplitude, mean_value, gaussian_sd, 5, emax
        )
        curve_for_mean_no_hop = create_gaussian_curve(
            amp_scale * gaussian_amplitude_no_hop, mean_value_no_hop,
            gaussian_sd_no_hop, 5, emax
        )
        gaussian_curves.append(curve_for_mean)
        
        if use_envelope:
            
            summed_curve += (
                curve_for_mean[:, 1] *
                np.exp(-curve_for_mean[:, 0] / tau) *
                env_amp
            )
            summed_curve_no_hop += (
                curve_for_mean_no_hop[:, 1] *
                np.exp(-curve_for_mean_no_hop[:, 0] / tau) *
                env_amp
            )
        else:
            
            summed_curve += curve_for_mean[:, 1]
            summed_curve_no_hop += curve_for_mean_no_hop[:, 1]

    fig = plt.figure(figsize=(10,6))
    ax = fig.add_subplot(111)
    for curve in gaussian_curves:
        ax.plot(curve[:, 0], curve[:, 1],
                linestyle='dashed', linewidth=1)
        
    ax.plot(gaussian_curves[0][:, 0], summed_curve[:], linewidth=2,
            label="Full curve")
    ax.plot(gaussian_curves[0][:, 0], summed_curve_no_hop[:],
            linewidth=2, label="Full curve (no hopping)")
    
    ax.set_ylim([0, vertical_axis_max])
    
    ax.set_xlabel("Energy (meV)")
    ax.set_ylabel("INS Intensity")
    
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

    plt.show()


def visualize_ins_model():
    
    hopping_parameter = widgets.IntSlider(
        min=0, max=500, step=1, value=128,
        description="Hopping t:"
    )
    superexchange_parameter = widgets.FloatSlider(
        min=0.0, max=40.0, step=0.1,
        value=7.5, readout_format='.1f',
        description="Superexchange J:"
    )
    emax = widgets.IntText(
        30, description="Energy window max"
    )
    gaussian_sd = widgets.FloatText(
        6.5, description="Width"
    )
    gaussian_sd_no_hop = widgets.FloatText(
        4.75, description="Width (no hopping)"
    )
    gaussian_amplitude = widgets.FloatText(
        3, description="Amplitude"
    )
    gaussian_amplitude_no_hop = widgets.FloatText(
        4.125, description="Amplitude (no hopping)"
    )
    use_envelope = widgets.Checkbox(
        value=False, description="Add exp. envelope"
    )
    env_amp = widgets.FloatText(
        2, description="Env. Amplitude"
    )
    tau = widgets.FloatText(
        20, description="Env. half-life"
    )
    vertical_axis_max = widgets.FloatText(
        8, description="Vertical axis max"
    )

    interactive_window = widgets.interactive(
        update_plot,
        hopping_parameter=hopping_parameter,
        superexchange_parameter=superexchange_parameter,
        emax=emax,
        gaussian_sd=gaussian_sd,
        gaussian_sd_no_hop=gaussian_sd_no_hop,
        gaussian_amplitude=gaussian_amplitude,
        gaussian_amplitude_no_hop=gaussian_amplitude_no_hop,
        use_envelope=use_envelope,
        env_amp=env_amp,
        tau=tau,
        vertical_axis_max=vertical_axis_max
    )

    widget_column1 = widgets.VBox(
        [interactive_window.children[1],
         interactive_window.children[0]]
    )
    widget_column2 = widgets.VBox(
        [interactive_window.children[3],
         interactive_window.children[4],
         interactive_window.children[5],
         interactive_window.children[6],
         interactive_window.children[2]]
    )
    widget_column3 = widgets.VBox(
        [interactive_window.children[7],
         interactive_window.children[8],
         interactive_window.children[9],
         interactive_window.children[10]
        ]
    )
    widgets_box = widgets.HBox(
        [widget_column1, widget_column2, widget_column3]
    )
    widgets_full = widgets.VBox(
        [widgets_box, interactive_window.children[11]]
    )

    IPython.display.display(widgets_full)

    
visualize_ins_model()

The following figure illustrates the "best case" scenario for the model, yielding the qualitative features of the INS spectra and approximate quantitative agreement.