# Comment on
“Capturing Phase Behavior of Ternary Lipid Mixtures with a Refined Martini
Coarse-Grained Force Field”

Matti Javanainen<EMAIL_ADDRESS>Institute of Organic Chemistry
and Biochemistry, Czech Academy of Sciences, Flemingovo nám. 542/2, 160 00
Prague 6, Czech Republic  Balazs Fabian  Institute of Organic Chemistry and
Biochemistry, Czech Academy of Sciences, Flemingovo nám. 542/2, 160 00 Prague
6, Czech Republic  Hector Martinez-Seara  Institute of Organic Chemistry and
Biochemistry, Czech Academy of Sciences, Flemingovo nám. 542/2, 160 00 Prague
6, Czech Republic

(August 27, 2024)

###### Abstract

We report here on the pitfalls of the simulation model introduced in the
“Capturing Phase Behavior of Ternary Lipid Mixtures with a Refined Martini
Coarse-Grained Force Field” [Journal of Chemical Theory and Computation 2018,
14, 11, 6050–6062]. This refined Martini model was reported to reproduce
experimental phase diagrams for a ternary DOPC/DPPC/cholesterol mixture,
including the coexistence of two liquid phases. However, we demonstrate that
this coexistence only emerged due to an unfortunate choice of simulation
parameters, which leads to poor energy conservation. Specifically, the
constraints on the cholesterol model drained energy out from the membrane,
resulting in two coexisting phases at drastically different temperatures.
Using the simulation parameters recommended for the used cholesterol model,
this artefact is eliminated, yet so is phase coexistence, i.e. experimental
phase diagrams are no longer reproduced.

It is important to highlight that the present comment was submitted to
Chemical Theory and Computation. However, it was rejected without peer-review
by the Editor-in-Chief, who stated that the journal “rarely publishes such
material”.

At room temperature and at certain ratios, the canonical mixture of
dipalmitoylphosphatidylcholine (DPPC) with two saturated acyl chains,
dioleoylphosphatidylcholine (DOPC) with two monounsaturated acyl chains, and
cholesterol (CHOL) undergoes phase separation into liquid ordered (Lo) and
liquid disordered (Ld) phases Veatch and Keller (2003). It is well-documented
that the standard implementation of the coarse-grained (CG) Martini model
Marrink  _et al._ (2007) does not capture this behavior Carpenter  _et al._
(2018); Domański  _et al._ (2012); Davis  _et al._ (2013). Therefore, phase
separation and phase coexistence studies performed using the Martini model
replace DOPC by dilinoleoylphosphatidylcholine (DLiPC) with two diunsaturated
acyl chains as low transition temperature ($T_{\mathrm{m}}$) lipid Risselada
and Marrink (2008). Unfortunately, no experimental phase diagrams for the
DPPC/DLiPC/CHOL mixture exist to our knowledge, rendering the comparison
between simulations and experiments qualitative at best.

Recently, Carpenter  _et al._ tackled this outstanding issue by a careful
refinement of the Martini parameters for DOPC and DPPC. They adjusted the
bonded parameters of each bead separately to reproduce the bond length and
angle distributions extracted from atomistic simulation data Carpenter  _et
al._ (2018). This approach diverts from the building block concept of the
version 2 family of the Martini force field, where the number of bonded
parameters is minimized. With this extended freedom to fine-tuning of the
interactions, Carpenter  _et al._ achieved almost quantitative agreement with
the experimental phase diagrams for the mixture of DPPC, DOPC, and CHOL Veatch
and Keller (2003); Davis  _et al._ (2009). Considering the known limitations
of the CG models Alessandri  _et al._ (2019), this is an impressive milestone
on the way to accurate modeling of complex lipid mixtures. In their work,
Carpenter  _et al._ derived two different models; “Extensible parameters”
where DPPC and DOPC head groups share identical interaction parameters, and
“Optimal parametes”, where they are allowed to differ. The latter parameter
set is of more interest, as it provides the best agreement with the
experimental phase diagram Carpenter  _et al._ (2018).

Here, we demonstrate that the reported liquid–liquid coexistence in the model
by Carpenter  _et al._ (“Optimal parameters”) is an artifact caused by an
unfortunate choice of simulation parameters for GROMACS. Namely, the used
current cholesterol model by Melo  _et al._ Melo  _et al._ (2015) was
parametrized using a very conservative set of parameters for the LINCS
constraint algorithm to properly model the ring structure of cholesterol.
Unfortunately, these more conservative constraint options were not followed by
Carpenter  _et al._ Carpenter  _et al._ (2018), who instead followed the
parameter set generally recommended for the Martini model. When the
cholesterol model is simulated with these generally used parameters, i.e. 4th
order expansion of the constraint coupling matrix (lincs_order equal to 4) and
1 iteration of constraints per integration time step (lincs_iter equal to 1),
energy is not properly conserved due to the presence of virtual sites and
corresponding coupled constraints in this model. This effect, which is also
highlighted in theGROMACS manual Abraham  _et al._ (2015), manifests itself
especially at large integration time steps.

In the case of the study by Carpenter  _et al._ Carpenter  _et al._ (2018),
CHOL constraints drain out a substantial amount of energy from their
neighborhood, which in the case of the studied ternary lipid mixtures consists
preferentially of DPPC. Thus, this leads to the cooling down of CHOL and the
surrounding DPPC lipids which together form a sort of ordered phase. In
contrast — to maintain the target temperature of the thermostat — DOPC-rich
region gets warmer than the target temperature of the thermostat, thus forming
a very disorered phase. This Maxwell demon causes the temperatures of DPPC
(and CHOL) and DOPC to diverge exponentially as a function of the integration
time step from the target temperature of the thermostat, therefore resulting
in the apparent Lo phase being dozens of K cooler than the Ld phase.
Importantly, Lo/Ld coexistence regions reported by the authors and matching
the experimental phase diagram in the DPPC/DOPC/CHOL mixtures are only
recovered in simulations in the presence of this artifact.

We show that this build-up of temperature difference can be at least partially
prevented by either 1) following the suggested LINCS parameters for CHOL
simulations by Melo  _et al._ , 2) coupling the different lipid types to
separate thermostats, or 3) decreasing the simulation time step. All these
options, that improve energy conservation, cause the Lo/Ld phase coexistence
to vanish. Instead, an almost ideal mixing is instead observed.

Thus, when achieving a proper NPT ensemble, the model by Carpenter  _et al._
does not present an improvement over the standard Martini v2.2 lipid model in
capturing the phase coexistence in DOPC/DPPC/CHOL mixtures Davis  _et al._
(2013); Risselada and Marrink (2008); Carpenter  _et al._ (2018). Therefore,
further work is required before a direct comparison of experimental and
Martini-based coarse grained simulations on phase-separating lipid mixtures is
viable. Meanwhile, the DPPC/DLiPC/CHOL mixture in the standard Martini
implementation will—despite its obvious limitations—serve as the main model of
choice for studies on Lo/Ld phase coexistence.

##  I Results

To demonstrate the issues in the model by Carpenter  _et al._ , we first
followed the original simulation protocol reported in Ref. 3. We performed 15
µs simulations of 3/3/2 and 2/1/1 compositions of the DPPC/DOPC/CHOL mixture
(see Methods and Set 1 in Table 1). Both of these compositions fall in the
heart of the Lo/Ld coexistence region Carpenter  _et al._ (2018). The
simulations were performed at 298 K using time steps of 10, 15, 20, 25, 30,
35, and 40 fs. Importantly, 35 fs was used in the original work Carpenter  _et
al._ (2018). We characterized the degree of phase separation using the mean
contact fraction Domański  _et al._ (2012) ($f_{\mathrm{mix}}$, see Methods)
extracted from the last 1 µs simulation. For each value of the integration
time step, we extracted $f_{\mathrm{mix}}$ and the temperature of the lipids.
With blue markers in Fig. 1A), we demonstrate the obvious correlation of
$f_{\mathrm{mix}}$ of the function of lipid temperature, both of which depend
on the used integration time step. Notably, larger time steps lead to
significantly cooler membranes ($<$ 298 K) and smaller values of
$f_{\mathrm{mix}}$, i.e. higher degree of phase separation. With time steps
that maintain the overall temperature of the lipids reasonably close to the
target temperature of 298 K, both mixtures display nearly ideal mixing
instead. The fact that the target temperature of the thermostat is not
preserved even with small time steps is in line with a recent report on the
imperfect energy conservation with the New-RF simulation settings used
Benedetti and Loison (2018).

Figure 1: LINCS settings and temperature coupling affect the degree of phase
separation. A) The degree of phase separation as a function of the simulation
temperature of the lipids, both of which depend on the integration time step,
indicated by point labels. Only the labels that are multiples of 10 are shown
for clarity. The dotted gray line shows the target temperature of the
thermostat. Data for the two studied mixtures are shown in the two panels, and
the $f_{\mathrm{mix}}$ values corresponding to ideal mixing are highlighted by
dashed gray lines. The blue markers show data for simulations in which all
lipids are coupled Together to the thermostat (Set 1 in Table 1). Green
markers show data for simulations in which lipid types are coupled to Separate
thermostats (Set 2 in Table 1). Red markers show data for simulations in which
more conservative LINCS parameters are used (Set 3 in Table 1). B)
Temperatures of the lipid types and the solvent beads as a function of the
integration time step. Data are only shown for the 3/3/2 mixture, yet the
behavior in the 2/1/1 mixture is essentially identical. C) Temperatures of the
lipid types and the solvent beads as a function of temperature coupling and
LINCS options. Data are again shown for the 3/3/2 mixture only.

We also extracted the temperatures of the different lipid types. As this
information not available in the original trajectories as no velocities were
saved, we extended each simulation by 100 ns, during which we also saved the
instantaneous velocities (see Methods). The temperatures of each lipid type
are shown as a function of the integration time step in Fig. 1B). It is
evident that the increase in time step causes the low-$T_{\mathrm{m}}$ lipid
(DOPC) to heat up and the high-$T_{\mathrm{m}}$ (DPPC) lipid to cool down. The
behaviors of DPPC and CHOL are similar. For time steps up to 20 fs, this
effect is negligible, yet with the 35 fs time step used by Carpenter  _et al._
, the temperature difference between the lipid components is already more than
20 K with DPPC being $\sim$30 K below its $T_{\mathrm{m}}$ of 314 K. Thus, it
is clear that decreasing the integration time step improves energy
conservation, yet leads to a smaller level of phase separation.

