Symmetry reduction of tensor networks in many-body theory I. Automated symbolic evaluation of $SU(2)$ algebra

The ongoing progress in (nuclear) many-body theory is accompanied by an ever-rising increase in complexity of the underlying formalisms used to solve the stationary Schr\"odinger equation. The associated working equations at play in state-of-the-art ab initio nuclear many-body methods can be analytically reduced with respect to angular-momentum, i.e. $SU(2)$, quantum numbers whenever they are effectively employed in a symmetry-restricted context. The corresponding procedure constitutes a tedious and error-prone but yet an integral part of the implementation of those many-body frameworks. Indeed, this symmetry reduction is a key step to advance modern simulations to higher accuracy since the use of symmetry-adapted tensors can decrease the computational complexity by orders of magnitude. While attempts have been made in the past to automate the (anti-) commutation rules linked to Fermionic and Bosonic algebras at play in the derivation of the working equations, there is no systematic account to achieve the same goal for their symmetry reduction. In this work, the first version of an automated tool performing graph-theory-based angular-momentum reduction is presented. Taking the symmetry-unrestricted expressions of a generic tensor network as an input, the code provides their angular-momentum-reduced form in an error-safe way in a matter of seconds. Several state-of-the-art many-body methods serve as examples to demonstrate the generality of the approach and to highlight the potential impact on the many-body community.


Introduction
In recent years, ab initio nuclear many-body theory has undergone a major renewal. In this process, expansion methods have become prominent in large-scale applications to mid-mass nuclei. The success obtained within the last two decades is leading to the design of more and more advanced approaches to continuously refine the accuracy of the calculations and extend them systematically to an even larger arXiv:2002.05011v1 [nucl-th] 12 Feb 2020 2 portion of the nuclear chart. This rise in the degree of sophistication of state-of-the-art many-body expansion schemes is leading to an increase of the formal complexity that is now at the edge of what is humanly processable.
When following the ab initio philosophy to solve the stationary Schrödinger equation, quasi-exact approaches based on Monte Carlo techniques [1][2][3][4] or configuration interaction (CI) [5,6] are limited by their computational scaling to the lightest systems. Moving to the realm of medium-and heavymass nuclei involves the use of expansion many-body techniques building a wave-function parametrization on top of a conveniently chosen reference state. These methods display a polynomial scaling with system size, the degree of the polynomial increasing with the targeted accuracy, i.e., with the order at which the expansion is truncated. This computational advantage typically comes at the price of being restricted to working in a non-variational scheme. Examples of such approaches are many-body perturbation theory (MBPT) [7][8][9][10][11][12][13][14][15][16], coupled-cluster (CC) theory [17][18][19][20][21][22][23], self-consistent Green's function (SCGF) theory [24][25][26][27][28] or the in-medium similarity renormalization group (IMSRG) method [29][30][31][32][33][34][35][36][37], all of which provide a consistent description of (at least) groundstate observables in nuclear many-body systems. In quantum chemistry in particular, MBPT and CC theories have a long tradition and both frameworks have been derived and implemented at very high truncation orders [8]. Although every member of the aforementioned approaches can be applied to much higher masses and larger system sizes than exact methods, the truncation levels needed for high-accuracy calculations require substantial effort in the derivation of the formalisms and for their numerical implementations.
While in earlier works the working equations were derived by hand, the rising computational power and the development of computer-aided algebraic manipulation tools have facilitated the derivation of more advanced truncation schemes in modern many-body approaches, many of which have undergone their pioneering studies in quantum chemistry [38][39][40][41][42]. A shining example is the tensor contraction engine that was developed in close collaboration with computer scientists and has been one of the most powerful tools to extend quantum-chemistry calculations to higher accuracy by generating working equations and source code for large-scale distributed implementations [43].
Even though large progress has been made in the development of software supporting the formal developments at play in quantum many-body research, only few are directly dedicated to the nuclear many-body problem [44]. While sharing many formal similarities with its electronic counterpart the nuclear many-body problem differs in two key points requiring a dedicated attention (1) At the mean-field level, single-nucleon states carry good total angular momentum j = l + s, i.e., one must employ the so-called j-coupling scheme to define appropriate one-nucleon states. Contrarily, electrons carry a welldefined projection s z of the intrinsic spin and are, thus, best described on the basis of the so-called ls-coupling.
The main consequence is that nucleons orbit in energy shells characterized by a greater degree of degeneracy, thus leading to the large dominance of open-shell groundstates, i.e., degenerate systems, over the nuclear chart. (2) The inclusion of three-body forces in a realistic nuclear Hamiltonian is mandatory to ensure a quantitatively correct description of nuclear observables, i.e. one is bound to use where T is the kinetic energy operator whereas V and W are two-and three-body potentials, respectively.
While expansion many-body methods are first formulated in terms of a generic single-particle basis, their actual implementations typically exploit symmetry properties of the basis functions and of the targeted many-body state, e.g., with respect to angular-momentum or parity quantum numbers. The adaptation of the generic formalism to a specific symmetry group defines a symmetry reduction of the many-body formalism. The goal is to use reduced many-body tensors associated to irreducible representations (IRREPs) of the symmetry group to pre-process a subset of the summations at play in the tensor networks defining the working equations. A simple, yet representative, example is the pre-processing of spin summations in spin-restricted quantum chemistry calculations. The counterpart in nuclear structure theory relates to the exploitation of rotational invariance associated with the conservation of total angular momentum and encoded in terms of the S U(2) nonabelian Lie group. In this particular case the reduction scheme will be referred as the angular-momentum reduction (AMR).
Eventually, it turns out that the AMR poses a nontrivial problem requiring the same amount of effort that the derivation of the initial working equations. However, there exists a highly systematic and elegant way to deal with this task that is close in spirit to the use of Feynman's diagrams as a mnemonic device to represent physical processes.
Consequently, it is highly desirable to parallel the efforts done to automatize the generation of working equations by devising a framework that automatically performs the tedious symmetry reduction in an error-safe way. Currently, there is -to the best of our knowledge -no open-source library that can deal with the requirements imposed by nuclear structure many-body methods to perform symbolic manipulations of angular-momentum algebra. Typically, existing software is restricted to the numerical evaluation of coupling coefficients instead of performing symbolic manipulations including the simplification of complex tensor networks. There have been similar attempts for symbolic simplifications of 3 angular-momentum expressions before without formally connecting it to many-body theory [45].
Therefore, the goal of the source code accompanying the present document is to support the implementation of advanced many-body frameworks in nuclear structure in an error-free way. Of course, this does not resolve the problem of writing an efficient and error-free numerical implementation of the symmetry-reduced formalism itself. While the generation of the source code is envisioned, it is, however, beyond the scope of the present work.
The document is organized as follows. In Sec. 2 the notion of symmetry in the context of many-body theory is introduced using a group theory formulation. Section 3 focuses on the angular-momentum algebra and its relation to states and operators. In Sec. 4 the diagrammatic allowing for the handling of the S U(2) algebra is laid out and the simplification rules for the graph theory reformulation of tensor networks are presented. Section 5 discusses several state-of-the-art many-body approaches that serve as pedagogical examples to demonstrate the generality of the approach. Ultimately, an outlook is provided in Sec. 6.
2 Symmetries and many-body theory

