Supercomputers against strong coupling in gravity with curvature and torsion

Many theories of gravity are spoiled by strongly coupled modes: the high computational cost of Hamiltonian analysis can obstruct the identification of these modes. A computer algebra implementation of the Hamiltonian constraint algorithm for curvature and torsion theories is presented. These non-Riemannian or Poincaré gauge theories suffer notoriously from strong coupling. The implementation forms a package (the ‘Hamiltonian Gauge Gravity Surveyor’ – HiGGS) for the xAct tensor manipulation suite in Mathematica. Poisson brackets can be evaluated in parallel, meaning that Hamiltonian analysis can be done on silicon, and at scale. Accordingly HiGGS is designed to survey the whole Lagrangian space with high-performance computing resources (clusters and supercomputers). To demonstrate this, the space of ‘outlawed’ Poincaré gauge theories is surveyed, in which a massive parity-even/odd vector or parity-odd tensor torsion particle accompanies the usual graviton. The survey spans possible configurations of teleparallel-style multiplier fields which might be used to kill-off the strongly coupled modes, with the results to be analysed in subsequent work. All brackets between the known primary and secondary constraints of all theories are made available for future study. Demonstrations are also given for using HiGGS – on a desktop computer – to run the Dirac–Bergmann algorithm on specific theories, such as Einstein–Cartan theory and its minimal extensions.