To pinpoint the culprit of the poor energy conservation, we performed
additional simulations with the input parameters slightly modified. First, we
used the same integration time step of 35 fs as Carpenter  _et al._ but with
the LINCS parameters provided by Melo  _et al._ Melo  _et al._ (2015)—namely
two LINCS iterations per time step, and an 8th order expansion of the
constraint coupling matrix (Set 3 in Table 1). Data for these settings is
shown in red in Fig. 1A). It is evident that this more accurate handling of
the constraints leads to a much smaller degree of cooling of the lipids. With
these settings, phase separation also essentially vanishes. The temperatures
of each lipid type in the simulations with a time step of 35 fs and with
various simulation settings are shown in Fig. 1C). The more accurate LINCS
settings clearly lead to all the lipids being at almost the same temperature
which also closely matches the target one. This indicates that energy
conservation can be improved by at least these two ways—decreasing the time
step or increasing the number of LINCS iterations and the order of the LINCS
expansion. The fact that “LINCS” and “20” overlap in Fig. 1A) suggests that
the effect on phase separation is independent on the means of improving energy
conservation. However, both approaches lead to a disagreement between the
simulation and the experimental phase diagram, contrary to the results
reported in the paper by Carpenter  _et al._.

We performed an additional simulation, in which we attempted to improve energy
conservation by using a separate thermostat for each of the lipid types while
still using an integration time step of 35 fs (Set 2 in Table 1). As shown in
Fig. 1A) by green dots, these settings lead to an almost as cool membrane as
in the case when all lipids are coupled to the same thermostat. Still, the
degree of phase separation is significantly decreased. Based on Fig. 1C), the
used temperature coupling does not prevent CHOL and DPPC from cooling down,
yet it does not allow DOPC to heat up. Therefore, the average lipid
temperature decreases, whereas the temperature difference between DPPC and
DOPC is minor resulting in no phase separation.

We also plot the final structures of selected simulations on the top row of
Fig. 2. It is evident that proper phase separation takes place with a time
step of 35 fs and following the simulation options of Carpenter  _et al._
Carpenter  _et al._ (2018) (Set 1 in Table 1), whereas decreasing the time
step to 10 fs results in an almost ideal mixing of the lipid types. Similar
behavior is observed with the conservative LINCS parameters (Set 3 in Table
1). However, for the case in which all lipid types are coupled separately to a
thermostat (Set 2 in Table 1), some heterogeneity is still observed. The local
temperature maps, calculated from the 100 ns simulation during which
velocities were saved, are shown on the bottom row of Fig. 2. It is evident
that the observed phase separation with the model by Carpenter  _et al._ goes
hand in hand with the uneven temperature distribution in the membrane,
indicating that the systems is not sampling the proper equilibrium ensemble.
With the simulation setup of Carpenter  _et al._ Carpenter  _et al._ (2018),
the temperatures of the Lo and Ld phases can differ as much as 100 K between
regions.

Figure 2: Top row: final structures of chosen simulations with the 3/3/2
DPPC/DOPC/CHOL mixture after 15 µs. DPPC, DOPC, and CHOL are shown in red,
green, and yellow, respectively. Only the system where all lipids are coupled
to the same thermostat (“Together”) and a 35 fs integration time step is used
undergoes phase separation. Some heterogeneity is also seen in other systems
with DPPC and CHOL clustering together. Bottom row: the temperature maps
calculated from the 100 ns simulations during which velocities are saved (see
Methods). For the “Together, 35 fs” system, the two phases show almost 100 K
difference in temperature. When the lipid types are coupled to separate
thermostats (“Separate”), the DPPC-rich regions are slightly cooler than the
remainder of the membrane. For the other systems, the temperature is close to
the target on (298 K) in the entire membrane.

##  II Conclusions

From the data presented above, we can draw the following conclusions: 1) The
use of LINCS parameters that are incompatible with the used CHOL model,
combined with a large integration time step, results in poor energy
conservation. In the case of DPPC/DOPC/CHOL mixture considered by Carpenter
_et al._ Carpenter  _et al._ (2018), the virtual site model of CHOL Melo  _et
al._ (2015) drains out a significant amount of energy from the system, mainly
from DPPC which locates itself close to CHOL. 2) When all lipids are coupled
together to the thermostat, as was done by Carpenter  _et al._ Carpenter  _et
al._ (2018), the cooling of CHOL and DPPC leads to the heating up of DOPC.
This temperature difference drives the separation of the lipids into hot
(DOPC-rich) and cool (DPPC- and CHOL-rich) phases, whose temperatures differ
by significantly. This behavior is in obvious violation of the second law of
thermodynamics. 3) The poor energy conservation can be cured by either
decreasing the integration time step or using the LINCS parameters that were
used in the parametrization of the virtual site model of CHOL Melo  _et al._
(2015). This better energy conservations leads to the loss of phase
coexistence, and thus to poor agreement with experimental phase diagrams.
Additionally, by coupling all lipid types to separate thermostats, the
temperature difference between DPPC and DOPC is limited, and only a small
degree of lipid demixing is observed.

It is also worth mentioning that the standard implementation of the current
version 2.2 of the Martini model does not display phase separation of the
DPPC/DOPC/CHOL mixture even when the simulation settings of Carpenter  _et
al._ are used (set 4 in Table 1). This results from the fact that the standard
Martini DOPC model mixes fairly well with cholesterol. Indeed, a hybrid
DPPC/DOPC/CHOL mixture with the parameters of Carpenter  _et al._ used for
DPPC and the standard Martini parameters used for DOPC did not phase-separate.
However, with DOPC parameters adapted from Carpenter  _et al._ and DPPC
parameters from the standard Martini implementation, phase separation was
recovered. This indicates that the unfavorable mixing of DOPC and CHOL in the
model by Carpenter  _et al._ is central for the co-cooling of DPPC with CHOL,
which is required for the formation of a cool Lo phase.

All in all, while the model of Carpenter  _et al._ seemed to finally reproduce
the experimental phase diagram of the well-studied DPPC/DOPC/CHOL lipid
mixture, our results demonstrate that it does due to a simulation artifact.
This artifact relates to a poor energy conservation, which results from an
unfortunate combination of large integration time steps and not using the the
LINCS parameters required by the CHOL model. These choices lead to a clear
violation of target NPT ensemble. Unfortunately, this compromises the validity
of any findings obtained with the model by Carpenter  _et al._ Carpenter  _et
al._ (2018), and suggests that further work is still needed until a Martini-
like coarse-grained model reproduces the phase behavior of the canonical
DPPC/DOPC/CHOL mixture.

##  III Methods

###  III.1 Simulations

In the first set of simulations, we strictly followed the protocol of
Carpenter  _et al._ Carpenter  _et al._ (2018) and built lipid membranes with
dimensions of $30\times 30\times 15$ nm2, which had a total of 3040 lipids and
$\sim$77000 solvent beads, out of which 10% were modeled as antifreeze
particles. The systems were set up by the insane tool Wassenaar  _et al._
(2015). We considered two compositions in the heart of the Lo/Ld coexistence
region, namely DPPC/DOPC/CHOL ratios of 3/3/2 and 2/1/1. We applied restraints
with a force constant of 2 kJ/(mol$\times$nm2) to the phosphate beads of the
phospholipids in one leaflet. All in all, we used the same equilibration
protocol and simulation parameters as Carpenter  _et al._ in the commented
paper Carpenter  _et al._ (2018). Namely, the New-RF simulation parameters de
Jong  _et al._ (2016) were used, unless otherwise mentioned. To pinpoint the
violation of the NPT ensemble to the used CHOL model, we performed the
following simulations (see also Table 1):

First, following Carpenter  _et al._ , the lipids and the solvent were
separately coupled to a thermostat with a target temperature of 298 K, the 15
µs simulations with time steps of 10, 15, 20, 25, 30, 35, and 40 fs, saving
the trajectories every 1 ns. Secondly, we performed 15 µs simulations with a
35 fs time step but with all the lipid types (DPPC, DOPC, and CHOL) coupled
separately to a thermostat. Thirdly, we performed 15 µs simulations with two
iterations of the LINCS algorithm (lincs_iter = 2) and using the 8th order
expansion of the constraint coupling matrix (lincs_order = 8). Here, all
lipids were coupled together to the thermostat.

We also run additional 100 ns simulations, starting from the final structures
of the aforementioned simulations, and saved the velocities (.trr file) so
that the temperatures of each lipid type could be extracted using gmx traj. As
the used cholesterol model has constraints, we corrected for the missing
degrees of freedom.

Table 1: Simulated systems. Composition is given as DPPC/DOPC/CHOL. “$\Delta t$” is the integration time step in fs. “$T$-coupl.” refers to the way the lipid temperatures are coupled; the lipids are being coupled either together to one thermostat or separately to three thermostats.  Composition | $\Delta t$ |  $T$-coupl. | lincs_order | lincs_iter
---|---|---|---|---
Set 1. Lipids coupled Together
3/3/2 | 10 | Together | 4 | 1
3/3/2 | 15 | Together | 4 | 1
3/3/2 | 20 | Together | 4 | 1
3/3/2 | 25 | Together | 4 | 1
3/3/2 | 30 | Together | 4 | 1
3/3/2 | 35 | Together | 4 | 1
3/3/2 | 40 | Together | 4 | 1
2/1/1 | 10 | Together | 4 | 1
2/1/1 | 15 | Together | 4 | 1
2/1/1 | 20 | Together | 4 | 1
2/1/1 | 25 | Together | 4 | 1
2/1/1 | 30 | Together | 4 | 1
2/1/1 | 35 | Together | 4 | 1
2/1/1 | 40 | Together | 4 | 1
Set 2. Lipids coupled Separately
3/2/2 | 35 | Separately | 4 | 1
2/1/1 | 35 | Separately | 4 | 1
Set 3. Conservative LINCS settings
3/2/2 | 35 | Together | 8 | 2
2/1/1 | 35 | Together | 8 | 2
Set 4. Unmodified Martini v2.2
3/2/2 | 35 | Together | 4 | 1

###  III.2 Analyses

The contact fraction $f_{\mathrm{mix}}$, describing the level of lipid
demixing, was adapted from Ref. 4, and is defined as

$\displaystyle
f_{\mathrm{mix}}=\frac{c_{\mathrm{US-S}}}{c_{\mathrm{US-S}}+c_{\mathrm{US-
US}}},$ (1)

where for example $c_{\mathrm{US-S}}$ refers to contacts between unsaturated
(US, here DOPC) and saturated (S, here DPPC) lipids. Thus, the smaller the
value of $f_{\mathrm{mix}}$, the sharper the separation is. A contact is
registered if the phosphate beads of the lipids were within 1.1 nm.