Symmetry group
Physical symmetries impact many-body formalisms at various stages of their elaboration. The existence of symmetries in finite systems is intimately connected to conservation laws, e.g., the existence of U(1) global gauge symmetry corresponds to particle-number conservation while S U(2) symmetry corresponds to angular momentum conservation. Mathematically, the invariance of a quantum system, characterized by its Hamiltonian H, is encoded in terms of transformation properties imposed by a symmetry group G Ham whose action leaves the physical system invariant or, equivalently, the existence of a unitary linear representation U acting on the space of states such that which can be rewritten as Given the eigenstates of the Hamiltonian Eqs. (2)(3) stipulate that the transformed states are also eigenstates with the same eigenvalues.
In the case of discrete symmetries such as parity or time reversal the corresponding symmetry group is finite, e.g. Z 2 .
Contrarily, continuous symmetries correspond to Lie groups allowing for a continuous parametrization of the (infinite number of) group elements in terms of a finite set of parameters. The present focus is on the nonabelian S U(2) Lie group associated with rotational invariance of nuclear systems. Relevant details about this symmetry group are provided in Sec. 3.3.
Eventually, symmetries enter the formulation of (nuclear) quantum many-body methods at three different levels (1) the symmetry group of the Hamiltonian G Ham specifying the invariance of the physical system under a given set of transformations along with the symmetry quantum numbers carried by its many-body eigenstates, (2) the symmetry group of the single-particle basis G bas specifying the symmetry properties of the computational basis, (3) the symmetry group of the reference state G ref employed in expansion methods specifying the symmetries of the auxiliary many-body problem that is solved to construct the reference state.
While the symmetry group of the Hamiltonian is fixed by the physical system under consideration, the symmetry properties of the single-particle basis and the reference state result from a choice such that various combinations of G bas and G ref can be employed.

Symmetries of the single-particle basis
Given H and its symmetry group G Ham , there is infinitely many different single-particle bases spanning the one-body Hilbert space H 1 that can be used to represent the operator in second-quantized form. The single-particle basis functions are typically obtained as eigenstates of an auxiliary one-body Hamiltonian H bas whose symmetries are characterized by [H bas , U(g)] = 0 (∀g ∈ G bas ).
When choosing G bas = S U(2), one-body basis states are eigenstates of the squared total angular-momentum operator 1 where J x , J y and J z denote the Cartesian components of the total angular-momentum vector. In most ab initio nuclear structure applications such a one-body basis is indeed employed, e.g., the eigenbasis of the three-dimensional spherical harmonic oscillator (sHO) Hamiltonian such that the one-body eigenstates of H sHO are proportional to spherical Harmonics. In other frameworks, e.g. nuclear energy density functional calculations, the single-particle basis is possibly taken as eigenfunctions of the axially deformed HO Hamiltonian that breaks rotational invariance and, thus, displays a smaller symmetry group G bas than H sHO .

Symmetries of the reference state
The rationale of expansion methods relies on the definition of a conveniently chosen A-body reference state |Φ that serves as starting point for the correlation expansion. Acting on the vacuum, the wave operator W yields the exact, e.g., ground state The wave operator is expanded and truncated according to a given many-body scheme, e.g., in MBPT, SCGF or CC theory. The resulting equations are symmetry-unrestricted and therefore make no use of symmetry properties of manybody operators. In practice, the reference state is typically obtained as the ground state of an 'unperturbed' Hamiltonian H ref capturing the average behavior of the system's dynamics and characterized by a symmetry group G ref such that |Φ typically belongs to the trivial IRREP of G ref .
In the following, the reference state |Φ G ref , thus, carries a subscript specifying the symmetry group of the Hamiltonian it is the ground state of. In the simplest case, the vacuum is chosen to be a Slater determinant |Φ G Ham obtained from a symmetry-restricted Hartree-Fock mean-field calculation, i.e.
In nuclear systems, |Φ G Ham typically belongs to the trivial IRREP of S U(2) and U(1), i.e., it carries good angular momentum J = 0 and a fixed number of particles. Dynamic correlations are introduced via the action of the wave operator that generates summations over elementary particle-hole excitations.
In open-shell systems, the above reference state is improper due to the partial filling of the last occupied shell. This leads to a degeneracy with respect to particle-hole excitations, thus, signalling the existence of a Goldstone mode More importantly, this lowering is accompanied by a lifting of the degeneracy of |Φ G Ham with respect to elementary excitations such that W can be expanded safely. In this case, however, the wave operator must not only capture dynamical correlations but also restore the symmetry G Ham associated with the exact eigenstates of H. Because of the necessary truncation, a standard expansion of W is not capable of restoring the symmetry such that the symmetry contamination needs to be retrieved by the explicit inclusion of a symmetry projector in the definition of W [21-23].

Reduction schemes and groups
While the symmetry group of the Hamiltonian is fixed from the outset, the choice of the single-particle basis and reference states leaves tremendous freedom to adapt G bas and G ref in order to deal with a specific situation. A case of particular interest arises when the symmetry groups of the Hamiltonian, the single-particle basis and the reference state coincide, i.e., In this setting, the common algebraic structure can be exploited to simplify the many-body formalism by expressing all working equations in terms of G sym -reduced tensors, thus, potentially providing a tremendous gain in the required runtime and memory resources. In the present paper, this situation is exploited relative to the S U(2) group (independently of the treatment of other symmetries such as U(1)).

Tensors and tensor networks
Due to the large variety of expansion schemes built to retrieve the solution of the many-body problem, it is desirable to introduce a unifying language for the various frameworks. This common ground is provided by the language of tensors and tensor networks. A mode-k symmetry-unrestricted tensor (SU-T) is a multi-variate data array carrying k indices with (possibly different) index ranges I 1 , ..., I k . Tensors constitute the basic 5 building blocks of many-body expansion methods. Given a set of SU-T's A, B, C, ... a contraction is defined as a summation over a common index, e.g., where the ellipses indicate indices that are not summed over 2 .
A symmetry-unrestricted tensor network (SU-TN) denotes a set of SU-T's combined according to a given contraction scheme specifying the way the tensors are contracted with each other. Furthermore, a SU-TN is said to be closed if all tensor indices are summed over and is said to be open otherwise.
In many-body applications tensors typically appear in two broad classes (1) input tensors that are known prior to addressing the actual solution of the Schrödinger equation in a given manybody framework, (2) output tensors that are specific to a given many-body approach and are typically the objects being solved for.
Examples for input tensors are matrix elements of many-body operators like the Hamiltonian whereas examples of output tensors are CC amplitudes or dressed propagators in SCGF theory. In most non-perturbative many-body frameworks, like CC, IMSRG or SCGF, open TN's specify the working equations required to determine the unknown output tensors while the calculation of observables, e.g. the energy, relates to the evaluation of closed TNs.

