--- # vim: spell tw=79 # Document meta data: title: New Projection Code subtitle: Architecture, Design, Manual author: - Martin Ueding () date: \today abstract: > The [existing projection code](https://github.com/HISKP-LQCD/sLapH-projection) only works for the $\rho$ resonance. We want to extend it to multiple particles and arbitrary isospin. This note contains the architecture and a survey of possible programming languages and libraries for the implementation. Also the design and implementation is covered here, including documentation for all functions in the Mathematica package. # LaTeX presentation options: urlcolor: blue toc: true numbersections: true documentclass: scrreprt #mainfont: Libertinus Serif sansfont: Noto Sans #monofont: Noto Mono #mathfont: Libertinus Math ... # Specification The [sLapH contraction code](https://github.com/HISKP-LQCD/sLapH-contractions) can only work with single particle operators multiplied together. These operators are just Dirac-$\Gamma$ structure and integer three-momentum $p$. The spin and isospin projection code needs to take these and give a correlator matrix for given total momentum $P^2$ and irrep $\Gamma$. We have to decide whether and when we want to include baryons. As discussed with Carsten we will stick with the single cover of the octahedral group and therefore just mesons for the time being. The result that we want from the projection is a data frame with these columns: - Total momentum $P^2$ - Irrep $\Gamma$ - Irrep row $\alpha$ - Correlator matrix row index (source operator ID) - Correlator matrix column index (sink operator ID) - HDF5 dataset name - Boolean indicating whether the diagram needs to be complex conjugated - Complex weighting factor With this we can take the output of the sLapH contraction code and project into the irreps. In order to get this data frame, we will do the first four of these steps: 1. Group theory 2. Isospin projection 3. Spin projection 4. Analytic Wick contractions 5. Actual projection of the correlators from the HDF5 files. # Architecture Isospin and spin projection should be independent of each other. The current projection code does both concepts at the same time and exploits certain simplifications that come with the $\rho$-channel. In this new code we want to refrain from any such premature simplifications to use it for any physical process to come. ## Group theory We will need to have the following of the octahedral group and each of the relevant little groups: - Irreducible representation matrices $D^\Gamma$ - Irreducible continuum matrices $D^J$, the Wigner-D matrices - Cartesian 3D representation matrices $R_g$ There must be correspondence between the various representations, otherwise we could not sensibly multiply different representations together. It still needs to be decided whether we want to have all the little groups corresponding to every momentum or whether we do rotations to a reference momentum ourselves. ## Isospin The isospin projection can be done by hand for the cases we have done so far. For the $I = 0$, $I_3 = 0$ case for instance we construct a two-pion operator as $$ \pi\pi(\vec p_1, \vec p_2) = \pi^+(\vec p_1) \; \pi^-(\vec p_2) + \pi^-(\vec p_1) \; \pi^+(\vec p_2) - \frac{1}{\sqrt{3}} \pi^0(\vec p_1) \; \pi^0(\vec p_2) \,. $$ For $I = 1$, $I_3 = 0$ the operator is $$ \pi\pi(\vec p_1, \vec p_2) = \frac{1}{\sqrt{2}} \left(\pi^+(\vec p_1) \; \pi^-(\vec p_2) - \pi^-(\vec p_1) \; \pi^+(\vec p_2) \right) \,. $$ For $I = 3$, $I_3 = 3$ the operator is just $$ \pi\pi\pi(\vec p_1, \vec p_2, \vec p_3) = \pi^+(\vec p_1) \; \pi^+(\vec p_2) \; \pi^+(\vec p_2) \,. $$ For scalar particles this is very easy. For vector mesons like the $\rho$ we need to build the various spin states $J_3$. As pointed out by Marcus this can also be done automatically. As I plan this to be independent of the remainder, it does not matter so much whether and when we do this automatically. ## Spin For the spin part we consult Markus' “projection” chapter and use Equation (3.34). For the generalization to $N_\text{p}$ particles we have the following expression: \begin{align*} O_\Gamma^{\alpha}(\vec p_\text{cm})^\dagger &= \sum_{M=-J}^J \phi_M \left[ \prod_{i=1}^{N_\text{p}} \sum_{M_i=-{J_i}}^{J_i} \right] \langle J, M | J_1, M_1, \ldots, J_{N_\text P}, M_{N_\text P} \rangle \\&\quad\times \sum_{\beta=1}^{\mathop{\mathrm{nrow}(\Gamma)}} \tilde\phi_\beta \sum_{g \in \mathop{\mathrm{LG}}(\vec p_\text{cm})} D_{\alpha\beta}^\Gamma(g)^* \\&\quad\times \prod_{i=1}^{N_\text{p}} \sum_{M_i'=-J_i}^{J_i} \sum_{M_i''=-J_i}^{J_i} D_{M_i' M_i}^{J_i}(g) \; D_{M_i'' M_i'}^{J_i}(\tilde g) \; O_{i M_i''}^{J_i}(R_{\tilde g}^{-1} R_g R_{\tilde g} \vec p_i)^\dagger \,. \end{align*} The variables $J$ and $J_i$ are determined by the physical process that one investigates. ## Analytic Wick contraction The analytic Wick contraction can be done with the “Quark Contraction Tool” in Mathematica. For the $I = 3$, $I_3 = 3$ channel for zero momentum with source at `x` and sink at `y` we have summands like the following, which we identify as `C6cD`: ```mathematica trace[Gamma^5.DE[{up, up}, {x, y}].Gamma^5.DE[{dn, dn}, {y, x}]]^3 ``` Also we have the `C6cCD`: ```mathematica trace[Gamma^5.DE[{up, up}, {x, y}].Gamma^5.DE[{dn, dn}, {y, x}]] trace[Gamma^5.DE[{dn, dn}, {y, x}].Gamma^5.DE[{up, up}, {x, y}]. Gamma^5.DE[{dn, dn}, {y, x}].Gamma^5.DE[{up, up}, {x, y}]] ``` And also the `C6cC`: ```mathematica trace[Gamma^5.DE[{dn, dn}, {y, x}].Gamma^5.DE[{up, up}, {x, y}]. Gamma^5.DE[{dn, dn}, {y, x}].Gamma^5.DE[{up, up}, {x, y}]. Gamma^5.DE[{dn, dn}, {y, x}].Gamma^5.DE[{up, up}, {x, y}]] ``` With the sLapH method at zero momentum we have it quite easy as the three sources and three sinks are indistinguishable. But with general momenta we cannot simplify it as much. The expression from the contraction code with sources `x1`, `x2`, `x3` and sinks `x4`, `x5`, `x6` contains `C6cD` diagrams like this one: ```mathematica trace[Gamma^5.DE[{up, up}, {x1, x4}].Gamma^5.DE[{dn, dn}, {x4, x1}]] trace[Gamma^5.DE[{up, up}, {x2, x6}].Gamma^5.DE[{dn, dn}, {x6, x2}]] trace[Gamma^5.DE[{up, up}, {x3, x5}].Gamma^5.DE[{dn, dn}, {x5, x3}]] ``` But also we have this one here: ```mathematica trace[Gamma^5.DE[{up, up}, {x1, x4}].Gamma^5.DE[{dn, dn}, {x4, x1}]] trace[Gamma^5.DE[{up, up}, {x2, x5}].Gamma^5.DE[{dn, dn}, {x5, x2}]] trace[Gamma^5.DE[{up, up}, {x3, x6}].Gamma^5.DE[{dn, dn}, {x6, x3}]] ``` The difference is just that the exchange $x_5 \leftrightarrow x_6$ was done. The output from the Quark Contraction Tool is basically a prescription to a linear combination of diagrams from the contraction code. The positions are equivalent to momentum labels. From this output we learn the way that we have to combine the different diagrams to yield a particular isospin. We have to convert this output into a data frame with the following columns: - Correlator matrix row index (source operator ID) - Correlator matrix column index (sink operator ID) - HDF5 data set name like `C2c_uu_p010.d000.g5_p010.d000.g5` - Boolean flag whether the correlator needs to be conjugated to compensate for reversal of quark flow direction - Complex weight factor Ideally we also apply certain symmetry simplifications where we can apply say the permutation $(13)(24)$ on a `C4cD` diagram and get the same result. ## Spin and isospin We have prescriptions for both spin and isospin. They both need to be resolved to yield the wanted linear combination of computed correlation functions. We have already expressed the spin projection in terms of multi-particle operators. The isospin projected operators can therefore be resolved there. From there we can spin project onto given momentum and irrep and so forth. The result will be multiple-particle operators with definite spin and isospin. Taking these we can do the Wick contractions. How do we get the dataset names? Just as we would have to parse the Mathematica output into HDF5 dataset name *patterns* after doing the Wick contractions right after the isospin projection, we now need to extract the full dataset name from the symbolic expression that contains the concrete momenta. That expression will likely have very many summands. Parsing of it needs to be done as part of the program. Either with the PyParsing library in Python or within Mathematica. # Choice of languages and libraries We have a natural interface between the group theory part that provides us with a table of matrices in various irreps and the projection which contains the physical process. The projection shall touch the actual numeric correlators as late as possible. I want to have it compute the spin and isospin projection operator as a data frame and only then apply it to the data. For the choice of programming languages and libraries keep in mind that programming languages are not tools but the very material from which we build things. They stick around once the project is finished. ## Group theory For the group theory there are a few options: Copying tables ~ There are various books and papers which list all the things that we need. Copying them is just error-prone and we would have to check for consistency ourselves. Apparently the `Bethe` Maple library also has it copied from a paper. GAP ~ A commenter on my [Stack Exchange question](https://mathematica.stackexchange.com/q/188767/1507) mentioned [GAP](https://www.gap-system.org/). There is seems to be rather easy to generate the groups. It is built in or can be created from the generators. Then it has the character tables and the conjugacy classes available right away. ``` Oh := SmallGroup(48,48); Oh := Group( (1,3,4,2)(5,7,8,6), (2,5,3)(4,6,7), (1,8)(2,7)(3,6)(4,5) ); Display(CharacterTable(Oh)); Display(ConjugacyClasses(Oh)); ``` Maple ~ This is bad because neither the university nor the institute has have a license for it. We could ask the institute to buy one, though I do not believe that the institute can buy the student version for around 150 EUR. The `Bethe` library has been used by Markus for the $\rho$-project. Perhaps we just ask him to export all we need and hope that we never have to touch it again. Mathematica ~ The university has a Mathematica campus license, so although that is license riddled as well, at least *we* can use it. The octahedral group is not directly available in Mathematica but I managed to get it via the representation as a permutation group. I can also construct the little groups (also known as stabilizers) but not the irreps. There also is the [GTPack](https://www.frontiersin.org/articles/10.3389/fphy.2018.00086/full) library which might add what we need. We have not looked into this in detail yet. Sage ~ I have been looking into alternatives a bit yesterday and found that Sage has the group in some form that one can at least get the multiplication and character table, see [this snippet](https://sagecell.sagemath.org/?z=eJyNjrFOwzAQhvdIeQdLDLWFY0pYGApLJRgQ4gGqNjrcwzVy4mBfI_z2xIGAQB3QLb_-T9_d2bb3gVgEg8oEf-yjaoGCfW_MGF9sZwldagx2GIBwXxZ5HtRquGU3_9fUOmnnybdW31l0e34tygLGDUM17K7YWXwLxOuxs7nb1WXxpCnzx2nrfT7BN7y6PLfiopbfgfGfpsphK9nmswO5lGwp-dTDVuS_-2A74ov1AQJowsAInh0uxEzyUaVn2kyUi18qJIfptDehk5LvXo8GdGLaQYwY_5ozb774qH8A9juBiQ==&lang=sage): ```python import sage.groups.matrix_gps.finitely_generated K. = sage.groups.matrix_gps.finitely_generated.CyclotomicField(8) a = v - v^3 i = v^2 Octa = MatrixGroup([(-1+i)/2, (-1+i)/2, (1+i)/2, (-1-i)/2], [(1+i)/a, 0, 0, (1-i)/a]) print('Character table') print(Octa.character_table()) ``` SciPy ~ This has Clebsch-Gordan coefficients and Wigner-3j symbols, but not really much more. We now read in the tables that are generated by Markus. ## Projection operator construction Also we need Clebsch-Gordan coefficients, but that should be covered in all the packages that we consider. There are very likely to be cancellations between different terms. The rotation matrices contain factors like $1/\sqrt{2}$ and therefore we likely want to do symbolic algebra. We have these options: Maple ~ Not a good idea as we do not want to depend on it more than strictly needed. Mathematica ~ We have a campus license, it is a very powerful language. However it is a proprietary tool and our group has limited experience with it. I personally do not want to invest too much learning into a proprietary tool. Wigner-D matrices are available as `WignerD`. Sage ~ For what we need it will likely be overpowered. As it is a pain to install, we should only do that when strictly needed. SymPy ~ This would be my favorite as I and others know Python and it is free software. Wigner-D matrices are available as `physics.quantum.spin.Rotation` We chose Mathematica for this as we also do the Wick contractions with Mathematica. This way we have just one code. ## Wick contraction Cadabra ~ While asking about Wick contractions in SymPy, the author of the [Cadabra](https://cadabra.science/) CAS system reached out to me offering help implementing Wick contractions in that system. It is programmable in Python. I have described this project to the author and he told me that he was going to try it a bit, but I never heard back from him. Mathematica ~ We have the [Quark Contraction Tool](https://arxiv.org/abs/1603.01576) already available. SymPy ~ There is `physics.secondquant.wicks`, which seems to be able to do Wick contractions. Also there is `physics.secondquant.contraction` which does a contraction between two fermionic operators. I have spend the afternoon of 2019-01-29 looking into it, but the definition of the fermion creation and annihilation operators seem to not support a tensorial structure, at least not with the Wick contractions provided. I have asked a [question on Stack Overflow](https://stackoverflow.com/q/54426941/653152) about this. We have chosen Mathematica for this step due to the great “Quark Contraction Tool”. ## Projecting the data Once we have the data frame with the desired linear combinations of the correlators from the contraction code, we just have to load the HDF5 data and project it. This could be done in various languages. As HDF5 reading seems to be very slow, we might want to consider adding a re-packaging step just as Markus does by converting the HDF5 files from the sLapH contraction code into Pandas HDF5 files. C++ ~ We do not really need much for projecting the data, except a HDF5 library. Perhaps the performance will be better if we do this in C++. But the most likely bottleneck is the HDF5 reading, therefore it does not matter which language we use. Declaring the HDF5 data types in C or C++ is rather bothersome. We can of course just copy them from the contraction code, but having them dynamically is likely easier. Python ~ The current projection code uses Python and Pandas. It can work with HDF5 files. I do like Python as a language, but knowledge of that language is not widespread in our work group. R ~ As the analysis is in R one can consider R for doing the actual projection. We have HDF5 routines there as well. The advantage would be that the people who do not know Python could work on this part of the code. I choose R for the implementation. This will also allow bundling it closer with the Lüscher-Analysis package. # Design and implementation (Wolfram Language) Most of the code will be done in Mathematica. This section contains the design and implementation of the different tasks along the way. As the Wolfram Language is very functional and I know some Haskell, I made great use of this. The following table gives you pointers to some of the operators used and their correspondence in other languages. Mathematica functions are in general hyperlinked directly to the official documentation. You can also just type `?x` to get the help for `x`. | Long | Infix | Python | R | | --- | --- | --- | --- | --- | | `f[x]` | `f @ x` or `x // f` | `f(x)` | `f(x)` | | [`Apply`]`[f, a]` | `f @@ a` | `f(*a)`[^kwargs] | `do.call(f, a)` | | [`Map`]`[f, x]` | `f /@ x` | `map(f, x)` | `lapply(x, f)` | | [`Composition`]`[f, g]` | `f @* g` | —[^pycomp] | [`purrr::compose`](https://purrr.tidyverse.org/reference/compose.html) | | [`Function`]`[x, x^2]` | `#^2 &` | `lambda` | `function` | | [`Part`]`[xs, i]` | `xs[[i]]` | `xs[i-1]` | `xs[[i]]` | | [`Rule`] | `->` | — | — | | [`RuleDelayed`] | `:>` | — | — | | [`Association`]`[a -> b]` | `<| a -> b |>` | `{a: b}` | `list(a = b)`[^rlist] | | [`Dot`]`[a, b]` | `a . b` | `a @ b` | `a %*% b` | | [`NonCommutativeMultiply`] | `**` | — | — | | [`ReplaceAll`] | `/.` | — | — | | [`StringJoin`] | `<>` | `+` | `paste0` | [^kwargs]: In the Wolfram Language the [`List`] instance `a` could also contain [`Rule`] instances, therefore making a function call with positional and named arguments possible. R's `do.call` also uses `names(a)` as the names for named arguments. In Python one cannot mix a list and a dictionary, therefore one would have to call `f(*a, **b)` where `b` was a dictionary with string keys containing the named arguments. [^pycomp]: Not directly available, but [there are ways](https://stackoverflow.com/a/24047214/653152). [^rlist]: This is not entirely correct, though. In R, the names of the list have to be strings. For a given `list` instance `l`, we find that `typeof(names(l))` is `character` (homogenious string vector). This means that we can only have a mapping from strings to arbitrary types, but not from arbitrary types to arbitrary types. The Wolfram Language does not have this limitation. Also Python only requires the dictionary key to be hashable, otherwise there are no constraints. Note that in R one can create user-defined infix operators, so writing the following would allow writing R code similar to infix Wolfram Language code: ```r `%.%` <- purrr::compose `%/@%` <- function (f, l) lapply(l, f) `%@@%` <- do.call ``` Common patterns in functional programming are [*tacit programming*](https://en.wikipedia.org/wiki/Tacit_programming) and expressing loops in terms of list comprehensions or *maps*. I use some of the former and plenty of the latter. This might make the code harder to read for a programmer who is not familiar with the Wolfram Language. From what I gathered the built-in documentation mechanism in Mathematica is not as formal as say Doxygen in C++ or Rd in R. The formal one which is used for the official documentation is too much overhead for my little package. Therefore I chose to explain the functions from my `sLapHProjection` package towards the end of each section in this document. The notation of the function arguments will be mixed text and math. Sometimes it is more expressive to give a mathematical expression instead of just the name that the argument has in code. Optional parameters can be recognized by the default value behind a colon, which is the Wolfram Language notation for default values. Since the Wolfram Language is a functional one, the most straightforward way to define a constant is as a function with zero arguments. In Haskell there is no distinction between constants and argument less functions at all. You need to properly install the Quark Contraction Tool. For this download the `qct.m` from their [GitHub repository](https://github.com/djukanovic/qct) and then follow [Wolfram's installation instructions](http://support.wolfram.com/kb/5648). --- - **`MonitoredMap`**(function, list, label) A version of [`Map`] that just maps the *function* over the *list*. It shows a progress bar together with the given *label*. With sensible labels it can therefore be nested. - **`DropEmpty`**(container) Drops all elements from the container which are of length 0. ## Momenta For this code we use the momentum parametrization for $n$ particles with a total momentum $\vec d_\text{tot}$ and $n - 1$ relative momenta $\{ \vec q_i | i = 1, \ldots, n-1\}$. These are converted into individual momenta using the relations $$ \vec p_1 = \vec d_\text{tot} - \sum_{j = 1}^{n - 1} \vec q_j \qquad \text{and} \qquad \vec p_i = \vec q_{i - 1} \,. $$ In this section we describe some functions that work with momenta only and do not require group theory. For those functions see the section on the [Cartesian representation]. --- - **`MomentumToString`**($\vec p$) Converts an integer momentum vector like `{1, 0, -1}` into the string `"10-1"`. - **`RelativeToIndividualMomenta`**($\vec d_\text{tot}$, $\{ \vec q_i \}$) Converts the given total momentum $\vec d_\text{tot}$ and $n - 1$ relative momenta $\{ \vec q_i \}$ to a list of $n$ particle momenta $\{ \vec p_i \}$. - **`AllRelativeMomenta`**($\vec d_\text{tot}$, $\{ \vec q_i \}$, cutoff) If any of the momenta $\{ \vec p_i \}$ have a magnitude greater than *cutoff*, the result is an empty list. - **`UniqueTotalMomenta`**($\vec d^2$) Gives a list of unique total momenta corresponding to the given total momentum magnitude. For $\vec d^2 = 1, 2, 3, 4$, these are the following: {0, 0, 1}, {0, 0, -1}, {1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0} {1, 1, 0}, {1, -1, 0}, {-1, 1, 0}, {-1, -1, 0}, {0, 1, 1}, {0, 1, -1}, {0, -1, 1}, {0, -1, -1}, {1, 0, 1}, {1, 0, -1}, {-1, 0, -1}, {-1, 0, 1} {1, 1, 1}, {1, -1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, -1}, {-1, 1, 1}, {-1, -1, -1} {0, 0, 2}, {0, 0, -2}, {2, 0, 0}, {-2, 0, 0}, {0, 2, 0}, {0, -2, 0} - **`RelMomentaFromIndividual`**($\{ \vec p_i \}$) Takes a list of $n$ individual momenta and returns $n-1$ relative momenta. Due to our parametrization this is just taking every element except the last one. - **`RelMomentaRef`**($\vec d_\text{tot}$, $\{ \vec q_i \}$) Takes relative momenta in the lattice frame and rotates them back into the reference system using the inverse `MatrixRGTilde` operation. - **`RelMomentaRefFromIndividual`**($\{ \vec p_i \}$) Convenience function using `RelMomentaFromIndividual` and `RelMomentaRef`. - **`RelMomentaRefLabelFromIndividual`**($\{ \vec p_i \}$) Takes individual momenta and creates a label like `"000,001"` from this. This is a convenience function using `RelMomentaRefFromIndividual` and `MomentumToString`. - **`MomentaMaxNorm`**($\{ \vec p_i \}$) Computes $\max_i \| \vec p_i \|^2$. ## Group theory ### Data format for group elements and irreps In `sciebo/Lattice/Representations` we have a bunch of text files that contain the group elements and all the irreducible representations for the full group and subgroup. In the file `Oh-elements.txt` we find a CSV table with elements like the following: | Operator | $\alpha/\pi$ | $\beta/\pi$ | $\gamma/\pi$ | | :--- | ---: | ---: | ---: | | `E` | 0 | 0 | 0 | | `C2x` | 0 | 1 | 1 | | `C2y` | 0 | 1 | 0 | | `C2z` | 0 | 0 | 1 | | `C31+` | 0 | 0.5 | 0.5 | These are the group elements and their action parametrized in terms of the three Euler angles $\alpha$, $\beta$ and $\gamma$. Then there is `Little-group-elements.txt`, a TSV table that contains the subset of operators for each subgroup associated with a momentum vector $d$. The entries look like this: | $\vec{d}$ | $\mathrm{LG}(\vec{d})$ | Operator | $\alpha/\pi$ | $\beta/\pi$ | $\gamma/\pi$ | | --- | --- | --- | ---: | ---: | ---: | | (0, 0, 1) | `C4v` | `E` | 0 | 0 | 0 | | (0, 0, 1) | `C4v` | `C4z+` | 0 | 0 | $0.5$ | | (0, 0, 1) | `C4v` | `C4z-` | 0 | 0 | $-0.5$ | | (0, 0, 1) | C4v | C2z | 0 | 0 | $1$ | | (0, 0, 1) | `C4v` | `sigma_x` | 0 | 1 | $1$ | | (0, 0, 1) | `C4v` | `sigma_y` | 0 | 1 | $0$ | | (0, 0, 1) | `C4v` | `sigma_d1` | 0 | 1 | $0.5$ | | (0, 0, 1) | `C4v` | `sigma_d2` | 0 | 1 | $-0.5$ | The representations for the individual subgroups are in files like `C2v-(-1,-1,0)-representations.txt`. There also is `Oh-(0,0,0)-representations.txt` which contains all group elements because the little group is the group itself: | Irrep | Oh | Lg | $\alpha/\pi$ | $\beta/\pi$ | $\gamma/\pi$ | row | col | matrix element | | --- | --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | | `A1g` | `E` | `E` | 0 | 0 | 0 | 1 | 1 | 1. | | `A1g` | `C2x` | `C2x` | 0 | 1 | 1 | 1 | 1 | 1. | | `A1g` | `C2y` | `C2y` | 0 | 1 | 0 | 1 | 1 | 1. | | `A1g` | `C2z` | `C2z` | 0 | 0 | 1 | 1 | 1 | 1. | | `Eg` | `E` | `E` | 0 | 0 | 0 | 1 | 1 | 1. | | `Eg` | `E` | `E` | 0 | 0 | 0 | 1 | 2 | 0. | | `Eg` | `E` | `E` | 0 | 0 | 0 | 2 | 1 | 0. | | `Eg` | `E` | `E` | 0 | 0 | 0 | 2 | 2 | 1. | | `Eg` | `C2x` | `C2x` | 0 | 1 | 1 | 1 | 1 | 1. | | `Eg` | `C2x` | `C2x` | 0 | 1 | 1 | 1 | 2 | 0. | | `Eg` | `C2x` | `C2x` | 0 | 1 | 1 | 2 | 1 | 0. | | `Eg` | `C2x` | `C2x` | 0 | 1 | 1 | 2 | 2 | 1. | The “Irrep” column denotes the irrep name, “Oh” denotes the name of the group element as part of the $O_h$ group and “Lg” the group element name it has within the little group. Little groups belonging to the momenta that are related by a lattice rotation are equivalent to each other but contain different members of the original group. All group theoretical results can be found in human-readable form in the Mathematica notebook `analytic_projection/Summary_of_Matrices.nb`. In case it exists see the PDF file with the same name for an exported version. ### Lattice irreps The irrep matrices are represented in the *long format*. For instance that `C2x` group element in the `Eg` irrep is a $2 \times 2$ matrix that is written out as $$ D^{E^+}(C_{2x}) = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \,. $$ There are two ways that one can represent tensors: 1. A high-dimensional array of the actual values. Each dimension corresponds to one index. Contractions are performed by looping over all indices that stay and then sum up the elements along the contracted dimension. The tracking of the index labels it a bit tedious. This also does not scale well when indices are added later as the dimensionality of the arrays has to be increased. Without proper labelling, this will quickly spiral out of control. Memory usage is good, though. 2. A data frame in the *long data format*. Contractions are done by *group-by* and *summarize* operations. To contract an index one groups by all the indices that are to stay and summarizes the sum of the values. This scheme easily scales with arbitrary many indices and adding indices later on is not a problem. Memory usage is higher, though zeros can be omitted easily. It seems that Mathematica has the notion of [`Dataset`], which means that one can just `Import` a CSV or TSV file and has a data frame as with R or Pandas. Also the `Dataset` supports having a complex matrix as column type, which would allow for a hybrid approach if we wanted to. 3. There is another way, using nested [`Association`] instances. They are just heterogeneous mappings from one thing to another. This corresponds to the `dict` in Python. The `std::map` in C++ is similar, but requires the same type among the keys and the same time among the values. Using this for our implementation, we could make associations that map from some named group element to an association that maps from indices to a c-number. As associations can be called like functions, we can really nicely work with them. So let's make up a nonsense association like this: ```mathematica assoc = Association[foo -> bar, 1 -> "one", x^2 -> Association] ``` When you now call `assoc[x^2]`, you get the function `Association`. Or you call it with `assoc[1]` and get `"one"`. Missing keys like `bar` are some sort of error condition via say `Missing["KeyAbsent", bar]`. Markus was immediately sold on this approach 3, and I quite like it as well. This is now implemented. --- - **`ReadDataframe`**(filename) Reads a CSV formatted table with a header and generates a `Dataset` from it. - **`ConvertMatrixElement`**(dataset) The matrix elements of the irreducible representations, $D^\Gamma_{\alpha\beta}(g)$ are stored as strings in the CSV file. Using [`ToExpression`] these are parsed into Mathematica symbols. This function takes the dataset and replaces the “matrix element” column with the parsed one. - **`ReadIrreps`**(filename) The composition of `ConvertMatrixElement` and `ReadDataframe`. - **`MatrixElementsToAssociation`**(dataset) Converts the “row”, “col” and “matrix element” into mappings of the form $(\alpha, \beta) \to D^\Gamma_{\alpha\beta}(g)$. The result is a [`Dataset`] with just one row per irrep and little group element and an [`Association`] which contains all the matrix elements. - **`IrrepsToAssociation`**(matrixElementsAssociations) Continuing with the result from `MatrixElementsToAssociation` this function then builds a more deeply nested [`Association`] which has the form $$ \Gamma \to g \to (\alpha, \beta) \to D_{\alpha\beta}^\Gamma(g) \,. $$ - **`DatasetToAssociations`**(dataset) Composition of `IrrepsToAssociation` and `MatrixElementsToAssociation`. - **`ExtractMomentumFromFilename`**(filename) Given a filename like `C2v-(0,1,1)-representations.txt` it extracts the three-momentum vector $\vec p_\text{cm}$. - **`IrrepDGammaAssoc`**() The $D_{\alpha\beta}^\Gamma(g)$ are represented as an association of this form: $$ \vec d \to \Gamma \to g \to (\alpha, \beta) \to D_{\alpha\beta}^\Gamma(g) \,. $$ You can think of this as a function taking a momentum vector $\vec d$ and giving you a new function which accepts an irrep $\Gamma$ and gives you a function which …. - **`IrrepSize`**($\vec d_\text{tot}$, $\Gamma$) Gives the number of rows (equal to number of columns) that a given irrep has. ### Cartesian representation In the file `Oh-elements.txt` we just have the names of the group elements and the Euler angles. The Cartesian representation comes for free, there is [`EulerMatrix`] which just takes such a vector. This will give us the mapping $g \to R_g$ that we need. Great, we're done here if the definition of the Euler angles matches. --- - **`ReadEulerAngles`**(filename) Reads the file with Euler angles and produces a handy association of the form $$ g \to (\alpha, \beta, \gamma) \,.$$ - **`EulerAnglesAssoc`**() All the Euler angles already read in available as a global constant. - **`MomentumRefScalar`**($|\vec p_\text{cm}|^2$) Gives the reference momentum. These are hand-chosen to be: $$ \begin{aligned} 0 &\to (0, 0, 0) \\ 1 &\to (0, 0, 1) \\ 2 &\to (1, 1, 0) \\ 3 &\to (1, 1, 1) \\ 4 &\to (0, 0, 2) \,. \end{aligned} $$ This only works only up to $|\vec p_\text{cm}|^2 \leq 8$ because after that the mapping is not unique any more, one can write 9 as $9 = 0^2 + 0^2 + 3^2$ but also as $9 = 2^2 + 2^2 + 1^2$. - **`MomentumRef`**($\vec p_\text{cm}$) Gives the reference momentum $\vec p_\text{ref}$ for given center-of-mass momentum $\vec p_\text{cm}$. As this function is implemented via `MomentumRefScalar`, the same restrictions apply. - **`EulerGTilde`**($\vec p_\text{cm}$) Finds a set of Euler angles $(\psi, \theta, \phi)$ such that the rotation matrix converts the reference momentum into center-of-mass momentum, i.e. $$\vec p_\text{cm} = R(\psi, \theta, \phi) \; \vec p_\text{ref} \,.$$ This transformation is not unique, but that *should* not be a problem as we sum over the little group elements anyway. - **`MatrixRGTilde`**($\vec p_\text{cm}$) Simply computes the rotation matrix corresponding to the Euler angles from \texttt{EulerGTilde}. - **`MomentumTransform`**($\vec p$, $\vec p_\text{cm}$, $(\psi, \theta, \phi)$) The momentum transformation that is needed in the argument of the single-particle operator: $$ R_{\tilde g}^{-1} R_g R_{\tilde g} \,. $$ - **`RotateMomenta`**($g$, $\{ \vec p_i \}$) Takes the group element name as a string and a set of momenta. Then it rotates all the momenta using the cartesian representation of the given group element. - **`MomentaOrbit`**($\{ \vec p_i \}$) Computes the whole orbit of the given momentum set under the rotations according to the little group corresponding to $\vec P = \sum_i \vec p_i$. The result is a list of momentum sets. For instance we can take $p_1 = (0, 0, 1)$ and $p_2 = (0, 0, -1)$. Then the according stabilizer subgroup will the the whole of the octahedral group. We can then run this: ```mathematica momentaOrbit[{{0, 0, 1}, {0, 0, -1}}] ``` And get the following result back: ```mathematica {{{-1, 0, 0}, {1, 0, 0}}, {{0, -1, 0}, {0, 1, 0}}, {{0, 0, -1}, {0, 0, 1}}, {{0, 0, 1}, {0, 0, -1}}, {{0, 1, 0}, {0, -1, 0}}, {{1, 0, 0}, {-1, 0, 0}}} ``` - **`RemoveRedundantMomenta`**($\{\{ \vec p_i \}\}$) Takes a list of sets of momenta and removes the ones that can be obtained via little group rotations. All these combinations must have the same total momentum, otherwise the results are likely strange. The picked set of unique momenta is deterministic but likely different from the hand-picked set that Markus Werner uses. - **`ParityEulerMatrix`**($\{ \alpha, \beta, \gamma \}$, $P$) Multiplies the Euler matrix associated with the three Euler angles with the additional parity factor $P$. - **`GetParity`**($\vec d_\text{tot}$, $\{ \alpha, \beta, \gamma \}$) Determines the parity eigenvalue $P$ for a given rotation by solving $$ P R(\alpha, \beta, \gamma) \vec d_\text{tot} = \vec d_\text{tot} $$ for $P$. $R$ is the cartesian rotation matrix corresponding to the three Euler angles. - **`Parities`**() An association which maps group element names to either a parity eigenvalue of either `+1` or `-1`. These values are determined from the $A_1^-$ irrep. ### Spin representation The spin-$J$ representations are for free as well, we just use [`WignerD`] which is a mapping $$ (j, m_1, m_2, \psi, \theta, \phi) \to D^J(g) \,. $$ ### Clebsch-Gordan coefficients The Clebsch-Gordan coefficients are implemented as [`ClebschGordan`]. We do need to implement higher Clebsch-Gordan coefficients ourselves, though. This is the relation that I came up with by figuring that we first couple two spins and then couple that to the next one. \begin{multline*} \langle J, M | j_1, m_1, j_2, m_2, j_3, m_3, \ldots \rangle \\= \sum_{\tilde J = |j_2 - j_3|}^{j_2 + j_3} \sum_{\tilde M = - \tilde J}^{\tilde J} \langle J, M | j_1, m_1, \tilde J, \tilde M \rangle \langle \tilde J, \tilde M | j_2, m_2, j_3, m_3, \ldots \rangle \,. \end{multline*} If the “$\ldots$” are empty, then we have the usual Clebsch-Gordan coefficients that are already available. This is the stop of the recursion. Creating a recursion relation is not easy at is seems. I have [asked on Physics Stack Exchange](https://physics.stackexchange.com/questions/459408/clebsch-gordan-coefficients-for-more-than-2-particles), and in one answer I was told that this is correct though not unique. I was pointed to the [Racah $W$ coefficient](https://en.wikipedia.org/wiki/Racah_W-coefficient). These do not seem to be implemented in Mathematica, though. For coupling pions these factors will all be 1 anyway, so we could consider not bothering with them for just now. --- - **`HigherClebschGordan`**($\{ j_i \}$, $\{ m_i \}$, $(J, M)$) Computes the Clebsch-Gordan coefficient for coupling all the spins $(j_i, m_i)$ together to $(J, M)$. Note that the API is different from [`ClebschGordan`] as here the $j_i$ and $m_i$ are grouped among each other instead of grouped by $i$. In contrast to [`ClebschGordan`] this function does not emit any warnings in case the coupling is unphysical. ## Spin We implement spin independent of isospin here. The spin projection will just create an expression that contains `SingleOperator`, which can then later be replaced with the actual operators used. But isospin projection gives us one big multi-particle operator. I let it compute with the single operators and then massage the expression such that the products of operators are directly adjacent. Then I use a pattern replacement to replace them all with one multi-particle operator. Markus has chosen the phase vector $\tilde \phi_\beta$ as $\delta_{\beta 1}$, effectively fixing the column of the irrep. In my implementation I just expose it as a parameter to the user and remove the sum over the $\beta$. The formula is implemented from the right to left. We have `MakeSingleOperator` applying the Wigner-$D$-matrices to the single operators. Then `MakeMultiOperator` creates a product of the single operators. `MakeGroupSum` sums over these operator products, `MakeMagneticSum` sums over the $M$ and $M_i$ quantum numbers. The last step is not finished as the higher Clebsch-Gordan coefficients are not finished. From the spin projected operator we extract the actual momenta --- - **`MakeSingleOperator`**($\vec p_i$, $\vec p_\text{cm}$, $\vec \Psi_g$, $J_i$, $M_i$, $i$) Creates an appropriate linear combination of the single particle operator $O$ from $$ \sum_{M_i'=-J_i}^{J_i} \sum_{M_i''=-J_i}^{J_i} D_{M_i' M_i}^{J_i}(g) \; D_{M_i'' M_i'}^{J_i}(\tilde g) \; O_{i M_i''}^{J_i}(R_{\tilde g}^{-1} R_g R_{\tilde g} \vec p_i)^\dagger \,. $$ The group element $g$ is specified via its Euler angles $\vec \Psi_g$. For instance taking $\vec p_i = (1, 1, 1)$, $\vec p_\text{cm} = (1, 1, 1)$, $\vec \Psi_g = (\pi/2, 1, 1)$, $J_i = 1$, $M_i = 1$, $i = 1$ gives us the following expression: ```mathematica I ConjugateTranspose[SingleOperator[1, 1, -1, {1, 1, -1}]] ``` The function `SingleOperator` is just a placeholder defined next. We can also take a more interesting case where the Wigner $D^J$ matrices are not so simple: ```mathematica MakeSingleOperator[{0, 0, 1}, {1, 1, 1}, Pi {1/2, 1/2, 0}, 2, 0, 1] ``` We then get the following linear combination of operators: ```mathematica -(1/2) Sqrt[3/2] ConjugateTranspose[SingleOperator[1, 2, -2, {0, 1, 0}]] - 1/2 ConjugateTranspose[SingleOperator[1, 2, 0, {0, 1, 0}]] - 1/2 Sqrt[3/2] ConjugateTranspose[SingleOperator[1, 2, 2, {0, 1, 0}]] ``` - **`SingleOperator`**($i$, $J_i$, $M_i$, $\vec p$) Just a placeholder without a definition. It holds the single particle operator $O_{i M_i}^{J_i}(\vec p)$. These will later be replaced with the isospin projected operators. - **`MakeMultiOperator`**($\{ \vec p_i \}$, $\vec \Psi_g$, $\{ J_i \}$, $\{ M_i \}$) Creates the multi particle operator as a product of single particle operators. The center-of-mass momentum $\vec p_\text{cm}$ is taken to be the sum of the individual momenta $\vec p_i$. It corresponds to this expression: $$ \prod_{i=1}^{N_\text{p}} \mathtt{SingleOperator}(i, J_i, M_i, \vec p) \,. $$ For example we use the following code: ```mathematica MakeMultiOperator[{{1, 1, 1} - {1, 0, 0}, {1, 0, 0}}, {1, 0, Pi}, {1, 1}, {1, 0}] ``` This then gives us the following product of two single operators: ``` (-E^I ConjugateTranspose[ SingleOperator[1, 1, 1, {Sin[1], -Cos[1], 1}]]) ** ConjugateTranspose[SingleOperator[2, 1, 0, {-Cos[1], -Sin[1], 0}]] ``` - **`MakeGroupSum`**($\Gamma$, $\alpha$, $\beta$, $\{ \vec p_i \}$, $\{ J_i \}$, $\{ M_i \}$) This computes the sum over the group elements: $$ \sum_{g \in \mathop{\mathrm{LG}}(\vec p_\text{cm})} D_{\alpha\beta}^\Gamma(g)^* \; \mathtt{MakeMultiOperator}(\{ \vec p_i \}, \vec \Psi_g, \{ J_i \}, \{ M_i \}) $$ The center-of-mass momentum is computed as the sum of the individual momenta. - **`MakeMagneticSum`**($\Gamma$, $\alpha$, $\beta$, $\{ \vec p_i \}$, $\{ J_i \}$, $\phi$) This function performs the sum over all the magnetic quantum numbers $M$ and $\{ M_i \}$. The resulting expression is \begin{align*} O_\Gamma^{J\alpha\beta}(\vec p_\text{cm})^\dagger &= \sum_{M=-J}^J \phi_M \left[ \prod_{i=1}^{N_\text{p}} \sum_{M_i=-J_i}^{J_i} \right] \langle J, M | J_1, M_1, \ldots, J_{N_\text P}, M_{N_\text P} \rangle \\&\quad\times \texttt{MakeGroupSum}(\Gamma, \alpha, \beta, \{ \vec p_i \}, \{ J_i \}, \{ M_i \}) \,. \end{align*} - **`ExtractMomenta`**(expr) Extracts all the momenta from the `SingleOperator`s used. The resulting expression contains summand of this form: ```mathematica``` (-(1/2) + (I Sqrt[3])/2) DTMomenta[{1, 0, 1}, {0, 1, 0}] ``` By replacing the `DTMomenta` with something else, one can extract just the momenta from this. One has to be a bit careful as Mathematica will try to distribute the common factor onto the elements in the lists if one replaces it with something simple like [`List`]. Therefore we use something which Mathematica does not know and therefore does not try to simplify. - **`momentaToRules`**(momenta, location) --- [not exported] Takes a product of momentum [`List`]s in a [`NonCommutativeMultiply`] and converts them into an [`Association`]. *location* shall be either `"so"` or `"si"` Calling this function with ```mathematica MomentaToRules[{0, 1, 1} ** {1, 0, 0}, "si"] ``` gives: ```mathematica <|"psi1" -> "011", "psi2" -> "100"|> ``` - **`DTMomenta`** An undefined symbol which acts as a data type. Used with [`ReplaceAll`]. - **`MomentaToAssoc`**($\{ \vec p_i \}$, location) Converts the momentum terms with [`Association`]s containing terms like ```mathematica DTMomenta[{0, 1, 1}, {1, 0, 0}] ``` to an association like this: ```mathematica DTMomentaAssoc[<|"psi1" -> "011", "psi2" -> "100"|>] ``` The name is build up from `p`, then the *location* and consecutive numbers. It is wrapped in `DTMomentaAssoc` such that Mathematica does not try to distribute numeric factors to the arguments. - **`DTMomentaAssoc`** An undefined symbol which acts as a data type. Used with [`ReplaceAll`]. - **`MomentaToAssocSourceSink`**(source, sink) Takes two expressions containing `DTMomenta` for the source and sink. They get converted into `DTMomentaAssoc` expressions like the following: ```mathematica DTMomentaAssoc[<|"psi1" -> "011", "psi2" -> "100", "pso1" -> "011", "pso2" -> "100"|>] ``` ## Isospin Isospin and the Wick contractions are independent of the spin. Therefore we can do the contractions on the isospin operators and then pattern match that to the spin expressions. This allows us to go a bit further before connecting isospin and spin. We define a $\pi^+$ operator that has two spin, one color and one position/momentum index as follows: ```mathematica \[Pi]Plus[s1_, s2_, c1_, x1_] := -FieldB[up, c1, s1, x1] ** (Gamma^5)[SI[{s1, s2}]] ** Field[dn, c1, s2, x1]; ``` Similarly we define the other two pion operators. From these we build for instance an $I = 2$ operator. ```mathematica \[Pi]\[Pi]I2[s1_, s2_, s3_, s4_, c1_, c2_, x1_, x2_] := \[Pi]Plus[s1, s2, c1, x1] ** \[Pi]Plus[s3, s4, c2, x2]; ``` For reasons not entirely clear to me we need to define the “bar” version of that operator ourselves. This likely is because the QCT does not assume anything about the Dirac algebra behind the `Gamma` and therefore cannot really do anything about that. Anyway, the second operator is this: ```mathematica \[Pi]\[Pi]I2Bar[s1_, s2_, s3_, s4_, c1_, c2_, x1_, x2_] := \[Pi]Minus[s3, s4, c2, x2] ** \[Pi]Minus[s1, s2, c1, x1]; ``` A few pion operators are defined as part of the library such that they do not have to be created for every projection projection from scratch. --- - **`\[Pi]Plus`**($s_1$, $s_2$, $c_1$, $x_1$) $$ [\pi^+]_{s_1 s_2}^{c_1}(x_1) = - \bar u_{s_1}^{c_1}(x_1) \; \gamma^5_{s_1 s_2} \; d_{s_2}^{c_1}(x_1) $$ - **`\[Pi]Minus`**($s_1$, $s_2$, $c_1$, $x_1$) $$ [\pi^+]_{s_1 s_2}^{c_1}(x_1) = \bar d_{s_1}^{c_1}(x_1) \; \gamma^5_{s_1 s_2} \; u_{s_2}^{c_1}(x_1) $$ - **`\[Pi]Zero`**($s_1$, $s_2$, $c_1$, $x_1$) $$ [\pi^0]_{s_1 s_2}^{c_1}(x_1) = \frac{1}{\sqrt 2} \left( \bar u_{s_1}^{c_1}(x_1) \; \gamma^5_{s_1 s_2} \; u_{s_2}^{c_1}(x_1) - \bar d_{s_1}^{c_1}(x_1) \; \gamma^5_{s_1 s_2} \; d_{s_2}^{c_1}(x_1) \right) $$ - **`\[Pi]\[Pi]I1`**($s_1$, $s_2$, $s_3$, $s_4$, $c_1$, $c_2$, $x_1$, $x_2$) $$ [\pi\pi^{I=1}]_{s_1 s_2 s_3 s_4}^{c_1 c_2}(x_1, x_2) = [\pi^+]_{s_1 s_2}^{c_1}(x_1) \; [\pi^-]_{s_3 s_4}^{c_2}(x_2) - [\pi^+]_{s_3 s_4}^{c_2}(x_2) \; [\pi^-]_{s_1 s_2}^{c_1}(x_1) $$ - **`\[Pi]\[Pi]I1Bar`**($s_1$, $s_2$, $s_3$, $s_4$, $c_1$, $c_2$, $x_1$, $x_2$) $$ [\bar{\pi\pi}^{I=1}]_{s_1 s_2 s_3 s_4}^{c_1 c_2}(x_1, x_2) = [\pi^-]_{s_3 s_4}^{c_2}(x_2) \; [\pi^+]_{s_1 s_2}^{c_1}(x_1) - [\pi^-]_{s_1 s_2}^{c_1}(x_1) \; [\pi^+]_{s_3 s_4}^{c_2}(x_2) $$ - **`\[Pi]\[Pi]I2`**($s_1$, $s_2$, $s_3$, $s_4$, $c_1$, $c_2$, $x_1$, $x_2$) $$ [\pi\pi^{I=2}]_{s_1 s_2 s_3 s_4}^{c_1 c_2}(x_1, x_2) = [\pi^+]_{s_1 s_2}^{c_1}(x_1) \; [\pi^+]_{s_3 s_4}^{c_2}(x_2) $$ - **`\[Pi]\[Pi]I2Bar`**($s_1$, $s_2$, $s_3$, $s_4$, $c_1$, $c_2$, $x_1$, $x_2$) $$ [\bar{\pi\pi}^{I=2}]_{s_1 s_2 s_3 s_4}^{c_1 c_2}(x_1, x_2) = [\pi^-]_{s_3 s_4}^{c_2}(x_2) \; [\pi^-]_{s_1 s_2}^{c_1}(x_1) $$ ## Wick contractions With the appropriate isospin operators defined, we can perform the Wick contraction using the “quark contraction tool”: ```mathematica wc = WickContract[\[Pi]\[Pi]I2Bar[s1, s2, s3, s4, c1, c2, x1, x2] ** \[Pi]\[Pi]I2[s5, s6, s7, s8, c5, c6, x3, x4]]; ``` We can then perform the quark contractions on this: ```mathematica qc = QuarkContract @ wc ``` When we take a look at the output (using `TraditionalForm`) we get something that looks like this:[^vimregex] \begin{align*} &\text{tr}\left(\gamma^5\;S^{\text{up}}(p_4,p_1)\;\gamma^5\;S^{\text{dn}}( p_1,p_4)\right) \text{tr}\left(\gamma^5\;S^{\text{up}}(p_3,p_2)\;\gamma^5\;S^{\text{dn }}(p_2,p_3)\right) \\&\quad +\text{tr}\left(\gamma^5\;S^{\text{up}}(p_3,p_1)\;\gamma^5\;S^{\text{dn}}(p_1,p_3)\right) \text{tr}\left(\gamma^5\;S^{\text{up}}(p_4,p_2)\;\gamma^5\;S^{\text{dn }}(p_2,p_4)\right) \\&\quad -\text{tr}\left(\gamma^5\;S^{\text{dn}}(p_2,p_3)\;\gamma^5\;S^{\text{up}}(p_4,p_2)\;\gamma^5\;S^{\text{dn}}(p_1, p_4)\;\gamma^5\;S^{\text{up}}(p_3,p_1)\right) \\&\quad -\text{tr}\left(\gamma^5\;S^{\text{up}}(p_4,p_1)\;\gamma^5\;S^{\text{dn}}(p_2,p_4)\;\gamma^5\;S^{\text{up}}(p_3,p_2)\;\gamma^5\;S^{\text{dn}}(p_1,p_3)\right) \end{align*} [^vimregex]: Unfortunately Mathematica does not give this output directly like that, it needs `TeXForm` and these replacements in Vim to get it as above: ```vim s/\\text{Gamma}/\\gamma/g s/trace/tr/g s/\\text{x\(\d\)}/p_\1/g s/\./\\;/g ``` The result is something that we will recognize as “C4cD” and “C4cC”. If we are at zero momentum and there is no distinction between them, this can be simplified down further to $$ 2 \cdot \text{C4cD} - 2 \cdot \text{C4cC} \,. $$ In this general case we cannot do that, though. ### Normalizing trace expressions In the contraction code in [`src/DiagramSpec.cpp`](https://github.com/HISKP-LQCD/sLapH-contractions/blob/8433f69825f5d5b097f6f30c045af037b06997b9/src/DiagramSpec.cpp) we have documentation describing the diagrams that can be computed: Box: $$ C_\text{C4cB} = \langle \Gamma_\mathtt{Op0} D_\mathtt{Q1}^{-1}(t|t') \Gamma_\mathtt{Op1} \gamma_5 D_\mathtt{Q2}^{-1}(t'|t')^\dagger \gamma_5 \Gamma_\mathtt{Op2} D_\mathtt{Q3}^{-1}(t'|t) \Gamma_\mathtt{Op3} \gamma_5 D_\mathtt{Q0}^{-1}(t|t)^\dagger \gamma_5 \rangle $$ Cross: $$ C_\text{C4cC} = \langle \Gamma_\mathtt{Op0} D_\mathtt{Q1}^{-1}(t|t') \Gamma_\mathtt{Op1} \gamma_5 D_\mathtt{Q2}^{-1}(t|t')^\dagger \gamma_5 \Gamma_\mathtt{Op2} D_\mathtt{Q3}^{-1}(t|t') \Gamma_\mathtt{Op3} \gamma_5 D_\mathtt{Q0}^{-1}(t|t')^\dagger \gamma_5 \rangle $$ Direct: $$ C_\text{C4cD} = \langle \Gamma_\mathtt{Op0} D_\mathtt{Q1}^{-1}(t|t') \Gamma_\mathtt{Op1} \gamma_5 D_\mathtt{Q0}^{-1}(t|t')^\dagger \gamma_5 \rangle \langle \Gamma_\mathtt{Op2} D_\mathtt{Q3}^{-1}(t|t') \Gamma_\mathtt{Op3} \gamma_5 D_\mathtt{Q2}^{-1}(t|t')^\dagger \gamma_5 \rangle $$ Vacuum: $$ C_\text{C4cV} = \langle \Gamma_\mathtt{Op0} D_\mathtt{Q1}^{-1}(t|t) \Gamma_\mathtt{Op1} \gamma_5 D_\mathtt{Q0}^{-1}(t|t)^\dagger \gamma_5 \rangle \langle \Gamma_\mathtt{Op2} D_\mathtt{Q3}^{-1}(t'|t') \Gamma_\mathtt{Op3} \gamma_5 D_\mathtt{Q2}^{-1}(t'|t')^\dagger \gamma_5 \rangle $$ In order to translate the expressions that we get from the Quark Contraction Tool we need to massage them until they match a diagram from the contraction code. There are three things that need to be done: - The trace is cyclic. Therefore similar expressions can be equivalent. The trace expression used, `trace` from `qct`, does not let Mathematica know that it is cyclic. Therefore we manually cycle them until the up-type propagator starting from the lowest vertex (to be defined) is the second element. The first element shall be the $\Gamma$-structure. The first vertex is the source vertex with the lowest number. If there is no source vertex, then the lowest sink vertex is to be used. The key observation is that [expressions can be used like lists](https://reference.wolfram.com/language/tutorial/PartsOfExpressions.html). With this we can access the different parts of expressions. The arguments of a function are subexpressions. - The direction of quark flow might be in the wrong direction. By adding a transposition to the argument of the trace one can reverse all the propagators. This might change the value of the trace if the Dirac structures give some signs under transposition. - In all the charged/conjugated diagrams the position of the down-type propagators is fixed. If from the Wick contraction we happen to have them in other positions we need to exchange all flavors in order to get to the form that is implemented. This is done by applying a conjugate transpose to the argument of the trace. The resulting diagram will need to be conjugated. --- - **`RotateGammaToFront`**(expression) Takes an expression like ```mathematica Gamma^5.DE[{"dn", "dn"}, {so[2], si[1]}]. Gamma^5.DE[{"up", "up"}, {si[2], so[2]}] ``` and cyclicly permutes the expression if it does not start with a `Gamma`. The resulting expression always starts with `Gamma`. - **`Starts`**(traceContent) Gives the starting times of the propagators in the trace content expression. The result could look like this: ```mathematica {si[1], so[2], si[2], so[1]} ``` - **`StartScore`**(propagator) Gives a sorting score for propagators based on their flavor and time slices. *up* and *down* get weighted with 1000 and 2000, respectively. The time slice index at source gets added with a scaling of 1, the time slice index at sink gets added with a scaling of 100. - **`IndexOfFirst`**(traceContent) Determines which propagator should be the first one in the expression. If the result is $n$, then the whole expression should be rotated left by $2(n - 1)$ iterations. It uses the `StartScore` to figure out the ordering. - **`NormalizeTrace`**(traceContent) Cyclicly permutes the trace such that the Dirac structure corresponding to the first vertex is at the front. - **`NormalizeTraceRecursive`**(expression) Recursively traverses an expression and applies `NormalizeTrace` on every `trace` expression encountered. The recursion can encounter three cases: - The length of the expression is zero, it is just a simple number or so. In that case the value is just the expression given. - The expresssion matches `trace[_]`, where `_` is a placeholder. In this case the argument is run through `NormalizeTrace`. - Else this function is mapped over all parts of the expression. - **`ReplacePropagators`**(expression) Replaces the QCT propagators with a simpler expression. For some reason the propagators contain the flavor twice redundantly and also contain sink and source time. For the later replacement with diagrams it is sufficient to have the flavor once and the source time only. This function does this replacement: ```mathematica qct`DE[{f_, f_}, {_, t_}] :> prop[f, t] ``` - **`FlowReversalRules`**() A list of rules that perform reversal of quark flow by taking the transpose of the content of the trace. We would need to know what the transpose of the Dirac structure is, which depends on the basis, which is not implemented yet. Therefore this is limited to the case of $\gamma_5$ which will be sufficient as long as we scatter just pions. - **`FlavorSwitchRules`**() A list of rules that exchange the flavors and might also change the direction of quark flow. Similar to the `FlowReversalRules` this currently only supports the $\gamma_5$ Dirac structure. ### Contractions into dataset name templates At the stage of the Wick contractions we still know which labels correspond to source and sink vertices. We can track this information by adding the undefined functions `so` and `si` which label the source and sink vertices instead of `x1` and so on. The contractions are then defined with operators evaluated like so: ```mathematica wc = WickContract[ \[Pi]\[Pi]I2Bar[s1, s2, s3, s4, c1, c2, so[1], so[2]] ** \[Pi]\[Pi]I2[s5, s6, s7, s8, c5, c6, si[1], si[2]]]; ``` In the contracted form we have propagator expressions like these: ```mathematica prop["up", so[2]] ``` These have to be mapped to HDF5 dataset names. Actually they have to be mapped to templates ([`StringTemplate`]) that can then be used with [`TemplateApply`] to insert the momenta. The following replacement can be used to replace an expression with a `C4cD` diagram: ```mathematica qct`trace[qct`Gamma^g1_ . prop["up", so[so1_]]. qct`Gamma^g2_ . prop["dn", si[si2_]]] qct`trace[qct`Gamma^g3_ . prop["up", so[so3_]]. qct`Gamma^g4_ . prop["dn", si[si4_]]] :> TemplateApply[ "C4cD_uuuu_" <> MakeTemplate[4], <|"g1" -> g1, "g2" -> g2, "g3" -> g3, "g4" -> g4, "x1" -> "`pso" <> ToString @ so1 <> "`", "x2" -> "`psi" <> ToString @ si2 <> "`", "x3" -> "`pso" <> ToString @ so3 <> "`", "x4" -> "`psi" <> ToString @ si4 <> "`"|>], ``` The resulting expression after applying these rules is a sum of strings, which is exactly what we want. In a later stage we can insert the momenta corresponding to the spin projected operators. --- - **`MakeTemplate`**($n$) Creates a template for an $n$-particle correlation function consisting of expressions like ```mathematica "g`g1`.p`x1`.d000" ``` separated by underscores and with increasing numbers for `g1` and `p1`. - **`DatasetNameRules`**() List of delayed rules ([`RuleDelayed`]) that can be used with with [`ReplaceAll`] to replace propagator expression with HDF5 dataset name templates ([`StringTemplate`]). The replacements are expressions like these: ```mathematica "C2cD_uu_g5.p`pso1`.d000_g5.p`psi1`.d000" ``` The following mappings are performed: ## Wick contraction and spin We now have the needed linear combination of contraction code diagrams. And using `MomentaToAssocSourceSink` we have a prescription for inserting the momenta. Inserting every set of momenta into every term of the diagram expression gives us the desired result, a linear combination of actual HDF5 dataset names. --- - **`CombineIsospinAndSpin`**(correlator templates, momenta associations) Takes an expression of *correlator templates* containing strings like ```mathematica "C4cD_uuuu_g5.p`pso1`.d000_g5.p`psi1`.d000_g5.p`pso2`.d000_g5.p`psi2`.d000" ``` and an expression containing *momenta associations* like ```mathematica DTMomentaAssoc[<|"pso1" -> "011", "pso2" -> "100", "psi1" -> "011", "psi2" -> "100"|>] ``` and gives an expression containing strings like these: ```mathematica "C4cD_uuuu_g5.p001.d000_g5.p011.d000_g5.p110.d000_g5.p100.d000" ``` This expression then contains everything that is needed to project actual data. ### Export to data frame We want to export the expression in a form like this: | datasetname | conjugate | re | im | | --- | --- | ---:| ---: | | `C4cB_uuuu_g5.p001.d000_…` | False | 1. | 0. | | `C4cB_uuuu_g5.p001.d000_…` | False | -0.5 | -0.8660254037844386 | | `C4cB_uuuu_g5.p001.d000_…` | False | 0.5 | -0.8660254037844386 | | `C4cB_uuuu_g5.p001.d000_…` | False | -0.5 | 0.8660254037844386 | In CSV format it will be this: ```csv datasetname,conjugate,re,im "C4cB_uuuu_g5.p001.d000_…","False",1.,0. "C4cB_uuuu_g5.p001.d000_…","False",-0.5,-0.8660254037844386 "C4cB_uuuu_g5.p001.d000_…","False",0.5,-0.8660254037844386 "C4cB_uuuu_g5.p001.d000_…","False",-0.5,0.8660254037844386 ``` --- - **`StringExpressionToAssociation`**(expression) Takes an expression containing a linear combination of strings and converts them into an association where the strings are the keys and the c-number valued factors are the values. It uses [`Coefficient`] to extract the leading coefficient in front of the strings. Also it prefixes the string `conj:` in front of diagram names that need to be conjugated. This function is based on a [a post by Kuba](https://mathematica.stackexchange.com/a/191718/1507) with an implicit MIT license. - **`NeedsConjugation`**(name) Checks whether the *name* starts with `conj:` and returns the list (bare name, True/False), where “bare name” is the name without the `conj:`. - **`DatasetnameAssocToCSV`**(association, filename) Takes the *association* which maps HDF5 dataset names to complex numbers (like the result from `StringExpressionToAssociation`) and converts that into a CSV structure with these columns: - `datasetname` - `conjugate` - `re` - `im` This function is based on a [a post by Kuba](https://mathematica.stackexchange.com/a/191718/1507) with an implicit MIT license. ## Creating a full correlator matrix ### Interface to numerical code So far we only have functions that generate a list of HDF5 dataset names that will yield a single correlation function. We of course want to build a full correlator matrix (with multiple correlators) and then many of these correlator matrices. This means that we will have the following independent variables describing the correlators: - All total momenta $\vec d$ for given $\vec d^2$ (and not just total momentum $\vec d^2$) - Irrep - Irrep rows - Correlator matrix row and column (operator ids & relative momenta $\vec q_i$) The nature of the operator and the relative momenta shall be passed down into the analysis. This can be done with the meta data notation that Markus has implemented in his projection code. Although it would be nice to use the same parser in R, writing it out in a machine readable format would be even better. Therefore I favor writing out the meta data in the YAML format. The information per correlator (`datasetname`, `conjugate`, `re`, `im`) can be written out as CSV with the current implementation. But how will we do multiple correlators? I see these options: 1. We can just write out one CSV file per correlator matrix element. The information about $\vec d^2$, $\Gamma$, $\alpha$ and the $\{ \vec q_i \}$ would then be encoded in the file name. This would be a classic thing, but these days I really dislike parsing filenames for information. 2. Have one huge CSV file, where there are columns for everything. Then use *group by* and *summarize* operations to boil it down. The result would then also be a correlator in the *long data format*. This might be too inflexible to work with. 3. Use a nested structure with a YAML representation. This way one can inject arbitrary amounts of meta data at every point in the structure. Also it could be parsed quite easily in R. But can it be generated from Wolfram Language? Apparently there are some third-party implementations, but the Wolfram Language can export JSON. I am perfectly fine with that, converting it to YAML is trivial with some Python or R script. As usual with these listings, I prefer the last option. The interface will be a an association of the following form: $$ d_\text{tot} \to \Gamma \to \beta \to \alpha \to O_i \to O_j \to \{ (C, \text{Re}, \text{Im}, \dagger) \} \,. $$ A typical JSON file looks like this: ```js { "001": { "E": { "1": { "1": { "0-10": { "0-10": [ { "conj": false, "datasetname": "C4cB_uuuu_p0-10.d000.g5_p0-10.d000.g5_p01-1.d000.g5_p011.d000.g5", "im": 0.0, "re": -0.0078125 }, { "conj": false, "datasetname": "C4cB_uuuu_p0-10.d000.g5_p010.d000.g5_p0-1-1.d000.g5_p011.d000.g5", "im": 0.0, "re": 0.0078125 }, // … ``` The levels are: 1. Total momentum (as string) 2. Irrep name 3. Irrep column 4. Irrep row 5. Label for the correlator matrix row, currently just a comma separated list of the relative momenta that have been used. In the future where we have other operators than just the pion, these will become more sophisticated labels. These are for human consumption anyway, so for the point of this code they are just strings. 6. Label for the correlator matrix column 7. List with summands that contribute. Each list element contains these four fields: conj ~ Whether the correlator is supposed to be complex conjugated before the scaling factor is applied. datasetname ~ Name of the HDF5 dataset that is to be taken from the contraction code. re, im ~ Real and imaginary part of the weight factor that is to be applied. ### Generating the structure First we need to iterate through all the momenta $\vec d$ that we want to use. For this we take the $\vec d_\text{ref}$ for each $\vec d^2$ and apply all rotations from the quotient group $O_h / \mathrm{LG}(\vec d_\text{ref})$. This way we get each unique $\vec d$ with the same magnitude. But in order to make this easier, we just apply the whole of the octahedral group and just use [`DeleteDuplicates`] afterwards. This works fine for $\vec d^2 \leq 3$. With $\vec d = (0, 0, 2)$ we have the same little group as for $\vec d = (0, 0, 1)$, but the $\vec d$ is different. This needs to be addressed. As we are interested in low momenta at this point, this corner is cut for now. We just use `Keys[IrrepDGammaAssoc[]]` to get all the momenta that we have little groups for. We just iterate through all $\vec d_\text{tot}$, then we iterate through all irreps and then iterate through all relative momenta. Using our `MakeGroupSum` function we construct the operators that will make up the correlator matrix. Later on we need to take the outer product in order to obtain the actual correlator matrix. From the given operators we need to extract the momenta and then perform the outer product of source and sink momenta. Then we need to apply the isospin templates in order to get the actual HDF5 dataset names. From there we massage the terms until we eventually get a JSON representation. This interface is described in the next chapter. --- - **`MultiGroupSum`**(irrep, $\{\{ \vec p_i \}\}$, hold : [`Identity`]) Calls `MakeGroupSum` for all the momentum sets given. Currently irrep row and column are fixed to $\alpha = \beta = 1$. Duplicates are automatically removed. Also this is limited to scalar particles with all $J_i = 0$ and $M_i = 0$. It can make sense to insert a [`Hold`] right before the `MakeGroupSum` such that you get a list of things that would be computed. As `MakeGroupSum` takes like 30 seconds per call, this function `MultiGroupSum` can take several hours to complete. If you pass [`Hold`] for *hold*, the actual spin projection will not be performed. Instead you will need to call `Map[expr, ReleaseHold, Infinity]` on the generated expression `expr` to evaluate all the `MakeGroupSum` calls. - **`GroupSumWholeIrrep`**($\vec d_\text{tot}$, irrep, $\{ \vec q_i \}$, cutoff, hold : [`Identity`]) Performs the `MultiGroupSum` for all individual momenta that give the given total momentum. - **`GroupSumWholeTotalMomentum`**($\vec d_\text{tot}$, $\{ \vec q_i \}$, cutoff, hold : [`Identity`]) Performs the `MultiGroupSum` for all irreps that are available with the given total momentum. Return value is an association from irrep to a list of spin projected operators. - **`PrescriptionToNumeric`**(prescription) Takes a prescription as an association with the fields - `"datasetname"` - `"re"` - `"im"` - `"conj"` and converts the symbolic expressions for real and imaginary part into floating point numbers such that it can be exported as JSON. - **`DatasetnameToObject`**(value, key) The key is a HDF5 dataset name string which might start with `conj:`. The value is a symbolic complex number. These are converted into an association with the fields - `"datasetname"` - `"re"` - `"im"` - `"conj"`. # Design and implementation (R) ## Projecting the computed correlators Taking the tables from the preceding step we must read in the prescribed HDF5 data sets and combine them given the factors. The interface that we get is described in [Interface to numerical code]. ### Iteration order The contraction code generates one HDF5 file per diagram type and per gauge configuration with names like `C4cC_cnfg4824.h5`. We have to decide how we iterate through all this data. In the end we want one file per correlator matrix element containing all the observations from all configurations. These iteration orders come to mind: 1. Treat the configurations independent of each other. Open all the diagram files for a particular configuration (like 4825) and then build all the correlator matrix elements from this. *Advantages*: Only one HDF5 file has to be opened at any one time. Also the file can be consumed completely in one go, amortizing the expensive indexing operation after opening the file. Parallelization over observations is trivial this way. *Disadvantages*: The intermediate result are lots of correlator matrix elements but with only one observation each. These have to be combined later on. Depending on the data format it is just a concatenation. One has to be careful not to create millions of files for intermediate output. 2. Focus on one correlator matrix element and build it for all observations at the same time. *Advantages*: There is no intermediate result, the whole statistics for a given correlator matrix element is created directly. *Disadvantages*: A lot of HDF5 files has to be opened, and only few datasets are going to be extracted. [We know that](https://github.com/HISKP-LQCD/hadron/issues/25) `h5ls` is about as expensive as `h5read`. Therefore this would make it rather costly to iterate this way. Also we would need a lot of RAM to store everything. As discussed with Markus, the first option seems to be the more sensible one. This decouples (numeric) projection and aggregation of the whole statistics. We just need to think about the intermediate data format. As we want to further process the data with R, there is no danger to use the `Rdata` format for serialization. For each configuration the whole correlator matrix will be in exactly the same outer nested list structure, just with a numeric vector as payload instead of the mapping between HDF5 dataset names and coefficients. Combining these is simple: Just load them one at a time and do `rbind` on the vectors. This will then be the data portion of an unsymmetrized `hadron::cf` correlation function. With this setup there will be no text IO for numeric data at all! --- Tasks: - Think about iteration order - Iterate through files - Open needed HDF5 files - Combine various numeric correlators # Using the programs A typical usage would be specifying the physical process, computing all the projection coefficients and then projecting the computed correlators. The directory and file structure works as follows. The git repository contains all the source files and is called *source directory* or `SRC`. The user must create a *work directory*. It is recommended that this is outside of the source directory such that none of the generated files appear to git. ## Projection coefficients The computation of the projection coefficients is parallel across the different total momenta $\vec d$ and the irreps $\Gamma$. This means that each of them can run in parallel at the cost of re-doing the Wick contractions. But as these are very cheap and the group theoretical sums are expensive this does not matter that much. Instead of using a Mathematica Notebook one uses a Wolfram Language script like the one found in `SRC/analytic_projection/driver`. This then calls the library functions that need. It takes command line arguments specifying the total momentum and the irrep. Using the `SRC/jobscript_generator/make-projection-jobs` script one generates an array of SLURM job scripts in `WORK/jobscripts` that each generate the projection coefficients for a single irrep. These jobs are then submitted to SLURM with `SRC` as the working directory. Using Mathematica 11.3 they take between an hour and around a day. With Mathematica 12 we face a severe performance regression and it is basically unusable. The projection coefficients will be a bunch of JSON files in `WORK/prescriptions`. ## Projecting numerical data Numeric projection needs the HDF5 files from the contraction code. These are to be placed into `WORK/correlators` and named like `C4cB_cnfg2552.h5`, just as the contraction code generates them. For the comparison to Markus Werner's data one needs the reference data for the same ensemble placed in `WORK/reference` in the form of TSV files like `rho_p2_B2_op8_gevp3.0.tsv`. For the actual numerical projection the code is prepared as a script in SRC/numeric_projection/number_crunching.R and takes the same command line arguments. The job script generator creates files for that into `WORK/jobscripts` as well, so you can submit these jobs also. The projected data will be placed into `WORK/projected` and comparison plots are created into `WORK/comparison`. # Tests The following are tests that we can perform to increase the confidence in $H_0$ which says that the code works perfectly fine. ## Comparison to Markus's data The comparison to Markus's data is performed in a separate Rmarkdown file: numeric_projection/Comparison_with_Markus_Werner.Rmd #### Momentum cutoff For $\vec d = (0, 0, 0)$ in the T1u we are missing something in the HDF5 files: ``` Object C4cB_uuuu_p0-10.d000.g5_p-1-1-1.d000.g5_p111.d000.g5_p010.d000.g5 does not exist in this HDF5 file. ``` Also for $\vec d = (0, 0, 0)$ in the A2u we have this: ``` Object C4cB_uuuu_p-1-1-1.d000.g5_p-1-1-1.d000.g5_p111.d000.g5_p111.d000.g5 does not exist in this HDF5 file. ``` This is because Markus uses $d^2 \le 2$ for the rest frame. This is a some additional rule that one has to keep in mind. #### Multiple irrep rows With some $d^2 = 1$ in the E we have this: ``` Error: nrow(filtered) == 1 is not TRUE ``` The issue here is that Markus has persisted multiple irrep rows. I need to compare to all the ones that I have or just fix it to some particular one. This also only happens in the E irreps, the T1u just contains one. #### Operators without coupling Sometimes there are no operator indices: ``` cannot open file 'B35.32/rho_p3_A2_operator-indices.tsv': No such file or directory ``` This likely is just because there is no coupling into that channel. My JSON files also show that there are matrix elements but none of them contain any actual summands. Markus has not created files in case that there is no coupling. ### Comparison results We currently (2019-08-08) see that the following: - For $d^2 = 0$ in the T1u we get agreement up to a factor of 2. We haven't resolved it yet. - Almost all other irreps (except the E) show perfect agreement. - Then for some reason all the E irreps in both $d^2 = 1$ and $d^2 = 3$ are off. They look as if they give somewhat the same effective masses, but they do not agree perfectly. The E are special because they contain complex irrep coefficients whereas the T1u contains only purely real or purely imaginary ones, all other irreps contain just real coefficients. ## Physical non-coupling There are a bunch of channels which should not couple with specific momenta and irreps, we can check for these. ## Cross irrep Correlation functions with operators from different irreps should just vanish. We can just pick a few examples and see how that works out. [`Apply`]: http://reference.wolfram.com/language/ref/Apply.html [`Association`]: http://reference.wolfram.com/language/ref/Association.html [`ClebschGordan`]: http://reference.wolfram.com/language/ref/ClebschGordan.html [`Coefficient`]: http://reference.wolfram.com/language/ref/Coefficient.html [`Composition`]: http://reference.wolfram.com/language/ref/Composition.html [`Dataset`]: http://reference.wolfram.com/language/ref/Dataset.html [`DeleteDuplicates`]: http://reference.wolfram.com/language/ref/DeleteDuplicates.html [`Dot`]: http://reference.wolfram.com/language/ref/Dot.html [`EulerMatrix`]: http://reference.wolfram.com/language/ref/EulerMatrix.html [`Function`]: http://reference.wolfram.com/language/ref/Function.html [`Hold`]: http://reference.wolfram.com/language/ref/Hold.html [`Identity`]: http://reference.wolfram.com/language/ref/Identity.html [`List`]: http://reference.wolfram.com/language/ref/List.html [`Map`]: http://reference.wolfram.com/language/ref/Map.html [`NonCommutativeMultiply`]: http://reference.wolfram.com/language/ref/NonCommutativeMultiply.html [`Part`]: http://reference.wolfram.com/language/ref/Part.html [`Prefix`]: http://reference.wolfram.com/language/ref/Prefix.html [`ReplaceAll`]: http://reference.wolfram.com/language/ref/ReplaceAll.html [`Rule`]: http://reference.wolfram.com/language/ref/Rule.html [`RuleDelayed`]: http://reference.wolfram.com/language/ref/RuleDelayed.html [`StringJoin`]: http://reference.wolfram.com/language/ref/StringJoin.html [`StringTemplate`]: http://reference.wolfram.com/language/ref/StringTemplate.html [`TemplateApply`]: http://reference.wolfram.com/language/ref/TemplateApply.html [`ToExpression`]: http://reference.wolfram.com/language/ref/ToExpression.html [`WignerD`]: http://reference.wolfram.com/language/ref/WignerD.html