The temperatures of the lipids were extracted using gmx energy. When the lipid
types were coupled separately, their temperatures were extracted separately
and the averaging was performed over the degrees of freedom. For CHOL, this
was less than $3N$ due to the constraints. The temperatures of each lipid type
were extracted from the 100 ns trajectories containing velocities (see
Simulations Methods above) using gmx traj. The values were corrected for the
missing degrees of freedom of CHOL.

The heat maps of temperature distribution were calculated by projecting the
lipid center of mass onto the macroscopic plane of the membrane. While
performing the binning, each point was weighted by the instantaneous
temperature of the corresponding molecule. For DPPC and DOPC, the number of
degrees of freedom in each lipid was taken as $3N$, while in the case of CHOL
the decrease of degrees of freedom in the presence of constraints was
accounted for.

###### Acknowledgements.

M.J., H.M.-S., and B.F. acknowledge support from the Czech Science Foundation
(EXPRO grant 19-26854X). M.J. thanks the Emil Aaltonen Foundation for funding.

## References

  * Veatch and Keller (2003) S. L. Veatch and S. L. Keller, Biophys. J. 85, 3074 (2003).
  * Marrink  _et al._ (2007) S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. De Vries, J. Phys. Chem. B 111, 7812 (2007).
  * Carpenter  _et al._ (2018) T. S. Carpenter, C. A. López, C. Neale, C. Montour, H. I. Ingólfsson, F. Di Natale, F. C. Lightstone, and S. Gnanakaran, J. Chem. Theory Comput. 14, 6050 (2018).
  * Domański  _et al._ (2012) J. Domański, S. J. Marrink, and L. V. Schäfer, BBA – Biomembranes 1818, 984 (2012).
  * Davis  _et al._ (2013) R. S. Davis, P. Sunil Kumar, M. M. Sperotto, and M. Laradji, J. Phys. Chem. B 117, 4072 (2013).
  * Risselada and Marrink (2008) H. J. Risselada and S. J. Marrink, Proc. Natl. Acad. Sci. USA 105, 17367 (2008).
  * Davis  _et al._ (2009) J. H. Davis, J. J. Clair, and J. Juhasz, Biophys. J. 96, 521 (2009).
  * Alessandri  _et al._ (2019) R. Alessandri, P. C. Souza, S. Thallmair, M. N. Melo, A. H. De Vries, and S. J. Marrink, J. Chem. Theory Comput. 15, 5448 (2019).
  * Melo  _et al._ (2015) M. Melo, H. Ingólfsson, and S. Marrink, J. Chem. Phys. 143, 12B637_1 (2015).
  * Abraham  _et al._ (2015) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1, 19 (2015).
  * Benedetti and Loison (2018) F. Benedetti and C. Loison, Comput. Phys. Commun. 228, 146 (2018).
  * Wassenaar  _et al._ (2015) T. A. Wassenaar, H. I. Ingólfsson, R. A. Böckmann, D. P. Tieleman, and S. J. Marrink, J. Chem. Theory Comput. 11, 2144 (2015).
  * de Jong  _et al._ (2016) D. H. de Jong, S. Baoukina, H. I. Ingólfsson, and S. J. Marrink, Comput. Phys. Commun. 199, 1 (2016).

##  IV Data Availability

The inputs and outputs for our simulations with the model by Carpenter  _et
al._ are available as follows:

  1. 1.

Long simulations with different time steps:

     * •

DOI: 10.5281/zenodo.3956709 (3/3/2 mixture)

     * •

DOI: 10.5281/zenodo.3956797 (2/1/1 mixture)

  2. 2.

Short simulations with different time steps and with velocities saved:

     * •

DOI: 10.5281/zenodo.3956761 (3/3/2 mixture)

     * •

DOI: 10.5281/zenodo.3956812 (2/1/1 mixture)

  3. 3.

Additional simulations of the 3/3/2 mixture with different temperature
coupling or LINCS settings:

     * •

DOI: 10.5281/zenodo.3956775 (3/3/2 mixture)

     * •

DOI: 10.5281/zenodo.3956814 (2/1/1 mixture)


Supersymmetry with Cadabra

Behzad Mansouri*

QuantGates Ltd, London, EC2A 4NE, UK

Abstract

> This article provides a quick guide to the implementation of supersymmetry
> with Cadabra, a symbolic computer algebra system. Details are provided on
> the implementation of Grassmann variables, fermionic fields, supercharges,
> and superderivatives and the treatment of superfield expressions and
> supersymmetric Lagrangians. Finally, the automation of supersymmetric
> Lagrangian generation in the Cadabra programming framework is discussed.

††*<EMAIL_ADDRESS>
##  1 Introduction to Supersymmetry and Cadabra

Supersymmetry, strings and branes are conjectured to be key ingredients toward
establishing a unique unified theory of physics. In addition to the extension
of possible symmetry groups relating bosons to fermions, new structures within
the theory seem to provide approaches to treating some conceptual problems in
high-energy physics, including the hierarchy problem.

The underlying calculations in supersymmetric theories, especially in higher
dimensions, are technically long, and using a symbolic computer algebra system
for the implementation of such highly structured theories can facilitate
calculations in the model under study. Cadabra is a new computer algebra
system (CAS) that was developed specifically for solutions to problems that
have emerged in field theory. It was first developed to study the
supersymmetry of higher-derivative effective actions [1, 2]. The extensive
capabilities of Cadabra, such as the ability to deal with any type of tensor,
anticommuting variables, Clifford algebras and Fierz transformations, as well
as its many ready-to-run simplification algorithms and fully programmable
features, make it a powerful tool to efficiently study problems emerging in
quantum field theory and string theory. In [3, 4, 5], Cadabra was introduced,
and the various capabilities of this new computer algebra system, such as its
comprehensive functionality for large-scale tensor computations, were
discussed. In [6, 7, 8, 9, 10], Cadabra was extensively used to verify or
derive results in the supergravity context.

The aim of this paper is to illustrate how Cadabra is used for supersymmetric
computations and provide the source code for a variety of explicit
calculations, paving the way for better and faster routes to study
supersymmetric aspects of string theory and M-theory in more depth via
Cadabra.

###  1.1 Indices, Properties and Algorithms in Cadabra

In Cadabra, all objects in the form of subscripts or superscripts are
considered to be $indices$. Cadabra recognizes free and dummy indices, and it
has an automatic index renaming feature. Cadabra uses LaTeX syntax.

    Ψex1:= \Gamma^{i j r} \Gamma^{k l m n} \Gamma_{r};