Symmetry-reduced tensor networks
The goal of this work is to transform an initial SU-TN into a symmetry-reduced tensor network (SR-TN) encapsulating the symmetry reduction according to the associated symmetry group. To do so, the SU-T's must be replaced by their symmetry-reduced counterparts. Given an initial SU-T, the corresponding symmetry-reduced tensor (SR-T) is obtained from a transformation f G sym mediating the symmetry reduction related to the group G sym . Here, the symbol λ denotes the relevant IRREP labels of the symmetry-reduced tensor. In the following, quantities with a tilde indicate symmetry-reduced objects. Note that the content of the indices themselves change, such that the set of quantum numbers labelling a SU-T and its SR-T counterpart are different. Thus, the SR-TN denotes the end product obtained via the replacement of the SU-T's by their SR-T 3 Angular-momentum algebra

Rationale
While the discussion on symmetry-reduction and SR-TN's has been generic so far, the present paper focus on the S U(2) group. The goal is, thus, to obtain angular-momentum-reduced tensor networks (AMR-TN's) from SU-TN ones. The procedure requires to (1) replace all the SU-T's by their AMR-T counterparts according to the transformation f S U(2) , (2) constrain the contraction pattern to only be left with summations over the reduced set of quantum numbers.
In practice, step (1) involves a set of substitution rules for every many-body tensor at play that specify how the symmetry reduction is performed. The resulting SR-TN-and its computational complexity-may strongly depend on the choice made to perform this initial step 3 . From this point of view at least a minimum level of human input (and experience) is necessary to come up with the most convenient choice. This does not pose a severe limitation in any of the examples discussed below.

Other symmetries
While presently focusing on rotational symmetry, other symmetries can be exploited in the same way. A key example relates to intrinsic spin in quantum chemistry that is analogous to the total angular-momentum when using a ls-coupling scheme. The spin projection being only two-fold degenerate, i.e. m s = ± 1 2 , spin-restricted many-body theories benefit less from the symmetry reduction than in the j-coupling scheme. Still, pre-processing the sums over spin projections is an important tool to reduce the computational cost and advance state-of-the-art expansion methods in strongly correlated electronic systems. Finite symmetry groups, e.g., the dihedral groups D n , may also arise in quantum molecules whereas cubic groups play an important role in the computation of homogeneous matter, e.g., the infinite electron gas or infinite nuclear matter, since periodic boundary conditions are employed to facilitate the calculation. In solid-state physics, symmetry properties of the many-body systems, e.g., helical 6 symmetries in nano tubes, can also be exploited to reduce computational complexity.
All the aforementioned examples correspond to a reduction of exact symmetries of a many-body system. In recent years, exploiting emergent approximate symmetries has also been shown to be highly beneficial, in particular in the context of nuclear CI-based approaches. In this case, the symmetry group of the configuration basis G bas is larger than the actual symmetry group of the Hamiltonian, thus, exploiting algebraic properties that are not strictly realized in nature. A prime example is the symplectic symmetry group S p(3, R) that is not an exact symmetry of the nuclear Hamiltonian but of the kinetic energy operator. In the symmetry-adapted no-core shell model (SA-NCSM) an Abody configuration basis is constructed from the Casimir operators of the approximate symmetry group S U(3) ⊂ S p(3, R). The use of symplectic algebra was shown to provide an efficient selection of many-body basis states, thus, yielding computational savings in the diagonalization of the many-body Hamiltonian at the price of a more involved handling of many-body operators [46].

S U(2) group
In order to move closer to a concrete implementation of the above procedure, let us introduce details about the nonabelian compact S U(2) ≡ {R(Ω), Ω ∈ D S U(2) } Lie group associated with the rotation of a A-body fermion system characterized by an integer or a half-integer angular momentum. The group is parametrized by three Euler angles Ω ≡ (α, β, γ) whose domain of definition is As S U(2) is considered to be a symmetry group of H, the commutation relations hold for Ω ∈ D S U(2) . Subsequently, the unitary representation of S U(2) on Fock space is utilized The components of the total angular-momentum vector make up the Lie algebra where i jk denotes the Levi-Civita tensor. The Casimir operator of the group built from the infinitesimal generators through a non-degenerate invariant bilinear form is the total angular momentum Matrix elements of the irreducible representations (IR-REPs) of S U(2) are given by the so-called Wigner D functions [47] where |ξJM is an eigenstate of J 2 and J z The index ξ collects all quantum numbers but J and M. The (2J+1)-dimensional IRREPs are labelled by J and are spanned by the set of states {|ξJM } for fixed J and ξ. An irreducible tensor operator T J of rank J is made of 2J + 1 operators T J K transforming under rotation as or, equivalently, fulfilling where J ± = J x ± iJ y denotes the usual raising (lowering) operators. The nuclear Hamiltonian is an example of a spherical tensor operator of rank zero. Such operators are denoted as scalar.
A powerful tool to treat spherical tensor operators is the celebrated Wigner-Eckart theorem (WET) where ≡ 2 j + 1. The theorem states that matrix elements of a given component of a spherical tensor in the basis spanning the IRREPs can be written as a product of a geometric part independent of the spherical tensor at play and of a reduced matrix element independent on the particular component of the spherical tensor and the members of the IRREPs under consideration [48].
In the special case of a scalar operator one has such that both the initial and reduced matrix elements are independent of any projection quantum number. In this particular case, the notion of reduced matrix element is thus irrelevant. 7

Fermionic algebra
One of the building blocks of quantum many-body theory are second-quantized operators where c † (c) denote single-particle creation (annihilation) operators associated with a basis B 1 of the one-body Hilbert space that is the tensor product of a spatial part, a spin part and an isospin part. Anti-symmetrized matrix elementsō k 1 ...k i+ j carrying (i + j) one-body indices constitute a mode-(i + j) SU-T. Creation and annihilation operators are assumed to fulfil the canonical anti-commutation rules defining the Fermionic algebra 4 . Processing many-body matrix elements of strings of such operators via various forms of Wick's theorem is at the core of quantum many-body methods and gives rise to the multitude of TN's at play in expansion methods.

