Modifying the Equations of Motion
Once we have made a decision regarding which operators to disregard during the evolution of our system, we still need to decide how to modify the equations of motion for the remaining operators. This is necessary because even if we have decided to neglect the time evolution of the "unimportant" operators, it is still possible for them to appear in the equations of motion. For example, one possible term in our original equations of motion might be
\[ \frac{\mathrm{d}}{\mathrm{d}t} \langle \sigma_1^x \sigma_2^y \rangle = \langle \sigma_1^x \sigma_2^y \sigma_3^z \rangle + \ldots \]
If we have made our approximation such that the operator \( \sigma_{1}^{x}\sigma_{2}^{y} \) is kept, but the operator \( \sigma_{1}^{x}\sigma_{2}^{y}\sigma_{3}^{z} \) is discarded, then we need to decide how to alter the right side of this equation in the absence of the discarded operator.
In general, there are many ways such a modification could be made, with some methods being more applicable than others depending on the particular problem at hand. Raqet currently supports two methods for modifying the equations of motion:
Truncate
: The simplest possible method of accounting for the omitted operators, in which we simply drop them from the equations of motion, effectively setting their expectation values to zero for the entire simulation.Ursell
: A method which uses the concept of an "Ursell function" in order to express the omitted expectation value in terms of the expectation values which are kept. This is explained in more detail in the next section.
The Ursell Method
The idea behind the Ursell method is to replace the expectation value of an operator with a sum over all possible partitions of that expectation value. An example of such an approximation might look like
\[ \langle \sigma_{1}^{x}\sigma_{2}^{y}\sigma_{3}^{z} \rangle \to p \langle \sigma_{1}^{x}\sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle + q \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \sigma_{3}^{z} \rangle + r \langle \sigma_{2}^{y} \rangle \langle \sigma_{1}^{x} \sigma_{3}^{z} \rangle + s \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle \]
for some coefficients p, q, r, and s. This "factoring" reflects the idea that the operator \( \sigma_{1}^{x}\sigma_{2}^{y}\sigma_{3}^{z} \) is "too big," and thus it can be approximated by "breaking it up" into smaller pieces.
To completely specify the method, we must determine these coefficients. The specific choice Raqet makes when determining these coefficients is based upon the notion of the joint cumulant, or Ursell function, of a collection of single-particle operators. If we have a collection of operators \( \{ \sigma_{i} \} \) on different sites, we define the joint cumulant of these operators to be
\[ \kappa \left ( \sigma_{1}, \dots, \sigma_{n} \right ) = \sum_{\pi} (-1)^{|\pi|-1} \left ( |\pi| - 1 \right )! \prod_{B \in \pi} \langle \prod_{i\in B}\sigma_{i} \rangle \]
In the above expression, the sum over \( \pi \) is a sum over all possible partitions of the collection of sites, \( |\pi| \) is the number of parts in the partition, and \( B \) represents a part inside of a partition. As a concrete example, we might have
\[ \kappa \left ( \sigma_{1}^{x}, \sigma_{2}^{y}, \sigma_{3}^{z} \right ) = \langle \sigma_{1}^{x} \sigma_{2}^{y} \sigma_{3}^{z} \rangle - \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \sigma_{3}^{z} \rangle - \langle \sigma_{1}^{x} \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle - \langle \sigma_{2}^{y} \rangle \langle \sigma_{1}^{x} \sigma_{3}^{z} \rangle + 2 \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle \]
We can use this definition to define our factoring formula. If we have an operator we wish to "factor" which involves single-particle operators on several sites, we find the appropriate factoring formula by setting the joint cumulant of these to zero,
\[ \kappa \left ( \sigma_{1}, \dots, \sigma_{n} \right ) = 0. \]
For our example above, this would correspond to
\[ \langle \sigma_{1}^{x} \sigma_{2}^{y} \sigma_{3}^{z} \rangle - \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \sigma_{3}^{z} \rangle - \langle \sigma_{1}^{x} \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle - \langle \sigma_{2}^{y} \rangle \langle \sigma_{1}^{x} \sigma_{3}^{z} \rangle + 2 \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle = 0 \]
This then implies
\[ \langle \sigma_{1}^{x} \sigma_{2}^{y} \sigma_{3}^{z} \rangle = \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \sigma_{3}^{z} \rangle + \langle \sigma_{1}^{x} \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle + \langle \sigma_{2}^{y} \rangle \langle \sigma_{1}^{x} \sigma_{3}^{z} \rangle - 2 \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle. \]
In some case, it may be necessary to apply this condition recursively. For example, if we had chosen our approximation such that the operator \( \sigma_{1}^{x} \sigma_{3}^{z} \) is also deemed "unimportant," then it should not appear in the equations of motion either. We can eliminate this operator by applying a similar factoring to it as well, so that
\[ \langle \sigma_{1}^{x} \sigma_{3}^{z} \rangle \to \langle \sigma_{1}^{x} \rangle \langle \sigma_{3}^{z} \rangle \]
Plugging this back into the first approximation and simplifying, we ultimately have
\[ \langle \sigma_{1}^{x} \sigma_{2}^{y} \sigma_{3}^{z} \rangle = \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \sigma_{3}^{z} \rangle + \langle \sigma_{1}^{x} \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle - \langle \sigma_{1}^{x} \rangle \langle \sigma_{2}^{y} \rangle \langle \sigma_{3}^{z} \rangle. \]
We emphasize here that setting a joint cumulant to zero is a way to find an approximate expression for an expectation value which we have discarded. It is not, however, a method for determining which operators to discard. In Raqet, the decision to discard an operator is always based on its diameter, and setting cumulants to zero is used to determine factoring expressions only after the list of important operators has been constructed. In this way, the factoring method differs from a standard cumulant expansion.