We can assign mathematical properties to symbols:

    \Gamma_{#}::GammaMatrix.
    \partial{#}::PartialDerivative.

With the above declaration, the properties of `GammaMatrix` and
`PartialDerivative` are attached to the symbols $\Gamma$ and $\partial$,
respectively.

After writing the mathematical expressions, computations can be executed on
them by applying Cadabra’s algorithms. There are many built-in algorithms for
substitution and variation, index manipulations, sorting and simplifications,
as well as those that have been designed to act on spinors and much more. In
the following sections, we demonstrate the use of a variety of Cadabra’s
algorithms to perform different calculations.

###  1.2 Gamma Matrices and Fierz transformations

In this section, we show how to derive useful gamma matrices and Fierz
identities within Cadabra. We define the indices and denote the underlying
dimension and attach gamma matrix property to the $\Gamma$ symbol, and we also
associate a metric with it:

    {i,j,k,l,m,n}::Indices(vector).
    {i,j,k,l,m,n}::Integer(0..7).
    \Gamma_{#}::GammaMatrix(metric=\delta).
    \delta_{m n}::KroneckerDelta.

In Cadabra, simplification steps are contained in the function `post_process`,
which is executed on every new input. Based on our needs, we update it to
consider `sort_product`, `eliminate_kronecker` and `canonicalise`:

    def post_process(ex):
        sort_product(ex)
        eliminate_kronecker(ex)
        canonicalise(ex)
        collect_terms(ex)

We aim to find an identity for the below expression:

    ex1:= \Gamma^{i j m} \Gamma^{k l} \Gamma_{m};

$\begin{array}[]{l}{}\Gamma^{ijm}\Gamma^{kl}\Gamma_{m}\end{array}$ (1)

By executing `join_gamma` on the expression, which joins two fully
antisymmetrized gamma matrix products, and distributing factors over sums by
`distribute`, we arrive at the below result:

    join_gamma(_);
    distribute(_);

$\begin{array}[]{l}\Gamma^{ijklm}\Gamma_{m}+\Gamma^{ijl}\Gamma^{k}-\Gamma^{ijk}\Gamma^{l}+\Gamma^{ilm}\Gamma_{m}\delta^{jk}-\Gamma^{ikm}\Gamma_{m}\delta^{jl}-\Gamma^{jlm}\Gamma_{m}\delta^{ik}+\Gamma^{jkm}\Gamma_{m}\delta^{il}+\\\
\Gamma^{i}\Gamma^{k}\delta^{jl}-\Gamma^{i}\Gamma^{l}\delta^{jk}-\Gamma^{m}\Gamma_{m}\delta^{ik}\delta^{jl}+\Gamma^{m}\Gamma_{m}\delta^{il}\delta^{jk}+\Gamma^{j}\Gamma^{l}\delta^{ik}-\Gamma^{j}\Gamma^{k}\delta^{il}\end{array}$
(2)

By repeating the above procedure once more:

    join_gamma(_);
    distribute(_);

$\begin{array}[]{l}2\Gamma^{ijkl}+4\Gamma^{il}\delta^{jk}-4\Gamma^{jl}\delta^{ik}-4\Gamma^{ik}\delta^{jl}+4\Gamma^{jk}\delta^{il}-6\delta^{ik}\delta^{jl}+6\delta^{il}\delta^{jk}\end{array}$
(3)

Now, we perform a substitution for the term $\Gamma^{ijkl}$:

    substitute(_, $\Gamma^{i j k l} -> \Gamma^{k l}\Gamma^{i j}
    +\Gamma^{j k}\delta^{i l}-\Gamma^{i k}\delta^{j l}
    -\Gamma^{j l}\delta^{i k}+\Gamma^{i l}\delta^{j k}
    -\delta^{i l}\delta^{j k}+\delta^{i k}\delta^{j l}$ );

Finally, one can arrive at the below useful identity:

$\begin{array}[]{l}{}\Gamma^{ijm}\Gamma^{kl}\Gamma_{m}=2\Gamma^{kl}\Gamma^{ij}+6\Gamma^{jk}\delta^{il}-6\Gamma^{ik}\delta^{jl}-6\Gamma^{jl}\delta^{ik}+6\Gamma^{il}\delta^{jk}+4\delta^{il}\delta^{jk}-4\delta^{ik}\delta^{jl}\end{array}$
(4)

The above substitution is obtained using the relation from the below
manipulation:

    ex2:= \Gamma^{k l}\Gamma^{i j};

$\begin{array}[]{l}\Gamma^{kl}\Gamma^{ij}\end{array}$ (5)

    join_gamma(_);

$\begin{array}[]{l}\Gamma^{ijkl}-\Gamma^{jk}\delta^{il}+\Gamma^{ik}\delta^{jl}+\Gamma^{jl}\delta^{ik}-\Gamma^{il}\delta^{jk}+\delta^{il}\delta^{jk}-\delta^{ik}\delta^{jl}\end{array}$
(6)

On a product of four spinors, one can apply a Fierz transformation. Here, in
eleven-dimensional Clifford algebra, we define our spinors as follows:

    {m,n,p,q,r,s}::Indices;
    {m,n,p,q,r,s}::Integer(0..10);
    \Gamma{#}::GammaMatrix;
    \bar{#}::DiracBar;
    {\epsilon, \lambda, \chi, \psi}::Spinor;

Consider the following expression of the spinor product:

    ex:=\bar{\epsilon} \Gamma_{m} \psi \bar{\chi} \Gamma^{m} \lambda;

$\begin{array}[]{l}{}\bar{\epsilon}\Gamma_{m}\psi\bar{\chi}\Gamma^{m}\lambda\end{array}$
(7)

Now, we apply a Fierz transformation, with the desired order as follows:

    fierz(_, $\epsilon, \lambda, \chi, \psi$);

$\begin{array}[]{l}{}-\frac{1}{32}\bar{\epsilon}\Gamma_{m}\Gamma^{m}\lambda\bar{\chi}\psi-\frac{1}{32}\bar{\epsilon}\Gamma_{m}\Gamma^{n}\Gamma^{m}\lambda\bar{\chi}\Gamma^{n}\psi-\frac{1}{64}\bar{\epsilon}\Gamma_{m}\Gamma^{np}\Gamma^{m}\lambda\bar{\chi}\Gamma^{pn}\psi\\\
-\frac{1}{192}\bar{\epsilon}\Gamma_{m}\Gamma^{npq}\Gamma^{m}\lambda\bar{\chi}\Gamma^{qpn}\psi-\frac{1}{768}\bar{\epsilon}\Gamma_{m}\Gamma^{npqr}\Gamma^{m}\lambda\bar{\chi}\Gamma^{rqpn}\psi-\frac{1}{3840}\bar{\epsilon}\Gamma_{m}\Gamma^{npqrs}\Gamma^{m}\lambda\bar{\chi}\Gamma^{srqpn}\psi\end{array}$
(8)

##  2 Superalgebra and Superfield Implementations

Here, we first look at super-Poincaré algebra. We establish the vector and
spinorial indices, and then, by declaring the antisymmetric and noncommuting
features available in super-Poincaré generators, we simplify a commutation
relation in this regard. Additionally, in the following subsection, by
introducing the superspace formalism and supercharges in Cadabra, we shed
light on the implementation of superfields, including scalar, chiral and
vector superfields specifically for $N=1$ and $D=4$, to derive supersymmetric
transformations of component fields.

First, we define the vector and spinor indices for the superalgebra generators
as follows:

    {a,b,c}::Indices(spinor);
    {\mu,\nu,\rho,\sigma,\lambda,\kappa,\alpha,\beta,
    \gamma,\xi}::Indices(vector);
    {\mu,\nu,\rho,\sigma,\lambda,\kappa,\alpha,\beta,
    \gamma,\xi}::Integer(0..3);
    \delta{#}::KroneckerDelta;
    \epsilon_{\mu\nu\lambda\rho}::EpsilonTensor(delta=\delta);

At this stage, we declare the antisymmetric and noncommuting nature of our
generators by attaching the properties `AntiSymmetric`, `SelfNonCommuting` and
`NonCommuting` to them.

    {J_{\mu\nu},\sigma^{4}_{\mu\nu}} ::AntiSymmetric;
    J_{\mu\nu}::SelfNonCommuting;
    { J_{\mu\nu}, P_{\mu}, W_{\mu} }::NonCommuting;
    { J_{\mu\nu}, W_{\mu},Q_{a}}::NonCommuting;
    {J_{\mu\nu}, P_{\mu}, W_{\mu},Q_{a}}::Depends(\commutator{#});

$\begin{array}[]{l}\text{Property Indices(position=free) attached
to~{}}{[}{a,~{}b,~{}c}{]}.\\\ \text{Property Indices(position=free) attached
to~{}}{[}{\mu,~{}\nu,~{}\rho,~{}\sigma,~{}\lambda,~{}\kappa,~{}\alpha,~{}\beta,~{}\gamma,~{}\xi}{]}.\\\
\text{Property Integer attached
to~{}}{[}{\mu,~{}\nu,~{}\rho,~{}\sigma,~{}\lambda,~{}\kappa,~{}\alpha,~{}\beta,~{}\gamma,~{}\xi}{]}.\\\
\text{Property KroneckerDelta attached to~{}}\delta{(}{\\#}{)}.\\\
\text{Property EpsilonTensor attached to~{}}\epsilon_{\mu\nu\lambda\rho}.\\\
\text{Property AntiSymmetric attached
to~{}}\left[J_{\mu\nu},~{}\sigma^{4}\,_{\mu\nu}\right].\\\ \text{Property
SelfNonCommuting attached to~{}}J_{\mu\nu}.\\\ \text{Property NonCommuting
attached to~{}}{[}{J_{\mu\nu},~{}P_{\mu},~{}W_{\mu}}{]}.\\\ \text{Property
NonCommuting attached to~{}}{[}{J_{\mu\nu},~{}W_{\mu},~{}Q_{a}}{]}.\\\
\text{Property Depends attached
to~{}}{[}{J_{\mu\nu},~{}P_{\mu},~{}W_{\mu},~{}Q_{a}}{]}.\end{array}$ (9)

As mentioned in the previous section, simplification treatments are contained
in the `post_process` function, which is executed on every new input. Here, we
update it to consider `unwrap(ex)`, `eliminate_kronecker(ex)`,
`canonicalise(ex)` and `rename_dummies(ex)`:

    def post_process(ex):
        unwrap(ex)
        eliminate_kronecker(ex)
        canonicalise(ex)
        rename_dummies(ex)
        collect_terms(ex)

Now, we establish the commutator relations for $J_{\mu\nu},Q_{a}$ and
$P_{\mu}$ generators.

    superpoincare:= { \commutator{J_{\mu\nu}}{Q_{a}} ->
    -\sigma^{4}_{\mu \nu} Q_{b},\commutator{P_{\mu}}{Q_{a}} -> 0 };

$\begin{array}[]{l}{}\left[{}\left[J_{\mu\nu},Q_{a}\right]{}\rightarrow-\sigma^{4}\,_{\mu\nu}Q_{b},~{}{}\left[P_{\mu},Q_{a}\right]{}\rightarrow
0\right]\end{array}$ (10)

Below, we aim to compute the commutator between the square of the
Pauli–Lubanski vector and supercharge $Q_{a}$:

    Wsq:= \commutator{W_\mu W_\mu}{Q_{a}};

$\begin{array}[]{l}{}{}{[}{W_{\mu}W_{\mu},Q_{a}}{]{}}\end{array}$ (11)

We substitute the explicit definition of the Pauli–Lubanski vector into the
commutator:

    substitute(_, $W_\mu ->
    1/2 \epsilon_{\mu\nu\lambda\rho} J_{\nu\lambda} P_\rho $);

$\begin{array}[]{l}{}\frac{1}{4}\epsilon_{\alpha\mu\nu\kappa}\epsilon_{\alpha\lambda\rho\sigma}{}{[}{J_{\mu\nu}P_{\lambda}J_{\rho\sigma}P_{\kappa},Q_{a}}{]{}}\end{array}$
(12)

One can replace the product of two epsilon tensors with a generalized delta by
the below command:

    epsilon_to_delta(_);

$\begin{array}[]{l}{}\frac{3}{2}\delta_{\mu\rho\nu\sigma\lambda\kappa}{}{[}{J_{\mu\nu}P_{\lambda}J_{\rho\sigma}P_{\kappa},Q_{a}}{]{}}\end{array}$
(13)

Now, we expand the generalized delta into standard two-index Kronecker deltas:

    expand_delta(_);

$\begin{array}[]{l}\frac{3}{2}\left(\frac{1}{6}\delta_{\lambda\kappa}\delta_{\mu\rho}\delta_{\nu\sigma}-\frac{1}{6}\delta_{\lambda\nu}\delta_{\kappa\sigma}\delta_{\mu\rho}-\frac{1}{6}\delta_{\lambda\kappa}\delta_{\mu\sigma}\delta_{\nu\rho}+\frac{1}{6}\delta_{\lambda\mu}\delta_{\kappa\sigma}\delta_{\nu\rho}+\frac{1}{6}\delta_{\lambda\nu}\delta_{\kappa\rho}\delta_{\mu\sigma}-\frac{1}{6}\delta_{\lambda\mu}\delta_{\kappa\rho}\delta_{\nu\sigma}\right)\\\
{}\left[J_{\mu\nu}P_{\lambda}J_{\rho\sigma}P_{\kappa},Q_{a}\right]\end{array}$
(14)

By applying `product_rule`, we can simplify the commutator in the above:

    product_rule(_);

$\begin{array}[]{l}{}\frac{3}{2}\left(\frac{1}{6}\delta_{\lambda\kappa}\delta_{\mu\rho}\delta_{\nu\sigma}-\frac{1}{6}\delta_{\lambda\sigma}\delta_{\kappa\nu}\delta_{\mu\rho}-\frac{1}{6}\delta_{\lambda\kappa}\delta_{\mu\sigma}\delta_{\nu\rho}+\frac{1}{6}\delta_{\lambda\sigma}\delta_{\kappa\mu}\delta_{\nu\rho}+\frac{1}{6}\delta_{\lambda\rho}\delta_{\kappa\nu}\delta_{\mu\sigma}-\frac{1}{6}\delta_{\lambda\rho}\delta_{\kappa\mu}\delta_{\nu\sigma}\right)\\\
\left({}\left[J_{\mu\nu},Q_{a}\right]{}P_{\lambda}J_{\rho\sigma}P_{\kappa}+J_{\mu\nu}{}\left[P_{\kappa},Q_{a}\right]{}J_{\rho\sigma}P_{\lambda}+J_{\mu\nu}P_{\lambda}{}\left[J_{\rho\sigma},Q_{a}\right]{}P_{\kappa}+J_{\mu\nu}P_{\kappa}J_{\rho\sigma}{}\left[P_{\lambda},Q_{a}\right]{}\right)\par\end{array}$
(15)

At this level, one can insert the super-Poincaré commutators defined
previously into the expression and obtain:

    substitute(_, superpoincare);

$\begin{array}[]{l}{}\frac{3}{2}\left(\frac{1}{6}\delta_{\lambda\kappa}\delta_{\mu\rho}\delta_{\nu\sigma}-\frac{1}{6}\delta_{\lambda\nu}\delta_{\kappa\sigma}\delta_{\mu\rho}-\frac{1}{6}\delta_{\lambda\kappa}\delta_{\mu\sigma}\delta_{\nu\rho}+\frac{1}{6}\delta_{\lambda\nu}\delta_{\kappa\rho}\delta_{\mu\sigma}+\frac{1}{6}\delta_{\lambda\mu}\delta_{\kappa\sigma}\delta_{\nu\rho}-\frac{1}{6}\delta_{\lambda\mu}\delta_{\kappa\rho}\delta_{\nu\sigma}\right)\\\
\left(-\sigma^{4}\,_{\rho\sigma}Q_{b}P_{\lambda}J_{\mu\nu}P_{\kappa}-J_{\rho\sigma}P_{\lambda}\sigma^{4}\,_{\mu\nu}Q_{b}P_{\kappa}\right)\end{array}$
(16)

Finally, we can distribute factors over sums and obtain a non-vanishing
result:

    distribute(_);

$\begin{array}[]{l}{}-\frac{1}{2}\sigma^{4}\,_{\mu\nu}Q_{b}P_{\rho}J_{\mu\nu}P_{\rho}-\frac{1}{2}J_{\mu\nu}P_{\rho}\sigma^{4}\,_{\mu\nu}Q_{b}P_{\rho}-\sigma^{4}\,_{\rho\mu}Q_{b}P_{\rho}J_{\mu\nu}P_{\nu}-J_{\mu\nu}P_{\mu}\sigma^{4}\,_{\nu\rho}Q_{b}P_{\rho}\end{array}$
(17)

###  2.1 Scalar, Chiral and Vector Superfield Implementations

In this section, we first define a general Lorentz scalar superfield in terms
of its power series expansion. Then, by defining the derivative and
antiderivative operators (for Grassmann variables) and establishing the
related anticommuting properties, we aim to derive supersymmetric
transformations of component fields. We achieve this by first defining the
infinitesimal supersymmetry variation and applying it to the superfield with
`substitute` and further simplifications via `distribute`, `product_rule` and
`unwrap`. At the final stage, we sort the coefficients based on the same
powers of $\theta$ in the superfield expansion to read the SUSY transformation
of the component field. We carry this out via hard-coded substitutions using
`take_match` and `replace_match`. In a broader design, one can develop sorting
rules to simplify general expressions containing Grassmann variables. A
similar approach can be applied for the SUSY transformation of chiral and
vector superfields.

Here, we start by defining a general scalar superfield `Phi`:

    \tbar{#}::LaTeXForm("\bar\theta").
    \cbar{#}::LaTeXForm("\bar\chi").
    \lambar{#}::LaTeXForm("\bar\lambda").
    \aldot{#}::LaTeXForm("\dot\alpha").
    \betdot{#}::LaTeXForm("\dot\beta").

    Superfield:= Phi-> (f(x)+\theta^{\beta} \phi_{\beta}+ \tbar_{\betdot}
    \cbar^{\betdot}+\indexbracket{\theta\theta}m(x)
    +\indexbracket{\tbar\tbar}n(x)
    +\theta^{\beta}\sigma_{\beta \betdot}^{\nu}\tbar^{\betdot} V_{\nu}
    +\indexbracket{\theta\theta}\tbar_{\betdot}\lambar^{\betdot}
    +\indexbracket{\tbar\tbar}\theta^{\beta}\psi_{\beta}
    +\indexbracket{\theta\theta}\indexbracket{\tbar\tbar}d(x));

$\begin{array}[]{l}Phi\rightarrow
f{(}{x}{)}+\theta^{\beta}\phi_{\beta}+\bar{\theta}_{\dot{\beta}}\bar{\chi}^{\dot{\beta}}+{(}{\theta\theta}{)}m{(}{x}{)}+{(}{\bar{\theta}\bar{\theta}}{)}n{(}{x}{)}+\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}V_{\nu}+{(}{\theta\theta}{)}\bar{\theta}_{\dot{\beta}}\bar{\lambda}^{\dot{\beta}}+{(}{\bar{\theta}\bar{\theta}}{)}\theta^{\beta}\psi_{\beta}\\\
+{(}{\theta\theta}{)}{(}{\bar{\theta}\bar{\theta}}{)}d{(}{x}{)}\end{array}$
(18)

    \parbar{#}::LaTeXForm("\bar\partial").
    \partialmu{#}::LaTeXForm("\hat\partial").
    \epsbar{#}::LaTeXForm("\bar\epsilon").

    {\partial{#},\parbar{#},\partialmu{#}}::PartialDerivative;

    {\theta^{#},\tbar^{#},\epsilon^{#},\epsbar_{#}}::AntiCommuting;
    {\theta^{\alpha},\tbar^{\aldot}}::SelfAntiCommuting;
    \partial{#}::ImplicitIndex(\partial{#}_{\theta^{#}});
    \parbar{#}::ImplicitIndex(\parbar{#}_{\tbar^{#}});

    susy:=\epsilon^{\alpha} \partial_{\alpha}{Phi}+\epsbar_{\aldot}
    \parbar^{\aldot}{Phi}
    +\indexbracket{i\theta\sigma^{\mu}\epsbar} \partialmu_{\mu}{Phi}
    -\indexbracket{i\epsilon\sigma^{\mu}\tbar} \partialmu_{\mu}{Phi};

    \delta_{\alpha}^{\beta}::KroneckerDelta.
    \delta^{\aldot}_{\betdot}::KroneckerDelta.
    superderivative:={\partial_{\alpha}{\theta^{\beta}}->\delta_{\alpha}^{\beta},
    Ψ\parbar^{\aldot}{\tbar_{\betdot}}->\delta^{\aldot}_{\betdot},
    Ψ\parbar^{\aldot}{\tbar^{\betdot}}->-e^{\aldot \betdot},
    Ψ\partial_{\alpha}{\indexbracket{\theta\theta}}->2\theta_{\alpha},
    Ψ\parbar^{\aldot}{\indexbracket{\tbar\tbar}}->2\tbar^{\aldot}};

$\begin{array}[]{l}\text{Property PartialDerivative attached
to~{}}{[}\partial{\\#},\bar{\partial}{\\#},\hat{\partial}{\\#}{]}.\\\
{}\text{Property AntiCommuting attached
to~{}}{[}{\theta^{\\#},~{}\bar{\theta}^{\\#}}{,\epsilon^{\\#},\bar{\epsilon}_{\\#}}].\\\
{}\text{Property SelfAntiCommuting attached
to~{}}{[}{\theta^{\alpha},~{}\bar{\theta}^{\dot{\alpha}}}{]}.\\\
{}\text{Property ImplicitIndex attached to~{}}\partial{(}{\\#}{)}.\\\
{}\text{Property ImplicitIndex attached to~{}}\bar{\partial}{(}{\\#}{)}.\\\
\epsilon^{\alpha}\partial_{\alpha}{(}{Phi}{)}+\bar{\epsilon}_{\dot{\alpha}}\bar{\partial}^{\dot{\alpha}}{(}{Phi}{)}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\hat{\partial}_{\mu}{(}{Phi}{)}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\hat{\partial}_{\mu}{(}{Phi}{)}\\\
{}{[}{\partial_{\alpha}{(}{\theta^{\beta}}{)}\rightarrow\delta_{\alpha}\,^{\beta},~{}\bar{\partial}^{\dot{\alpha}}{(}{\bar{\theta}_{\dot{\beta}}}{)}\rightarrow\delta^{\dot{\alpha}}\,_{\dot{\beta}},~{}\bar{\partial}^{\dot{\alpha}}{(}{\bar{\theta}^{\dot{\beta}}}{)}\rightarrow-e^{\dot{\alpha}\dot{\beta}},~{}\partial_{\alpha}{(}{{(}{\theta\theta}{)}}{)}\rightarrow
2\theta_{\alpha},~{}\bar{\partial}^{\dot{\alpha}}{(}{{(}{\bar{\theta}\bar{\theta}}{)}}{)}\rightarrow
2\bar{\theta}^{\dot{\alpha}}}{]}\end{array}$ (19)

    substitute(susy,Superfield);

    distribute(_)
    product_rule(_)

    {\theta^{#},\indexbracket{\theta\theta}}
    ::Depends(\partial{#});
    {\tbar_{#},\indexbracket{\tbar\tbar}}
    ::Depends(\parbar{#});
    {f(x),m(x),n(x),\phi_{\alpha},\cbar^{\aldot},V_{\nu},\lambar^{\betdot},
    Ψ\psi_{\beta}}::Depends(\partialmu{#});

$\begin{array}[]{l}\text{Property Depends attached
to~{}}{[}{\theta^{\\#},{(}{\theta\theta}{)}}{]}.\\\ \text{Property Depends
attached to~{}}{[}{\bar{\theta}_{\\#},{(}{\bar{\theta}\bar{\theta}}{)}}{]}.\\\
\text{Property Depends attached
to~{}}{[}{f{(}{x}{)},~{}m{(}{x}{)},~{}n{(}{x}{)},~{}\phi_{\alpha},~{}\bar{\chi}^{\dot{\alpha}},~{}V_{\nu},~{}\bar{\lambda}^{\dot{\beta}},~{}\psi_{\beta}}{]}.\end{array}$
(20)

    unwrap(_);

$\begin{array}[]{l}\epsilon^{\alpha}\partial_{\alpha}{\theta^{\beta}}\phi_{\beta}+\epsilon^{\alpha}\partial_{\alpha}{(}{{(}{\theta\theta}{)}}{)}m{(}{x}{)}+\epsilon^{\alpha}\partial_{\alpha}{\theta^{\beta}}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}V_{\nu}+\epsilon^{\alpha}\partial_{\alpha}{(}{{(}{\theta\theta}{)}}{)}\bar{\theta}_{\dot{\beta}}\bar{\lambda}^{\dot{\beta}}+\epsilon^{\alpha}{(}{\bar{\theta}\bar{\theta}}{)}\partial_{\alpha}{\theta^{\beta}}\psi_{\beta}+\\\
\epsilon^{\alpha}\partial_{\alpha}{(}{{(}{\theta\theta}{)}}{)}{(}{\bar{\theta}\bar{\theta}}{)}d{(}{x}{)}+\bar{\epsilon}_{\dot{\alpha}}\bar{\partial}^{\dot{\alpha}}{\bar{\theta}_{\dot{\beta}}}\bar{\chi}^{\dot{\beta}}+\bar{\epsilon}_{\dot{\alpha}}\bar{\partial}^{\dot{\alpha}}{(}{{(}{\bar{\theta}\bar{\theta}}{)}}{)}n{(}{x}{)}-\bar{\epsilon}_{\dot{\alpha}}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\partial}^{\dot{\alpha}}{\bar{\theta}^{\dot{\beta}}}V_{\nu}+\bar{\epsilon}_{\dot{\alpha}}{(}{\theta\theta}{)}\bar{\partial}^{\dot{\alpha}}{\bar{\theta}_{\dot{\beta}}}\bar{\lambda}^{\dot{\beta}}+\\\
\bar{\epsilon}_{\dot{\alpha}}\bar{\partial}^{\dot{\alpha}}{(}{{(}{\bar{\theta}\bar{\theta}}{)}}{)}\theta^{\beta}\psi_{\beta}+\bar{\epsilon}_{\dot{\alpha}}{(}{\theta\theta}{)}\bar{\partial}^{\dot{\alpha}}{(}{{(}{\bar{\theta}\bar{\theta}}{)}}{)}d{(}{x}{)}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\hat{\partial}_{\mu}{f{(}{x}{)}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\phi_{\beta}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\chi}^{\dot{\beta}}}+\\\
{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\theta\theta}{)}\hat{\partial}_{\mu}{m{(}{x}{)}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\hat{\partial}_{\mu}{n{(}{x}{)}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}\hat{\partial}_{\mu}{V_{\nu}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\theta\theta}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\lambda}^{\dot{\beta}}}\\\
+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\psi_{\beta}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\hat{\partial}_{\mu}{f{(}{x}{)}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\phi_{\beta}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\chi}^{\dot{\beta}}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}{(}{\theta\theta}{)}\hat{\partial}_{\mu}{m{(}{x}{)}}\\\
-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\hat{\partial}_{\mu}{n{(}{x}{)}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}\hat{\partial}_{\mu}{V_{\nu}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}{(}{\theta\theta}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\lambda}^{\dot{\beta}}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\psi_{\beta}}\end{array}$
(21)

    take_match(_, $\indexbracket{i\theta\sigma^{\mu}\epsbar}
    \indexbracket{\theta\theta} Q??$)
    substitute(_, $A??+ B?? -> 0$)
    replace_match(_)

    take_match(_, $\indexbracket{i\epsilon\sigma^{\mu}\tbar}
    \indexbracket{\tbar\tbar} Q??$)
    substitute(_, $A??+ B?? -> 0$)
    replace_match(_)

    substitute(_,superderivative)
    eliminate_kronecker(_);

$\begin{array}[]{l}\epsilon^{\beta}\phi_{\beta}+2\epsilon^{\alpha}\theta_{\alpha}m{(}{x}{)}+\epsilon^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}V_{\nu}+2\epsilon^{\alpha}\theta_{\alpha}\bar{\theta}_{\dot{\beta}}\bar{\lambda}^{\dot{\beta}}+\epsilon^{\beta}{(}{\bar{\theta}\bar{\theta}}{)}\psi_{\beta}+2\epsilon^{\alpha}\theta_{\alpha}{(}{\bar{\theta}\bar{\theta}}{)}d{(}{x}{)}+\bar{\epsilon}_{\dot{\beta}}\bar{\chi}^{\dot{\beta}}\\\
+2\bar{\epsilon}_{\dot{\alpha}}\bar{\theta}^{\dot{\alpha}}n{(}{x}{)}+\bar{\epsilon}_{\dot{\alpha}}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}e^{\dot{\alpha}\dot{\beta}}V_{\nu}+\bar{\epsilon}_{\dot{\beta}}{(}{\theta\theta}{)}\bar{\lambda}^{\dot{\beta}}+2\bar{\epsilon}_{\dot{\alpha}}\bar{\theta}^{\dot{\alpha}}\theta^{\beta}\psi_{\beta}+2\bar{\epsilon}_{\dot{\alpha}}{(}{\theta\theta}{)}\bar{\theta}^{\dot{\alpha}}d{(}{x}{)}\\\
+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\hat{\partial}_{\mu}{f{(}{x}{)}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\phi_{\beta}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\chi}^{\dot{\beta}}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\hat{\partial}_{\mu}{n{(}{x}{)}}\\\
+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}\hat{\partial}_{\mu}{V_{\nu}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\psi_{\beta}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\hat{\partial}_{\mu}{f{(}{x}{)}}\\\
-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\theta^{\beta}\hat{\partial}_{\mu}{\phi_{\beta}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\chi}^{\dot{\beta}}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}{(}{\theta\theta}{)}\hat{\partial}_{\mu}{m{(}{x}{)}}\\\
-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}\hat{\partial}_{\mu}{V_{\nu}}-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}{(}{\theta\theta}{)}\bar{\theta}_{\dot{\beta}}\hat{\partial}_{\mu}{\bar{\lambda}^{\dot{\beta}}}\end{array}$
(22)

    take_match(_, $\theta^{\beta}\partialmu_{\mu}{\phi_{\beta}} Q??$)
    take_match(_, $\indexbracket{i\theta\sigma^{\mu}\epsbar} Q??$)
    substitute(_,$A?? B?? C??->\indexbracket{\theta\theta}
    \indexbracket{-i/2 \partialmu_{\nu}{\phi}\sigma^{\nu}\epsbar}$)
    replace_match(_)
    replace_match(_)

    take_match(_, $\tbar_{\betdot}\partialmu_{\mu}{\cbar^{\betdot}} Q??$)
    take_match(_, $\indexbracket{i\epsilon\sigma^{\mu}\tbar} Q??$)
    substitute(_,$A?? B?? C??-> -\indexbracket{\tbar\tbar}
    \indexbracket{i/2 \epsilon\sigma^{\nu}\partialmu_{\nu}{\cbar}}$)
    replace_match(_)
    replace_match(_)

    take_match(_, $\indexbracket{i\theta\sigma^{\mu}\epsbar}
    \partialmu_{\mu}{V_{\nu}} Q??$)
    substitute(_,$A?? B?? C?? D?? E??->
    i/2 \indexbracket{\theta\theta}\indexbracket{\tbar\epsbar}
    \partialmu^{\mu}{V_{\mu}}$)
    replace_match(_)

    take_match(_, $\partialmu_{\mu}{m(x)} Q??$)
    substitute(_,$\indexbracket{i\epsilon\sigma^{\mu}\tbar}
    \indexbracket{\theta\theta}->-i \indexbracket{\theta\theta}\tbar_{\betdot}
    \indexbracket{\epsilon\sigma^{\mu}e}^{\betdot}}$)
    replace_match(_)

    take_match(_, $\partialmu_{\mu}{V_{\nu}} Q??$)
    substitute(_,$K?? L?? M?? N?? O??->
    i/2 \indexbracket{\tbar\tbar}\indexbracket{\theta\epsilon}
    \partialmu^{\mu}{V_{\mu}}$)
    replace_match(_)

    take_match(_, $\partialmu_{\mu}{\psi_{\beta}} Q??$)
    substitute(_,$K?? L?? M?? N??  ->i/2 \indexbracket{\theta\theta}
    \indexbracket{\tbar\tbar}\partialmu_{\nu}{\psi}\sigma^{\nu}\epsbar$)
    replace_match(_)

    take_match(_, $\partialmu_{\mu}{\lambar^{\betdot}} Q??$)
    substitute(_,$K?? L?? M?? N??  ->
    -i/2 \indexbracket{\theta\theta}\indexbracket{\tbar\tbar}
    \epsilon\sigma^{\nu}\partialmu_{\nu}{\lambar}$)
    replace_match(_)

    take_match(_, $\tbar_{\betdot}\lambar^{\betdot} Q??$)
    substitute(_,$K?? L?? M?? N?? ->
    1/2 \indexbracket{\theta\sigma^{\mu}\tbar}
    \indexbracket{\epsilon\sigma_{\mu}\lambar}$)
    replace_match(_)

    take_match(_, $\theta^{\beta}\psi_{\beta} Q??$)
    substitute(_,$K?? L?? M?? N?? ->
    1/2 \indexbracket{\theta\sigma^{\mu}\tbar}
    \indexbracket{\psi\sigma_{\mu}\epsbar}$)
    replace_match(_)

    take_match(_, $\partialmu_{\mu}{\phi_{\beta}} Q??$)
    substitute(_,$K?? L?? M?? ->
    -i/2 \indexbracket{\theta\sigma^{\nu}\tbar}
    \indexbracket{\epsilon\partialmu_{\nu}{\phi}$)
    replace_match(_)

    take_match(_, $\tbar_{\betdot}\partialmu_{\mu}{\cbar^{\beta}} Q??$)
    substitute(_,$K?? L?? M?? ->
    -i/2 \indexbracket{\theta\sigma_{\nu}\tbar}
    \indexbracket{\partialmu^{\nu}{\cbar}\epsbar}$)
    replace_match(_);

$\begin{array}[]{l}\epsilon^{\beta}\phi_{\beta}+2\epsilon^{\alpha}\theta_{\alpha}m{(}{x}{)}+\epsilon^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}\bar{\theta}^{\dot{\beta}}V_{\nu}+\epsilon^{\beta}{(}{\bar{\theta}\bar{\theta}}{)}\psi_{\beta}+2\epsilon^{\alpha}\theta_{\alpha}{(}{\bar{\theta}\bar{\theta}}{)}d{(}{x}{)}+\bar{\epsilon}_{\dot{\beta}}\bar{\chi}^{\dot{\beta}}+2\bar{\epsilon}_{\dot{\alpha}}\bar{\theta}^{\dot{\alpha}}n{(}{x}{)}\\\
+\bar{\epsilon}_{\dot{\alpha}}\theta^{\beta}\sigma_{\beta\dot{\beta}}\,^{\nu}e^{\dot{\alpha}\dot{\beta}}V_{\nu}+\bar{\epsilon}_{\dot{\beta}}{(}{\theta\theta}{)}\bar{\lambda}^{\dot{\beta}}+2\bar{\epsilon}_{\dot{\alpha}}{(}{\theta\theta}{)}\bar{\theta}^{\dot{\alpha}}d{(}{x}{)}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}\hat{\partial}_{\mu}{f{(}{x}{)}}+{(}{i\theta\sigma^{\mu}\bar{\epsilon}}{)}{(}{\bar{\theta}\bar{\theta}}{)}\hat{\partial}_{\mu}{n{(}{x}{)}}\\\
-{(}{i\epsilon\sigma^{\mu}\bar{\theta}}{)}\hat{\partial}_{\mu}{f{(}{x}{)}}+{(}{\theta\theta}{)}{(}{-\frac{1}{2}i\hat{\partial}_{\nu}{\phi}\sigma^{\nu}\bar{\epsilon}}{)}+{(}{\bar{\theta}\bar{\theta}}{)}{(}{\frac{1}{2}i\epsilon\sigma^{\nu}\hat{\partial}_{\nu}{\bar{\chi}}}{)}+\frac{1}{2}i{(}{\theta\theta}{)}{(}{\bar{\theta}\bar{\epsilon}}{)}\hat{\partial}^{\mu}{V_{\mu}}\\\
+i{(}{\theta\theta}{)}\bar{\theta}_{\dot{\beta}}{(}{\epsilon\sigma^{\mu}e}{)}\,^{\dot{\beta}}\hat{\partial}_{\mu}{m{(}{x}{)}}-\frac{1}{2}i{(}{\bar{\theta}\bar{\theta}}{)}{(}{\theta\epsilon}{)}\hat{\partial}^{\mu}{V_{\mu}}+\frac{1}{2}i{(}{\theta\theta}{)}{(}{\bar{\theta}\bar{\theta}}{)}\hat{\partial}_{\nu}{\psi}\sigma^{\nu}\bar{\epsilon}\\\
+\frac{1}{2}i{(}{\theta\theta}{)}{(}{\bar{\theta}\bar{\theta}}{)}\epsilon\sigma^{\nu}\hat{\partial}_{\nu}{\bar{\lambda}}+{(}{\theta\sigma^{\mu}\bar{\theta}}{)}{(}{\epsilon\sigma_{\mu}\bar{\lambda}}{)}+{(}{\theta\sigma^{\mu}\bar{\theta}}{)}{(}{\psi\sigma_{\mu}\bar{\epsilon}}{)}\\\
+\frac{1}{2}i{(}{\theta\sigma^{\nu}\bar{\theta}}{)}{(}{\epsilon\hat{\partial}_{\nu}{\phi}}{)}-\frac{1}{2}i{(}{\theta\sigma_{\nu}\bar{\theta}}{)}{(}{\hat{\partial}^{\nu}{\bar{\chi}}\bar{\epsilon}}{)}\end{array}$
(23)

Now, one can easily read off the supersymmetry transformation of the component
fields from the above expression.

One may find that the various hard-coded substitution routines at specific
terms that we used to achieve the canonical form for $\theta$ expansions make
the computation procedure more difficult. As mentioned previously, a broader
algorithm would be more feasible to sort these Grassmann expressions at least
to a better stopping point. For example, in the explicit indices, we can have:

    {\tbar^{#}, \theta^{#}}::SelfAntiCommuting;
    {\tbar^{#},\epsbar^{#}, \theta^{#},\epsilon^{#},\psi_{#}}::AntiCommuting;
    {\tbar^{#},\theta^{#},\epsilon^{#}}::SortOrder;
    \partial{#}::PartialDerivative;
    ex:=\epsilon^{\alpha} M^{\mu}_{\alpha \aldot} \tbar^{\aldot}
    \theta^{\beta} M^{\nu}_{\beta \betdot} \tbar^{\betdot}
    \partial_{\mu}{V_{\nu}} +
    \theta^{\beta} M^{\mu}_{\beta \aldot} \epsbar^{\aldot}
    \tbar_{\betdot}\tbar^{\betdot} \theta^{\alpha} \partial_{\mu}{\psi_{\alpha}};

$\begin{array}[]{l}{}\text{Property SelfAntiCommuting attached
to~{}}{[}{\bar{\theta}^{\\#},~{}\theta^{\\#}}{]}.\\\ {}\text{Property
AntiCommuting attached
to~{}}{[}{\bar{\theta}^{\\#},~{}\bar{\epsilon}^{\\#},~{}\theta^{\\#},~{}\epsilon^{\\#},~{}\psi_{\\#}}{]}.\\\
{}\text{Property SortOrder attached
to~{}}{[}{\bar{\theta}^{\\#},~{}\theta^{\\#},~{}\epsilon^{\\#}}{]}.\\\
{}\text{Property PartialDerivative attached to~{}}\partial{\\#}.\\\
{}\epsilon^{\alpha}M^{\mu}\,_{\alpha\dot{\alpha}}\bar{\theta}^{\dot{\alpha}}\theta^{\beta}M^{\nu}\,_{\beta\dot{\beta}}\bar{\theta}^{\dot{\beta}}\partial_{\mu}{V_{\nu}}+\theta^{\beta}M^{\mu}\,_{\beta\dot{\alpha}}\bar{\epsilon}^{\dot{\alpha}}\bar{\theta}_{\dot{\beta}}\bar{\theta}^{\dot{\beta}}\theta^{\alpha}\partial_{\mu}{\psi_{\alpha}\par}\end{array}$
(24)

Here, $M^{\mu}$ is the Pauli matrices and the anticommutativity of our spinors
is indicated using the `AntiCommuting` property. We also define our preferred
order using `SortOrder`. Now, we apply `sort_product` to sort factors in our
Grassmann expression defined above:

    sort_product(_);

$\begin{array}[]{l}{}M^{\mu}\,_{\alpha\dot{\alpha}}M^{\nu}\,_{\beta\dot{\beta}}\bar{\theta}^{\dot{\alpha}}\bar{\theta}^{\dot{\beta}}\theta^{\beta}\epsilon^{\alpha}\partial_{\mu}{V_{\nu}}+M^{\mu}\,_{\beta\dot{\alpha}}\bar{\epsilon}^{\dot{\alpha}}\partial_{\mu}{\psi_{\alpha}}\bar{\theta}_{\dot{\beta}}\bar{\theta}^{\dot{\beta}}\theta^{\alpha}\theta^{\beta}\end{array}$
(25)

After this stage, one can use the relation among the Grassmann variables for
further simplifications. Finally, working with $\theta$-expansions can be
cumbersome. One can handle things in a covariant way, using projection
techniques to find component field expansions of supersymmetric actions. The
implementation of projection operators via supercovariant derivatives in the
current stage of Cadabra is a little tricky, but it is still possible to
perform computations on them.

At the end, chiral and vector superfields can be implemented in Cadabra in the
same way as general scalar superfields. Here, we quickly review chiral and
vector superfields. We use these multiplets in the next section for the
automation of supersymmetric Lagrangians.

A superfield $\Phi$ satisfying the constraint $\bar{D}\Phi=0$ is called a left
chiral superfield, where $\bar{D}$ is a covariant derivative for the
superfield. In terms of new variable
$y^{\mu}=x^{\mu}+i\theta\sigma^{\mu}\bar{\theta}$, the superfield has the
following power series expansion in $\theta$:

$\begin{array}[]{l}\Phi({y^{\mu}},{\theta^{\alpha}})=A(y)+\sqrt{2}\theta\psi(y)+\theta\theta
F(y)\end{array}$ (26)

where A(y) and F(y) are complex scalar fields, and $\psi(y)$ is a left-handed
Weyl spinor. By expanding the component fields in terms of variables
$x,\theta,\bar{\theta}$, we have

$\begin{array}[]{l}\Phi({x^{\mu}},{\theta^{\alpha}},{{\bar{\theta}}^{\dot{\alpha}}})=\phi(x)+\sqrt{2}\theta\psi(x)+\theta\theta
F(x)+i\theta{\sigma^{\mu}}\bar{\theta}{\partial_{\mu}}\phi(x)\\\
-\frac{i}{{\sqrt{2}}}(\theta\theta){\partial_{\mu}}\psi(x){\sigma^{\mu}}\bar{\theta}-\frac{1}{4}(\theta\theta)(\bar{\theta}\bar{\theta}){\partial_{\mu}}{\partial^{\mu}}\phi(x)\end{array}$
(27)

Under supersymmetry transformation, we have

$\delta\Phi=i(\varepsilon Q+\bar{\varepsilon}\bar{Q})\Phi$ (28)

As we defined the differential operator representations of these supercharges,
one can easily compute the SUSY transformation for the chiral field using the
same approach as we used for a general scalar superfield.

A vector superfield
$V({x^{\mu}},{\theta^{\alpha}},{{\bar{\theta}}^{\dot{\alpha}}})$ is defined by
the reality condition $V=V^{\dagger}$. In the Wess–Zumino gauge, a vector
superfield is expressed as:

$V_{W.Z}=(\theta\sigma^{\mu}\bar{\theta})V_{\mu}+(\theta\theta)\bar{\theta}\bar{\lambda}(x)+(\bar{\theta}\bar{\theta})\theta\lambda(x)+(\theta\theta)(\bar{\theta}\bar{\theta})D(x)$
(29)

where $V_{\mu}(x)$ is a real vector field, $\lambda(x)$ is a complex spinor
field and D(x) is a real scalar field.

In section 3, we choose a vector superfield to define the invariant kinetic
part of the superspace Lagrangian density and to show how one can achieve
simple automation for SUSY Lagrangian generation.

##  3 Supersymmetric Lagrangian Implementations

In an automated way, one can construct supersymmetric Lagrangians and actions
and verify the various properties and structures of the related SUSY model
directly within Cadabra.

Here, we shed light on the potential capabilities of Cadabra in the automatic
generation of supersymmetric Lagrangians, which can speed up verification
tasks and facilitate the study of SUSY theories that are analytically
computation-intensive. We do not address all of the details of Cadabra
implementations and only elaborate on areas in which one can develop
Cadabara’s algorithms and functions to handle the automation.

Part of the supersymmetric Lagrangian describing the dynamics of different
component fields in chiral and vector supermultiplets can be automated. The
most general supersymmetry Lagrangian describing the dynamics between various
multiplets can be constructed based on invariant kinetic and superpotential
parts and gauge interactions:

${\mathcal{L}={\left.{{\Phi^{\dagger}}_{i}{e^{{q_{j}}{V^{j}}}}{\Phi_{i}}}\right|_{\theta\theta\bar{\theta}\bar{\theta}}}+{\left.{W[\Phi]}\right|_{\theta\theta}}+{\left.{\bar{W}[{\Phi^{\dagger}}]}\right|_{\bar{\theta}\bar{\theta}}}+\frac{1}{{16}}{\left.{W_{A}^{\alpha}W_{\alpha}^{A}}\right|_{\theta\theta}}+\frac{1}{{16}}{\left.{\bar{W}_{\dot{\alpha}}^{A}\bar{W}_{A}^{\dot{\alpha}}}\right|_{\bar{\theta}\bar{\theta}}}}$
(30)

where $W[\Phi]$ is a polynomial of the superfield $\Phi$ and is called the
superpotential. $W_{A}^{\alpha}$ and $\bar{W}_{\dot{\alpha}}^{A}$ are our
spinorial superfields. Based on the above structure, Lagrangian density
generation can be automated by defining three main functions in Cadabra as the
extractor of the ${\theta\theta\bar{\theta}\bar{\theta}}$, ${\theta\theta}$
and ${\bar{\theta}\bar{\theta}}$ components.

Here, for educational purposes, we show how to achieve simple automation for
the gauge covariant kinetic term ( ${\theta\theta\bar{\theta}\bar{\theta}}$
component extractor). Similarly, one can develop ${\theta\theta}$ and
${\bar{\theta}\bar{\theta}}$ component extractors out of scalar superfields of
$W_{A}^{\alpha}W_{\alpha}^{A}$ and
$\bar{W}_{\dot{\alpha}}^{A}\bar{W}_{A}^{\dot{\alpha}}$ and from
superpotentials $W[\Phi]$ and $\bar{W}[{\Phi^{\dagger}}]$.

We define a simple toy function called `T2TBar2` to derive the invariant
kinetic part of the SUSY Lagrangian from a list of chiral superfields in the
Wess–Zumino gauge. Here, for simplicity of the algorithm, we do not search for
the component of the vector superfield and use preconfigured components of
$V_{\mu}(x)$ , $\lambda(x)$ and $D(x)$.

    def T2TBar2(SF):
        F_terms:=0;
        Scalar_terms:=0;
        Spinor_lambda:=0;
        Spinors_psi:=0;
        Total_Terms:=0;

        for element in SF["F"]:
            F_term:= \mid \indexbracket{@(element)} \mid**2.
            F_terms += F_term

        for element in SF["A"]:
            indexa = element.indices().__next__()
            Scalar_term:= q_{@(indexa)}D(x) \mid
            \indexbracket{@(element)}\mid**2}
            + (\partial_{\mu}{@(element)}\partial^{\mu}{@(element)^{\star}}
            +1/4 q_{@(indexa)}**2 V_{\mu}(x) V^{\mu}(x) \mid
            \indexbracket{@(element)}\mid**2.
            Scalar_terms += Scalar_term

            for psi in SF["\psi"]:
                indexpsi = psi.indices().__next__()
                if str(indexpsi)==str(indexa):
                    Spinor_term1:= -1/\sqrt{2}} q_{@(indexa)}
                    \indexbracket{\lambda(x) @(psi)} @(element)^{\star}.
                    Spinor_term2:= -1/\sqrt{2} q_{@(indexa)}\bar{\lambda(x)}
                    \bar{@(psi)} @(element).
                else:
                    Spinor_term1:=0 ;
                    Spinor_term2:=0 ;

                Spinor_lambda +=Spinor_term1+Spinor_term2

        for psi in SF["\psi"]:
            indexpsi = psi.indices().__next__()
            Spinor:= i \partial_{\mu}{\bar{@(psi)}}\bar{\sigma}^{\mu}
            @(psi)-1/2 q_{@(indexpsi)} V(x)_{\mu}\bar{@(psi)}
            \bar{\sigma}^{\mu} @(psi).
            Spinors_psi += Spinor
        Total_Terms=  F_terms + Scalar_terms + Spinor_lambda + Spinors_psi
        return Total_Terms

The input is a simple list of field contents. For example, consider the case
of two scalar superfields (chiral):

    ΨSF:=[A_{1}(x),\psi_{1}(x),F_{1}(x),q_{1},A_{2}(x),\psi_{2}(x),
    ΨF_{2}(x),q_{2},V_{\mu}(x),\lambda(x),D(x)];

${[}{A_{1}{(}{x}{)},~{}\psi_{1}{(}{x}{)},~{}F_{1}{(}{x}{)},~{}q_{1},~{}A_{2}{(}{x}{)},~{}\psi_{2}{(}{x}{)},~{}F_{2}{(}{x}{)},~{}q_{2},~{}V_{\mu}{(}{x}{)},~{}\lambda{(}{x}{)},~{}D{(}{x}{)}}{]}$
(31)

    T2TBar2(SF);

$\begin{array}[]{l}{}\mid{(}{F_{1}{(}{x}{)}}{)}{\mid}^{2}+\mid{(}{F_{2}{(}{x}{)}}{)}{\mid}^{2}+q_{1}D{(}{x}{)}\mid{(}{A_{1}{(}{x}{)}}{)}{\mid}^{2}+\partial_{\mu}{[}{A_{1}{(}{x}{)}}{]}\partial^{\mu}{[}{A_{1}{(}{x}{)}\,^{\star}}{]}+\\\
\frac{1}{4}{q_{1}}^{2}V_{\mu}{(}{x}{)}V^{\mu}{(}{x}{)}\mid{(}{A_{1}{(}{x}{)}}{)}{\mid}^{2}+q_{2}D{(}{x}{)}\mid{(}{A_{2}{(}{x}{)}}{)}{\mid}^{2}+\partial_{\mu}{[}{A_{2}{(}{x}{)}}{]}\partial^{\mu}{[}{A_{2}{(}{x}{)}\,^{\star}}{]}+\\\
\frac{1}{4}{q_{2}}^{2}V_{\mu}{(}{x}{)}V^{\mu}{(}{x}{)}\mid{(}{A_{2}{(}{x}{)}}{)}{\mid}^{2}-{2}^{-\frac{1}{2}}q_{1}{(}{\lambda{(}{x}{)}{[}{\psi_{1}{(}{x}{)}}{]}}{)}{[}{A_{1}{(}{x}{)}\,^{\star}}{]}-{2}^{-\frac{1}{2}}q_{1}\bar{\lambda{(}{x}{)}}\bar{\psi_{1}{(}{x}{)}}{[}{A_{1}{(}{x}{)}}{]}-\\\
{2}^{-\frac{1}{2}}q_{2}{(}{\lambda{(}{x}{)}{[}{\psi_{2}{(}{x}{)}}{]}}{)}{[}{A_{2}{(}{x}{)}\,^{\star}}{]}-{2}^{-\frac{1}{2}}q_{2}\bar{\lambda{(}{x}{)}}\bar{\psi_{2}{(}{x}{)}}{[}{A_{2}{(}{x}{)}}{]}+i\partial_{\mu}{(}{\bar{\psi_{1}{(}{x}{)}}}{)}\bar{\sigma}\,^{\mu}{[}{\psi_{1}{(}{x}{)}}{]}-\\\
\frac{1}{2}q_{1}V{(}{x}{)}\,_{\mu}\bar{\psi_{1}{(}{x}{)}}\bar{\sigma}\,^{\mu}{[}{\psi_{1}{(}{x}{)}}{]}+i\partial_{\mu}{(}{\bar{\psi_{2}{(}{x}{)}}}{)}\bar{\sigma}\,^{\mu}{[}{\psi_{2}{(}{x}{)}}{]}-\frac{1}{2}q_{2}V{(}{x}{)}\,_{\mu}\bar{\psi_{2}{(}{x}{)}}\bar{\sigma}\,^{\mu}{[}{\psi_{2}{(}{x}{)}}{]}\end{array}$
(32)

The above example can provide the reader with simple guidance on how to
initiate programming in Cadabra specifically for supersymmetry use cases.
Obviously, the real value of such automation or other types of computations
will be much more noticeable when dealing with supersymmetric models in higher
dimensions. For example, many fields emerge in supergravity theories in ten
dimensions, usually restricting us to checking supersymmetry only with respect
to a subset of the transformation rules. With the help of Cadabra, heavy
computations, especially in the context of field theory, can be addressed more
easily.

## Acknowledgments

The author thanks Martin Cederwall and Daniel Butter for their helpful
comments.

## References

  * [1] K. Peeters, P. Vanhove, A. Westerberg, _Supersymmetric higherderivative actions in ten and eleven dimensions, the associated superalgebras and their formulation in superspace_ , Class.Quant.Grav. 18 (2001) 843-890, [arXiv:hep-th/0010167].
  * [2] K. Peeters, A. Westerberg, _The Ramond-Ramond sector of string theory beyond leading order_ , Class.Quant.Grav. 21 (2004) 1643-1666, [arXiv:hep-th/0307298].
  * [3] K. Peeters, _Introducing Cadabra: a symbolic computer algebra system for field theory problems_ , [arXiv:hep-th/0701238].
  * [4] K. Peeters, _A field-theory motivated approach to symbolic computer algebra_ , Comput. Phys.Commun. 176 (2007) 550 [arXiv:cs/0608005].
  * [5] L. Brewin, _A brief introduction to Cadabra: a tool for tensor computations in general relativity_ , Comput. Phys. Commun. 181 (3) (2010) 489–498 [arXiv:0903.2085].
  * [6] N. Ozdemir, M. Ozkan, O. Tunca, U. Zorba, _Three-dimensional extended Newtonian (super) gravity_ , JHEP 05 (2019) 130 [arXiv:1903.09377].
  * [7] K. Gubarev and E.T. Musaev, _Polyvector deformations in eleven-dimensional supergravity_ , Phys. Rev. D 103 (2021) 066021 [arXiv:2011.11424].
  * [8] D. Butter, J. Novak, M. Ozkan, Y. Pang and G. Tartaglino-Mazzucchelli, _Curvature squared invariants in six-dimensional N = (1, 0) supergravity_ , JHEP 04 (2019) 013 [arXiv:1808.00459].
  * [9] D. Butter, F. Ciceri and B. Sahoo, _N = 4 conformal supergravity: the complete actions_ , JHEP 01 (2020) 029 [arXiv:1910.11874].
  * [10] D. Butter, J. Novak and G. Tartaglino-Mazzucchelli, _The component structure of conformal supergravity invariants in six dimensions_ , JHEP 05 (2017) 133 [arXiv:1701.08163].