S U(2) symmetry and basis states
In the following, B 1 is taken to be the eigenbasis of a S U(2)invariant Hamiltonian H bas such that basis states are conveniently labeled as where n k denotes the radial quantum number, l k the orbital quantum number, j k the total angular-momentum quantum number, m j k its projection and t k the isospin projection. This constitutes a so-called j-coupled basis, i.e. it is not a direct product of bases of H r 1 and H s 1 but a coupled basis whose members are eigenstates of the total angular momentum j 2 . While the eigenstates of the aforementioned sHO Hamiltonian provide a example of practical interest, eigenbases 4 When breaking U(1) symmetry, one employs the quasi-particle algebra associated with the Bogoliubov transformation [49] such that operators are expressed in this basis and that the indices of the associated matrix elements relate to quasi-particles. This feature does not change fundamentally what follows regarding the handling of S U(2) symmetry.
of other one-body Hamiltonians characterized by rotational symmetry are equally valid. Later on, the AMR-T's employed throughout the symmetry reduction will carry reduced labelsk characterized bỹ where the angular-momentum projection, i.e. the magnetic quantum number m k , is explicitly excluded compared to the definition of k through the given of the basis in Eq. (33). The tensor product of two one-body states defines a basis state of the two-body Hilbert space H 2 Contrarily, in the coupled representation the two total angular momenta j k 1 and j k 2 are coupled to a total two-body angular momentum J with projection 5 M, where the vector space inner product denotes the Clebsch-Gordan (CG) coefficient mediating the transformation from the uncoupled to the coupled basis. The left-hand side of Eq. (36) defines two-body eigenstates of J 2 , the Casimir operator of the group. The inverse transformation of Eq. (36) is given by Along the same lines, the uncoupled three-body basis states of H 3 , i.e., the tensor product of three single-particle states is defined. Performing the angular-momentum coupling requires fixing the coupling order which is subsequently chosen to be i.e., the first two single-particle states are coupled to an intermediate two-body angular-momentum quantum number J 12 , which is further coupled to the third state to yield the overall (half-integer) three-body angular-momentum J. In analogy to the two-particle case, Eq. (40) defines three-body eigenstates of J 2 .
The choice of the coupling order for three-body states employed in Eq. (40) is arbitrary such that an alternative coupling scheme is given by where the second and third single-particle states are coupled to an intermediate angular momentum J 23 that is subsequently coupled with j k 1 to an overall J. Both coupling schemes enable for the construction of a basis of H 3 that is an eigenbasis of J 2 . The transformation between the two representations is given by where the Wigner 6 j-symbol was introduced.
Recursively, N-body states can be introduced for N ≥ 3, e.g. for the uncoupled representation Since in current ab initio implementations four-and higherbody operators play no dominant role yet, this extension is not discussed here.

Many-body matrix elements
With the operator O being the µ component of a spherical tensor of rank λ, its uncoupled matrix elements are defined bȳ where it is not assumed that the numbers of indices labelling the bra and the ket states coincide. By means of the transformations between uncoupled and coupled representation of the bra and ket states, coupled expressions for matrix elements can be derived. Focusing on a two-body operator characterized by uncoupled matrix elementsō k 1 k 2 k 3 k 4 , their angular-momentum-coupled counterparts are 6 6 The transformation can in principle be accompanied by an additional phase factor.
Analogously, coupled three-body matrix elements are obtained as Neither Eq. (45)

Diagrammatic method
Even though all manipulations necessary to simplify angularmomentum expressions can be performed solely in terms of the expressions introduced in Sec. 3 it is at the heart of this work to introduce a more convenient representation of the involved algebraic steps that, additionally, allows for computer-aided derivations. As Feynman or Goldstone diagrams are used to efficiently capture the results of cumbersome applications of Wick's theorem, diagrams can be introduced to restate complicated identities associated with angular momentum algebra [47]. A modern account of the underlying group-theoretic properties is provided in Ref. [45]. For completeness, let us mention that similar frameworks can be introduced to tackle other (more involved) symmetry groups 7 . The interested reader is referred to Ref. [50] for an extensive discussion.

Preliminaries
As seen in Sec. 3, CG coefficients constitute the basic building blocks of angular-momentum theory. However, CG coefficients are somewhat inconvenient due to their asymmetry with respect to the involved angular-momenta. A more symmetric representation can be obtained in terms of Wigner 3 jm-symbols 8 9 Wigner 3 jm-symbols are invariant under cyclic column permutations, whereas anti-cyclic permutations induce a phase factor Wigner 3 jm-symbols with opposite magnetic quantum numbers are related via the identity Furthermore, 3 jm-symbols with one vanishing ( j, m) pair simplify according to

Vertices
Wigner 3 j-symbols provide the building blocks of the diagrammatic formalism. They are represented by vertices in the so-called Yutsis graphs. 9 More specifically, a vertex carrying three outgoing lines, each labelled by a tuple ( j k , m k ), represents the 3 j-symbol The vertex sign denotes a convention specifying the column order that must be used to write the corresponding 3 jmsymbol, i.e., a plus (minus) sign stipulates that the lines and the associated angular-momentum labels must be read counterclockwise (clockwise). Furthermore, the vertex with one ingoing line represents Starting from the two above definitions, the vertices with two and three ingoing lines are obtained by applying the operation consisting of inverting the directions of all three lines at once. Starting for example from the vertex with three outgoing lines, one obtains the vertex with three ingoing lines 9 Named after Lithuanian Adolfas Jucys, these are also known as Jucys graphs. Since many of his works were initially published in Russian, the transliteration of the name from Russian, Yutsis graph, is the betterknown one. whose expression is given by as a testimony of Eq. (50) and where the magnetic quantum numbers have been added to the phase at no cost given that m 1 + m 2 + m 3 = 0 holds. Through this operation, the sign is not altered. Additionally, performing the operation twice does give back the original vertex thanks to the identity Changing the sign carried by the vertex can be performed at the price of the phase factor where the lower index 'ns' stipulates the node sign reversal. Indeed, moving from a clockwise to a counterclockwise (or vice versa) reading of the vertex corresponds to performing one column inversion in the 3 j-symbol whose effect is characterized by Eq. (49). Notice that changing the vertex sign is equivalent to moving one line across another one.

Yutsis graphs
The network of 3 jm-symbols generated via step (1) of the angular-momentum reduction of a SU-TN (see Sec. 3.1) is represented by a Yutsis graph. Those graphs are, thus, obtained by contracting a set of vertices through their edges in a way that consistently represent the network of 3 j-symbols.
Contracting the edges of two vertices is possible if both lines carry the same angular momentum quantum numbers ( j, m) and go in the same direction, i.e., one must be going out of the first vertex while the other one must be going into the second vertex Reading the vertices according to the definitions given previously, the algebraic expression resulting from the contraction reads as Given a Yutsis graph, the direction of an internal line carrying angular momentum j can be reversed at the price of accounting for the phase factor An example of practical interest relates to fully contracting the two vertices actually corresponding to the so-called Wigner 3 j-symbol j 1 j 2 j 3 , also called triangular delta 10 or triangular inequality. The corresponding algebraic expression is given by which vanishes unless the inequalities are satisfied.