I. INTRODUCTION
The modern frontier in the search for alternatives to general relativity (GR) is characterised by a very large number of competing models [1][2][3].These models are problematic by dint of their heterogeneity: differences in the mathematical formulation hinder comparison between different classes (e.g.mimetic [4] and MOND [5,6] theories).This seems intractable, but we expect the number of classes to grow only in proportion to the community.A more serious problem presents when a class contains a large number of parameters, e.g. the couplings in a Lagrangian.Critical theoretical properties, such as the number and health of propagating degrees of freedom (d.o.f) may be sensitive to these parameters in ways which are hard not only to characterise in general, but even to calculate in detail for a given parameter choice.
This scenario commonly arises in the Hamiltonian analysis, used routinely to determine the d.o.fs of a proposed theory [7][8][9][10][11].If the coupling parameters elimniate some velocity ̇ from the motivated Lagrangian , then the total Hamiltonian  T must be modified to express the constraint ∼ ̇  ≈ 0. That constraint is only preserved if ̇ ∼ { ,  T } ≈ 0, which either vanishes identically or constitutes a new constraint, ad infinitum.Each constraint subtracts d.o.fs from the countable fields naïvely present in .The full chain of constraints is systematically elucidated by means of the Dirac-Bergmann 1.On a supercomputer, 448 processors obtain simplified, covariant expressions for all nonlinear Poisson brackets, among all discovered primary and secondary constraints, for 192 modified gravity theories.These theories are based on Einstein-Cartan gravity with an extra, massive spin-parity 1 + , 1 − or 2 + torsion particle, in which various parts (2 3 × 2 3 choices) of the curvature and torsion are disabled by multipliers.Each colour is a different theory, black/gray is initialisation.The objective is eventually to discover whether multipliers can be used to suppress the strongly-coupled 1 − , 1 + or 2 − modes, which respectively ruin these theories.Ready access to the Poisson bracket structure is vital for this analysis.
algorithm [7,10,12], the fundamental computational unit of which is the Poisson bracket.In theories constructed from the higher-spin representations of the Lorentz group, including most tensor theories of gravity, a single covariant bracket can be surprisingly cumbersome for manual evaluation [13,14]; less so, as we will find, for computer assistance (see Fig. 1).
Through the algorithm, the d.o.f and symmetry structure depend not only on eliminated velocities, but also on all brackets between all constraints, which themselves may be contingent on the couplings in ways that are, ab initio, unknowable.
The strong coupling problem is commonly diagnosed in the Hamiltonian picture.It is usually preferred that two gravita-tional d.o.f propagate [15,16].Additional d.o.f may be tolerated unless ghostly or tachyonic 1 , but they must kept under close theoretical and phenomenological control (e.g. by large [19] -or small [20,21] -masses, screening [22] or other measures).Frequently we start with many fields (ten in the case of GR), so that reduction to two d.o.f is an achivement in the linearised Hamiltonian picture: strongly coupled modes may increase this number the nonlinear analysis [23].Strong coupling can be imagined as a finely-tuned suppression of the kinetic coefficients2 in a mode's linear wave equation [24].The nonlinear operators are generally still present, however, and describe a non-perturbative dynamics which may be unacceptable, having for example an elliptic or parabolic character.The linearisation in this case refers to some motivated exact vacuum solution.This solution need not be Minkowskian: for GR in higher (odd) dimensions, a fine-tuned admixture of the Gauss-Bonnet invariant strongly couples the whole graviton on a maximally symmetric but curved background [25].In many versions of Hořava gravity [26], the 'detailed balance' which defines the Lorentz-asymmetric theory has a similar effect [27,28].
There are some prominent scenarios where strong coupling seems desirable.The linearised, massive graviton of Fierz and Pauli [29] appears troubled even in the massless limit by its fifth d.o.f -the helicity-0 mode or van Dam-Veltman-Zakharov (vDVZ) scalar [30,31] -which couples to the matter trace.However, nonlinear completions [32,33] of Fierz-Pauli theory revealed that rather than persisting as a light d.o.f, the vDVZ scalar instead becomes strongly coupled, i.e.Vainshtein-screened [34,35].Unfortunately, these same theories were also plagued by a sixth d.o.f -the nonlinear Boulware-Deser ghost [32] -whose origin was apparently connected back to the same strong coupling effect [36].These matters remain strongly contested [37], and massive gravity continues to evolve [38,39].
Strong coupling in massive gravity had another (equally contested [37]) association with superluminality [40,41].Both properties are sometimes conflated with acausality, and illposedness of the Cauchy-Kovalevskaya problem.While systems might exist in which these are all mutually implicated pathologies, it is by no means general.In each case one should study the characteristic surfaces and sift for coordinate artifacts [42,43].Strongly coupled or singular surfaces in the Hamiltonian phase space, where the rank of the matrix of constraint brackets changes, have also been accused of causality violation in the Poincaré gauge theory (PGT) of gravity [23,44,45].PGT is the general theory which may be constructed from spacetime curvature and torsion.Whether or not the same caveats resolve the causal question here [13], the proclivity for strong coupling in this class remains, and has yet to be seriously addressed.To this end, we target the Poincaré gauge theory class in this paper.
The Poincaré gauge theory encompasses a large sector of the general category of non-Riemannian theories, whose current popularity warrants a quick introduction.Geometrically interpreted, we retain for non-Riemannian gravity the system of clocks and rulers in the gravitational metric potential .In GR, however, the connection Γ is fixed to the Levi-Civita form ≡ ( ) ≡ The multiplers 3and ̂ constrain the geometry to be torsion-free and metric, while the dynamics are those of curvature.At the time of writing, the attention of the community is drawn to the remarkable non-uniqueness of (1) as a realisation of GR in this broader, non-Riemannian context.By cycling multipliers onto different pairs in the curvature-torsionnon-metricity triad, and identifying suitable dynamical terms, one can construct the flat, metric teleparallel equivalent of GR (TEGR) from pure torsion [47], and the flat, torsion-free symmetric alternative (STEGR) from pure non-metricity [48].
With (1), these points in the space of non-Riemannian theories define the geometrical trinity of gravity [46].
There is now activity to determine a preferred vertex of the trinity, and the viability of the surrounding non-Riemannian landscape.Various considerations must be balanced in a very large parameter space [49][50][51][52].As we mentioned previously, the PGT may be thought of as the sector of the landscape in which is suppressed, e.g. by a multiplier as in (1).PGT is a convenient sector to study -without prejudice to the ultimate rôle of non-metricity in constructing a viable gravity theory -since it has a relatively self-contained history stemming from the Einstein-Cartan model [53], and the gauge-theoretic interpretation pioneered by many authors [8,[54][55][56], beginning with Kibble [57], Utiyama [58] and also Sciama [59].
The linear Hamiltonian structure of PGT is particularly well developed [8,[60][61][62].In the nonlinear structure, strong coupling phenomena are know to be abundant, and in this regard PGT is feared to be representative of the broader non-Riemannian landscape.Critically however, the challenge of the Hamiltonian analysis has prevented this structure from being mapped in any comprehensive detail.A few islands in the landscape were probed at the turn of the millenium [23]: minimal extensions to the Einstein-Cartan theory in which a single extra massive spin-parity ( ) 1 + , 1 − or 2 − torsion particle is present.In each case respectively, 1 − , 1 + and 2 + modes were suggested to be strongly coupled.Since then, the PGT sector of the landscape has been treated with a sense of 'hic sunt dracones' -i.e.avoided as potentially dangerous -with apparently the only safe configurations being exclusive activation of the 0 + or 0 − scalar torsion particles [45].In light of our comments above, one can even afford to remain agnostic on the pathology of this strong coupling 4 .It would seem, however, given the recent, promising developments in the non-Riemannian approach [46,[49][50][51][52], that the extent of the phenomenon should still be understood, and general tools be developed to that end.
In this paper we present a computer algebra implementation of the nonlinear Hamiltonian analysis for a generalisation of the full PGT, in which one may covariantly disable arbitrary irreps of the torsion and curvature by means of Lagrangian multiplier fields.This generalisation was put forward as an antistrong-coupling measure in [14], and its Hamiltonian structure is elucidated in the companion paper [64].In the conventions of (1), we may write generalised PGT as where the usual 'quadratic' PGT is spanned by the ten parameters ̂ 0 , { ̂ }, { ̂ }.The projections ̂ and ̂ extract the SO + (1, 3) field strength irreps.The core of the implementation is (version 1.0.0 of [65]) the Hamiltonian Gauge Gravity Surveyor (HiGGS), a Mathematica package grounded in the popular open-source xAct tensor manipulation suite [66][67][68][69][70][71].
The HiGGS package is suitable for targeted use on a desktop computer.Since it is parallelised over Poisson brackets, HiGGS also scales to clusters and supercomputers.Development in this direction is with the aim of surveying the constraint structure of the non-Riemannian landscape at scale.Modules from HiGGS, in particular those concerned with bracket evaluation, can be used as a back-end in searching the parameter space for desirable canonical features -as such features become better understood with time. 4It is important to understand that the main historical objection to the strongly coupled modes of the PGT has been their ghostly character, as inferred by an inspection of the signs of squared momenta in the Hamiltonian.These signs are fixed by the unitarity requirements of the desired, linearly active modes in each case, and in each case they are negative [23].This 'catch-22' is particular to the PGT, but, as cogently explained in [63], there are principled reasons to be suspicious of stongly coupled surfaces, which appeal to neither the non-perturbative dynamics, nor the causality arguments mentioned above.Generically, a background which is strongly coupled had better not be one which is also seen in nature, since it cannot have been reached by any smooth trajectory through the phase space.
For the moment, we perform the brute-force 'calibration' survey illustrated in Fig. 1.In this run, which takes a little over 1 h, the 1 + , 1 − and 2 − Einstein-Cartan extensions are modified with all possible configurations of curvature-and torsiondisabling multipliers: 192 generalised Poincaré gauge theories in total.During the run HiGGS obtains, for every theory, all the primary and secondary constraints which can be inferred from a knowledge of the literature, and then computes simple covariant expressions for the nonlinear Poisson brackets between all constraint pairs.All brackets identified in this survey can be found in the supplemental materials [72].
Examples are also provided of how to use HiGGS to calculate constraint velocities when implementing the Dirac-Bergmann algorithm.Focus is on the minimal Einstein-Cartan extensions with strong coupling, and the viable 0 + and 0 − extensions.The unmodified Einstein-Cartan theory is also studied, and HiGGS is used to show that it propagates only the two graviton polarisations.Finally, we will use the results of the initial survey to show how multipliers might conceivably be used to suppress strongly coupled fields.Note that the major undertaking of fully analysing the results, isolating and confirming any viable multiplier configurations, is left to future work.
The remainder of this paper is structured as follows.In Section I A we briefly set out our conventions for PGT, using the non-geometric, gauge-theoretic formulation.In Section II we describe the HiGGS implementation, including general tools for the canonical manipulation curvature and torsion, up to Poisson brackets and higher-level functionality for the theoryspecific calculation of constraints and velocities.Solutions for scaling the Hamiltonian analysis to high-performance computing (HPC) resources are also described.In Section III we present examples of the algorithm, and the results of our initial survey in Fig. 1. Conclusions follow in Section IV.
The geometric covariant derivative appearing in the definition ≡ ∇ , as it acts on a vector , is written In the Poincaré gauge theory we enforce = 0 by assumption, though a multiplier could also be used.The general non-Riemannian connection thus loses its disformation part ≡ where we define the (inverse of the) translational gauge field ℎ and the rotational gauge field ≡ [ ] -we will return to this derivative in Section II C.These fields guarantee invariance under the general coordinate transformations (GCTs or passive diffeomorphisms) on , and Lorentz rotations, and together they gauge the Poincaré group ℝ 1,3 ⋊ SO + (1,3).How to connect back to the geometric picture?The Minkowskian metric components can be recovered via the identities ℎ ≡ and ℎ ≡ .On the other hand, the usual system of clocks and rulers can be recovered using ≡ and ≡ ℎ ℎ .The covariant measures on  and  are respectively √ − , i.e. the conventional ≡ det , and ≡ ℎ −1 ≡ det .Finally, the geometric field strength tensors -referred to as the curvature and the torsion -are provided by the formulae The conversion of these back to the (numerical values of) the geometric components, is done by contraction with the translational gauge fields.The theory (2) meanwhile becomes where our metricity assumption dispenses with the need for the final term in (2), and we use the Planck mass rather than the Einstein constant p ≡ 1∕ √ .Once again, we reiterate that our use of the gauge theory setup over geometric alternatives (including those which use the tetrad and spin-connection over the metric and contorsion) is merely a matter of convenience in the present work.
In this article we will follow [77] by using the following syntax highlighting for code listings: keywords for Mathematica are typeset in green, for xAct in blue and for HiGGS in red.This paper uses the 'West coast' signature (+, −, −, −).

II. IMPLEMENTATION
In this section the implementation of the Hamiltonian analysis in the HiGGS package is described.Many aspects of our particular approach will be inefficient: the package is monolithic and expensive in terms of memory and maintenance.At the time of writing, however, HPC is a cheap resource.By carefully tuning only a few aspects of the implementation it is therefore possible to produce a product which surveys the theory space in a matter of hours.It is important to emphasise that these 'high-level' features of HiGGS are currently limited to the generalised Poincaré gauge theory in (7).The question of scaling to arbitrary theories will be addressed in Section IV, where we will find that such a scaling is actually likely to reduce the complexity of the package, so long as the PoissonBracket[] function is upgraded so as to be 'aware of' Leibniz's rule.Moreover, the 'low-level' functionality, which is not specific to a Lagrangian, will be of standalone utility.In what follows, we assume a basic familiarity with the xAct suite and Wolfram Language.

A. Geometric setup
We begin by describing the way in which the Riemann-Cartan geometry is implemented in HiGGS.Whilst xAct is perfectly capable of accommodating not only a Riemann-Cartan curvature, but also a torsion tensor, we prefer to adhere to the 'particle physics' picture of gravitational gauge theories [8,14,56,78], and set up all the physics on a flat spacetime.The following equivalent xAct commands are issued when the HiGGS environment is initially built (see Section III A) with the command BuildHiGGS[] This offers several advantages, and foremost among these is an easy comparison with the very substantial body of literature on the Hamiltonian structure as mentioned in Section I A. We will also not be limited by the fact that the non-Riemannian features of xAct are (presently) less flexible and comprehensive than those of a simple Minkowskian setup.This is to be expected, given that GR is the preferred effective theory of gravitation.An essential feature of this setup, which goes beyond the particle picture of the literature, is the conflation of holonomic and non-holonomic tangent spaces: a single collection of indices a, b, c, etc. is ascribed to TangentM4, and these represent both the Lorentz indices , , and coordinate indices , , .This practice might be anathema to the field of differential geometry, but it is acceptable for pragmatic, computational purposes.The coordinates are assumed implicitly to be Cartesian, so that the components of the flat metric are Minkowskian (equal to those of ), and moreover that a rotational gauge is chosen in which 0 = ̂ 0 , 1 = ̂ 1 , 2 = ̂ 2 and 3 = ̂ 3 .As a consequence of this gauge-fixed setup, the covariance of quantities is not guaranteed internally, as it would be if all the features of xAct (such as user-defined connections) were fully exploited in HiGGS.It is instead possible to check the covariance by visually inspecting the final results for explicit gauge fields: in practice this turns out to be very easy, and the elimination of bare gauge fields in HiGGS is reliable.
It is important to note that whilst the gauge is fixed to conveniently overload the indices of TangentM4, the tensor structure of PGT is wholly preserved.This is in contrast, for example, with the seminal paper [60] in which the time gauge imposes ℎ 0 = 0, thereby massively simplifying the various algebraic expressions.On the contrary, the algebraic expressions produced from HiGGS are valid for any gauge, once the various shared indices are interpreted as being either Greek or Roman.Care is taken in the definitions to avoid any ambiguity over this division of indices.
Following on from our discussion of indices, we note that the 1, 2, 3 spacelike indices , , and , , are also subsumed into a, b, c.These may be extracted by means of the projection operator G3[a,b], which represents , and lies at the heart of the Arnowitt-Deser-Misner (ADM) split.The ADM or 3 + 1 split uses a spacelike foliation which is characterised by timelike unit vector , defined as where we recall that the (gravitational) metric components are recovered by ≡ ℎ ℎ .Any vector with local Lorentz indices may then be decomposed into perpendicular and parallel components  =  ⟂ +  .An overbar is used to denote the parallel indices.There are then some identities ℎ ≡ and ℎ ≡ , which follow from (8).The gauge fields ℎ , and are denoted by H[-a,b], B[a,-b] and A[a,b,-c].The timelike vector is represented by V[a], and we see that a variety of identities are then implied within the built HiGGS environment In[ ]:= ToCanonical[quantity] These functions carry information about the part 0 of the translational gauge field (which, as we will shortly see, is nonphysical), wheres is independent of this quantity.In HiGGS the lapse is Lapse[] and we use also the spatial measure J[] for the quantity ≡ ∕ , where ≡ det .Some further identities are then These identities are also incorporated into the HiGGS environment, which prefers to extract -from all derivatives of quantities dependent on the translational gauge field -the form and so we find Now that the geometric setup is in place, we introduce the canonical setup [7,8,79].The canonical momenta are As with our earlier work [13,14,64], we will neglect the matter Lagrangian M .The field strengths in (6a) and (6b) are independent of the velocities for 0 and 0 , so (11) imply 10 primary constraints of the form From ( 12) we arrive at the result mentioned above, that the conjugate field 0 is non-physical; the same applies to 0 .By (≈), the weak equality is denoted, i.e. not an approximation.The constraints (12) are referred to as the 'sure' primary, first class (sPFC) constraints: they are a consequence of Poincaré symmetry, and they apply for all choices of the { ̂ }, { ̂ }, { ̄ }, { ̄ }.According to (11) The purpose of the parallel momenta is seen above: they are the physical parts not touched by the sPFCs, and defined according to ̂ ≡ and ̂ ≡ .We will see in Section II D that the specific internal rule PiPToPi should not often need to be used in practice, and has a more general alternative in the ToBasicForm[] command which is provided officially by the package.In general, parallel and perpendicular quantities can be accessed with some projections The flat manifold of course has vanishing Riemannian curvature = 0, and so we must be careful not to use the Instead, the field strengths  and  are given by their own tensors, and these expand to give the definitions in Eqs.(6a) and (6b) The final set of fields which we introduce are the multipliers and -these are precisely the same shape as  and  , and we write them It should be emphasised again that rules such as ExpandStrengths, PADMActivate and PiPToPi should not often be needed, and are subsumed under the ToBasicForm[] command.

B. Irreducible decompositions
The Hamiltonian analysis invites decomposition into irresducible representations of the Lorentz group SO + (1, 3), and (through the ADM split) the special orthogonal group SO(3).The former is useful also in the Lagrangian picture, since the Lagrangian formulation of a theory is typically Lorentzcovariant, it is useful to split field representations into blocks which transform only among themselves.In the Hamiltonian case, the choice of slicing introduces a 'preferred' timelike vector -this is not unique, so covariance is not ultimately lost -but the symmetry in the context of the slicing is reduced to the spatial rotations SO(3).In the case of rotations, the irreps of tensor fields correspond to states of definite spin and parity, allowing us to designate the parts as .

Lorentz group
We do not often encounter the Lorentz decomposition in the course of Hamiltonian calculations, but it is implemented across the higher-rank tensors, for example the 'humanreadable' decomposition into familiar irreps such as the Weyl the Ricci and the Ricci scalar and the tensor torsion, defined as the remainder of  after removing the vector (i.e.torsion contraction) and pseudovector components -see e.g.[23,45,80].These components may be accessed as follows Out []= 2ℛ ab c d More commonly, we might need access to the complete, orthonormal SO + (1, 3) operators which appear in (7).These may be associated with the couplings { ̄ } as follows, for the purposes of constructing a Lagrangian5 In[ ]: From this we see how to access the  and  projections.

Rotation group
The SO(3) projections are similarly defined, but of course refer to the vector , not just the Lorentzian metric .The par-allel momenta can be decomposed in this way: for the translational case where the second term in (16a) is the 1 − vector mode, and the terms in (16b) are respectively the 0 + scalar, skew-symmetric 1 + vector and symmetric-traceless 2 + tensor modes.The rotational case is similarly decomposed as where (17b) are the 0 + , 1 + and 2 + modes and (17c) are the 0 − , 1 − and 2 − modes.Accordingly, we access these as follows The parallel parts of the field strengths can also be decomposed, since they are canonical, but we neglect the perpendicular parts which depend on the unphysical fields 0 and 0 , and also on non-canonical velocities.The field strength decomposition is We notice that  and  ⟂ both share all six representations present in the rotational momentum.In the parallel case, these are denoted by the 0 + scalar , 1 + dual vector -just as with the translational momentum.
In HiGGS, only the canonical, parallel parts of these tensors are decomposed.The multipliers share the same tensor structure, and both their parallel and perpendicular parts are decomposed: this is because all multiplier fields are assumed to be canonical.The resultant decomposition is


Note that our typeset convention for the SO(3) multiplier irreps will be simply to recycle the field strength expressions above, substituting the symbols ,  → .This is a general rule unless stated otherwise: i.e., tensors whose (Lorentz-invariant) index structure are identical will use the same SO(3) irrep notation, but a different symbol.It is possible, if needed, to recover the explicit 'humanreadable' projection operators used in the above calculations, such as appear in ̂ ́ ≡  ́ ̂ and ̂ ́ ≡  ́ ̂ , where the sectors are represented by and indices according to our conventions in [64]; the ́ accent denotes variable indices, and these projections are defined as obtaining precisely the irreps in Eqs.(16a), (16b), (17b) and (17c).To reflect the various conventions used in the literature, these projections are sometimes defined differently for spin representations contained within parts of the field strengths and momenta whose tensor structures are identical -for example ̂ and  .We therefore use different symbol names for momenta and field strength projections (or those of multipliers), and find e.g. for the selection


Note above that it is sometimes necessary to repeatedly apply routines in order to recover the desired form.In prinicple, the HiGGS commands are constructed so as to be idempotent, but broadly in Mathematica and xAct repeated commands are sometimes helpful when multiple functions are nested.

C. Derivatives
While the gauge-fixed Minkowskian setup initially makes for easy development, it sometimes means that we have to 'reinvent the wheel' in order to access machinery for which there is already a very sophisticated implementation in xAct.A clear example of this is given by the way in which HiGGS handles gauge-covariant derivatives in the Poincaré gauge theory.Generalising from (5), we recall from [8,14,73,74] that for some matter field we have where Σ are the Lorentz group generators specific to the representation of .This construction could be implemented using the xAct command DefCovD[], in such a way that the Σ generators are automatically calculated for tensorial representations of the kind that arise in the Hamiltonian analysis of the matter-free theory.In particular, the derivative is geometrically interpreted as ∇ as it appears in (3).Rather than following this route, HiGGS defines a pair of first derivatives (broadly corresponding to the two definitions in ( 19)) for every canonical quantity which might be of interest.This process is not very efficient or flexible, but it is sufficient when explicit covariant derivatives are scarce.Indeed, the most common occurance of the gauge covariant derivative is through its commutator in the form of the field strengths: as we saw in Section II A, these have a separate implementation.Explicit gradients will tend to arise mostly at the end of the calculations for which HiGGS is designed.Gradients of the field strengths cannot usually be fed back into the algorithm anyway, since they would require an implementation of the second order Euler-Lagrange equations in order to be processed.
Based on the understanding that we are always dealing with arbitrarily-indexed (possibly mixed) tensorial representations, HiGGS' two preferred forms of the derivative are ́ ́ and The former is straightforward 6 , and the latter is the parallel projection of the gradient on all indices.Picking a couple of SO(3) irreps at random, we find for example One cannot generally pass from ́ ́ to a parallel derivative, unless only the spacelike indices are involved, e.g.́ ́ .Since HiGGS tries to achive parallel derivatives wherever possible, this serves as one of the internal checks on the canonical status of a quantity, making sure that any unacceptable time derivatives would appear explicitally in the output.Accordingly, continuing from above We see in the penultimate expression above the two parallel derivatives.If an expression containing gradients is covariant, velocity-independent and parallel, then the various other terms should cancel among themselves.As we have mentioned, covariance also requires that gauge fields do not explicitally appear, but are implicit in covariant quantities.Overall, the work of packaging an expression into a canonical, covariant form, parallel if possible and free from unphysical fields, is done by the ToNesterForm[] command, which we introduce in Section II D. As a final comment on the use of derivatives, we recall that the formulae Eqs.(6a) and (6b) imply a pair of Bianchi identities [8,80] as follows These can be verified using the tools already introduced, but the process will be easier with the ToBasicForm[] and ToNesterForm[] commands.

D. Low-level functions
Throughout Sections II A and II B we have used many rules, such as DpRPDeactivate, DPiPActivate, etc., whichwhilst they are not confined to xAct'HiGGS'Private' and so are available to the user -should not often be needed within a HiGGS science session.Instead these rules are wrapped into two 'official' functions: ToNesterForm[], which strives to collect expressions 7 , and ToBasicForm[], which expands them.Naturally ToNesterForm[] is the more complicated of the two.It is, in some sense, the extension of the ToCanonical[] command from xAct into HiGGS.

Module: 'ToBasicForm'
We begin by breaking some expressions which we know to be covariant, using ToBasicForm[] We see that the results are expressed in terms of the bare gauge fields, non-parallel momenta, independent ADM quantities and coordinate derivatives.If a Poisson bracket were to be manually evaluated between the fields in quantity, it would first be necessary to perform the expansion provided by ToBasicForm[], before calculating variational derivatives and multiplying.For this reason, the PoissonBracket[] command which we introduce shortly relies on ToBasicForm[] as an initial step.In the case above, all the brackets would be relatively straightforward to evaluate by hand: to see why the strong coupling problem demands an implementation such as HiGGS, let us break open a simple momentum gradient The gradient of a more complicated irrep such as DpPiPA2m[-i,-j,-k] would span many pages after an application of ToBasicForm[], rendering manual evaluation of brackets impractical, yet we will see in Section III that there is nothing to stop these terms arising in the analysis. 7The chosen naming refers to the fact that the irrep conventions and notation of collected expressions tends to align most closely with those of a collection of articles -very useful during development -for which Nester is a common author, see e.g.[23,44,45,54,76,[82][83][84][85][86][87].
Now is a convenient time to verify the second Bianchi identity in (21b).Setting up the derivative, and using the gaugefixed which follows from -epsilonG[a,b,c,d]we find Evidently, this quantity vanishes, and so what we are seeing is a side-effect of the HiGGS geometric setup: xAct does not know that CD[-i][] is the partial derivative, so we need a final step Out []= 0

Module: 'ToNesterForm'
Let us now undo the breaking, using ToNesterForm[].In general, much of the functionality of ToNesterForm[] has to do with fully incorporating the known primary and secondary constraints of the theory during the course of simplification -so as to leave no extra conditions for the user to worry about.Before a theory has been defined (see Section II E) using DefTheory[], we will have to suppress this activity by passing the option "ToShell"→False, yielding In[ ]:= quantity = ToNesterForm[quantity, " → ToShell" → False]


Note that whilst the results are covariant, we do not get back the same form from ToNesterForm[] as that which was provided to ToBasicForm[], instead we recover the SO(3) expansion.The reason for this is that reducible quantities such as T[-a,-b,-c] are usually quickly broken up during the course of the Hamiltonian analysis: the output of ToNesterForm[] is tuned so as to be useful in that context.Moreover, we note that this operation will not now work in reverse, since input such as TP2m[-a,-b,-c]//ToBasicForm will not return a broken expression.There is no special reason behind this: ToBasicForm[] is not a sophisticated function, it essentially imposes a list of internal rules on its argument, and the expansion of SO(3) field strength irreps, as with those of momenta, would be perfectly straightforward to implement in the source if and when needed.
When ToNesterForm[] is passed a non-covarint quantity, it is unable to return a covariant result.Nonetheless, it tries, returning for the innocuous gradient As with this case, it is usually easy to identify nonphysical expressions which indicate a human error, through the appearance of bare gauge fields.If needed, this covariance check is easy to automate through an output search with the Mathematica Head[] function for an unwanted HiGGS quantity, such as A or B. We can also observe in the output above part of the route taken by ToNesterForm[].Gauge field gradients are converted, where possible, to field strengths and covariant derivatives.The residual spin connection terms are then extensively manipulated in an attempt to cancel them.Leftover asymmetric derivatives of the translational gauge field can sometimes be ascribed to covariant quantities through (10), for example the following gradient cannot be expressed through the torsion alone Note that the above result could equally be written as the single term -DpV[-q,-y], but ToNesterForm[] takes the opportunity to separate out the antisymmetric part.Before moving on, we return to verify the first Bianchi identity in (21b).We can access the canonical (i.e.velocityindependent) part of this identity by projecting Eq. (21b) with ⟂ , or equivalently using in place of the foliation equivalent ⟂ , which is given in HiGGS by Eps[i,j,k], In[ ]:= quantity = ToNesterForm[quantity, " → ToShell" → False]; The output here is then equal to the SO(3) irreps contained within ⟂  , as expected.

Module: 'PoissonBracket'
The PoissonBracket[] command is the third 'official' function provided by the HiGGS package.The Poisson bracket appearing in this article is defined for general functionals  and  of the gravitational fields and their conjugate momenta with a natural extension of the formula when multiplier fields are admitted.The formula (22) may appear no more daunting than a commonplace action variation, but in practice  and  are frequently local tensors rather than nonlocal scalars.Locality signifies that the underlying functionals contain Dirac distribuions, themselves subject to the total derivatives of the generalised Euler-Lagrange equations.The full ramifications of covariantly removing these Dirac gradients are detailed in [14,64], and some special cases are discussed in electrodynamics on the lightcone [88] and noncritical string theory [89].Currently, HiGGS is able to accommodate the first order Euler-Lagrange formalism in PoissonBracket[].
First order brackets, evaluated by inserting the spatial depen-dence into (22), produce four terms of the form where the 1 , 2 , 3 and 4 can be determined by certain formulae.Note that by our conventions in [13,14,64], we will in future denote by 3 the equal-time Dirac function 3 ( 1 − 2 ).Without any special instructions, PoissonBracket[] returns a List of the four Dirac coefficients in (23).This behaviour is tied into calls to PoissonBracket[] from within Velocity[], which we introduce in Section II E. Note that PoissonBracket[] contains calls to ToNesterForm[], and so works to exhaust transformations which can be applied to the output by virtue of the known primary and secondary constrints.In this sense, it depends on the theory introduced by DefTheory[] and so for the time being we must again pass the option "ToShell"→False.We begin with a very simple bracket, whose output can be understood in terms of our previous "ToShell" → False]; Out []=  .
The bracket is not 'surficial', in the sense that the latter three entries vanish and the nonvanishing part of the bracket is a compact, covariant expression.If the latter three coefficients are nonvanishing however, explicit covariance in this expression will be lost, as we can see by modifying the previous example In such 'surficial' cases the default output of PoissonBracket[] will not be helpful for visual inspection.To resolve this we can recall, again from [14], that ( 23) can be alternatively expressed as The three-component List output corresponding to ( 24) can be produced by passing the option "Surficial"→True.Trying again with this option, we obtain → "ToShell" → False, "Surficial" → True → ]; We see that the result is indeed covariant, and fairly simple.We will come back to this 'surficial' case in Section III B, where we attempt to recover the historical results in [23,45].

E. High-level functions
As we mentioned in Section II D, much of the work done by ToNesterForm[] has to do with the imposition of the theory shell so as to simplify the argument.The particular shell used is not specifically that of primary vs secondary constraints, but it is restricted to the constraints of which we have prior knowledge from the literature, and does not include new constraints discovered in the course of a HiGGS session.In particular, we rely on the so-called if-constraint structure which was discovered by Blagojević and Nikolić in [60].Depending on the Lagrangian parameters in (7), the number and type of primary constraints may be radically different: these contingent primaries are called primary if-constraints (PiC).There are similarly secondary (SiC) and tertiary (TiC) quantities, etc. Returning to [64], the PiCs of the theory (7) take the form where we defined e.g.̂ ∥∥ ≡ ∑ ∥∥ ̂ , with the matrix In the handling of PiCs, it is extremely useful to refer to the functions and we note that the PiC functions defined in Eqs.(25a) and (25b) are only constrained when ( ̂ ⟂⟂ ) = 1 or The structure of Eqs.(25a) and ( 25b) is quite useful in that it allows momenta to substituted for parallel field strengths.These substitutions are among those performed every time ToNesterForm[] is called, unless the option "ToShell"→False is passed as above in Section II D.

Module: 'DefTheory'
In order to discover these PiCs, we must first set up the shell using DefTheory[].Let us consider the simple example of Einstein-Cartan theory, i.e. the substantial restriction of (7) to the simple Einstein-Hilbert term We implement the theory ( 27) by passing to the DefTheory[] command the system of equations which deactivates all the { ̂ }, { ̂ }, { ̄ }, { ̄ }, while leaving ̂ 0 untouched.We want to store our knowledge of the shell once it has been obtained, and so we pass the label "Export"→"EinsteinCartan", which will be used to construct a filename.The input is Out []= ** DefTheory: Found the following primary if-→ constraints: ** DefTheory: Found the following secondary → perpendicular if-constraints: ** DefTheory: Found the following secondary → parallel if-constraints: ** DefTheory: Found the following secondary → singular if-constraints: ** DefTheory: The super-Hamiltonian is: ** DefTheory: The linear super-momentum is: .
** DefTheory: Exporting the binary at svy/ → EinsteinCartan.thr.mxIn the output above we can see the listing of the PiCs.It is clear from (27) and also from the form of the constraints in Eqs.(25a) and (25b), why we get this specific list.All the coefficients associated with the SO + (1, 3) irreps of  and  are vanishing, and so all the PiC functions become constraints.Moreover, these constraints are very simple in their form: they are all pure momenta, with the single exception being in the case of the 0 + roton constraint, which contains a constant Following this listing of the PiCs, there are references to perpendicular, parallel and singular SiCs, none of which are present in the Einstein-Cartan theory.We will return to these in Section III C, but for now we note that they follow from the imposition of multipliers in (7), and their presence is fully understood in [64], just as the presence of the PiCs of the basic Poincaré gauge theory is understood in [60,62].This is the intended scope of DefTheory[]: to elucidate not only the primary constraints, but all of that part of the constraint structure which is already known from the literature.Every if-constraint identified by DefTheory[] is used during the DefTheory[] call in the construction of very large internal rule sets $StrengthPShellToStrengthPO3, $PiPShellToPiPPO3, $TheoryCDPiPToCDPiPO3 and $TheoryPiPToPiPO3.These rules are applied -by defaultduring ToNesterForm calls.Most of the if-constraints that can arise, as detailed in [64], are of the form 1 ̂ ́ + ⋯ ≈ 0 or 1 ̂ ́ + ⋯ ≈ 0, and so the shell is defined by replacing all possible instances of the momenta in favor of other quantities.There are also the 'parallel' SiCs, which deactivate irreps within the canonical parts of the field strength tensors in Eqs.(18a) and (18b) -these are straightforwardly implemented.
Following the listing of the if-constraints, the output above details the structure of the canonical Hamiltonian, which from [64] is written This Dirac form [61,90] is useful because it expresses part of the Hamiltonian as a linear combination of the nonphysical fields , and 0 .The physical (and canonical) coefficients of these undetermined fields are 8 and they form the 'sure' secondary FC constraints (sSFCs) We see that these 10 constraints are divided up under SO(3) to give the final four entries in the output above, with . The (linearised) values returned for  ⟂ and  seem sensible in the context of the constraint in (28) - 8 Note that ≡ 0 + 1 2 0 .
we can see from this that HiGGS has begun imposing the PiC shell.What about the angular super-momentum?Let us verify by hand that the answer is correct.The first term in (30c) vanishes on the PiC shell, because HiGGS has told us that For the second term, we find from (28) after a few lines that from which we determine These are indeed the 1 + and 1 − parts of the torsion, as HiGGS is claiming.
The linearised velocity of (e.g.) some PiC ́ is calculated using the formula where from (29) we need evaluate only the commutator with  ⟂ , which is re-expressed in [64] on the PiC shell as A particular problem now is that ( 35) has a quadratic structure.If we apply PoissonBracket[] to evaluate a velocity for some PiC { ́ ,  ⟂ }, we will actually be calling ToBasicForm[] directly on  ⟂ .This will generally result in a very expensive computation, which we avoid in Velocity[] by manually applying the Leibniz rule.Accordingly in the HiGGS source there is a collection of prepared formulae for generic expressions such as which collectively act as a template for each velocity.For each SO(3) index shown in the sum in the second equality in (36), we need the four portions of a single Poisson bracket, as given in (23) and returned as a List by a call to PoissonBracket[].Accordingly, the whole of ( 35) is broken into blocks, each of four terms, which follow from a single bracket, and these are evaluated in parallel as discussed in Section II F. As we mentioned in Section II D 3, the formula ( 23) is not explicitally covariant, and hence the appearance in (36) of partial derivatives.Ultimately, a call to ToNesterForm[] is needed to restore explicit covariance to the quadratic expression (36).This setup is not so good: the use of block formulae restricts us to calculating velocities over the given  ⟂ -it is also susceptible to human error, and concludes with an expensive simplification process.We propose that any future iterations of HiGGS implement velocities based two extensions to

PoissonBracket[]:
1.The operands should be expressed as sums of products of covariant quantities, and the Leibniz rule be used to distribute the operation.2. The operation should return the covariant form given in (24), so that covariance is maintained at all steps.An improvement of PoissonBracket[] along these lines should be straightforward to implement, and would also be equally emenable to parallelisation: we defer its development to future work.

Module: 'ViewTheory'
The final module ViewTheory[] is used to produce a human-readable summary of the results from StudyTheory[].Since we already have a summary of the literature constraint structure from the DefTheory[] call in Section II E 1, we pass for brevity our earlier theory name with the option "Literature"→False, producing In[ ]:= ViewTheory["EinsteinCartanVelocities", " → Literature" → False]; Out []= ** DefTheory: Incorporating the binary at svy → /EinsteinCartanVelocities.thr.mx We see how the results are stored in the theory binary at What about the remaining PiCs P and T ?These 0 − and 2 − irreps 9 are not represented in the translational sector.Previous analyses tell us that their relevant commutators will be with their own secondaries [60].Accordingly, we turn to the velocities, which are also provided in the above output.Unlike the if-constraints, the linearised sSFC Hamiltonian constraints are not automatically implemented when "ToShell"→True is passed: we therefore note that since  ♭ ≈  ♭ ⟂ ≈ 0, from the output of DefTheory[] above, both rotational and translational (Hamiltonian) multipliers for the 0 + and 1 − sectors will also vanish at linear order.The same is true of the translational 1 + and 2 + Hamiltonian multipliers, but the rotational counterparts will be linearly proportional to Returning to the 'lonely' 0 + and 2 − sectors, we can go right ahead and calculate the expected brackets mentioned above.We first find → ToShell" → True]; As expected, this commutator survives at linear order: one could determine also the 0 − multiplier, but it is enough to notice that both P and its linearised secondary10 P ♭ ≡  P ♭ ≈ 0 also become SC.Moving on, we find The same, then, is true of T and its linearised secondary Launch sub-kernels The parallelisation in HiGGS is not sophisticated, but it allows for pragmatic scaling of the analysis to surveys.A large batch of theories can be broken via SBATCH directives over computation nodes, with a reasonable group size of one theory per core per node.Each node is controlled by a master kernel, which coordinates its portion of the survey through a single call to StudyTheory[].Three phases of work are then delegated to the sub-kernels.First, the constraint shell is partially reconstructed from the literature knowledge of the theory.The commutators between all constraints of all theories are then evaluated in parallel.Finally, the velocities of the constraints are decomposed into ther constituent Poisson brackets: these are evaluated and simplified in parallel before returning to the master kernel for synthesis of the results.number and class (all SC) of the if-constraints.The final d.o.f count is thus i.e. the massless graviton.

F. Parallelisation and HPC
The current version of HiGGS is not only designed for the local or desktop operations conducted in Sections II A to II E. The intended use is the evaluation of quantities useful to the Hamiltonian analysis, in an environment where multiple parallel kernels are available.As mentioned in Section I, there are some generic features of the Hamiltonian analysis which lend themselves very well to parallelisation.
Naïvely, we notice that once the primary constraints have been identified, the resulting constraint chains (i.e. the recursive consistency conditions of each primary) can notionally be evaluated in parallel.This seems a natural route down which to parallelise, but in reality chains might seem to continue indefinitely in a kernel with knowledge only of the primary constraint shell.Velocities are the most expensive quantities, and so it seems prudent to pause after each is calculated, to see if the corresponding acceleration is strictly necessary.However, when we parallelise over chains we also find that the velocities, accelerations and jerks rapidly desynchronise: this is to be expected given the varying complexities of the SO(3) irreps which underlie each chain.Accordingly, pre-velocity checks would require e.g. a centralised knowledge of the complete shell to be kept on the master kernel, to which sub-kernels could refer as needed.We could even imagine a setup where one sub-kernel is interrupted in mid-evaluation on the basis of a report from another.
While it is tempting to develop something sophisticated along these lines, the route turns out not to be practical.In practice, the number of chains is typically very few per theory (and strictly not more than ten in the case of the Poincaré gauge theory without multipliers).The desynchronisation effect is then so severe that most kernels would sit idle whilst waiting for the highest-spin chain to complete.More importantly, the question of whether chain need continue based on a new constraint from chain , is not a trivial problem in computer algebra.In the case of the PiCs, each constraint has as one term a unique part of the momentum: as mentioned in Section II E 1 this lends itself to a replacement rule which reliably implements the PiC shell, but in more general cases it is less clear how to proceed.
Ultimately, the experience of [23,45] suggests that evaluation of the accelerations is not usually necessary anyway.Once the brackets between the PiCs are known, it is often possible to determine from a visual inspection whether any chains are worth continuing.Velocities are expensive because of the structure of the Hamiltonian:  ⟂ is a sum of terms which are quadratic in nontrivial covariant quantities such as the ́ and  ́  . Calculation of an overall bracket with  ⟂ thus requires very many sub-brackets.
The quadratic structure of  ⟂ underlies other inefficiencies in the current HiGGS implementation, affecting all terms except for the final − .As discussed already in Section II E 2, the way in which HiGGS calculates velocities is based on (22) rather than (23).As a result, the calculation of a bracket while finding the velocity of  ́ ( 2 ), involving a quadratic term  ́ ( 2 ) ́ ( 2 ) in  ⟂ producesin the first instance -a non-covariant quadratic expression which must be covariantised internally by expensive calls to ToNesterForm[].This could be readily improvable with development, and in many cases we find that the extra penalty incurred by the quadratic covariantisation remains comparable to the cost of the prerequisite brackets.
Since brackets turn out to be the fundamental unit of the analysis, it is over brackets that HiGGS is parallelised.Brackets are cumbersome for humans to evaluate and covariantise: these tasks are relatively easy for computers, now that the infrastructure of Sections II A to II E is in place.On the other hand, shells made from constraints whose format is arbitrary are less easy for computers to use.Once covariantised however, such constraints often comprise straightforward tensor equations which humans can readily manipulate.This suggests a pragmatical division of labour in which HiGGS returns an organised structure of covariantised brackets -formerly the primary Poisson matrix (PPM) [13,14,23] -and velocities for human inspection.If velocities are expensive becuase they require multiple brackets, we can break them up accordingly, along with the quadratic covariantisation steps: in this way we can reduce chain desynchronisation.
The final structure of a HiGGS survey is as follows; 1.A list of theories is passed to StudyTheory[] in the master kernel.This function runs DefTheory[], for each theory, in its own parallel sub-kernel.The literature knowledge of the constraint shell of each theory is cached as a binary.2. Within the same StudyTheory[] call, the master kernel imports the shells and prepares a combined List of brackets which need to be evaluated for the Poisson matrices.The List elements are delegated to all available parallel sub-kernels -as many of which are launched as necessary.Brackets are transferred to the next available sub-kernel using the Wolfram Symbolic Transfer Protocol (WSTP).The WSTP overhead is marginal, since each sub-kernel only needs to transfer data after its PoissonBracket[] call11 .3. Within the same StudyTheory[] call, the master kernel decomposes each linearised velocity into blocks dependent on a single bracket.For the velocity of some PiC ́ the blocks are; The blocks are delegated to the sub-kernels again using WSTP.After the PoissonBracket[] call, each subkernel pre-processes its block using ToNesterForm[] before passing the result back to the master kernel.4. Within the same StudyTheory[] call the blocks are recombined into velocities, and all results are cached in a final binary.
We illustrate this process in Fig. 2. In the parallel HiGGS environment of any node, one has at any one time a collection of theories, each of which has associated with it a complex and growing structure of constraints and commutators.This environment suggests an object-oriented approach.While Mathematica allows for object-oriented programming using e.g.Association[] (see also [91]), the scope of HiGGS is simple enough that we can retain a procedural approach.
In fact, the shell structure of any one theory is referenced fairly infrequently, matching the rate of PoissonBracket[] calls (see Fig. 1 for an illustration of this).This means that while WSTP is necessary for the scheduling of evaluations, the information needed to switch between theories within a sub-kernel can be sourced by using relay system of binary files.Consequently, very minimal use is made by HiGGS of DistributeDefinitions[] to transfer theory-specific data between kernels.In a more serious implementation of the Dirac-Bergmann algorithm, we would envisage more extensive use of WSTP and the object-oriented approach.How efficient is the above approach?Pending an implementation of the Leibniz rule, the efficacy of Fig. 2 as it applies to the whole constraint algorithm is hard to gauge.The full run is only implemented in Section III B for a handful of previously studied theories (without multipliers).Averaged over those cases, and for the non-fundamental reasons outlined above, a minority of the Velocity[] calls contain serial tasks which turn out to dominate the workload.The truly parallel fraction of the workload (including the initial evaluation of all constraint brackets) is then as low as ∼ 10 −2 .Modelling the whole implementation in Fig. 2 as an Amdahl task [92], the -core speedup would appear very limited.We get a fairer understanding, however, from the 'calibration' survey set out in Section III C.
In Table I we run a portion of this survey -spanning only 2 6 = 64 modified gravity theories -on a cluster with a vairable number of cores per node.We take this to be a task of 'reasonable size' when using HiGGS to learn about the canonical structure of a theory.Since the survey does not contain Velocity[] calls, ( ) will be more sensitive to the relevant trade-off between WSTP overhead and bracket evaluation.
There are relatively few brackets, and this portion of the survey is now dominated by the (fast but serial) BuildHiGGS[] and DefTheory[] calls: the benchmarking is thus quick to run, though it still gives an impression of low efficiency due to (39).The normalising time is based on five cores per node: roughly equal to the per-node number of theories.We expect ( ) ≡ ⟨ (5)⟩∕⟨ ( )⟩ to be concave and sublinear, but we do not see it peak in our use-case.No Hamiltonian analysis tools have previously been made that could provide a more meaningful benchmark 12 , so rather than obtaining more compre- hensive statistics we only want to observe here that, when the analysis is reduced to its constituent Poisson brackets, there is a clear benefit from parallelisation.

III. EXAMPLES
Having introduced the implementation in Section II, we now give some basic examples.These will include a reproduction of the analysis in [23,45] and a simple HPC survey which extends those same 'minimal' theories with the use of multipliers.

A. Preparing a science session
To prepare a science session, we begin by loading the package into a fresh Mathematica kernel In[ ]:= <<xAct'HiGGS'; If the HiGGS sources have been correctly placed with respect to the xAct directory tree, a copyright greeter should be displayed, indicating that the context xAct'HiGGS' and its dependancy contexts xAct'xTensor', xAct'xPerm', xAct'xCore', xAct'xTras' have been loaded.
However, loading the package does not yet introduce the physics: the kernel is still in a fresh xAct session, without even a differential manifold.To construct the very many physical definitions needed for science we must build the package using the command This begins an execution in xAct'HiGGS'Private' of xAct/HiGGS/HiGGS_sources.m, most of which is taken up with symbol definitions, and in particular calls to DefTensor[] and MakeRule[].To shorten the process, the most expensive definitions (including those for SO + (1, 3) projection operators) have been stored in binary files under xAct/HiGGS/bin/build/*.mx -those binaries are im-ported at this time 13 .In an active front end, the output cells displaying the progress of this build process are periodically deleted, and should finally be replaced by the following message (adjusting for memory) Out []= ** BuildHiGGS: The HiGGS environment → is now ready to use and is occupying → 63184600 bytes in RAM.The context xAct'HiGGS' should now be populated with many physical quantities, and as a result the session alone is very memory-intensive.We can view some of these quantities through the xAct variables $Tensors and $ConstantSymbols.From this point we may perform the various operations in Section II, or proceed directly to science.

B. Yo-Nester unit tests
As a first application, and to verify that HiGGS has been correctly calibrated, we recapitulate the historical analysis in [23,45].
The theory in which the 0 + mode is active, as considered in [45], is given by imposing the following constraints on ( 7) All the constraint couplings { ̄ } and { ̄ } are of course assumed to vanish, since they are only recently proposed in [14,64].No other 'simple' linear identities are assumed among the couplings, and this avoids collision with other theories.The other theory in [45], in which the 0 − mode is instead allowed to propagate, is defined by In [23] the analysis was extended to 'higher-spin' modes.The theory with only the 1 + mode propagating is The theory with only the 1 − mode propagating is and the theory with only the 2 − mode propagating is There is also tested in [23] a pair of more complex theories in which the 0 − and 2 − modes are both active at the same time.These are given respectively by and The seven configurations Eqs. ( 40) to ( 46) can of course be processed in parallel.We set up a list called JobsBatch, which stores the coupling conditions and assigns a string label to each theory In[ ]:= JobsBatch={}; In[ ]:= JobsBatch~AppendTo~{"spin_0p", {Alp1 == → 0, Alp2 == 0, Alp3 == 0, Alp4 == 0, Alp5 → == 0, 2Bet1 + Bet2 == 0, Bet1 + 2 Bet3 → == 0, cAlp1 == 0, cAlp2 == 0, cAlp3 == → 0, cAlp4 == 0, cAlp5 == 0, cAlp6 == 0, → cBet1 == 0, cBet2 == 0, cBet3 == 0}}; Note that definition in terms of equalities is sufficient: the many strict inequations fixing the absence of other special linear conditions are always assumed by HiGGS to be implicit, and signs associated with other inequalities do not affect the constraint structure we wish to probe (though they may well affect the unitarity).
Assuming we have similarly entered the parameters for all theories above, the major undertaking of evaluating all seven PPMs can be initiated with the command In[ ]:= JobsBatch~StudyTheory~("Import"→True); The option "Import"→True assumes that DefTheory[] has already been run on all seven cases.In that case binaries such as svy/spin_0p.thr.mxetc., will have been created to store a summary of our prior knowledge of each constraint chain.
Even in the case of the small batch in JobsBatch, the calculation is barely viable using sub-HPC resources.We perform the run on a dedicated data processing server with 129 GB memory and eight available 2.90 GHz Intel ® Xeon ® E5-2690 0 CPU cores, corresponding to a maximum of eight parallel Mathematica kernels.This computation lasts ∼ 14 h, but as mentioned in Section II F almost all of this time is spent on the inefficient evaluation of velocities.
Once the calculations of the StudyTheory[] call are complete, we recall from Section II E 3 that the results are displayed in human-readable form via the ViewTheory[] command.We will not provide full analysis here, but focus on the exemplar 1 + case.Skipping the breakdown of PiCs and neglecting the velocities, all the nonlinear brackets are In[ ]:= ViewTheory["simple__spin_1p", " → Literature" → False, "Velocities" → → False]; Out []= ** DefTheory: Incorporating the binary at svy → /simple__spin_1p.thr.mx This information (which by itself is obtained within the first minutes of the StudyTheory[] run) encodes the whole nonlinear primary Poisson matrix (PPM) of the theory: it may be compared with the expressions carefully obtained in [23].
We see that the brackets } are strictly nonlinear: they vanish in the linear theory for which the 1 + mode is moving.The effect of this discontinuity in the constraint structure, when moving from the linear to the nonlinear theory, is well decribed in [23].In short, a counting analogous to that performed in (37) shows that the three extra d.o.f of the massive 1 − mode are nonlinearly activated: HiGGS thus tells us that a vector torsion mode is strongly coupled.
The results of the remaining cases in JobsBatch are provided in the supplemental materials [72].We find that, up to terms which vanish due to multi-term symmetries (i.e. which cannot be eliminated in the xAct architecture), the nonlinear brackets agree with those found in [23,45].There is a possible exception in the case of the general bracket { ∧ , ⟂ }, which we find to be 'surficial' in ways illustrated already by the bracket in Section II D 3. This subtlety appears only to touch the scalar mode results in [45]: it seems worthwhile to investigate -in future work -if and how this might affect the final d.o.f counting.For brevity we will omit the velocities here 14 , since even the covariantised, linear expressions can be very cumbersome.Examples are provided in Appendix A for the minimal 1 + case above: these illustrate how gradients of the field strengths can arise as the algorithm progresses, necessitating a future extension of HiGGS to the second-order Euler-Lagrange formalism.

C. Simple HPC survey
The purpose of this section is to further extend our knowledge of the theories examined in Section III B, by running HiGGS on basic HPC resources, so as to chart the effects of introducing multiplier fields.The findings of this survey, and any viable multiplier configurations, will be presented in future work.The starting point will be the 'minimal' Poincaré gauge theories set out in Eqs.(40) to (46).Of these, we know that (40) and ( 41) are the traditional cases in which massive 0 + and 0 − scalars propagate safely alongside the usual 2 + graviton.Of the cases with active higher-spin modes, we know that (45) and (46) propagate both 0 − and 2 − modes, and may be safe 15 .
The minimal 'problem' cases are Eqs.(42) to (44).These are ostensibly simple linearised theories which propagate extra 1 + , 1 − and 2 − modes respectively, but which appear to suffer from mode activation in the fully nonliear regime.This mode activation is identified by counting the canonical degrees of freedom in the Hamiltonian analysis: its association with a 14 The few velocities which are known from [45] corroborate our results, however we note that ToNesterForm[] actually fails to fully covariantise the velocity of ⟂ in the simple 0 − case.This 'bug' is not known to occur elsewhere, and certainly not during the evaluation of nonlinear brackets: the main feature of HiGGS.We suggest that a fully nonlinear 'Leibniz rule' implementation of velocities, as recommended above, would greatly reduce the risk of such effects by demanding less of ToNesterForm[] in each call.particular sector is not shown explicitly, but inferred by examining the unconstrained PiC functions in each case.These are as follows; • In the minimal 1 + theory ( 42 are unconstrained, the 2 + mode is thought to be activated.In [64] a potential avenue was outlined for preventing mode activation by means of the multiplier fields set out in (7).For any conventional Poincaré gauge theory, there are 2 6 × 2 3 possible multiplier configurations.Any pair of configurations differs, if not by the constrained sectors, then by the 'singular' or 'parallel' nature of a constrained sector.These differences may have quite unpredictable consequences in the analysis.In order to efficiently explore the space of multiplier configurations, it is therefore important to have all the commutators between the known constraints to hand.

Scope of survey
By adding various multipliers, there are 3 × 2 6 × 2 3 = 1536 separate theories which can be constructed from the minimal PGTs.However, some of these can be ruled out immediately on phenomenological grounds.In [64] we focussed on the modification of ( 42) by allowing ̄ 2 ≠ 0. In the Lagrangian picture, this maps to a pair of constraints on the 1 − part of the torsion tensor, suppressing velocities or equating them to gradients.The main phenomenological constraint on the minimal extensions is the requirement that they contain the Einstein-Cartan dynamics, for which the gravitational field still is described by the Riemann-Cartan curvature.For this reason, we suspect that it will generally be safer to impose the { ̄ } than the { ̄ }.How dangerous is it to suppress parts of the rotational sector?To get an idea, we consider how the two d.o.fs in Einstein-Cartan gravity might be manifest in a gravitational wave-type solution.Rather than studying the Einstein-Cartan waves directly, we instead introduce novel solutions to the special (purely quadratic, i.e. ̂ 0 = 0) PGT from [13,14,74,75] which describe null pp-waves on the Minkowski background [94].The wave solutions are formulated in the Brinkmann gauge, rather than the more popular transversetraceless (TT) setup [95].To introduce the Brinkmann gauge we start with a Cartesian coordinate system, and the rotation gauge is first chosen as in the HiGGS environment, so that the local Lorentz and coordinate bases are aligned.We then define 'perpendicular' and 'null' vectors as We will denote the 'wave coordinate' as ≡ − ; the wave amplitude is taken in all solutions to be a smooth and compact scalar function  ≡ ( ).The Brinkmann gauge is complementary to the TT gauge in the sense that it confines waves to the time and longitudinal components of the metric perturbation.The wave will therefore alter the definition of the unit timelike vector , which defines the foliation 16 , and unit spacelike vector which defines the direction of travel, but not the polarisation vector .Restricting to the case of weak waves, we will then have and ≡ cos( )( ) + sin( )( ) .While the Brinkmann gauge is opposed to the TT gauge at the level of the metric, this turns out to be somewhat reversed at the level of the field strengths.Accordingly, it is useful to define the 'TT symmetric-traceless' operation on the indices of a general TT tensor as These components turn out to be identical [94] to those of the Riemann tensor in the presence of the vacuum pp-waves known from GR. Recall that Eqs. ( 42) to (44) were initially set up as a modifications of Einstein-Cartan theory.Since the Einstein-Cartan theory differs from GR only by a contact torsion interaction, and since the linear 1 + , 1 − and 2 + modes are massive, it would seem strange if the the null pp-wave solution Eqs. (48a) to (48d) is not also mandatory in the minimal extensions.Referring back to [64], we find that we should then always take ̄ 1 = 0, since we would otherwise encounter the Lagrangian constraints which both force the wave amplitude  → 0. This already halves the volume of our parameter space, and so we make no attempt exhaustively determine the various other phenomenological constraints.Other limitations on the allowed multipliers come from the linearised particle spectrum.For the 1 + theory with ̄ 1 = 0 only the group { ̄ 3 , ̄ 4 , ̄ 6 , ̄ 1 , ̄ 2 , ̄ 3 } does not immediately constrain ∧ ̂ ⟂ .Similarly for the 1 − theory the space 16 For this reason the Brinkmann gauge will not be a natural choice for the canonical analysis. is For the 2 − theory, we consider the group { ̄ 3 , ̄ 5 , ̄ 6 , ̄ 1 , ̄ 2 , ̄ 3 }.Rightly, we ought to allow ̄ 4 ≠ 0 in this case, especially since the momentum ∼ ̂ ⟂ of the (anticipated) strongly coupled 2 + mode would then be disabled: just for this initial survey, however, the complexity of the 2 − calculations is such that restricting to ̄ 4 = 0 results in a significant economy.
The aim is then to obtain covariant expressions for all the possible commutators among all known if-constraints for all 3×2 3 ×2 3 = 192 theories stipulated above.The current HiGGS setup is supposed to be able to do this, and produce binaries of the results suitable for use in a database.The requisite calculations would take years on a single desktop computer core, so we distribute them over 14 nodes of the Peta4 supercomputerthe CPU cluster component of the heterogeneous CSD3 facility.Each node has two 2.60 GHz Intel ® Xeon ® Skylake 6142 CPUs, each having 16 cores with 6 GB of memory apiece, amounting to 448 processors.As described in Section II F, HiGGS runs on a master kernel within each node.The master kernel delegates the analysis of a (randomly allocated) batch of theories, and has a total of 32 parallel sub-kernels at its disposal.An illustration of the survey was shown already in Fig. 1.The node count of 14 is a service-level restriction on Peta4 rather than limitation of the implementation.We can see from a visual inspection that not all theories are equally expensive: some multiplier configurations engender more constraints and more brackets, while others more heavily involve the higherspin sectors.An examination of the stack trace from each node suggests that the most time-consuming cases are those for which a multiplier is used to constrain a sector which was thought to be nonlinearly activated.This is an interesting observation in the context of the strong coupling considerations.In particular, delays occur when new 'parallel' or 'singular' SiCs of the form are introduced.It is known from manual calculations in [14] that these secondaries (specifically their field strength terms) fail to commute with many other constraints, and can produce brackets of surprising complexity.Physics binaries (HiGGS extension *.thr.mxfiles) of all nonlinear brackets identified in this survey, along with plaintext stack traces and clocking times, can be found in the supplemental materials [72].

Results and prospects for new physics
In this final section, we make some preliminary observations on the new data provided by the 'calibration' survey [72] -the bulk of this analysis being reserved for future work.We focus on the 1 + and 1 − modes.
The investigation in [23] seems to identify the source of strong coupling as follows.In the linear theory without multipliers, rotational PiCs ́ may be paired off with conjugate PiCs ́ , so as to determine ́ and ́ and so terminate both chains.Unpaired PiCs generate SiCs, whose consistency conditions fix the original PiC Hamiltonian multipliers -all within the same sector.As the theory becomes nonlinear, the PiC structure does not, of course, change.However there are generally more PiC-PiC commutators which emerge, and which do not respect the law of conjugate pairs.Consider a theory where the collective { ́ } span d.o.f, and the { ́ } span < d.o.f.Then the emergence of nonlinear commutators can in principle alter the PPM rank so that the rotational consistencies determine up to rotational Hamiltonian multipliers { ́ }, overflowing into − SiCs (denoted [ − ] in [23]).The translational PiCs and − SiCs generically fail to commute (again due to nonlinearity) with the rotational PiCs, whose Hamiltonian multipliers { ́ } they then determine.
The (reduced) number − of SiCs dictates the number of strongly coupled modes in the d.o.f counting.Being unconstrained by the conjugacy, the number − may not map on to the 2 + 1 integer spin multiplicities, or it may seem to contribute a half-integer d.o.f.These matters should then be clarified by closer study, and identification of FC combinations to restore integer d.o.f.Critically, the division of translational and rotational PiCs is efficacious because the rotational sector commutes with itself nonlinearly.
Geometric multipliers contribute new primaries by direct analogy to (12) i.e. the multiplier momenta and (defined as in (11) to be conjugate to and ), which are also resistant to nonlinear commutators.Irreps of the fields and appear in other if-constraints, but these produce predictable, linear commutators with the conjugate sectors of the and .For the case of translational multipliers (as discussed in Section III C 1, we suspect these to be safer), an approach might be to use the consistencies of the constrained irreps to solve for a large number of { ́ } in the linear theory.This would force the rotational sector into developing the maximum number of SiCs ab initio.Assuming nonlinear commutators between the SiCs and rotational PiCs proliferate as usual, one then expects to solve for all the remaining Hamiltonian multipliers as before, with a generally SC system whose SiCs are not contingent on nonlinear effects.
The 'calibration' survey covers the simple extension of (42) by the condition ̄ 2 ≠ 0, and we used the resultant brack-ets when exploring this mechanism in [64].That theory does not appear to be strongly coupled, but it is also unlikely to be unitary.We close this section by suggesting another option.
In [23] the basic PGT with only ̂ 0 ≠ 0 and ̂ 5 ≠ 0 was briefly considered, in which both 1 + and 1 − modes were strongly coupled.This theory contains all the translational PiCs17 : , 1 ≠ 0, we can obtain all the translational Hamiltonian multipliers except for in the 0 + sectorwhich remains under traditional conjugacy.In this way, an extra ∼ ⟂ is forced in both regimes, and its natural conjugacy (comparing to the experience of the 0 + in [64]) will be where solid arrows indicate that a Hamiltonian multiplier is determined, and dashed arrows indicate that a secondary must be constructed.Four-step consistency chains such as in (52) are not seen in the original PGT.If the overall SC structure is indeed preserved in the nonlinear theory, it would seem that only two d.o.f overall propagate.This result seems suspicious, and must be very carefully tested -for example with HiGGS.If it is true, then the resulting theory might conceivably introduce extra contact interactions to the Einstein-Cartan model, and the overall phenomenological differences with GR would be very interesting to study.

IV. CONCLUSIONS
In this paper we have presented the package HiGGS, written for the tensor manipulation suite xAct and the computer alge-bra software Mathematica.The HiGGS package performs calculations -such as Poisson brackets -which frequently arise in the Hamiltonian (canonical) analysis of modified gravity theories with curvature and/or torsion.The current iteration of the package is tailored to the generalised Poincaré gauge theory in Eq. ( 7): given an action of that class, HiGGS can identify primary and secondary constraints, and calculate arbitrary field velocities.For more original torsionful actions, HiGGS may still be used to evaluate brackets, and canonicalise expressions by irreducible decomposition.
Parallelisation is a core feature of HiGGS.We have argued that those aspects of the Dirac-Bergmann Hamiltonian constraint algorithm which can become cumbersome during manual evaluation, are actually very well suited to parallel computing.On a laptop or desktop computer, HiGGS can take advantage of available cores to shorten the analysis of a given theory.Parallelisation is done within the HiGGS environment.This capability meshes well with the problem of surveying large numbers of action configurations using high-performance computing (HPC), since it avoids the need to educate general-purpose job scheduling tools (such as SLURM [96] or TORQUE [97]) about the physical details of the Hamiltonian constraint structure.In our example HPC survey, all job-scheduling decisions were made from within the HiGGS environment, to whole-node granularity.
The HPC survey we have performed here -whose results will be discussed further in future work -targets minimal extensions to the Einstein-Cartan theory in which a single extra massive spin-parity 1 + , 1 − or 2 − torsion particle is present, i.e. extending the two usual graviton polarisations by three or five extra degrees of freedom (d.o.f).These minimal extensions were believed [98][99][100] to be ghost free according to the linearised analysis.Subsequent nonlinear Hamiltonian analysis [23] suggested that the linear regime strongly couples extra modes of spin-parity 1 − , 1 + or 2 + respectively, so that these modes spoil the theories' viability.Based on the hypothesis that multiplier fields could selectively suppress these modes -in the style of teleparallel gravity -our survey lists all the commutators among all the known primary and secondary 'ifconstraints' (i.e.constraints which may arise due to a choice of Lagrangian couplings), for the various possible multiplier configurations.Our level of analysis at least matches that of mode activation 18 in [23], where the primary Poisson matrices of the minimal Einstein-Cartan extensions are obtained.We note that some of these results (which are available in the supplemental materials [72]) have already been used in [64].
While the HiGGS package may be of use to researchers in the ways described above, we must observe some of its many limitations; • The PoissonBracket[] module automatically expands its operands, rather than first taking advantage of the Leibniz rule wherever those operands are products of covariant factors.This can result in costly at- 18 We do not, however, perform the extra step in [23] of considering the problem of constraint bifurcation using the Poisson matrix pseudodeterminant.
tempts by ToNesterForm[] to covariantise arbitrarily complex brackets.The lack (to our knowledge) of a general formula in Poincaré gauge theory for the second-order bracket was also a limiting factor for our previous analysis in [13].• In general, HiGGS relies quite heavily on the gaugecovariant derivatives or  .The xAct suite already has a very sophisticated functionality to accommodate the definition of such derivatives, which is not exploited by the HiGGS implementation.
• In general, HiGGS relies very heavily on the SO + (1,3) and (particularly) the SO(3) decompositions of tensorial fields.However, these are manually defined in the implementation, for each decomposed field.A clear case is made for a general decomposition tool, more closely integrated with the existing functionality in xAct'SymManipulator'.

• Notwithstanding the Lagrangian structure implied by
Velocity[], inclusion of new dynamical variables constitutes a different problem.Extension to the metric affine gauge theory (MAGT) structure would be a straightforward, if time-consuming, exercise.More work would likely be needed for the inclusion of (e.g.fermionic) matter fields.
• The background assumed by HiGGS when the option "ToOrder"→0 or "ToOrder"→1 is passed to modules such as ToNesterForm[], is Minkowski spacetime with vanishing background torsion.There are, however, obvious motivations for considering curved backgrounds, such as de Sitter or Schwarzschild.Moreover in [13,14,74,75] it is argued that non-minimal gravitational gauge theories, which are detatched from the geometrical trinity, only become viable in an attractor background of constant axial torsion.Whether these limitations are best addressed by improving the HiGGS implementation, or beginning ab initio with an improved understanding of the challenges posed by canonical computer algebra, remains to be seen.For the moment, we hope at least to have shown that modified gravity is now ready for computer algebra assistance at scale.In this sense we build on the recent work of Lin, Hobson and Lasenby [81,101], who used computer algebra to systematically obtain all ghost and tachyon-free cases of the linearised Poincaré gauge theory, and did so with far more limited resources than are brought to bear in this paper.A future is then suggested in which vague concerns -viz, a potential for strong coupling demonstrated among a handful of cases -are no longer valid grounds upon which to dismiss a rich class of theories.These expressions are expressed on the PiC shell, but not on the sSFC shell.We note in particular the general presence of momentum gradients, which arise at the end of the HiGGS run.
In non-minimal theories (i.e.those with PiCs which depend on the field strengths), we can expect gradients of the field strengths to also appear.These are second derivative quantities, even in the first-order formulation of gravity that is PGT: they cannot be further processed by HiGGS, which uses a order Euler-Lagrange implementation.

⟂
svy/EinsteinCartan.thr.mx.There are four nonvanishing Poisson brackets, all between the translational and rotational pairs of the same SO(3) irrep, vs ⟂ , .These are the conjugate pairs, whose commutators are are proportional to mass parameters in the theory[60] -here mediated by the Einstein-Hilbert term ̂ 0 .The involved PiCs are second class (SC).
and so the algorithm terminates.What can HiGGS tell us about the physics of Einstein-Cartan theory?The PGT contains 2 × (24 + 16) = 80 naïve d.o.fs in its gauge fields.The non-physicality of the lapse and shift (the Poincaré gauge symmetry), remove 2 × 10 d.o.fs through the sure, primary (sPFC) constraints; a further 20 d.o.f are removed via the sSFCs.We further learn from HiGGS about the

1 − 1 +
mode is thought to be activated.•In the minimal 1 − theory (43), , mode is thought to be activated.•In the minimal 2 − theory (44), ∼ ⟂ and T solution which concerns us describes a wave in the Riemann-Cartan curvature, with vanishing torsion and two d.o.fs quantified by a polarisation vector.It has the components

.
The rotational PiCs are only conjugate in the 0 + and 2 + sectors, so the linear theory produces SiCs ∧ , nonlinear case, all 12 d.o.f in the rotational and translational sectors are assumed to solve exactly for each others' Hamiltonian multipliers: no SiCs are produced and the total d.o.f rises by two massive vectors 1 2 (1+5+3+3) = 3+3.However by imposing ̄ develops in the 2 − sector.Denoting consistency conditions with arrows, we arrive at We now move away from this geometric picture, to a 'particle physics' setup where the underlying manifold is always flat Minkowski space .The metric associated with (curvillinear) coorinate tangents is then ≡ ⋅ , and this metric is flat.The coordinate basis is accompanied by a Lorentz basis, whose dot products give the Minkowskian metric components ≡ ̂ ⋅ ̂ .The Lorentz basis can rotate locally under the proper, orthochronous Lorentz rotations, and is non-holonomic.Following on from (3), a vector  with Lorentz indices has a covariant derivative 1 2 − ( | | ) = 0, but still conveys torsion through the contorsion d.o.fs ≡ 1 2 + ( | | ) .Aside from the contorsion, the metric d.o.fs which source the Levi-Civita part are defined, as usual, by tangents to the coordinate curves ≡ ⋅ .

TABLE I .
(42)oximate benchmarks for parallel brackets.The spinparity 1 + theory in(42)is augmented with all configurations of the multipliers { ̄ 3 , ̄ 4 , ̄ 6 , ̄ 1 , ̄ 2 , ̄ 3 } -see (7) -and all the Poisson brackets are obtained.Either four or five theories are allocated per node, wallclock time ⟨ ( )⟩ is averaged over all 14 nodes, but apparent remains uniformly high (i.e.suspicious) due to a paucity of expensive brackets.For survey structure and specifications of the Peta4 cluster, see Section III C.
(7))e Velocity[] module is directly associated with the canonical Hamiltonian in(35)of the generalised Poincaré gauge theory in(7).This could easily be avoided (i.e. in favour of user-defined Hamiltonia) with Leibniz rule functionality in PoissonBracket[].•The PoissonBracket[] module is incapable of processing the second-order Euler-Lagrange formulation.