Unfactorizable graphs
Wigner 3n j-symbols provide relevant examples of Yutsis graphs that cannot be simplified via factorization rules. The first example is the Wigner 6 j-symbol that is graphically represented as a tetrahedral structure 10 Some texts call the triangular delta a 3 j-symbol, as it technically is the first of the 3n j-symbols, while naming the coupling coefficients 3 jm-symbols. This nomenclature in is not adopted in order to avoid confusion. Independently of which of the two diagrams is used, the corresponding algebraic expression is The case n = 3 yields the Wigner 9 j-symbol whose algebraic expression can be represented by the Yutsis graph given by the following hexagon involving six vertices and nine lines by inverting the signs of the magnetic quantum numbers in the last three 3 jm-symbols. While higher-order 3n j-symbols only rarely arise in nuclear many-body theory, they can be equally represented by an unfactorizable Yutsis graph. They do in fact naturally enter in the partial-wave decomposition of nuclear k-body Hamiltonians for k ≥ 4. In practice, Wigner 3n j-symbols play an important role given that they can be pre-calculated and stored in cache in large-scale applications. This is typically done for 6 jsymbols and if necessary for (a subset of) 9 j-symbols. Since the number of 9 j-symbols is very large for a selected model space it is often useful to re-express 9 j-symbols as sums of products of 6 j-symbols and resort to much smaller 6 j-caches if the structure of the angular-momentum networks supports such a strategy.

From tensor networks to Yutsis graphs
The crucial first step consists of extracting the Yutsis graph associated with the SU-TN of interest. Following step (1) in Sec. 3.1, this is achieved by expressing the original SU-T's in terms of SR-T's and a set of CG coefficients that are consecutively replaced by their 3 jm-symbol equivalents. The next step consists of splitting each involved summation according to k → n k l k j k t k m k → k m k . (60) In doing so, one can isolate the networks of 3 jm-symbols along with the sums over the magnetic quantum numbers. This corresponds to extracting the associated Yutsis graph.

Factorization rules
Having the Yutsis graph at hand, the goal is to simplify it a much as possible. This corresponds to identifying specific subparts in the graph that can be reduced via the application of identities satisfied by appropriate (sub)sets of 3 jmsymbols. Once this is completed, one is left with an expression involving irreducible Wigner 3n j-symbols (see Sec. 4.4) and no magnetic quantum number dependence anymore. The benefit of using Yutsis graphs is that the search for reducible parts can be automated while their actual reduction can be realized by applying systematic factorization rules on the graph. The rules are characterized by the length of the cycles involved in the factorization process. Below, the factorization rules are introduced one after another with increasing degree of complexity, i.e., cycle length. For the proofs of the factorization formula, the reader is referred to Ref. [45]. A more extensive list of angular-momentum-algebra identities that can be used to define factorization rules can be found in Ref. [47].

Cycles of length two
The next simplest factorization corresponds to the reduction of a 2-cycle. Algebraically, the corresponding identity is the orthogonality relation Figure 1 provides the diagrammatic representation of the identity stated in Eq. (61). Thus, the 2-cycle rule replaces two vertices connected by two lines by a single line in a Yutsis graph.

Cycles of length three
The simplest factorization rule leading to a non-trivial Wigner 3n j-symbol corresponds to the factorization of a 3-cycle as displayed in Fig. 2. Algebraically, the factorization corresponds to the identity YG 3c ≡ Equation (62) allows one to factorize a topology involving three vertices into an irreducible part, i.e. the 6 j-symbol, and a single vertex. Therefore, the resulting graph contains two vertices and three lines less than the initial one.

Cycles of length four
The most involved factorization rule employed in this work corresponds to a cycle of length four as displayed in Fig. 3. The underlying algebraic identity is given by Equation (63) allows to factorize the topology involving four vertices into an irreducible part made of two 6 j-symbols and two vertices. Therefore, the resulting graph contains two vertices and three lines less than the initial one. Note that the Yutsis graph in Fig. 3 has the symmetry of a square: rotations by multiples of 90 • leave it invariant. The rotation is equivalent to relabeling the edges, and leads to a different equivalent factorization for rotation angles of 90 • and 270 • . The handling of cycles of length four is thus nonunique.

Cycles beyond length four
While the present code supports factorizations involving up to cycles of length four, there exist topologies, sketched in Fig. 4, which cannot be simplified through the above stated rules but require more involved identities. In principle, this restricts the range of applicability to topologies that do not contain cycles of length five or higher. The smallest cubic graph involving a cycle of length five is the so-called Peterson graph containing ten vertices and 15 edges. Consequently, the simplest many-body diagram potentially leading to this topology must contain at least five two-body vertices, e.g., corresponding to a fifth-order MBPT diagram or a CC diagram with T 3 amplitudes.
In the testing phase of the current version of the code, the factorization rules were applied to hundreds of many-body diagrams including topologies that are far beyond current state-of-the-art applications. In none of these test cases a Yutsis graph involving a cycle of length five or higher appeared. Future versions of the program will be extended along these lines by including factorizations of more complex topologies, or including more elaborate techniques such as the interchange rule [51].

Applications
A number of different many-body formalisms are now used to exemplify the steps at play in the symmetry reduction process of SU-TNs. The emphasis is on the angular-momentum reduction and details of the formalisms themselves are not within the scope of the present work. All considered examples are typical of state-of-the-art nuclear structure applications.
The formulae derived below are not in the computationally most optimized form. For instance, below parityand isospin-conservation can be further exploited to yield more efficient implementations. However, processing S U(2)symmetry yields by far the highest computational benefit due to larger dimensionality of the associated IRREPs.

Many-body perturbation theory
In many-body perturbation theory (MBPT) an infinite power series is taken as an ansatz for the exact ground-state energy and wave function [8] where the lower index k enumerates excited states in the spectrum and |Φ ≡ |Ψ   interet. Since the following is exclusively concerned with the description of nuclear ground states, i.e., k = 0, the subscript is dropped for simplicity. The starting point is given by the definition of a splitting of the full Hamiltonian into an unperturbed part H 0 and a perturbation H 1 such that the reference energy is given by The first contribution to the correlation energy is, thus, obtained at second order. The simplest choice is to take |Φ as a Slater determinant, typically obtained as the solution of a S U(2)-restricted Hartree-Fock (HF) calculation.
In recent years, more sophisticated vacua have been used in order to account for so-called static correlations in open-shell systems. Both multi-configurational reference states obtained from a configuration interaction (CI) diagonalization in a small model space [13] and particle-number-broken HFB vacua [12] have shown to provide computationally cheap benchmarks without loss in accuracy when employing softened chiral potentials. For a recent review, see Ref. [15].
Presently, the simplest single-reference case of low-order canonical HF-MBPT is discussed 12 . Examples are worked out in detail to enable a deeper understanding of each of the individual algorithmic steps.

Second-order energy correction
The second-order energy correction reads as 13 where a, b and i, j denote particle and hole states, respectively, i.e., states that are unoccupied and occupied in the reference Slater determinant, respectively. Additionally, a short-hand notation for the energy denominator is used where k denotes HF single-particle energies. According to the previous definitions Eq. (68) provides a closed SU-TN involving two mode-4 tensors, i.e. H i jab and ab i j . Expressing the two involved tensors in Eq. (68) in terms of their AMR-T counterparts according to (the inverse of) Eq. (45) yields where the fact was used that ab i j = ãb ij is already an AMR-T given that the single-particle energies are m-independent, i.e., 14 k = k . The tensor network in Eq. (70) is split into an S U(2)invariant part that does not depend on single-particle angularmomentum projection quantum numbers and a part carrying the full dependence of magnetic quantum numbers that will be subsequently simplified. In a first step, CG coefficients are converted into 3 jm-symbols yielding where each phase factor gives in fact Focusing on the 3 jm-symbols network in Eq. (71), the second step consists of reversing all m quantum numbers in the second and fourth 3 jm-symbols at the price of an extra phase factor, where the magnetic quantum numbers have been added to the phase at no cost given that m i + m j − M 1 = 0 and m a + m b − M 2 = 0 hold and that M 1 and M 2 are integers. The expression in Eq. (73) is now in the proper form to allow for its identification with an appropriate Yutsis graph Now that the working graph has been built, the next step consists in simplifying it via the application of appropriate factorization rules. The application of the 2-cycle rule, see Fig. 1, requires the direction of the edges carrying j a and j b to be reversed, thus bringing the phase Φ lr = (−1) 2 j a (−1) 2 j b = (−1) 2 = 1 and yielding the diagram where the red box indicates the subpart of the diagram that is factorizated in the next step. Factorizing the 2-cycle provides the intermediate factor and leaves the diagram In the last step, the 3 j-symbol is identified after reversing the orientation of the edges carrying j i and j j leading to the additional phase Φ lr = (−1) 2 j i +2 j j = (−1) 2 = 1 and providing the overall result Replacing the m-dependent part of Eq. (71) by Eq. (75) finally provides the AMR form of the second-order energy correction where triangular inequalities coming from 3 j-symbols are assumed. While the initial SU-TN is of N 4 complexity, the AMR-TN is ofÑ 4 · (J max + 1) complexity, whereÑ is the number of reduced basis statesk and J max corresponds to the maximum number of channels (i.e. allowed values) of the two-body angular-momentum given the maximum one-body angular momentum retained in the (truncated) basis B 1 . For large model spaces the difference in runtime is improved by several orders of magnitude even for this very simple example.

Third-order energy correction
A more elaborate example is given by the third-order energy correction to the ground-state binding energy Following the same procedure as for E (2) 0 , one obtains the Yutsis graph associated with the particle-particle contribution (i.e. E (3) pp ) is given by and, similarly, the one extracted from the hole-hole contribution (i.e. E (3) hh ) which are topologically identical. Due to the presence of one more Hamiltonian matrix element compared to the secondorder energy correction, the number of 3 jm-symbols, i.e. the number of nodes, is increased by two. In both cases, the red boxes indicate the subgraphs that are factorized by the application of the 2-cycle rule. Applying it twice and identifying the resulting Yutsis graph as a 3 j-symbol leads to the result Considering theĴ 2 1Ĵ 2 2Ĵ 2 3 factor coming from the prior conversion of CG coefficients into 3 jm-symbols, the final AMR form of the two contributions is which can be read as simple matrix-matrix products within each J channel. The symmetry reduction of the particle-hole term (i.e. E (3) ph ) is more involved such that, following the same steps, the associated Yutsis graph is and can be re-arranged in a more convenient way as which is nothing but a 9 j-symbol. Consequently, the AMC form of the particle-hole term leads to the algebraic expression In practical applications, Eq. (81) is conveniently re-written by expressing the 9 j-symbols as a sum of products of three 6 j-symbols which can be obtained graphically by a successive application of the 4-cycle rule, the 3-cycle rule and finally the identification of a redundant 3 j-symbol. Based on the introduction of so-called Pandya-transformed matrix elements [48] 16 Eq. (82) can eventually be written as which reads as the trace of a two-fold matrix-matrix product of Panya-transformed Hamiltonian matrix elements. Equation (84) clearly shows the computational benefit of an appropriate choice of the coupling order which in practice is not at all obvious.

Coupled-cluster theory
Contrary to a simple power series ansatz, coupled-cluster theory aims at a non-perturbative resummation of large classes of MBPT diagrams.

General formalism
The starting point in the CC framework is an exponential ansatz to parameterize the exact ground state [8], in terms of the connected cluster operator T defined as The second-quantized form of the individual terms in Eq. 86 is given by T n ≡ 1 (n!) 2 a 1 ...a n i 1 ...i n t a 1 ...a n i 1 ...i n c † a 1 · · · c † a n c i n · · · c i 1 , with t a 1 ...a n i 1 ...i n the n-tuple cluster amplitudes characterizing a mode-2n tensor. Thanks to the exponential form for the wave operator, the CC approach is manifestly size-extensize. In actual applications, T is truncated at a fixed truncation level defining a particular CC model, e.g., where the lower index c stipulates the connected character of the expansion.

Energy equation
In the absence of three-body operators in the input Hamiltonian, the correlation energy is given for arbitrary CC truncations by Equation (90) defines a closed TN involving at most four internal contractions. Note that higher-order amplitudes affect the energy only implicitly by relaxing T 1 and T 2 without entering the energy equation explicitly. Contrary to canonical MBPT, the CC energy equation involves mode-2 tensors associated with one-body operators, i.e. the T 1 amplitudes and the matrix elements f pq of the Fock operator. The Fock operator is S U(2)-invariant as long as the mean-field calculation is performed in a symmetry-restricted way. As the reference Slater determinant is presently characterized by J = 0, cluster amplitudes are irreducible S U(2) tensors of rank J = 0 such that a similarity-transformed operatorŌ has the same irreducible S U(2) tensor rank as its non-transformed counterpart O. Hence, Wigner-Eckart's theorem trivially enables the introduction of reduced matrix elements Inserting such forms into the first contribution to the CC the energy yields from which the m-dependent part can be extracted to yield the Yutsis graph Reversing the direction of the j a edge (Φ lr = (−1) 2 j a = −1) together with changing the sign of the leftmost node (Φ ns = 17 (−1) j a + j j ) allows to make use of the 2-cycle rule which leads to −(−1) j a + j j j a j j 0 = −(−1) j a + j j δ j a j j = δ j a j j , such that one obtains the final AMR form ai f ia t ai = ξ a ξ i j a (ξ i j a |F|ξ a j a )(ξ a j a |T 1 |ξ i j a ) .
For the second term of the energy equation, one has The detailed derivation of the last contribution to the energy equation is omitted given that it is formally identical to the derivation of the second-order MBPT correction, i.e. the appropriate Yutsis graph is the one displayed in Sec. 5

Amplitude equations
The unknown cluster amplitudes are obtained by solving a set of CC amplitude equations Equations (99) constitute a set of coupled non-linear equations that must be solved iteratively for every external index combination. They also provide typical examples of open SU-TNs containing external indices that are not summed over 14 .
In order to perform the symmetry reduction, one must sum over all magnetic quantum numbers, and in particular the external ones. This will lead to a closed Yutsis graph. To do so, an external coupling order has to be fixed. The coupling is used in the case of the T 1 amplitude equations and is such that can be used equally well. However, it turns out that the resulting equation will be much simpler when the first option is employed since the coupling order is consistent with the coupling order used for the Hamiltonian matrix elements. This is an example where prior experience provides a strong guidance for the proper choice of the angular-momentum coupling scheme even though ultimately all choices yield equivalent results.
To exemplify the coupling of open SU-TNs, one particular contribution to the CCSD doubles amplitude equation is chosen where (k, l, c, d) denote internal indices that are summed over while (a, b, i, j) characterize the external indices. The construction of the angular-momentum network originating from the application of the external coupling, defined in Eq. (102), to Eq. (106) requires to sum over a product of (i) two 3 jm-symbols coming from the external coupling of a, b and i, j, (ii) two 3 jm-symbols with zero-edges coming from the application of Wigner-Eckart theorem to the T 1 amplitudes, (iii) four 3 jm-symbols coming from the coupling of H and The red box indicates a subgraph to which the 4-cycle factorization rule can be applied. However, first applying twice the zero-line rule to the leftmost and rightmost nodes ( 1 j a δ j a j k and 1 j j δ j j j d ) directly yields a Yutsis graph that is topologically equivalent to the one of the third-order particle-hole contribution in MBPT, i.e., which corresponds to a 9 j-symbol. The final expression reads as × (ξ d j j |T 1 |ξ j j j )(ξ a j a |T 1 |ξ k j a ) t J 2 cbĩl × j i j b K j a j j J j c j l K j a j j J 1 j i j l J 2 j b j c K .

In-medium similarity renormalization group
As a final example, the IMSRG approach is considered providing a non-perturbative alternative to CC theory. Throughout the last decade, IMSRG has been successfully applied to various nuclear observables, including low-lying excited states and electromagnetic transitions whose treatments were pioneered and applied to mid-mass closed-shell nuclei in Ref. [53]. Without the use of angular-momentum reduction, such studies in the mid-mass regime would have been impossible from a computational point of view. Thus, non-scalar operators associated with, e.g., electromagnetic transitions constitute an excellent playground to yet extend the application of our automated treatment of angular-momentum reduction.

General formalism
The IMSRG formalism is based on a unitary transformation U(s) of operators parametrized by a continuous variable s ∈ [0, ∞) such that Equation (108) can be recast into a first-order ordinary differential equation (ODE) involving an anti-Hermitian generator η that can be chosen conveniently to obtain a desired decoupling pattern. A standard choice is given by the Wegner generator defined as the commutator of the suitably chosen 'diagonal' and 'off-diagonal' parts of H, the end result being that H od (s) is eventually driven to zero. Even though the initial operator may contain up to two-body parts only, the evaluation of the commutator in Eq. (109) increases the particle rank of the operator, thus, inducing many-body operators up to the A-body level. In practice, the IMSRG(2) truncation is typically employed in which operators of higher rank than two-body operators are discarded. As discussed in Ref. [54], the IMSRG(2) approximation is exact up to third order in MBPT for the ground-state energy while resumming large classes of higher-order diagrams.

Evolution of non-scalar operators
The form to Eq. (109) is completely generic and valid for an arbitrary Hermitian operator O, independently of its transformation properties with respect to S U(2) symmetry. For practical applications, the specific tensorial properties of O need however to be taken into account. The evaluation of the ground-state energy provides the simplest case since both O = H and the generator are scalar operators in this case.
In the general case where O is a spherical tensor operator of rank λ, the evaluation of the AMR form of the commutator appearing in Eq. (109) is key. The associated form can be generically written as where S λ 1 µ 1 and T λ 2 µ 2 are spherical tensor operator of rank λ 1 and λ 2 , respectively, which are subsequently coupled to give a tensor of rank λ. This coupling is obtained via spherical tensor product defined through where the left-hand-side is indeed a spherical tensor operator of rank λ.
While the complete list of contributions can be found in Ref. [53], the so-called particle-particle contribution to the two-body part of the evolved operator is considered as an example. The associated SU-TN expression is given by has been used. Equation (119) can be rewritten in terms of angular-momentum-coupled matrix elements (Eq. (45)) instead of reduced matrix elements giving

Conclusions
In the present work, an automated tool to perform symbolic angular-momentum algebra operations has been designed. This tool relates to the fact that the working equations, i.e. the symmetry-unrestricted tensor networks, at play in stateof-the-art nuclear many-body methods can be analytically reduced with respect to angular-momentum quantum numbers whenever they are effectively employed in a symmetryrestricted context. The corresponding time-consuming and error-prone derivation of the angular-momentum-reduced form of the tensor networks is thus performed in a matter of seconds. The design of the tool is based on the use of Yutsis graph representing networks of Wigner 3 jm-symbols and fulfilling sets of factorization rules whose repeated application eventually provides the angular-momentum-reduced form of the equations. While examples of applications have been provided for many-body perturbation theory, coupled cluster theory and the in-medium similarity renormalization group method, the code can be interfaced with any many-body formalism of interest. While the present paper focuses on S U(2) symmetry, extensions are envisioned for the future, e.g. to the subgroup of S U(2) at play in axially deformed nuclei, or to other symmetry groups.
In view of obtaining the error prone, fast and numerically optimized implementation of involved many-body formalisms, the present code serves as the missing link between an automated tool used to generate the initial symmetryunrestricted equations and and an automated tool used to produce the efficient source code dedicated to numerical applications.

Command-Line Interface and Input Files
For simple usage of the code, the amc program is provided. The amc program is a command-line interface to the code that can be used to reduce a set of equations and output the reduced equations to a LaTeX document. The unreduced equations are supplied as an AMC file, in a domain-specific language described in Sec. 7.2.

Command-Line Options
There are a few options, which can be passed to amc, that modify the behavior of the program: -o OUTPUT, --output OUTPUT Write the resulting LaTeX document to OUTPUT. By default the code strips the extension from the input file, adds a .tex extension, and creates a file of that name in the same directory as the input file.
--collect-ninejs Try to reconstruct Wigner 9 j symbols from products of 6 j symbols in the reduced expressions. This results in shorter expressions, but might obscure opportunities to identify intermediates, e.g., when some of the 6 j symbols only depend on the indices of single tensors.
--keep-trideltas Print triangular inequalities generated during the reduction process. Most often, these inequality constraints are implicitly contained in tensor variables, so removing them does not generate any loss of information. Constraints that are implicit in 6 j or 9 j symbols are never shown.
--wet-convention CONVENTION Switch the convention used for reduced matrix elements. Currently, the code supports two conventions: the wigner convention which has been adopted in the body of the paper, and the sakurai convention The wigner convention is used by default.

Angular-Momentum Coupling Language
To facilitate the use of the code the AMC language is provided to specifiy the tensors and their coupling schemes as well as to enter the equations to be coupled. The basic building blocks of the language are integers and fractions (i / j), strings mode The number of indices of the tensor. The mode is either specified as an even integer, e.g. 2 for a one-body operator or 4 for a two-body operator, or as a tuple of two numbers (x,y) to specify x creator and y annihilator indices. scalar A boolean indicating that the tensor is scalar (rank 0). The code exploits additional angular-momentum constraints for scalar tensors, and uses the unreduced matrix elements by default. reduce A boolean indicating that this scalar tensor uses Wigner-Eckart reduced matrix elements. This key is ignored for nonscalar tensors, which always use reduced matrix elements. The default values is false , so unreduced matrix elements are assumed. diagonal A boolean value that specifies whether the tensor is diagonal. Diagonal tensors have half the number of indices, i.e., a mode-2 diagonal tensor has one index. scheme The coupling scheme of the tensor. There are multiple ways to couple the angular momenta of the tensor indices. By default, the angular momenta of the first two creator indices are coupled, the resulting angular momentum is coupled with the third, etc., until all angular momenta have been coupled, and the process is repeated for the annihilator indices. This key accepts nested tuples that specify the coupling order of the tensor indices. Creator indices are numbered from 1 to x, annihilator indices from x + 1 to x + y. The elements of each tuple are either tuples themselves or index numbers. Index numbers may be negated to request coupling of the time-reversed state. latex The LaTeX command used to typeset the tensor in the program output. By default, the name of the tensor is used.
To give an example, d e c l a r e X { mode = ( 2 , 2 ) , scheme = ( ( 1 , − 4 ) , ( 3 , − 2 ) ) , s c a l a r = t r u e } declares a scalar tensor X with two creator and two annihilator indices, whose m-scheme, i.e. S U(2) uncoupled, matrix elements can be recovered via X pqrs = (−1) j s −m s + j q −m q × JM j p j s J m p −m s M j r j q J m r −m q M X J pqrs .
Equations are declared as variable = expression ; The variable on the left-hand side is a declared tensor with index subscripts, such as X_pqrs. Indices consisting of more than one character can be used by enclosing the subscript with braces and separating the indices with spaces, like in X_{k1 k2 k3 k4}. Index names can consist of letters, numbers, and underscore characters. The expression on the righthand side consists of sums of products of tensor variables, denoted by + and * operators. Two special operators are available: sum and P. The sum operator lists the indices to be summed over. It is used in the following way: Z_abcd = sum_pq ( X_abpq * Y_pqcd ) ; The subscript lists the indices. The same rules apply as for tensor variables. All indices have to be mentioned exactly once either on the left-hand side of the equation or in the subscript of the sum operator. The other operator is the permutation operator P. It supports two modes of operation: used as P(i j), it permutes indices i and j in the expression to its right. Used as P(i 1 . . . i m / j 1 . . . j n / . . . /k 1 . . . k p ) , it generates all distinct permutations between the index sets separated by slashes. Concretely, P(i/ j) = 1−P(i j), P(i j/k) = 1 − P(ik) − P( jk), and P(i/ j/k) = 1 − P(i j) − P(ik) − P( jk) + P(i j)P( jk) + P(ik)P( jk).
As an example, an equation arising from the three-body part C3 of the commutator of a normal-ordered two-body operator A2 with a three-body operator B3, needed for the IMSRG(3), can be entered like this: C 3 _ p q r s t u = 1 / 2 * sum_ab ( ( n b a r _ a * n b a r _ b − n_a * n_b ) * ( P ( pq / r ) * A2_pqab * B 3 _ a b r s t u − P ( s t / u ) * B3_pqrabu * A 2 _ a b s t ) + ( n b a r _ a * n_b − n_a * n b a r _ b ) * P ( pq / r ) * P ( s t / u ) * B 3 _ p q a s t b * A2_brau ) ; The tensors n and nbar are diagonal one-body tensors containing occupation numbers.

Organization of the code
The AMC code is organized into five modules: ast, output, parser, reduction, and yutsis. The ast module defines classes whose instances make up the abstract syntax trees that are processed by the package, along with some helper classes that simplify working with the trees themselves. The output module contains functions that turn abstract syntax trees back into other formats. Currently, it only contains a 22 module for LaTeX output. The parser module provides a parser based on the PLY parser generator [55] that produces abstract syntax trees from AMC files. The reduction module contains functions to perform the angular-momentum reduction itself. Finally, the yutsis module contains classes and functions for building and manipulating Yutsis graphs, as well as for simplifying the resulting expressions. The AMC package is directly executable, and the installer creates a wrapper named amc for convenience. Executing the package provides a command-line interface for parsing an AMC file, reduction of the contained expressions, and output of a LaTeX file.
The program flow is the following: First, the AMC file is parsed into an abstract syntax tree. The tree of each equation is expanded until it consists of a sum of products. For each term, a Yutsis graph is constructed according to the coupling schemes of the mentioned tensors. The reduction procedure looks for 2-, 3-and 4-cycles in the graph and applies the rules discussed in Sec. 4.6, iteratively factorizing the graph until it is completely expressed in terms of Kronecker deltas, triangular deltas, and 6 j-symbols. If enabled, a post-processing step tries to reconstruct 9 j-symbols by combining sets of three 6 j-symbols. The resulting abstract syntax tree is constructed by replacing all tensor variables with reduced ones and adding the objects resulting from the reduction of the Yutsis graph. This syntax tree represents the reduced equation, and is subsequently converted to a LaTeX expression and written to the output document.

Testing files
The AMC package contains 7 example input files along with the outputs generated by the amc program. The examples cover all applications discussed in section 5. Additionally, a more complex example is provided in the form of commutators of three-body operators that appear in IMSRG (3), and a file showing how to derive the Pandya transform of a scalar and a non-scalar tensor in a few lines of AMC code.

Methods
In this section only a pointer to the central methods is provided. See the API documentation that accompanies the package for more information.
parser.Parser.parse (instance method) Parse a string into an abstract syntax tree according to the AMC language grammar.

reduction.reduce_equation
Reduce an equation, given as an abstract syntax tree, to symmetry-restricted form.
output.latex.equations_to_document Turn a list of equations into a LaTeX document. The equations can be in symmetry-reduced or unreduced form.