# Symmetric Mass Generation in the 1+1 Dimensional Chiral Fermion 3-4-5-0
Model

Meng Zeng  Department of Physics, University of California San Diego, La
Jolla, California 92093, USA  Zheng Zhu  Kavli Institute for Theoretical
Sciences, University of Chinese Academy of Sciences, Beijing 100190, China
CAS Center for Excellence in Topological Quantum Computation, University of
Chinese Academy of Sciences, Beijing 100190, China  Juven Wang  Center of
Mathematical Sciences and Applications, Harvard University, MA 02138, USA  Yi-
Zhuang You  Department of Physics, University of California San Diego, La
Jolla, California 92093, USA

###### Abstract

Lattice regularization of chiral fermions has been a long-standing problem in
physics. In this work, we present the density matrix renormalization group
(DMRG) simulation of the 3-4-5-0 model of (1+1)D chiral fermions with an
anomaly-free chiral U(1) symmetry, which contains two left-moving and two
right-moving fermions carrying U(1) charges 3,4 and 5,0, respectively.
Following the Wang-Wen chiral fermion model, we realize the chiral fermions
and their mirror partners on the opposite boundaries of a thin strip of (2+1)D
lattice model of multi-layer Chern insulator, whose finite-width implies the
quantum system is effectively (1+1)D. By introducing carefully designed two
sets of six-fermion local interactions to the mirror sector only, we
demonstrate that the mirror fermions can be gapped out by the interaction
above a critical strength without breaking the chiral U(1) symmetry, via the
symmetric mass generation (SMG) mechanism. We show that the interaction-driven
gapping transition is in the Berezinskii-Kosterlitz-Thouless (BKT)
universality class. We determine the evolution of Luttinger parameters before
the transition, which confirms that the transition happens exactly at the
point when the interaction term becomes marginal. As the mirror sector is
gapped after the transition, we check that the fermions in the light chiral
fermion sector remain gapless, which provides the desired lattice
regularization of chiral fermions.

_Introduction._ — It has been a long-standing issue to regularize chiral gauge
theories (e.g., the weak interaction in the Standard Model) on the lattice due
to the Nielsen-Ninomiya no-go theorem [1], which asserts that any free fermion
lattice model in even-dimensional spacetime with locally-realized chiral
symmetry will necessarily give rise to equal numbers of left-handed and right-
handed fermion fields at low energy, hence rendering the theory vector-like.
Over the past few decades, much effort [2, 3, 4, 5, 6, 7] has been devoted to
circumvent the fermion doubling problem by lifting different assumptions of
the no-go theorem.

In particular, the no-go theorem assumes the fermion theory to be infrared
free, i.e., fermion interactions, if there are any, must be perturbatively
irrelevant under the renormalization group (RG) flow. Lifting this assumption
by introducing non-perturbative (strong enough) fermion interactions could
potentially circumvent the problem. Efforts along this line are generally
referred to as the mirror fermion approach, which dates back to Eichten and
Preskill [8]. The basic idea is to start with a vector-like theory containing
both chiral fermions and their mirror fermion partners, which can be put on a
lattice without any issue. Then one attempts to generate a mass gap in the
mirror sector by introducing interactions among mirror fermions, such that the
remaining light (chiral fermion) sector survives in the low-energy spectrum,
providing the basis for lattice realizations of chiral gauge theories.
However, early numerical tests [9, 10, 11, 12, 13, 14, 15, 16, 17, 18]
appeared to invalidate the mirror fermion approach, as strong fermion
interactions typically result in the condensation of fermion bilinear mass at
low energy, which spontaneously breaks the chiral symmetry and gaps out the
light sector together with the mirror sector.

In recent years, a series of developments [19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33, 34] in the many-body quantum matter community have
significantly deepened our understanding. It is realized that in order to gap
out the mirror sector by interactions without breaking the chiral symmetry,
two conditions must be satisfied: (i) the mirror fermions must be _anomaly
free_ under the full spacetime-internal symmetry, (ii) the interaction must be
appropriately designed to satisfy certain _consistent gapping conditions_ [26,
33]. Along this line, recent numerical studies [35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49, 50, 51] have successfully demonstrated
examples of interaction-driven fermion mass generation without spontaneous
symmetry breaking in various spacetime dimensions. The phenomenon is known as
the _symmetric mass generation_ (SMG) [52, 53, 54, 55, 56]. Therefore, solving
the chiral fermion problem boils down to achieving the SMG for mirror fermions
in even spacetime dimensions.

Nevertheless, most numerical works realizing SMG in even spacetime dimensions
have been focused on vector-like lattice models [39, 40, 41, 44, 45, 48, 49,
51], which still have some distance from the goal of regularizing chiral
fermions. Recently, Catterall [50] studied the SMG of a chiral fermion lattice
model with a chiral discrete $\mathbb{Z}_{4}$ symmetry. In this work, we
demonstrate the SMG in the 3-4-5-0 model of (1+1)D chiral fermions that
cancels the $\mathbb{Z}$-class _perturbative_ local anomaly of the chiral
continuous $\mathrm{U}(1)$ symmetry, which is closer to the situation of
perturbative chiral anomaly cancellation in the (3+1)D Standard Model (such as
the chiral U(1)Y electroweak hypercharge). We propose a lattice model of
interacting fermions, and investigate the model using the density matrix
renormalization group (DMRG) numerical method [57, 58]. Our numerical results
provide clear evidence for the SMG in the mirror sector, successfully
achieving our goal of regularizing chiral fermions in the 3-4-5-0 model on a
lattice.

_The 3-4-5-0 Model._ — The 3-4-5-0 model describes four gapless complex
fermions in (1+1)D,

$S=\int\mathrm{d}t\;\mathrm{d}x\;\sum_{I=1}^{4}\psi_{I}^{\dagger}(\mathrm{i}\partial_{t}+\mathrm{i}v_{I}\partial_{x})\psi_{I},$
(1)

with two left-moving modes $\psi_{1},\psi_{2}$ (of $v_{1}=v_{2}=+1$) and two
right-moving modes $\psi_{3},\psi_{4}$ (of $v_{3}=v_{4}=-1$). The fermions are
charged under a chiral $\mathrm{U}(1)$ symmetry:
$\psi_{I}\to\mathrm{e}^{\mathrm{i}q_{I}\theta}\psi_{I}$, with the charge
assignment $(q_{1},q_{2},q_{3},q_{4})=(3,4,5,0)$ (hence the name “3-4-5-0”).
This seemly peculiar charge assignment is designed to cancel the
$\mathrm{U}(1)$ symmetry’s ’t Hooft anomaly, which is an $\mathbb{Z}$-class
perturbative local anomaly. The anomaly index is given by
$\sum_{I}v_{I}q_{I}^{2}=3^{2}+4^{4}-5^{2}-0^{2}=0$, which vanishes for the
charge assignment of the 3-4-5-0 model. The model is also free of the
gravitational anomaly. As the field theory is anomaly-free, it should admit a
lattice regularization in (1+1)D spacetime.

Figure 1: (Color online.) (a) The fermion hopping pattern on the two-leg
ladder lattice for the first layer. Arrow link:
$t_{1}\mathrm{e}^{\mathrm{i}\pi/4}$ (along the arrow direction), solid link:
$t_{2}$, dashed link: $-t_{2}$. This (2+1)D thin strip is effectively the same
as (1+1)D by regarding the finite-width dimension as internal degrees of
freedom of the (1+1)D system. (b) Energy dispersion for $t_{1}=1,t_{2}=0.5$.
Gapless edge modes are strictly localized on the two boundaries of the ladder.
(c) Schematic diagram showing the configuration of the four flavors of chiral
fermions on the edges.

Following Wang-Wen’s chiral fermion model [26, 33], the (1+1)D chiral fermions
and their mirror partners can be viewed as the chiral edge modes on the
opposite boundaries of a (2+1)D multi-layer Chern insulator [59], each layer
with a Chern number $\pm 1$. To construct the chiral fermions on a lattice, we
start with four layers of Chern insulators on a two-leg ladder as shown in
Fig. 1(a). On each lattice site $i$, we introduce four complex fermions,
described by the annihilation operators $\psi_{i,I}$ (with $I=1,2,3,4$ being
the layer/flavor index). The fermion hopping is governed by the lattice
Hamiltonian

$H_{\text{free}}=\sum_{I=1}^{4}\sum_{i,j}(t_{I,ij}\psi_{I,i}^{\dagger}\psi_{I,j}+\text{h.c.}),$
(2)

where the hopping parameters $t_{I,ij}$ are non-zero only on the nearest and
next-nearest neighbor links. For the first two layers $I=1,2$, the nearest
neighbor hoppings are purely imaginary with
$t_{I,ij}=\mathrm{e}^{\mathrm{i}\pi/4}t_{1}$ if $j\to i$ follows the link
direction, and the next-nearest neighbor hoppings are real with
$t_{I,ij}=t_{2}$ (or $-t_{2}$) on the solid (or dashed) links, as shown in
Fig. 1(a). We fix $t_{1}=1$ and $t_{2}=0.5$. This hopping pattern ensures a
$\pi$ Berry flux through each square plaquette, realizing a minimal model of
Chern insulator in each layer. For the last two layers $I=3,4$, the hopping
parameters are complex conjugated, such that the band Chern numbers in the
last two layers are opposite to those of the first two layers.

The lattice model has a four-site unit cell that repeats along the ladder
direction, hence the lattice momentum $k$ along the ladder direction is a good
quantum number, and the system is effectively (1+1)D. In each layer, the
single-particle energy dispersion (band structure) is shown in Fig. 1(b),
which includes two gapped bulk bands together and two gapless edge modes of
opposite velocities (localized separately on the two boundaries). Stacking all
layers together, the lattice model realizes four chiral fermions (as two pairs
of counter-propagating modes) on each edge, as illustrated in Fig. 1(c). Since
the four layers of fermions are decoupled at the free fermion level, we are
free to assign them with the 3,4,5,0 chiral $\mathrm{U}(1)$ charges
respectively, such that the low-energy edge modes realize the 3-4-5-0 chiral
fermions and their mirror partners. We treat the edge A as the light (chiral
fermion) sector, and the edge B as the mirror sector (to be gapped out). If we
can generate a mass gap for the edge B fermions only without breaking the
chiral $\mathrm{U}(1)$ symmetry, we will succeed in achieving a lattice
regularization of the 3-4-5-0 field theory Eq. (1) in this (1+1)D system in
terms of the gapless edge A fermions.

The fact that the $\mathrm{U}(1)$ ’t Hooft anomaly vanishes for the 3-4-5-0
model indicates that it should be possible to gap out the edge B fermions
trivially without breaking the chiral $\mathrm{U}(1)$ symmetry. However, the
chiral $\mathrm{U}(1)$ symmetry is restrictive enough to prevent the gapping
to happen on the free-fermion level. Therefore, we resort to the idea of
gapping out the mirror fermions by interactions, which has been previously
explored by Chen, Giedt and Poppitz (CGP) [18] in the 3-4-5-0 lattice model,
where all $\mathrm{U}(1)$ symmetry allowed interactions are included. However,
the CGP result shows a singular non-local behavior for the gauge field
polarization tensor in the mirror sector, which indicates the mirror sector
still has surviving gapless modes charged under the gauge field. The reason
could be that the CGP approach introduces too many interaction terms, and some
of them are harmful. In order to achieve the SMG, the fermion interaction must
be carefully selected to satisfy the gapping condition (i.e. the interaction
operators must be self-bosonic and mutual-bosonic in terms of the operator
braiding statistics [60, 61, 62, 63]), as elaborated in recent works [26, 33].
It turns out that the lowest order interactions that satisfy the gapping
condition are the following six-fermion local interactions [26],

$\begin{split}H_{\text{int}}=\sum_{i\in\text{B}}g_{1}&(\psi_{1,i}\psi_{2,i}^{\dagger}\psi_{2,i+1}^{\dagger}\psi_{3,i}\psi_{4,i}\psi_{4,i+1}+\text{h.c.})\\\
+g_{2}&(\psi_{1,i}\psi_{1,i+1}\psi_{2,i}\psi_{3,i}^{\dagger}\psi_{3,i+1}^{\dagger}\psi_{4,i}+\text{h.c.}),\end{split}$
(3)

which are seemly irrelevant dimension-5 operators in the perturbative RG
around the gapless free fermion fixed point. The interaction respects the
chiral $\mathrm{U}(1)$ symmetry, and is only applied to sites on the B edge
(denoted as $i\in\text{B}$), with $i+1$ being the next site of $i$ along the B
edge. The interaction looks highly irrelevant in the free-fermion limit.
However, strong enough interaction (strong in the sense that the interaction
energy scale $E_{\rm int}$ is large but still in the same order of magnitude
as the kinetic energy $E_{\rm free}$, thus $E_{\rm int}/E_{\rm
free}\simeq\mathcal{O}(1)$ is nonperturbative) may still generate non-
perturbative effect that gaps out the edge B fermions. Our central goal is to
numerically verify that the proposed interaction Eq. (3) indeed drives the SMG
in and only in the mirror sector.

Figure 2: (Color online.) Ground state (GS) energy per unit cell as a function
of interaction strength $g$. The inset shows the first-order derivative of the
GS energy with respect to $g$. The features around $g\approx 5.7$ (indicated
by the grey dashed line) signal a quantum phase transition.

_DMRG results._ — We study the lattice model
$H=H_{\text{free}}+H_{\text{int}}$ by the DMRG method [57] using the ITensor
software library [64]. For simplicity, we set $g_{1}=g_{2}=g$ as the only
interaction parameter. The simulation is performed on a two-leg ladder of 20
unit cells, where three different matrix product state (MPS) bond dimensions
$\mathcal{D}=6000,7000,8000$ are used (the code to reproduce our result is
available at [65]). Computed physical quantities are then extrapolated to the
$\mathcal{D}\to\infty$ limit assuming a $1/\mathcal{D}$ scaling. Figure 2
shows the ground state energy $E_{\text{GS}}$ (of the full Hamiltonian $H$)
per unit cell as a function of the interaction strength $g$, where the inset
shows its first-order derivative $\partial_{g}E_{\text{GS}}$. The onset of a
non-zero $\partial_{g}E_{\text{GS}}=g^{-1}\langle H_{\text{int}}\rangle$
around $g_{c}\approx 5.7$ signifies the development of the $\langle
H_{\text{int}}\rangle\neq 0$ condensation across the SMG transition. The
smooth kink of $\partial_{g}E_{\text{GS}}$ indicates a (high-order) continuous
transition (in fact, as to be shown soon, this should be an infinite-order
transition with no singularity in $E_{\text{GS}}$ and its derivatives to any
order).

To further confirm the existence of the critical point $g_{c}$, we calculate
the fermion correlation functions
$C_{\psi}(r)\equiv\langle\psi_{I,i+r}^{\dagger}\psi_{I,i}\rangle$ on both edge
A and edge B across the transition. It turns out that the behavior of
$C_{\psi}$ is the same for all $I=1,2,3,4$, such that it is sufficient to show
one of the four flavors. Figures 3 (a,c) and Figs. 3 (b,d) show the
correlation functions for each edge before and after the transition,
respectively. We observe that edge A is always gapless with power-law
correlations. In contrast, edge B is gapless when $g<g_{c}$ but becomes gapped
with an exponential-decay correlation when $g>g_{c}$. The two qualitatively
different behaviors must be separated by a quantum phase transition.

Figure 3: (Color online.) Correlations on both edges before and after
transition. Linear fit (red line) is performed for intermediate distances from
$r=2$ to $r=6$ in each case, in order to faithfully extract the low energy
physics while avoiding the artifacts due to the gap caused by finite bond
dimension in the matrix product state representation. (a) $g=5.0<g_{c}$ for
edge A. The log-log plot shows a power-law decay for intermediate distances.
(b) $g=7.0>g_{c}$ for edge A. The log-log plot again shows a power-law decay.
(c) $g=5.0<g_{c}$ for edge B. The log-log plot shows a power-law decay. (d)
$g=7.0>g_{c}$ for edge B. The semi-log plot indicates an exponential decay,
i.e. edge B becomes gapped. Figure 4: (Color online.) (a) The evolution of
fermion scaling dimension $\Delta_{\psi}$ on both edges as the interaction
strength $g$ approaches the critical point. The scaling dimension is obtained
from the power-law fitting as in Fig. 3. The horizontal dashed line indicates
the free fermion limit. The gray stripe shows the estimated critical
interaction strength $g_{c}$ with some uncertainty. (b) The solved scaling
dimension for the interaction terms on edge B based on the scaling dimensions
of multiple operators (Refer to Appendix C for details). The horizontal dashed
line indicates the marginal value 2 of $\Delta_{\text{int}}$, across which the
phase transition is expected to happen.

By fitting the power-law correlation function $C_{\psi}(r)\sim
1/r^{2\Delta_{\psi}}$ before the transition, we can extract the fermion
scaling dimension $\Delta_{\psi}$ on both edges. The result is shown in Fig.
4(a). In the free-fermion limit ($g=0$), the fermion scaling dimension is
$\Delta_{\psi}=\frac{1}{2}$ on both edges. The finite-size effect tends to
reduce the scaling dimension slightly. A finite-size scaling of the scaling
dimension in the free fermion limit is performed in Appendix A, confirming
that our result converges to the long-distance limit correctly. As $g$
increases towards $g_{c}$, the fermion scaling dimension on the edge B
increases continuously from $\frac{1}{2}$ to about $0.67$ (near $g_{c}$),
indicating that fermion operators get renormalized by the interaction
significantly. For $g>g_{c}$, the correlation on the edge B becomes short-
ranged, such that the fermion scaling dimension is no longer defined (although
the power-law fitting on the finite-size data will continue give some
estimated exponent that extrapolates beyond the critical point before the
correlation length shrinks below the system size). However, on the edge A, the
fermion scaling dimension, while experiencing some fluctuations near the
critical point, generally stays close to the free fermion limit regardless of
the interaction strength. The scaling dimension remains stable even after $g$
goes across the transition point $g_{c}$ by a significant amount. This implies
that the edge A remains gapless and almost free, as the edge B interaction can
only induce a perturbative interaction on the edge A through the proximity
effect.

To verify that the chiral $\mathrm{U}(1)$ symmetry is not broken spontaneously
by the condensation of fermion bilinear masses, we measure correlation
functions of Dirac and Majorana mass operators on the B edge, i.e.
$C_{\psi_{I}^{\dagger}\psi_{J}}(r)\equiv\langle\psi_{J,i+r}^{\dagger}\psi_{I,i+r}\psi_{I,i}^{\dagger}\psi_{J,i}\rangle$
and
$C_{\psi_{I}\psi_{J}}(r)\equiv\langle\psi_{J,i+r}^{\dagger}\psi_{I,i+r}^{\dagger}\psi_{I,i}\psi_{J,i}\rangle$.
Fig. 6 in the Appendix B shows the correlations for all the eight mass terms
are short-ranged (exponential decay) along the B edge in the strong coupling
gapped phase ($g>g_{c}$), which confirms that the mirror fermions on the B
edge are gapped by the SMG mechanism without long-range ordering of bilinear
masses. Therefore, the remaining gapless fermions on the A edge successfully
realize the lattice regularization of chiral fermions in the 3-4-5-0 model
preserving the chiral $\mathrm{U}(1)$ symmetry.

_Luttinger liquid RG analysis._ — To better understand the nature of the SMG
transition at $g_{c}$, we perform the Luttinger liquid RG analysis for the
edge B fermions. We first bosonize the mirror fermions by
$\psi_{I}\sim\mathrm{e}^{\mathrm{i}\varphi_{I}}$. Then the (1+1)D interacting
mirror fermions can be described by the Luttinger liquid effective field
theory in terms of the
$\varphi=(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})^{\intercal}$ fields

$\mathcal{L}=\frac{1}{4\pi}(\partial_{t}\varphi^{\intercal}K\partial_{x}\varphi+\partial_{x}\varphi^{\intercal}V\partial_{x}\varphi)+\sum_{\alpha=1,2}g_{\alpha}\cos(l_{\alpha}^{\intercal}\varphi),$
(4)

where $K=\sigma^{30}$ and $V=\sigma^{00}$ (in the UV limit) are $4\times 4$
matrixes (where $\sigma^{\mu\nu}=\sigma^{\mu}\otimes\sigma^{\nu}$ denotes the
tensor product of Pauli matrices). The two interaction terms $g_{1},g_{2}$ in
Eq. (3) correspond to the cosine-terms in Eq. (4) specified by the vectors
$l_{1}=(1,-2,1,2)^{\intercal}$ and $l_{2}=(2,1,-2,1)^{\intercal}$,
respectively. The RG flow with respect to the log-energy-scale
$\ell=-\ln\Lambda$ is given by [66]

$\begin{split}\frac{\mathrm{d}g_{\alpha}}{\mathrm{d}\ell}&=(2-\Delta_{l_{\alpha}})g_{\alpha}-\frac{1}{2}\sum_{l_{\beta}\pm
l_{\gamma}=l_{\alpha}}g_{\beta}g_{\gamma},\\\
\frac{\mathrm{d}V^{-1}}{\mathrm{d}\ell}&=\frac{1}{2}\sum_{\alpha}g_{\alpha}^{2}(K^{-1}l_{\alpha}l_{\alpha}^{\intercal}K^{-1}-V^{-1}l_{\alpha}l_{\alpha}^{\intercal}V^{-1}),\end{split}$
(5)

where $\Delta_{l}=\frac{1}{2}l^{\intercal}V^{-1}l$ denotes the scaling
dimension of the vertex operator
$\mathrm{e}^{\mathrm{i}l^{\intercal}\varphi}$. Under the RG flow, the $V$
matrix gets renormalized to the general form

$V=\sqrt{1+y_{1}^{2}+y_{2}^{2}}\;\sigma^{00}-y_{1}\,\sigma^{10}-y_{2}\,\sigma^{22},$
(6)

where $y_{1},y_{2}$ are Luttinger parameters that depend on the RG scale
$\ell$. Tab. 1 concludes the analytical expressions of the scaling dimensions
$\Delta_{l}$ of the fermion, Dirac mass and Majorana mass operators given
their corresponding charge vectors $l$. We numerically determine the scaling
dimensions of these operators before the transition $(g<g_{c})$, by fitting
the power-law exponents of their correlation functions (see Appendix. C for
details).

Table 1: Operator scaling dimensions and Luttinger parameters at the free-fermion limit ($g=0$) and at the critical point ($g=g_{c}$). | $l^{\intercal}$ | $\Delta_{l}$ | free | critical
---|---|---|---|---
$\psi_{1}$ | $(1,0,0,0)$ | $\frac{1}{2}\sqrt{1+y_{1}^{2}+y_{2}^{2}}$ | $\frac{1}{2}$ | $0.67\pm 0.07$
$\psi_{1}^{\dagger}\psi_{3}$ | $(-1,0,1,0)$ | $\sqrt{1+y_{1}^{2}+y_{2}^{2}}-y_{1}$ | $1$ | $0.76\pm 0.05$
$\psi_{1}\psi_{4}$ | $(1,0,0,1)$ | $\sqrt{1+y_{1}^{2}+y_{2}^{2}}-y_{2}$ | $1$ | $0.73\pm 0.01$
|  | $y_{1}$ | 0 | $0.72\pm 0.18$
|  | $y_{2}$ | 0 | $0.70\pm 0.15$

From the scaling dimensions, we infer the Luttinger parameters $y_{1},y_{2}$,
and calculate the scaling dimension of the interaction operator
$\Delta_{\text{int}}:=\Delta_{l_{1}}=\Delta_{l_{2}}=\sqrt{1+y_{1}^{2}+y_{2}^{2}}-3y_{1}-4y_{2}$.
The evolution of $\Delta_{\text{int}}$ is shown in Fig. 4(b), which drops
continuously from $\Delta_{\text{int}}=5$ at the free-fermion limit ($g=0$) to
$2.17\pm 0.27$ at the SMG transition ($g=g_{c}$). Although the interaction is
perturbatively irrelevant at the free-fermion fixed point, finite strength of
the interaction can renormalize the Luttinger parameter, which reduces its own
scaling dimension. Our numerical result indicates that the SMG transition is
triggered exactly when the interaction scaling dimension is reduced to
marginal $\Delta_{\text{int}}=2$, which matches the mechanism of the
Berezinskii-Kosterlitz-Thouless (BKT) transition. This scenario was also
proposed by Tong in a recent theoretical study [55]. Our numerical study
provides more detailed RG analysis and more solid evidence in support of the
BKT transition scenario.

_Conclusion and Discussions._ — We numerically demonstrate the lattice
regularization of (1+1)D chiral fermions in the 3-4-5-0 model. This is
achieved by gapping out the anomaly-free mirror sector using properly designed
interactions via the SMG mechanism, leaving the light sector gapless. By
simulating the lattice model with the DMRG method, we identify the SMG
transition point $g_{c}$. In the strong coupling phase ($g=g_{c}$), we show
that the mirror fermions are gapped without breaking the chiral symmetry, and
the light fermions remain gapless. We numerically determine the scaling
dimension of the interaction operator before the transition, which evolves
continuously from irrelevant to marginal. This behavior clearly indicates the
BKT nature of the SMG transition in our model. Once the anomaly-free U(1)
symmetry is dynamically gauged, we expect to obtain a (1+1)D lattice chiral
gauge theory coupled to chiral fermions, which could potentially be simulated
by the quantum Monte Carlo method [67], as our proposed six-fermion
interaction in Eq. (3) admits the following Yukawa decomposition (with site
indices omitted for brevity)

$\begin{split}H_{\text{Yuk}}=&(\phi_{1}^{2}\psi_{1}\psi_{3}+\phi_{1}^{\dagger}\psi_{2}^{\dagger}\psi_{4}+\text{h.c.})+\frac{1}{\tilde{g}_{1}}\phi_{1}^{\dagger}\phi_{1}\\\
+&(\phi_{2}^{2}\psi_{2}\psi_{4}+\phi_{2}^{\dagger}\psi_{1}\psi_{3}^{\dagger}+\text{h.c.})+\frac{1}{\tilde{g}_{2}}\phi_{2}^{\dagger}\phi_{2},\end{split}$
(7)

such that integrating out the Yukawa bosons $\phi_{\alpha}$ reproduces our
interaction at the leading order of $g_{\alpha}\sim\tilde{g}_{\alpha}^{2}$.
Based on the equivalence between the U(1) anomaly-free and gapping conditions
in (1+1)D [26, 33], hopefully our work can prompt future simulations on other
(1+1)D lattice chiral fermion/gauge theory models.

Acknowledgement. — MZ would like to thank Miles Stoudenmire for helpful
discussions on using the ITensor package. ZZ is supported by the National
Natural Science Foundation of China (Grant No. 12074375) and the start-up
funding of KITS at UCAS (Grant No.118900M026). JW is supported by Harvard
CMSA. MZ and YZY are supported by a start-up funding of UCSD.

## References

  * Nielsen and Ninomiya [1981] H. B. Nielsen and M. Ninomiya,  _No-go theorum for regularizing chiral fermions_ , Tech. Rep. (Science Research Council, 1981).
  * Ginsparg and Wilson [1982] P. H. Ginsparg and K. G. Wilson, Physical Review D 25, 2649 (1982).
  * Swift [1984] P. Swift, Physics Letters B 145, 256 (1984).
  * Smit [1986] J. Smit, Acta Phys. Polon. B 17, 531 (1986).
  * Banks and Dabholkar [1992] T. Banks and A. Dabholkar, Phys. Rev. D 46, 4016 (1992), arXiv:hep-lat/9204017 [hep-lat] .
  * Kaplan [1992] D. B. Kaplan, Physics Letters B 288, 342 (1992).
  * Shamir [1993] Y. Shamir, Nuclear Physics B 406, 90 (1993).
  * Eichten and Preskill [1986] E. Eichten and J. Preskill, Nuclear Physics B 268, 179 (1986).
  * Bock and De [1990] W. Bock and A. K. De, Physics Letters B 245, 207 (1990).
  * Bock  _et al._ [1992] W. Bock, A. K. De, and J. Smit, Nuclear Physics B 388, 243 (1992).
  * Bock  _et al._ [1993] W. Bock, A. K. De, E. Focht, and J. Smit, Nuclear Physics B 401, 481 (1993), arXiv:hep-lat/9210022 [hep-lat] .
  * Bock  _et al._ [1994] W. Bock, J. Smit, and J. C. Vink, Nuclear Physics B 414, 73 (1994), arXiv:hep-lat/9306012 [hep-lat] .
  * Hasenfratz  _et al._ [1990] A. Hasenfratz, W. Liu, and T. Neuhaus, Physics Letters B 236, 339 (1990).
  * Hasenfratz  _et al._ [1991] A. Hasenfratz, P. Hasenfratz, K. Jansen, J. Kuti, and Y. Shen, Nuclear Physics B 365, 79 (1991).
  * Lee  _et al._ [1990] I.-H. Lee, J. Shigemitsu, and R. E. Shrock, Nuclear Physics B 334, 265 (1990).
  * Montvay [1992] I. Montvay, Nuclear Physics B Proceedings Supplements 29, 159 (1992), arXiv:hep-lat/9205023 [hep-lat] .
  * Golterman  _et al._ [1993] M. F. L. Golterman, D. N. Petcher, and E. Rivas, Nuclear Physics B 395, 596 (1993), arXiv:hep-lat/9206010 [hep-lat] .
  * Chen  _et al._ [2013] C. Chen, J. Giedt, and E. Poppitz, Journal of High Energy Physics 2013, 1 (2013).
  * Fidkowski and Kitaev [2010] L. Fidkowski and A. Kitaev, Phys. Rev. B 81, 134509 (2010), arXiv:0904.2197 [cond-mat.str-el] .
  * Fidkowski and Kitaev [2011] L. Fidkowski and A. Kitaev, Phys. Rev. B 83, 075103 (2011), arXiv:1008.4138 [cond-mat.str-el] .
  * Turner  _et al._ [2011] A. M. Turner, F. Pollmann, and E. Berg, Phys. Rev. B 83, 075102 (2011), arXiv:1008.4346 [cond-mat.str-el] .
  * Ryu and Zhang [2012] S. Ryu and S.-C. Zhang, Phys. Rev. B 85, 245132 (2012), arXiv:1202.4484 [cond-mat.str-el] .
  * Qi [2013] X.-L. Qi, New Journal of Physics 15, 065002 (2013), arXiv:1202.3983 [cond-mat.str-el] .
  * Yao and Ryu [2013] H. Yao and S. Ryu, Phys. Rev. B 88, 064507 (2013), arXiv:1202.5805 [cond-mat.str-el] .
  * Wen [2013] X.-G. Wen, Chin. Phys. Lett. 30, 111101 (2013), arXiv:1305.1045 [hep-lat] .
  * Wang and Wen [2013] J. Wang and X.-G. Wen, arXiv e-prints , arXiv:1307.7480 (2013), arXiv:1307.7480 [hep-lat] .
  * Wang and Senthil [2014] C. Wang and T. Senthil, Phys. Rev. B 89, 195124 (2014), arXiv:1401.1142 [cond-mat.str-el] .
  * Gu and Levin [2014] Z.-C. Gu and M. Levin, Phys. Rev. B 89, 201113 (2014), arXiv:1304.4569 [cond-mat.str-el] .
  * Metlitski  _et al._ [2014] M. A. Metlitski, L. Fidkowski, X. Chen, and A. Vishwanath, arXiv e-prints , arXiv:1406.3032 (2014), arXiv:1406.3032 [cond-mat.str-el] .
  * You and Xu [2014] Y.-Z. You and C. Xu, Phys. Rev. B 90, 245120 (2014), arXiv:1409.0168 [cond-mat.str-el] .
  * Witten [2016] E. Witten, Phys. Rev. B 94, 195150 (2016), arXiv:1605.02391 [hep-th] .
  * Kikukawa [2019] Y. Kikukawa, PTEP 2019, 073B02 (2019), arXiv:1710.11101 [hep-lat] .
  * Wang and Wen [2018] J. Wang and X.-G. Wen, Phys. Rev. D 99, 111501 (2018), arXiv:1807.05998 [hep-lat] .
  * Wang and Wen [2020] J. Wang and X.-G. Wen, Phys. Rev. Res. 2, 023356 (2020), arXiv:1809.11171 [hep-th] .
  * Ayyar and Chandrasekharan [2015] V. Ayyar and S. Chandrasekharan, Phys. Rev. D 91, 065035 (2015), arXiv:1410.6474 [hep-lat] .
  * Slagle  _et al._ [2015] K. Slagle, Y.-Z. You, and C. Xu, Phys. Rev. B 91, 115121 (2015), arXiv:1409.7401 [cond-mat.str-el] .
  * Catterall [2016] S. Catterall, Journal of High Energy Physics 1, 121 (2016), arXiv:1510.04153 [hep-lat] .
  * Ayyar and Chandrasekharan [2016a] V. Ayyar and S. Chandrasekharan, Phys. Rev. D 93, 081701 (2016a), arXiv:1511.09071 [hep-lat] .
  * Catterall and Schaich [2016] S. Catterall and D. Schaich, ArXiv e-prints (2016), arXiv:1609.08541 [hep-lat] .
  * Ayyar and Chandrasekharan [2016b] V. Ayyar and S. Chandrasekharan, Journal of High Energy Physics 10, 58 (2016b), arXiv:1606.06312 [hep-lat] .
  * Ayyar [2016] V. Ayyar, ArXiv e-prints (2016), arXiv:1611.00280 [hep-lat] .
  * He  _et al._ [2016] Y.-Y. He, H.-Q. Wu, Y.-Z. You, C. Xu, Z. Y. Meng, and Z.-Y. Lu, Phys. Rev. B 94, 241111 (2016), arXiv:1603.08376 [cond-mat.str-el] .
  * DeMarco and Wen [2017] M. DeMarco and X.-G. Wen, arXiv e-prints , arXiv:1706.04648 (2017), arXiv:1706.04648 [hep-lat] .
  * Ayyar and Chandrasekharan [2017] V. Ayyar and S. Chandrasekharan, Phys. Rev. D 96, 114506 (2017), arXiv:1709.06048 [hep-lat] .
  * Schaich and Catterall [2018] D. Schaich and S. Catterall, in  _European Physical Journal Web of Conferences_, European Physical Journal Web of Conferences, Vol. 175 (2018) p. 03004, arXiv:1710.08137 [hep-lat] .
  * Catterall and Butt [2018] S. Catterall and N. Butt, Phys. Rev. D 97, 094502 (2018), arXiv:1708.06715 [hep-lat] .
  * Butt and Catterall [2018] N. Butt and S. Catterall, in  _The 36th Annual International Symposium on Lattice Field Theory. 22-28 July_ (2018) p. 294, arXiv:1811.01015 [hep-lat] .
  * Butt  _et al._ [2018] N. Butt, S. Catterall, and D. Schaich, Phys. Rev. D 98, 114514 (2018), arXiv:1810.06117 [hep-lat] .
  * Catterall  _et al._ [2020] S. Catterall, N. Butt, and D. Schaich, arXiv e-prints , arXiv:2002.00034 (2020), arXiv:2002.00034 [hep-lat] .
  * Catterall [2021] S. Catterall, Phys. Rev. D 104, 014503 (2021), arXiv:2010.02290 [hep-lat] .
  * Butt  _et al._ [2021a] N. Butt, S. Catterall, and G. C. Toga, arXiv e-prints , arXiv:2111.01001 (2021a), arXiv:2111.01001 [hep-lat] .
  * You  _et al._ [2018a] Y.-Z. You, Y.-C. He, C. Xu, and A. Vishwanath, Physical Review X 8, 011026 (2018a), arXiv:1705.09313 [cond-mat.str-el] .
  * You  _et al._ [2018b] Y.-Z. You, Y.-C. He, A. Vishwanath, and C. Xu, Phys. Rev. B 97, 125112 (2018b), arXiv:1711.00863 [cond-mat.str-el] .
  * Razamat and Tong [2021] S. S. Razamat and D. Tong, Physical Review X 11, 011063 (2021), arXiv:2009.05037 [hep-th] .
  * Tong [2021] D. Tong, arXiv e-prints , arXiv:2104.03997 (2021), arXiv:2104.03997 [hep-th] .
  * Butt  _et al._ [2021b] N. Butt, S. Catterall, A. Pradhan, and G. C. Toga, Phys. Rev. D 104, 094504 (2021b), arXiv:2101.01026 [hep-th] .
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck [2005] U. Schollwöck, Reviews of Modern Physics 77, 259 (2005), arXiv:cond-mat/0409292 [cond-mat.str-el] .
  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Haldane [1995] F. D. M. Haldane, Phys. Rev. Lett. 74, 2090 (1995), arXiv:cond-mat/9501007 .
  * Kapustin and Saulina [2011] A. Kapustin and N. Saulina, Nucl. Phys. B 845, 393 (2011), arXiv:1008.0654 [hep-th] .
  * Wang and Wen [2015] J. Wang and X.-G. Wen, Phys. Rev. B 91, 125124 (2015), arXiv:1212.4863 [cond-mat.str-el] .
  * Levin [2013] M. Levin, Phys. Rev. X 3, 021009 (2013), arXiv:1301.7355 [cond-mat.str-el] .
  * Fishman  _et al._ [2020] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations (2020), arXiv:2007.14822 .
  * [65] https://github.com/meng-zeng/julia-code-for-chiral-fermion-dmrg.
  * Katzir  _et al._ [2020] B. A. Katzir, A. Stern, E. Berg, and N. H. Lindner, arXiv e-prints , arXiv:2011.13950 (2020), arXiv:2011.13950 [cond-mat.mes-hall] .
  * Chen  _et al._ [2013] C. Chen, J. Giedt, and E. Poppitz, Journal of High Energy Physics 2013, 131 (2013), arXiv:1211.6947 [hep-lat] .

##  Appendix A Scaling dimensions in free fermion limit

In the free fermion limit, correlation functions and consequently scaling
dimensions can be calculated analytically. Here we work out the free fermion
case on the lattice as a benchmark to the DMRG calculation at $g=0$.

###  A.1 Fermion scaling dimension

For the free fermion-fermion correlation, we can just pick the fermion flavor
$f_{3}$ without loss of generality. The two-point correlation is given by
(with the flavor index ignored)

$C_{\psi}(r)=\langle\psi_{i+r}^{\dagger}\psi_{i}\rangle\sim r^{-\nu}.$ (8)

On a finite lattice with real space Hamiltonian $\mathcal{H}$, we can do a
change basis from the fermionic operators in real space fermionic operators in
energy eigenspace, so that the correlation at half-filling can be more
conveniently calculated. Assuming $\mathcal{H}$ has eigenstates
$|\epsilon_{n}\rangle$ with corresponding eigenenergies $\epsilon_{n}$, then
the change of basis is given by
$\psi^{\dagger}(i)=\sum_{n}\langle\epsilon_{n}|i\rangle\psi_{n}^{\dagger}$.
Eventually the correlation becomes

$C_{\psi}(r)=\sum_{n,m}\langle\epsilon_{n}|i\rangle\langle
j|\epsilon_{m}\rangle\langle\psi_{n}^{\dagger}\psi_{m}\rangle.$ (9)

Since the ground state is half-filled, the above summation in $n$ or $m$ is
only over the lower half of the energy spectrum. The system we have in the
main text consists of 20 unit cells. Exact diagonalization can be done in the
free fermion limit and the correlation can then be calculated using Eq. (9).
With open boundary condition (same as the DMRG setup), the correlation with
log-log scale is plotted in Fig. 5(a) together with a linear fit. The power
law exponent obtained is less than 1, which explains the deviation of the
scaling dimension $\Delta_{\psi}$ from 0.5 for $g=0$ in Fig. 4. This deviation
is mainly a finite-size effect, demonstrated in Fig. 5(b). In the large-
system-size limit, the free fermion scaling dimension of 0.5 is recovered.

Figure 5: (a) Linear fit for the correlation function on a log-log scale when
the system has 20 unit cells. The power law exponent obtained is around 0.936,
which is smaller than 1. (b) Finite-size scaling for the exponent using a
polynomial function for system sizes $L=20,40,80,160,320$. The extrapolation
to $L=\infty$ recovers the ideal $\nu=1$ limit.

###  A.2 Mass term scaling dimensions

Other than the fermion-fermion correlation, correlations of various mass terms
on a finite lattice can also be calculated analytically in the free fermion
limit. In this case, the expected scaling dimension for a bilinear in the
large-system-size limit is 1. Taking the fermion bilinear
$\psi_{1}\psi_{3}^{\dagger}$ as an example, the correlation on the lattice is
given by (using Wick’s theorem)

$\begin{split}C_{\psi_{1}^{\dagger}\psi_{3}}(r)&=\langle\psi_{1,i+r}^{\dagger}\psi_{3,i+r}\psi_{3,i}^{\dagger}\psi_{1,i}\rangle\\\
&=\langle\psi_{1,i+r}^{\dagger}\psi_{1,i}\rangle\langle\psi_{3,i+r}\psi_{3,i}^{\dagger}\rangle\\\
&\sim r^{-2\nu},\end{split}$ (10)

i.e., the fermion bilinear exponent is simply double of that for the single
fermion, as expected. Therefore, the finite-size behavior should also be the
same.

##  Appendix B Mass term correlations in the gapped phase for edge B

In this section, we present the correlations for the mass terms on edge B
after the gapping transition to demonstrate that the U(1) chiral symmetry is
preserved in the gapped phase. The correlations for the Dirac masses are shown
in Fig. 6(a) and the correlations for the Majorana masses are shown in Fig.
6(b). We see clear evidence for exponential decays for all the eight different
mass terms at relatively shorter length scales $r\lesssim 5$, where the
correlation is expected to be dominated by the SMG gap. The non-monotonic
behavior of the correlation function for larger distance $r\gtrsim 5$ is an
artifact arising from the finite MPS bond dimension, and should not be
trusted. Different interaction strengths are chosen for the different mass
terms in order to better demonstrate the exponential decay features. Thus, we
conclude that the chiral U(1) symmetry is preserved in the gapped phase on
edge B.

Figure 6: Semi-log correlations for the various mass terms in the gapped
phase. (a) Dirac masses $\psi_{I}\psi_{J}^{\dagger}$. The interaction
strengths chosen for the four mass terms $\psi_{1}\psi_{3}^{\dagger}$,
$\psi_{1}\psi_{4}^{\dagger}$, $\psi_{2}\psi_{3}^{\dagger}$ and
$\psi_{2}\psi_{4}^{\dagger}$ in the gapped phase are $6.5$, $6.5$, $6.6$ and
$6.4$ respectively. (b) Majorana masses $\psi_{I}\psi_{J}$. The interaction
strengths chosen for the four mass terms $\psi_{1}\psi_{3}$,
$\psi_{1}\psi_{4}$, $\psi_{2}\psi_{3}$ and $\psi_{2}\psi_{4}$ in the gapped
phase are $6.4$, $6.3$, $6.3$ and $7.2$ respectively.

##  Appendix C Scaling dimensions under RG

The scaling dimensions of some of the interesting operators are summarized in
Table 2.

Operators | $l$ | $\Delta_{l}$
---|---|---
$\psi_{1}$ | $(1,0,0,0)$ | $\frac{1}{2}\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{1}\psi_{3}^{\dagger}$ | $(1,0,-1,0)$ | $-y_{1}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{1}\psi_{4}^{\dagger}$ | $(1,0,0,-1)$ | $y_{2}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{2}\psi_{3}^{\dagger}$ | $(0,1,-1,0)$ | $-y_{2}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{2}\psi_{4}^{\dagger}$ | $(0,1,0,-1)$ | $-y_{1}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{1}\psi_{3}$ | $(1,0,1,0)$ | $y_{1}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{1}\psi_{4}$ | $(1,0,0,1)$ | $-y_{2}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{2}\psi_{3}$ | $(0,1,1,0)$ | $y_{2}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{2}\psi_{4}$ | $(0,1,0,1)$ | $y_{1}+\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{1}\psi_{3}\psi_{2}^{\dagger}\partial_{x}\psi_{2}^{\dagger}\psi_{4}\partial_{x}\psi_{4}$ | $(1,-2,1,2)$ | $-3y_{1}-4y_{2}+5\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
$\psi_{1}\partial_{x}\psi_{1}\psi_{3}^{\dagger}\partial_{x}\psi_{3}^{\dagger}\psi_{2}\psi_{4}$ | $(2,1,-2,1)$ | $-3y_{1}-4y_{2}+5\sqrt{1+y_{1}^{2}+y_{2}^{2}}$
Table 2: Scaling dimensions of the various operators, some of which are used
to solve the two parameters $y_{1}$ and $y_{2}$.

These scaling dimensions can in principle be measured by calculating the
power-law correlation using DMRG. Here we are interested in the scaling
dimensions of the gapping terms under RG before the Luttinger liquid becomes
gapped. However, for large scaling dimensions, the power-law decay is too fast
to be measured accurately. Therefore, in order to solve for the parameters
$y_{1}$ and $y_{2}$, we only use the correlations for operators with
relatively smaller scaling dimensions, which, based on Table 2, include
$\psi_{1},\psi_{1}\psi_{3}^{\dagger},\psi_{2}\psi_{3}^{\dagger},\psi_{2}\psi_{4}^{\dagger}$
and $\psi_{1}\psi_{4}$. Fig. 7 shows the evolution of the scaling dimensions
for some of the mass terms under RG flow.

Figure 7: Evolution of scaling dimensions for some of the mass terms under RG
on edge B. The fermion scaling dimension $\Delta_{\psi}$ is also shown for
comparison. Figure 8: Solutions to the Luttinger parameters $(y_{1},y_{2})$
based on RG equation and some of the measured scaling dimensions using DMRG.

Instead of using two of the different scaling dimensions to solve for the two
parameters, we try to make use of all the five aforementioned operators to
have a faithful representation of the available DMRG data. For a particular
interaction strength $g$, denote the five scaling dimensions from RG
calculation by $\Delta_{i}^{(\text{RG})}(y_{1}(g),y_{2}(g))$ and the DMRG
counterparts by $\Delta_{i}^{(\text{DMRG})}(g)$ with errors
$\delta\Delta_{i}^{(\text{DMRG})}(g)$. Then $y_{1}$ and $y_{2}$ are solved by
minimizing the following error function at each $g$:

$f_{\text{error}}(g)=\sum_{i}\frac{\left(\Delta_{i}^{(\text{RG})}(y_{1}(g),y_{2}(g))-\Delta_{i}^{(\text{DMRG})}(g)\right)^{2}}{\delta\Delta_{i}^{(\text{DMRG})}(g)^{2}},$
(11)

where $\frac{1}{\delta\Delta_{i}^{(\text{DMRG})}(g)^{2}}$ can be considered as
the weight of each contribution to the total error function. The solutions
$(y_{1},y_{2})$ with error bars are obtained in the following way. At each
$g$, a numerical value for the DMRG scaling dimension is randomly drawn from
the interval
$[\Delta_{i}^{(\text{DMRG})}-\delta\Delta_{i}^{(\text{DMRG})},\Delta_{i}^{(\text{DMRG})}+\delta\Delta_{i}^{(\text{DMRG})}]$,
then the error function $f_{\text{error}}$ is minimized to find
$(y_{1},y_{2})$, with the weight of this particular solution given by
$\frac{1}{f_{\text{error}}^{2}}$. This process is repeated 100 times for each
$g$, after which weighted average is taken to obtain the mean of
$(y_{1},y_{2})$ with the errors given by the weighted uncertainty. The
evolution of the two parameters with RG flow is shown in Fig. 8. The solved
$y_{1}$ and $y_{2}$ can then be used to calculate the scaling dimension of the
gapping terms $\Delta_{\text{int}}$ plotted in Fig. 4(b) in the main text.


# A Systematic Construction Approach for All $4\times 4$ Involutory MDS
Matrices

Yogesh Kumar111Scientific Analysis Group, DRDO, Metcalfe House Complex,
Delhi-110054
email<EMAIL_ADDRESS>P.R.Mishra222Scientific Analysis Group, DRDO,
Metcalfe House Complex, Delhi-110054
email<EMAIL_ADDRESS>Susanta Samanta333R. C. Bose Centre for
Cryptology and Security, Indian Statistical Institute, Kolkata-700108
email<EMAIL_ADDRESS>Atul Gaur444Department of Mathematics,
University of Delhi, Delhi-110007
email<EMAIL_ADDRESS>
###### Abstract

Maximum distance separable (MDS) matrices play a crucial role not only in
coding theory but also in the design of block ciphers and hash functions. Of
particular interest are involutory MDS matrices, which facilitate the use of a
single circuit for both encryption and decryption in hardware implementations.
In this article, we present several characterizations of involutory MDS
matrices of even order. Additionally, we introduce a new matrix form for
obtaining all involutory MDS matrices of even order and compare it with other
matrix forms available in the literature. We then propose a technique to
systematically construct all $4\times 4$ involutory MDS matrices over a finite
field $\mathbb{F}_{2^{m}}$. This method significantly reduces the search space
by focusing on involutory MDS class representative matrices, leading to the
generation of all such matrices within a substantially smaller set compared to
considering all $4\times 4$ involutory matrices. Specifically, our approach
involves searching for these representative matrices within a set of
cardinality $(2^{m}-1)^{5}$. Through this method, we provide an explicit
enumeration of the total number of $4\times 4$ involutory MDS matrices over
$\mathbb{F}_{2^{m}}$ for $m=3,4,\ldots,8$.

Keywords: Diffusion Layer, MDS matrix, Involutory matrix, Finite field

##  1 Introduction

The concepts of confusion and diffusion [19] play a crucial role in the design
of symmetric key cryptographic primitives. In general, the linear
transformation is used to achieve the diffusion property, which can be
represented as a matrix. Additionally, it is common practice to include
maximum distance separable (MDS) matrices in the linear transformations to
achieve optimal diffusion. One application is in the MixColumns operation of
AES [1]. Moreover, MDS matrices find extensive application in numerous other
cryptographic primitives, including stream ciphers such as MUGI [22] and hash
functions like Maelstrom [2], Gr$\phi$stl [3], and PHOTON [4]. Consequently,
the efficacy of MDS matrices within diffusion layers is widely acknowledged,
leading to the proposal of various techniques for their design.

Generally, three techniques are employed for generating MDS matrices. The
first technique, referred to as direct construction, utilizes algebraic
methods to produce an MDS matrix of any order without the necessity of search
procedures. The second method involves search-based construction. The third
approach integrates the strategies mentioned above into a hybrid methodology.
Direct constructions typically derive from Cauchy and Vandermonde matrices,
along with their extensions. Examples from the literature include references
[5, 6, 10, 14, 15].

Although direct construction methods allow for the possibility of obtaining
MDS matrices of any order, there is no assurance of achieving a matrix with
the most efficient hardware area, even for smaller sizes. The second
technique, referred to as the search-based technique, is currently the only
method that can give us the best possible MDS matrix in terms of area.
Nevertheless, this methodology is only viable in situations where the field
size is not too large and the matrix size is small. When it comes to search
techniques, there are different matrix structures to consider. These include
circulant, left-circulant, Hadamard, and Toeplitz matrices. A considerable
amount of progress has been made in this area, as demonstrated in [7, 11, 13,
18, 20].

In the third method, referred to as the hybrid approach, a representative MDS
matrix is usually obtained using a search method. Afterwards, this matrix is
used to generate several MDS matrices. Following a concise overview of prior
research in this field, we will now discuss the motivations of this paper. The
authors in [13] introduced a hybrid method to effectively produce MDS matrices
of size $2^{n}\times 2^{n}$. Their method presents the GHadamard matrix, a new
matrix form that may be employed to create new (involutory) MDS matrices from
a Hadamard (involutory) MDS matrix as a representative matrix. In [16], a
method was presented to produce new $n\times n$ MDS matrices that are
isomorphic to pre-existing ones. The authors in [9] demonstrate a hybrid
approach for generating all MDS matrices, which results in a search space of
$2^{9m}$ to generate all $4\times 4$ involutory MDS matrices. Whereas,
according to [23], the search space for generating all $4\times 4$ involutory
MDS matrices amounts to $2^{8m}$. Furthermore, a recent study [21] proposes a
hybrid method to generate all $4\times 4$ involutory MDS matrices over
$\mathbb{F}_{2^{m}}$, reducing the search space to $(2^{m}-1)^{8}$ to identify
all $4\times 4$ representative involutory MDS matrices. However, the paper
lacks the mathematical background necessary to define a representative
involutory matrix form, as well as the discussion on the existence of a unique
representative involutory matrix for every involutory MDS matrix.

Also, it is important to note that the enumeration of all $4\times 4$
involutory MDS matrices in these papers is limited to finite fields of small
order. This motivates us to explore the theoretical aspects of effectively
generating all $4\times 4$ involutory MDS matrices over $\mathbb{F}_{2^{m}}$
and reducing the search space to $(2^{m}-1)^{5}$. We present a hybrid method
to construct all $4\times 4$ involutory MDS matrices over $\mathbb{F}_{2^{m}}$
and demonstrate how our method significantly reduces the search space.
Initially, we partition the set of all $4\times 4$ involutory MDS matrices
into equivalence classes and identify representative matrices for each class.
Utilizing these representative matrices, we efficiently generate all
involutory MDS matrices within the same class through direct computation. We
illustrate that these involutory MDS class representative matrices of order
$4$ can be found by examining only a small subset of all $4\times 4$
involutory matrices over $\mathbb{F}_{2^{m}}$.

Additionally, we provide some characterizations of involutory MDS matrices of
even order. Furthermore, we propose a new matrix form for obtaining all
involutory MDS matrices of even order and compare it with other matrix forms
available in the literature.

The paper is structured as follows. Section 2 briefly discusses the
mathematical background and notations used in the article. In Section 3, we
provide some characterizations of involutory MDS matrices of even order.
Section 4 presents the theory underlying the generation of all $4\times 4$
involutory MDS matrices over $\mathbb{F}_{2^{m}}$. The algorithm to generate
all $4\times 4$ involutory MDS matrices and counting all $4\times 4$
involutory MDS matrices over some finite fields is detailed in Section 5.
Finally, we conclude in Section 6.

##  2 Mathematical Preliminaries

In this section, we discuss some definitions and mathematical preliminaries
that are important in our context. Let $\mathbb{F}_{2^{m}}$ be a finite field
of order $2^{m}$, where $m$ is a positive integer. Let
$\mathbb{F}_{2^{m}}^{*}$ denote the multiplicative group of the finite field
$\mathbb{F}_{2^{m}}$.

Let $M$ be any $n\times n$ matrix over $\mathbb{F}_{2^{m}}$ and let $|M|$
denote its determinant. A matrix $D$ of order $n$ is diagonal if $(D){i,j}=0$
for $i\neq j$. Using $d_{i}=(D){i,i}$, we represent the diagonal matrix $D$ as
$\operatorname{Diag}(d_{1},d_{2},\ldots,d_{n})$. The determinant of $D$ is
$|D|=\prod_{i=1}^{n}{d_{i}}$. Hence, $D$ is non-singular over
$\mathbb{F}_{2^{m}}$ if and only if $d_{i}\neq 0$ for $1\leq i\leq n$.

In this paper, we denote the identity matrix of order $n$ as $I_{n}$. Let
$\mathcal{L}_{2n}$ be the set of $2n\times 2n$ matrices over
$\mathbb{F}_{2^{m}}$ of type $\begin{pmatrix}A&B\\\ C&D\end{pmatrix}$, where
$A,B,C,D$ are $n\times n$ non-singular matrices.

An MDS matrix is utilized in various cryptographic primitives as a diffusion
layer. The notion of an MDS matrix originates from coding theory, particularly
from the domain of maximum distance separable (MDS) codes. An $[n,k,d]$ code
is considered MDS if it satisfies the singleton bound $d=n-k+1$.

######  Theorem 1.

[12, page 321] An $[n,k,d]$ code $C$ with generator matrix $G=[I~{}|~{}M]$,
where $M$ is a $k\times(n-k)$ matrix, is MDS if and only if every square sub-
matrix (formed from any $i$ rows and any $i$ columns, for any
$i=1,2,\ldots,min\\{k,n-k\\}$) of $M$ is non-singular.

######  Definition 1.

A matrix $M$ of order $n$ is said to be an MDS matrix if $[I~{}|~{}M]$ is a
generator matrix of a $[2n,n]$ MDS code.

Another way to define an MDS matrix is as follows:

######  Fact 1.

A square matrix $M$ is an MDS matrix if and only if every square sub-matrices
of $M$ are non-singular.

Implementing involutory diffusion matrices is more advantageous since it
enables using a single module for both the encryption and decryption
processes.

######  Definition 2.

An involutory matrix is a square matrix $M$ of order $n$ that is self-
invertible, or equivalently, satisfies $M^{2}=I_{n}$.

One of the fundamental row operations on matrices involves multiplying a row
or column of a matrix by a non-zero scalar. The MDS property remains unchanged
under such operations. Therefore, we have the following result concerning MDS
matrices.

######  Lemma 1.

Let $M$ be an involutory MDS matrix, then for any non-singular diagonal matrix
$D$, $D^{-1}MD$ will also be an involutory MDS matrix.

##  3 Construction of Involutory MDS Matrices

In this section, we explore some characterizations of involutory MDS matrices
for even orders. There exists a direct correspondence between the set of
$n\times n$ nilpotent matrices of index $2$ over $\mathbb{F}_{2^{m}}$ and the
set of $n\times n$ involutory matrices. It is straightforward to demonstrate
that each mapping from $M\mapsto M+I_{n}$ establishes a bijection between
these sets. This correspondence yields valuable insights into identifying
involutory matrices.

Rather than focusing on involutory matrices directly, we can examine nilpotent
matrices of index $2$. The relationship between involutory matrices and
nilpotent matrices of index $2$ is detailed in the following proposition.

######  Proposition 1.

Let $M$ be an $n\times n$ involutory matrix of even order over
$\mathbb{F}_{2^{m}}$, then $(M+I_{n})$ is a nilpotent matrix over
$\mathbb{F}_{2^{m}}$ with index $2$, and its rank can be at most
$\frac{n}{2}$.

###### Proof.

We have

$(M+I_{n})^{2}=(M^{2}+I_{n})=0.$

Thus, the index of $(M+I_{n})$ is 2. Let’s say $N=M+I_{n}$.

Next, we demonstrate that the rank of N cannot be more than $\frac{n}{2}$.
Let’s define the null space and range space of $M$ as $N(N)$ and $R(N)$,
respectively. Consequently, $rank(N)=\dim(R(N))$. Let $y\in R(N)$ be
arbitrary. Then a vector $x$ exists such that $y=Nx$. Now,
$Ny=N(Nx)=N^{2}x=0$. It means $y\in N(N)$. We now have $R(N)\subset N(N)$. It
implies

$\dim(R(N))\leq\dim(N(N)).$ (1)

From rank-nullity theorem

$\dim(R(N))+\dim(N(N))=n.$ (2)

(1) and (2) together imply that $rank(N)=\dim(R(N))\leq\frac{n}{2}$. ∎

If $M$ is an involutory MDS matrix in $\mathcal{L}_{2n}$ 555$\mathcal{L}_{2n}$
denotes the set of $2n\times 2n$ matrices over $\mathbb{F}_{2^{m}}$ of type
$\begin{pmatrix}A&B\\\ C&D\end{pmatrix}$, where $A,B,C,D$ are $n\times n$ non-
singular matrices., then, the next theorem demonstrates that the rank of
$M+I_{2n}$ is exactly half of the order of $M$. Here, $I_{2n}$ is a $2n\times
2n$ identity matrix.

######  Theorem 2.

Let $M$ be an involutory MDS matrix in $\mathcal{L}_{2n}$. Then,
rank$(M+I_{2n})=n$.

###### Proof.

Let $M=(m_{ij})$ and $M+I_{2n}=N=(n_{ij})$. Then we have.

$m_{ij}=\begin{cases}n_{ij}\text{ if }i\neq j\\\ n_{ij}+1\text{
otherwise}.\end{cases}$ (3)

According to proposition 1, $N$ is a nilpotent matrix with index $2$ and
rank$(N)\leq\frac{2n}{2}=n$. There are two possible outcomes: either
rank$(N)=n$ or rank$(N)<n$. If possible, assume that rank$(N))\neq n$. Then,
we must have rank$(N)<n$. In this scenario, all $n\times n$ minors of $N$ are
zero.

Take the $n\times n$ sub-matrix $W=(w_{ij})$ obtained from $N$ by the last $n$
rows and the first $n$ columns from it. The determinant of $W$ must be zero.
It is evident from (3) that $W$ is also a sub-matrix of $M$ whose determinant
is zero. This is a contradiction because $M$ is an MDS matrix. As a result,
rank$(N)=n$. ∎

The following theorem discusses a specific form of an involutory matrix $M$ in
$\mathcal{L}_{2n}$.

######  Theorem 3.

$M$ be an involutory matrix in $\mathcal{L}_{2n}$ if and only if $M$ can be
written as follows:

$\displaystyle M$ $\displaystyle=\begin{pmatrix}PC&PCP\\\
C&CP\end{pmatrix}+I_{2n},$ (4)

where $P$ and $C$ are $n\times n$ non-singular matrices over
$\mathbb{F}_{2^{m}}$.

###### Proof.

Let $M=N+I_{2n}$. So $N$ is nilpotent of index $2$. Since
$M\in\mathcal{L}_{2n}$, it has the form

$M=\begin{pmatrix}A^{\prime}&B^{\prime}\\\
C^{\prime}&D^{\prime}\end{pmatrix},$

where $A,B,C,D$ are $n\times n$ non-singular matrices. Then,

$N=\begin{pmatrix}A^{\prime}+I_{n}&B^{\prime}\\\
C^{\prime}&D^{\prime}+I_{n}\end{pmatrix}=\begin{pmatrix}A&B\\\
C&D\end{pmatrix},$

where $A=A^{\prime}+I_{n},B=B^{\prime},C=C^{\prime}$ and $D=D^{\prime}+I_{n}$.

We have $rank(N)\geq n$ since $C=C^{\prime}$ is non-singular. Furthermore,
since $N$ is nilpotent, its rank cannot exceed $n$. As a result, $rank(N)=n$
is required. Because of this, rows of $A$ can be expressed as a linear
combination of rows of $C$. In this instance, a $n\times n$ non-singular
matrix $P$ exists such that $A=PC$. Similarly, a linear combination of columns
of $C$ can be used to write columns of $D$. It implies that a non-singular
matrix $Q$ exists such that $D=CQ$. Then, there is

$N=\begin{pmatrix}PC&B\\\ C&CQ\end{pmatrix}.$ (5)

Let $O$ be a zero matrix of order $n$. Pre-multiply $N$ by
$\begin{pmatrix}P^{-1}&I_{n}\\\ O&I_{n}\end{pmatrix}$ and post-multiply by
$\begin{pmatrix}I_{n}&I_{n}\\\ O&Q^{-1}\end{pmatrix}$ get the following:

$\displaystyle\begin{pmatrix}P^{-1}&I_{n}\\\
O&I_{n}\end{pmatrix}N\begin{pmatrix}I_{n}&I_{n}\\\ O&Q^{-1}\end{pmatrix}$
$\displaystyle=$ $\displaystyle\begin{pmatrix}P^{-1}&I_{n}\\\
O&I_{n}\end{pmatrix}\begin{pmatrix}PC&B\\\
C&CQ\end{pmatrix}\begin{pmatrix}I_{n}&I_{n}\\\ O&Q^{-1}\end{pmatrix}$
$\displaystyle=$ $\displaystyle\begin{pmatrix}O&P^{-1}B+CQ\\\
C&CQ\end{pmatrix}\begin{pmatrix}I_{n}&I_{n}\\\ O&Q^{-1}\end{pmatrix}$
$\displaystyle=$ $\displaystyle\begin{pmatrix}O&P^{-1}BQ^{-1}+C\\\
C&O\end{pmatrix}.$

Given that $N$ has the rank $n$, we must have $P^{-1}BQ^{-1}+C=O$ either or
$B=PCQ$. By substituting the value of $B$ in (5) we get

$N=\begin{pmatrix}PC&PCQ\\\ C&CQ\end{pmatrix}.$

In addition, nilpotency of $N$ provides us $N^{2}=0$. To put it another way,
i.e.,

$\begin{pmatrix}PCPC+PCQC&PCPCQ+PCQCQ\\\ CPC+CQC&CPCQ+CQCQ\end{pmatrix}=O.$

Comparing both sides, we obtain the following:

$PCPC+PCQC=O\implies PC(P+Q)C=0\implies P=Q.$

When $Q$ is substituted in (5), we have $M$ in desired form as follows:

$M=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{2n}.$

Conversely, if $M$ has the following form:

$M=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{2n}.$

We have,

$\displaystyle M^{2}=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}^{2}+I_{2n}$
$\displaystyle\implies$ $\displaystyle
M^{2}=\begin{pmatrix}PCPC+PCPC&PCPPCP+PCPCP\\\
CPC+CPC&CPCP+CPCP\end{pmatrix}+I_{2n}$ $\displaystyle\implies$ $\displaystyle
M^{2}=I_{2n}.$

Therefore, $M$ is involutory. Thus, the proof is complete. ∎

Since each sub-matrix of an MDS matrix is non-singular, according to Theorem
3, we can find all involutory MDS matrices of even order $2n$ by iterating
over the entries of the matrices $P$ and $C$ of order $n$. Consequently, the
search space for finding all involutory MDS matrices of even order $2n$ is
reduced to $2^{2n^{2}m}$ over the field $\mathbb{F}_{2^{m}}$. Next, we will
explore the relationship between our proposed matrix form and other matrix
forms available in the literature for finding all involutory MDS matrices of
even order.

###  3.1 Relation to the Matrix Forms available in the literature

In [23], the authors propose a matrix form to find all involutory MDS matrices
of even order. The proposed matrix form is given by:

$\displaystyle M$ $\displaystyle=\begin{bmatrix}A_{1}&A_{2}\\\
A_{2}^{-1}(I_{n}+A_{1}^{2})&A_{2}^{-1}A_{1}A_{2}\end{bmatrix},$

where $A_{1}$ and $A_{2}$ are non-singular matrices over $\mathbb{F}_{2^{m}}$.
Comparing our proposed matrix form given in (4), we have:

$\displaystyle PC+I_{n}=A_{1},$ $\displaystyle PCP=A_{2}$ $\displaystyle
C=A_{2}^{-1}(I_{n}+A_{1}^{2}),$ $\displaystyle CP+I_{n}=A_{2}^{-1}A_{1}A_{2}.$

Note that since $M$ is MDS, we have $A_{2}^{-1}(I_{n}+A_{1}^{2})$ is non-
singular. Therefore, $(I_{n}+A_{1}^{2})$ must also be non-singular. Now,
$(I_{n}+A_{1}^{2})=(I_{n}+A_{1})^{2}$, implying that $(I_{n}+A_{1})$ is also
non-singular. Therefore, from the above equations, we can derive that:

$\displaystyle
P=(I_{n}+A_{1})^{-1}A_{2}~{}\text{and}~{}C=A_{2}^{-1}(I_{n}+A_{1}^{2}).$

The proposed matrix form in [17] for generating all involutory MDS matrices of
even order is given by:

$\displaystyle M$
$\displaystyle=\begin{bmatrix}A_{1}&(I_{n}+A_{1}^{2})A_{3}^{-1}\\\
A_{3}&A_{3}A_{1}A_{3}^{-1}\end{bmatrix},$

where $A_{1}$ and $A_{3}$ are non-singular matrices over $\mathbb{F}_{2^{m}}$.
Comparing our proposed matrix form given in (4), we have:

$\displaystyle PC+I_{n}=A_{1},$ $\displaystyle
PCP=(I_{n}+A_{1}^{2})A_{3}^{-1}$ $\displaystyle C=A_{3},$ $\displaystyle
CP+I_{n}=A_{3}A_{1}A_{3}^{-1}.$

Note that since $M$ is MDS, $A_{3}$ is non-singular. Additionally, since
$(I_{n}+A_{1}^{2})A_{3}^{-1}$ is non-singular and
$(I_{n}+A_{1}^{2})=(I_{n}+A_{1})^{2}$, implying that $(I_{n}+A_{1})$ is also
non-singular. Therefore, from the above equations, we can derive that:

$\displaystyle P=(I_{n}+A_{1})A_{3}^{-1}~{}\text{and}~{}C=A_{3}.$

##  4 Construction of $4\times 4$ Involutory MDS Matrices

This section introduces a hybrid technique for generating all $4\times 4$
involutory MDS matrices over $\mathbb{F}_{2^{m}}$. There are $(2^{m})^{16}$
matrices of order $4\times 4$ over $\mathbb{F}_{2^{m}}$ in existence, of which
$((2^{m})^{4}-1)((2^{m})^{4}-2^{m})((2^{m})^{4}-(2^{m})^{2})((2^{m})^{4}-(2^{m})^{3})$
are non-singular. It is not easy to filter out all $4\times 4$ involutory MDS
matrices from the collection of all $4\times 4$ non-singular matrices. By
resolving a small set of non-linear equations over $\mathbb{F}_{2^{m}}$, the
authors in [8] propose a matrix form to generate all $3\times 3$ involutory
MDS matrices over $\mathbb{F}_{2^{m}}$. However, it is not simple to
generalize the concept of [8] for orders of $4\times 4$ or greater.

For a brute force search for involutory MDS matrices over
$\mathbb{F}_{2^{m}}$, the search space is the entire set of non-singular
matrices over $\mathbb{F}_{2^{m}}$. Some obvious checks may be imposed on the
structure of the matrix for MDS search, such as choosing the entries of the
matrices to be non-zero. However, this would not result in a significant
reduction in the space and search time. Another idea for refining the search
is to search over involutory matrices only. This would provide a good amount
of reduction in search space but still offers a fairly large search space even
for smaller matrix dimensions. For this purpose, we utilize a systematic
construction that yields a search space much smaller than the set of all
involutory matrices over $\mathbb{F}_{2^{m}}$.

The main goal of our suggested approach is to first locate all involutory MDS
class representative matrices of order $4$ through search and then to obtain
all $4\times 4$ involutory MDS matrices by multiplying a diagonal matrix and
its inverse (pre and post) with these involutory MDS class representative
matrices.

Let $\mathcal{M}_{4}$ denote the matrix space containing all $4\times 4$
involutory MDS matrices. According to Lemma 1, if $A\in\mathcal{M}_{4}$, then
for any non-singular diagonal matrix $D$, $D^{-1}MD\in\mathcal{M}_{4}$. Now,
let’s define an equivalence relation on $\mathcal{M}_{4}$, such that matrices
$A$ and $B$ are related if there exists an invertible diagonal matrix
$D=\operatorname{Diag}(1,b_{1},b_{2},b_{3})$ over $\mathbb{F}_{2^{m}}^{*}$
such that $A=D^{-1}BD$. Thus, for $A\in\mathcal{M}_{4}$, if $\bar{A}$ denotes
the equivalence class of $A$, all the matrices in $\bar{A}$ can be obtained by
$D^{-1}AD$, varying the diagonal matrices $D=\text{Diag}(1,b_{1},b_{2},b_{3})$
over $\mathbb{F}_{2^{m}}^{*}$.

To find the class representative matrices of all $4\times 4$ involutory MDS
matrices over $\mathbb{F}_{2^{m}}^{*}$, we define the matrix $R$ as follows:

$R=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{4},$ (6)

where $I_{4}$ is $4\times 4$ identity matrix, and $P$ and $C$ are the $2\times
2$ matrices given by

$\displaystyle C$ $\displaystyle=c\begin{pmatrix}pq+r&p\\\ q&1\end{pmatrix}$
$\displaystyle P=\begin{pmatrix}d+1&d\\\ d&d+1\end{pmatrix},$

where $p,q,r,c,d\in\mathbb{F}_{2^{m}}^{*}$.

It can be observed that the sum of any rows or columns of $P$ is $1$.
Additionally, the $2\times 2$ non-singular matrices $C$ and $P$ over
$\mathbb{F}_{2^{m}}^{*}$ define the matrix form $R$ mentioned above for
determining all involutory MDS class representative matrices of order $4$.
Consequently, by searching within this reduced subset of all $4\times 4$
involutory matrices, the search space is now constrained to $(2^{m}-1)^{5}$.

Let $\mathcal{R}_{4}$ be the matrix space containing all representative
involutory matrices of the equivalence classes. Thus, we can write the set
$\mathcal{M}_{4}$ of all $4\times 4$ involutory MDS matrices as follows:

$\mathcal{M}_{4}=\underset{R\text{ is
MDS}}{\underset{R\in\mathcal{R}_{4}}{\bigcup}}\\{D^{-1}RD:~{}D=\operatorname{Diag}(1,b_{1},b_{2},b_{3}),~{}\text{where}~{}b_{1},b_{2},b_{3}\in\mathbb{F}_{2^{m}}^{*}\\}.$

We can associate elements of $\mathcal{R}_{4}$ with the tuple $(p,q,r,c,d)$
because $\mathcal{R}_{4}$ consists of five variables: $p,q,r,c$, and $d$ in
$\mathbb{F}_{2^{m}}^{*}$. Now, let’s demonstrate that two different tuples
correspond to two different elements of $\mathcal{R}_{4}$.

######  Proposition 2.

The cardinality of the set $\mathcal{R}_{4}$ is $(2^{m}-1)^{5}$.

###### Proof.

Consideration the two tuples, $(p_{1},q_{1},r_{1},c_{1},d_{1})$ and
$(p_{2},q_{2},r_{2},c_{2},d_{2})$, which are linked to the same elements
$R\in\mathcal{R}_{4}$. Let $P_{1},P_{2},C_{1},C_{2}$ as

$\displaystyle C_{1}=c_{1}\begin{pmatrix}p_{1}q_{1}+r_{1}&p_{1}\\\
q_{1}&1\end{pmatrix},C_{2}=c_{2}\begin{pmatrix}p_{2}q_{2}+r_{2}&p_{2}\\\
q_{2}&1\end{pmatrix}$ $\displaystyle P_{1}=\begin{pmatrix}d_{1}+1&d_{1}\\\
d_{1}&d_{1}+1\end{pmatrix}\text{ and }P_{2}=\begin{pmatrix}d_{2}+1&d_{2}\\\
d_{2}&d_{2}+1\end{pmatrix}.$

Now, based on the assumption, we have

$\displaystyle\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}+I_{4}=\begin{pmatrix}P_{2}C_{2}&P_{2}C_{2}P_{2}\\\
C_{2}&C_{2}P_{2}\end{pmatrix}+I_{4}$ $\displaystyle\implies$
$\displaystyle\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}=\begin{pmatrix}P_{2}C_{2}&P_{2}C_{2}P_{2}\\\
C_{2}&C_{2}P_{2}\end{pmatrix}$ $\displaystyle\implies$ $\displaystyle
C_{1}=C_{2}\text{ and }C_{1}P_{1}=C_{2}P_{2}$ $\displaystyle\implies$
$\displaystyle
c_{1}=c_{2},c_{1}p_{1}=c_{2}p_{2},c_{1}q_{1}=c_{2}q_{2}~{}\text{and}~{}d_{1}=d_{2}$
$\displaystyle\implies$
$\displaystyle(p_{1},q_{1},r_{1},c_{1},d_{1})=(p_{2},q_{2},r_{2},c_{2},d_{2}).$

Hence the proposition. ∎

######  Remark 1.

It is noteworthy that the authors of [23] demonstrate that the search space
for finding involutory MDS matrices of even order over $\mathbb{F}_{2^{m}}$
may be reduced to $2^{(mn^{2})/2}$, which is $2^{8m}$ when $n=4$. In contrast,
the search space in our method is $(2^{m}-1)^{5}$, significantly smaller than
that in [23].

Now, we will delve into the theoretical aspects behind the proposed matrix
form given in (6). For this discussion, we require the following significant
characteristics for involutory MDS matrices.

######  Theorem 4.

For any involutory MDS matrix $M\in\mathcal{L}_{4}$ given by

$\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{4},$

the matrices $C$ and $P$ have all non-zero entries.

###### Proof.

Let $M=(m_{ij}),C=\begin{pmatrix}c_{11}&c_{12}\\\ c_{21}&c_{22}\end{pmatrix}$
and $P=\begin{pmatrix}p_{11}&p_{12}\\\ p_{21}&p_{22}\end{pmatrix}$.

Since $M$ is an MDS matrix, we can say that all the entries in $C$ are non-
zero. Now, consider the following sub-matrices of $M$:

$M_{1}=\begin{pmatrix}m_{12}&m_{14}\\\
m_{32}&m_{34}\end{pmatrix},M_{2}=\begin{pmatrix}m_{12}&m_{13}\\\
m_{42}&m_{43}\end{pmatrix},M_{3}=\begin{pmatrix}m_{21}&m_{24}\\\
m_{31}&m_{34}\end{pmatrix},M_{4}=\begin{pmatrix}m_{21}&m_{23}\\\
m_{41}&m_{43}\end{pmatrix}.$

These sub-matrices can be represented in terms of the elements of $C$ and $P$
as follows:

$\displaystyle M_{1}$
$\displaystyle=\begin{pmatrix}c_{12}p_{11}+c_{22}p_{12}&c_{11}p_{11}p_{12}+c_{21}p_{12}^{2}+c_{12}p_{11}p_{22}+c_{22}p_{12}p_{22}\\\
c_{12}&c_{11}p_{12}+c_{12}p_{22}\end{pmatrix}$ $\displaystyle M_{2}$
$\displaystyle=\begin{pmatrix}c_{12}p_{11}+c_{22}p_{12}&c_{11}p_{11}^{2}+c_{21}p_{11}p_{12}+c_{12}p_{11}p_{21}+c_{22}p_{12}p_{21}\\\
c_{22}&c_{21}p_{11}+c_{22}p_{21}\end{pmatrix}$ $\displaystyle M_{3}$
$\displaystyle=\begin{pmatrix}c_{11}p_{21}+c_{21}p_{22}&c_{11}p_{12}p_{21}+c_{21}p_{12}p_{22}+c_{12}p_{21}p_{22}+c_{22}p_{22}^{2}\\\
c_{11}&c_{11}p_{12}+c_{12}p_{22}\end{pmatrix}$ $\displaystyle M_{4}$
$\displaystyle=\begin{pmatrix}c_{11}p_{21}+c_{21}p_{22}&c_{11}p_{11}p_{21}+c_{12}p_{21}^{2}+c_{21}p_{11}p_{22}+c_{22}p_{21}p_{22}\\\
c_{21}&c_{21}p_{11}+c_{22}p_{21}\end{pmatrix}.$

The following are the determinants of these sub-matrices.

$\displaystyle|M_{1}|=c_{12}c_{21}p_{12}^{2}+c_{11}c_{22}p_{12}^{2}=p_{12}^{2}\det{(C)}$
$\displaystyle|M_{2}|=c_{12}c_{21}p_{11}^{2}+c_{11}c_{22}p_{11}^{2}=p_{11}^{2}\det{(C)}$
$\displaystyle|M_{3}|=c_{12}c_{21}p_{22}^{2}+c_{11}c_{22}p_{22}^{2}=p_{22}^{2}\det{(C)}$
$\displaystyle|M_{4}|=c_{12}c_{21}p_{21}^{2}+c_{11}c_{22}p_{21}^{2}=p_{22}^{2}\det{(C)}.$

Given that $M$ is an MDS matrix, all the sub-matrices $M_{1},M_{2},M_{3}$ and
$M_{4}$ are non-singular. Therefore, $p_{11}$, $p_{12}$, $p_{21}$, and
$p_{22}$ must be non-zero. ∎

By Theorem 4, we know that neither the matrices $P$ nor $C$ can have zero
entries, leading to the following result.

######  Corollary 1.

Let $R=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{4}\in\mathcal{R}_{4}$ be
an MDS matrix, where $C=c\begin{pmatrix}pq+r&p\\\ q&1\end{pmatrix}$ and
$P=\begin{pmatrix}d+1&d\\\ d&d+1\end{pmatrix}$. Then we have

  1. 1.

$r\neq pq$,

  2. 2.

$d\neq 1$.

The following theorem states that we can generate only one $4\times 4$
involutory MDS matrix from a representative involutory MDS matrix of order $4$
by using the non-singular diagonal matrix
$D=\operatorname{Diag}(1,b_{1},b_{2},b_{3})$ over $\mathbb{F}_{2^{m}}^{*}$.

######  Theorem 5.

If a $4\times 4$ involutory MDS matrix $M\in[R]$, where $[R]$ denotes the
equivalence class of the representative matrix $R\in\mathcal{R}_{4}$, then
there exists a unique diagonal matrix
$D=\operatorname{Diag}(1,b_{1},b_{2},b_{3})$ over $\mathbb{F}_{2^{m}}^{*}$
such that $D^{-1}RD=M$.

###### Proof.

Let $M$ be any involutory MDS matrix of order $4$ over $\mathbb{F}_{2^{m}}$
and suppose that there exist a non-singular diagonal matrix
$D=\operatorname{Diag}(1,b_{1},b_{2},b_{3})$ and an involutory MDS class
representative matrix $R\in\mathcal{R}_{4}$ such that

$D^{-1}RD=M.$ (7)

Since $M$ is an involutory MDS matrix, by using Theorem 3, it can be expressed
as follows:

$M=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{4},$

where $C,P$ are $2\times 2$ non-singular matrices. Let
$C=\begin{pmatrix}c_{11}&c_{12}\\\ c_{21}&c_{22}\\\ \end{pmatrix}$ and
$P=\begin{pmatrix}p_{11}&p_{12}\\\ p_{21}&p_{22}\\\ \end{pmatrix}$.

We can write diagonal matrix $D$ as $\begin{pmatrix}B_{1}&O\\\
O&B_{2}\end{pmatrix}$, where $B_{1}=\begin{pmatrix}1&0\\\
0&b_{1}\end{pmatrix}$ and $B_{2}=\begin{pmatrix}b_{2}&0\\\
0&b_{3}\end{pmatrix}$. Given that $R\in\mathcal{R}_{4}$, so we can also write
$R$ as follows:

$R=\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}+I_{4},$

where $C_{1}$ and $P_{1}$ are still undetermined.

Now since $DMD^{-1}=R$, we have

$\begin{pmatrix}B_{1}&O\\\ O&B_{2}\end{pmatrix}\left[\begin{pmatrix}PC&PCP\\\
C&CP\end{pmatrix}+I_{4}\right]\begin{pmatrix}B_{1}^{-1}&O\\\
O&B_{2}^{-1}\end{pmatrix}=\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}+I_{4}$

or,

$\begin{pmatrix}B_{1}&O\\\ O&B_{2}\end{pmatrix}\begin{pmatrix}PC&PCP\\\
C&CP\end{pmatrix}\begin{pmatrix}B_{1}^{-1}&O\\\
O&B_{2}^{-1}\end{pmatrix}=\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}\\\ $

or,

$\begin{pmatrix}B_{1}PCB_{1}^{-1}&B_{1}PCPB_{2}^{-1}\\\
B_{2}CB_{1}^{-1}&B_{2}CPB_{2}^{-1}\end{pmatrix}=\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}.$ (8)

We have from (8)

$\displaystyle B_{1}PCB_{1}^{-1}=P_{1}C_{1}$ $\displaystyle
B_{1}PCPB_{2}^{-1}=P_{1}C_{1}P_{1}$ $\displaystyle
B_{2}CPB_{2}^{-1}=C_{1}P_{1}$ (9) $\displaystyle B_{2}CB_{1}^{-1}=C_{1}.$ (10)

Let $C_{1}=c\begin{pmatrix}s&p\\\ q&1\end{pmatrix}$. Then by using (10), we
obtain

$\displaystyle s$ $\displaystyle=$ $\displaystyle b_{2}c_{11}c^{-1}$
$\displaystyle p$ $\displaystyle=$ $\displaystyle b_{2}b_{1}^{-1}c_{12}c^{-1}$
$\displaystyle q$ $\displaystyle=$ $\displaystyle b_{3}c_{21}c^{-1}$
$\displaystyle c$ $\displaystyle=$ $\displaystyle b_{3}b_{1}^{-1}c_{22}.$

Since $c_{11},c_{12},c_{21}$, and $c_{22}$ are elements of the sub-matrix $C$
of the MDS matrix $M$, we can conclude that all of $c_{11},c_{12},c_{21}$, and
$c_{22}$ are non-zero. Additionally, since $b_{1},b_{2}$, and $b_{3}$ are non-
zero, we can infer that $s,p,q$, and $c$ are also non-zero. Thus, we have

$\displaystyle C_{1}$
$\displaystyle=\begin{pmatrix}b_{2}c_{11}&b_{2}b_{1}^{-1}c_{12}\\\
b_{3}c_{21}&b_{3}b_{1}^{-1}c_{22}\end{pmatrix}.$ (11)

Now,
$|C_{1}|=b_{2}b_{3}b_{1}^{-1}c_{11}c_{22}+b_{2}b_{3}b_{1}^{-1}c_{21}c_{12}=b_{2}b_{3}b_{1}^{-1}|C|$,
which indicates that $|C_{1}|\neq 0$. Since $C_{1}$ is non-singular, let
$\begin{vmatrix}s&p\\\ q&1\end{vmatrix}=r$, where $r\neq 0$. So
$r=pq+s\implies s=pq+r$. Thus, we have $C_{1}=c\begin{pmatrix}pq+r&p\\\
q&1\end{pmatrix}$. Also, according to our proposed representative matrix form,
we can write the matrix $P_{1}$ as $P_{1}=\begin{pmatrix}d+1&d\\\
d&d+1\end{pmatrix}$, where $d\in\mathbb{F}_{2^{m}}^{*}\setminus\set{1}$.

Now, from (9) and (10), we have

$\displaystyle(B_{2}CB_{1}^{-1})^{-1}B_{2}CPB_{2}^{-1}=C_{1}^{-1}C_{1}P_{1}$
(12) $\displaystyle\implies$ $\displaystyle
B_{1}C^{-1}B_{2}^{-1}B_{2}CPB_{2}^{-1}=C_{1}^{-1}C_{1}P_{1}$
$\displaystyle\implies$ $\displaystyle B_{1}PB_{2}^{-1}=P_{1}$
$\displaystyle\implies$ $\displaystyle\begin{pmatrix}1&0\\\
0&b_{1}\end{pmatrix}\begin{pmatrix}p_{11}&p_{12}\\\
p_{21}&p_{22}\end{pmatrix}\begin{pmatrix}b_{2}^{-1}&0\\\
0&b_{3}^{-1}\end{pmatrix}=P_{1}$ $\displaystyle\implies$
$\displaystyle\begin{pmatrix}p_{11}&p_{12}\\\
b_{1}p_{21}&b_{1}p_{22}\end{pmatrix}\begin{pmatrix}b_{2}^{-1}&0\\\
0&b_{3}^{-1}\end{pmatrix}=P_{1}$ $\displaystyle\implies$ $\displaystyle
P_{1}=\begin{pmatrix}p_{11}b_{2}^{-1}&p_{12}b_{3}^{-1}\\\
b_{1}p_{21}b_{2}^{-1}&b_{1}p_{22}b_{3}^{-1}\end{pmatrix}.$

Thus, from Equations (11) and (12), we observe that the matrices $P_{1}$ and
$C_{1}$, i.e., the representative matrix $R$, can be expressed in terms of the
entries of the sub-matrices $C$ and $P$ of the MDS matrix $M$.

Now, we will determine the values $b_{1},b_{2}$ and $b_{3}$ of the diagonal
matrix $D$ in terms of the entries of the sub-matrices $C$ and $P$. Since the
sum of each row or column of $P_{1}$ is equal to 1, we have

$\displaystyle p_{11}b_{2}^{-1}+p_{12}b_{3}^{-1}=1$ (13) $\displaystyle
b_{1}p_{21}b_{2}^{-1}+b_{1}p_{22}b_{3}^{-1}=1$ (14) $\displaystyle
p_{11}b_{2}^{-1}+b_{1}p_{21}b_{2}^{-1}=1$ (15) $\displaystyle
p_{12}b_{3}^{-1}+b_{1}p_{22}b_{3}^{-1}=1.$ (16)

Now, from (15) and (16), we have

$\displaystyle b_{2}=p_{11}+b_{1}p_{21}$ (17) $\displaystyle
b_{3}=p_{12}+b_{1}p_{22}.$ (18)

By solving the above equations, we have
$b_{1}=\sqrt{\frac{p_{11}p_{12}}{p_{21}p_{22}}}$. Thus, $b_{1}$ has been
uniquely determined by the terms of the entries of the sub-matrix $P$. We can
also uniquely determine $b_{2}$ and $b_{3}$ by substituting the value of
$b_{1}$ in (17) and (18). Thus, we have

$\displaystyle b_{1}$
$\displaystyle=\sqrt{\frac{p_{11}p_{12}}{p_{21}p_{22}}},~{}b_{2}=p_{11}+p_{21}\sqrt{\frac{p_{11}p_{12}}{p_{21}p_{22}}},~{}b_{3}=p_{12}+p_{22}\sqrt{\frac{p_{11}p_{12}}{p_{21}p_{22}}}.$
(19)

Now, according to Theorem 4, we have
$p_{11},p_{12},p_{21},p_{22}\in\mathbb{F}_{2^{m}}^{*}$ which implies
$b_{1}\in\mathbb{F}_{2^{m}}^{*}$.

If possible let $b_{2}=0$. This implies that

$\displaystyle p_{11}+b_{1}p_{21}=0$ $\displaystyle\implies$ $\displaystyle
b_{1}=\frac{p_{11}}{p_{21}}$ $\displaystyle\implies$
$\displaystyle\sqrt{\frac{p_{11}p_{12}}{p_{21}p_{22}}}=\frac{p_{11}}{p_{21}}$
$\displaystyle\implies$ $\displaystyle p_{11}p_{22}+p_{12}p_{21}=0.$

This leads to a contradiction since $P$ is non-singular. Therefore,
$b_{2}\in\mathbb{F}_{2^{m}}^{*}$. Similarly, we can demonstrate that
$b_{3}\in\mathbb{F}_{2^{m}}^{*}$. Consequently, we can uniquely determine the
values $b_{1},b_{2}$, and $b_{3}$ of the diagonal matrix $D$. In other words,
we can conclude that there exists a unique diagonal matrix
$D=\text{Diag}(1,b_{1},b_{2},b_{3})$ over $\mathbb{F}_{2^{m}}^{*}$ such that
$D^{-1}RD=M$. ∎

In the following example of an involutory MDS matrix $M$, taken from the
literature [13, Example 4], we illustrate the MDS class representative matrix
$R$ and the non-singular diagonal matrix $D$, which are used to create the
involutory MDS matrix $M$.

######  Example 1.

Let $x^{4}+x+1$ be the primitive polynomial that generates
$\mathbb{F}_{2^{4}}$. Let $\alpha$ be a primitive element that is a root of
this polynomial. Now, consider the $4\times 4$ involutory MDS matrix $M$ over
$\mathbb{F}_{2^{4}}$.

$M=\begin{pmatrix}1&1&1&1\\\ 1&\alpha&\alpha^{2}&\alpha^{5}\\\
\alpha^{7}&\alpha^{10}&\alpha&\alpha^{5}\\\
\alpha^{9}&\alpha^{2}&\alpha^{10}&1\end{pmatrix}.$

Theorem 3 allows $M$ to be expressed as

$M=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{4},$

where the matrices $C,P$ are $2\times 2$ non-singular. Let
$C=\begin{pmatrix}c_{11}&c_{12}\\\ c_{21}&c_{22}\end{pmatrix}$ and
$P=\begin{pmatrix}p_{11}&p_{12}\\\ p_{21}&p_{22}\end{pmatrix}$.

Here $C=\begin{pmatrix}\alpha^{7}&\alpha^{10}\\\ \alpha^{9}&\alpha^{2}\\\
\end{pmatrix}$ and $CP=\begin{pmatrix}\alpha+1&\alpha^{5}\\\
\alpha^{10}&0\end{pmatrix}=\begin{pmatrix}\alpha^{4}&\alpha^{5}\\\
\alpha^{10}&0\end{pmatrix}.$ Then, $P=\begin{pmatrix}\alpha^{10}&\alpha^{8}\\\
1&1\end{pmatrix}$.

From (19), we can find the unique diagonal matrix
$D=\operatorname{Diag}(1,b_{1},b_{2},b_{3})$ as

$\displaystyle b_{1}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{p_{11}p_{12}}{p_{21}p_{22}}}=\sqrt{\alpha^{10}.\alpha^{8}}=\alpha^{9}$
$\displaystyle b_{2}$ $\displaystyle=$ $\displaystyle
p_{11}+b_{1}p_{21}=\alpha^{10}+\alpha^{9}=\alpha^{13}$ $\displaystyle b_{3}$
$\displaystyle=$ $\displaystyle
p_{12}+b_{1}p_{22}=\alpha^{8}+\alpha^{9}=\alpha^{12}.$

We will now locate an involutory MDS class representative matrix $R$, to which
the involutory MDS matrix $M$ belongs. Theorem 3 also allows $R$ to be
expressed as

$R=\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}+I.$

From (11) and (12), we can find $C_{1}$ and $P_{1}$ as follows:

$C_{1}=\begin{pmatrix}b_{2}c_{11}&b_{2}b_{1}^{-1}c_{12}\\\
b_{3}c_{21}&b_{3}b_{1}^{-1}c_{22}\end{pmatrix}=\begin{pmatrix}\alpha^{5}&\alpha^{14}\\\
\alpha^{6}&\alpha^{5}\end{pmatrix},\text{and
}P_{1}=\begin{pmatrix}p_{11}b_{2}^{-1}&p_{12}b_{3}^{-1}\\\
b_{1}p_{21}b_{2}^{-1}&b_{1}p_{22}b_{3}^{-1}\end{pmatrix}=\begin{pmatrix}\alpha^{12}&\alpha^{11}\\\
\alpha^{11}&\alpha^{12}\end{pmatrix}.$

Now, we get

$P_{1}C_{1}=\begin{pmatrix}0&\alpha^{6}\\\
\alpha^{9}&\alpha^{4}\end{pmatrix},C_{1}P_{1}=\begin{pmatrix}\alpha^{4}&\alpha^{6}\\\
\alpha^{9}&0\end{pmatrix},\text{and
}P_{1}C_{1}P_{1}=\begin{pmatrix}\alpha^{2}&\alpha^{3}\\\
\alpha^{13}&\alpha^{2}\end{pmatrix}.$

Hence, we have

$R=\begin{pmatrix}P_{1}C_{1}&P_{1}C_{1}P_{1}\\\
C_{1}&C_{1}P_{1}\end{pmatrix}+I_{4}=\begin{pmatrix}1&\alpha^{6}&\alpha^{2}&\alpha^{3}\\\
\alpha^{9}&\alpha&\alpha^{13}&\alpha^{2}\\\
\alpha^{5}&\alpha^{14}&\alpha&\alpha^{6}\\\
\alpha^{6}&\alpha^{5}&\alpha^{9}&1\end{pmatrix}.$

In [21], the authors propose a representative matrix form to generate all
involutory MDS matrices of order $4$. In their representative matrix, the sum
of any rows (and columns) equals $1$. In the following theorem, we prove that
these conditions are also guaranteed for our representative form.

######  Theorem 6.

The sum of any rows (and columns) in the involutory matrix
$M=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{2n}$ in $\mathcal{L}_{2n}$ is
equal to 1 if and only if the sum of any rows (and columns) of $P$ is equal to
1.

###### Proof.

Let $M=N+I_{2n}$, where $N=\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}$. Assume
that the sum of any rows (and columns) in $M$ equals 1. Let $J$ be a
$n$-dimensional vector with all entries 1. Then, we have

$\displaystyle M\begin{pmatrix}J\\\ J\\\ \end{pmatrix}=\begin{pmatrix}J\\\
J\\\ \end{pmatrix}$ $\displaystyle\implies$
$\displaystyle(N+I)\begin{pmatrix}J\\\ J\\\ \end{pmatrix}=\begin{pmatrix}J\\\
J\\\ \end{pmatrix}$ $\displaystyle\implies$ $\displaystyle
N\begin{pmatrix}J\\\ J\\\ \end{pmatrix}=\begin{pmatrix}O\\\ O\\\
\end{pmatrix}$ $\displaystyle\implies$ $\displaystyle\begin{pmatrix}PC&PCP\\\
C&CP\end{pmatrix}\begin{pmatrix}J\\\ J\\\ \end{pmatrix}=\begin{pmatrix}O\\\
O\end{pmatrix},$

where $O$ is the $n$-dimensional zero-column vector. Similar to this, we have

$\begin{pmatrix}J^{T}&J^{T}\\\ \end{pmatrix}\begin{pmatrix}PC&PCP\\\
C&CP\end{pmatrix}=\begin{pmatrix}O^{T}&O^{T}\\\ \end{pmatrix}.$

The matrix equations above lead us to the conclusion that

$\displaystyle(PC+PCP)J=O$ $\displaystyle\implies$
$\displaystyle(I+P)J=O~{}\hskip 20.00003pt[\text{since P and C are non-
singular matrices}]$ $\displaystyle\implies$ $\displaystyle PJ=J.$

Likewise, we can show that $J^{T}P=J^{T}$. This demonstrates that the sum of
any rows (and columns) in $P$ is equal to 1.

Conversely, suppose that the sum of any rows (and columns) in $P$ equals 1.
Then, we have

$\displaystyle M\begin{pmatrix}J\\\ J\end{pmatrix}$ $\displaystyle=$
$\displaystyle(N+I_{2n})\begin{pmatrix}J\\\ J\end{pmatrix}$ $\displaystyle=$
$\displaystyle N\begin{pmatrix}J\\\ J\end{pmatrix}+\begin{pmatrix}J\\\
J\end{pmatrix}$ $\displaystyle=$ $\displaystyle\begin{pmatrix}PC&PCP\\\
C&CP\end{pmatrix}\begin{pmatrix}J\\\ J\end{pmatrix}+\begin{pmatrix}J\\\
J\end{pmatrix}$ $\displaystyle=$ $\displaystyle\begin{pmatrix}(PC+PCP)J\\\
(C+CP)J\end{pmatrix}+\begin{pmatrix}J\\\ J\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}PC(I+P)J\\\
C(I+P)J\end{pmatrix}+\begin{pmatrix}J\\\ J\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}PC(J+PJ)\\\
C(J+PJ)\end{pmatrix}+\begin{pmatrix}J\\\ J\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}O\\\ O\end{pmatrix}+\begin{pmatrix}J\\\
J\end{pmatrix}~{}\hskip 20.00003pt[\text{Since }PJ=J]$ $\displaystyle=$
$\displaystyle\begin{pmatrix}J\\\ J\end{pmatrix}.$

Similarly, we can demonstrate that $J^{T}M=J^{T}$. Hence, the sum of any rows
(and columns) in $M$ equals 1. ∎

##  5 Experimental Results

This section presents our algorithm for generating all $4\times 4$ involutory
MDS matrices over $\mathbb{F}_{2^{m}}$. Now, let us consider the collection of
diagonal matrices of order $4$ as follows:

$\mathcal{D}_{4}=\\{\operatorname{Diag}(1,b_{1},b_{2},b_{3})\mid
b_{1},b_{2},b_{3}\in\mathbb{F}_{2^{m}}^{*}\\}.$

We follow the steps below:

  1. 1.

First, we generate all MDS class representative matrices of order $4$.

  2. 2.

From these MDS class representative matrices, we generate all $4\times 4$
involutory MDS matrices.

######  Remark 2.

It should be noted that while generating an involutory MDS matrix, using the
involutory matrix form $R$ offers an additional benefit. Since $R$ is always
involutory, we have $|R|=1$. Also, we know that if entries of the inverse
(here $R^{-1}=R$) of a $4\times 4$ non-singular matrix are non-zero, then all
its $3\times 3$ sub-matrices are non-singular. Thus, we need to check only
sub-matrices of $R$ of order $1$ and $2$ to verify the MDS property of $R$.

Algorithm 1 Printing all $4\times 4$ involutory MDS matrices over
$\mathbb{F}_{2^{m}}$

Input: Finite Field $\mathbb{F}_{2^{m}}$.
Output: All Involutory MDS matrices of order $4\times 4$.

1:for $p,q,r,c\in\mathbb{F}_{2^{m}}^{*}$ and
$d\in\mathbb{F}_{2^{m}}^{*}\setminus\set{1}$ do

2: $C\leftarrow c\begin{pmatrix}pq+r&p\\\ q&1\end{pmatrix}$ and
$P\leftarrow\begin{pmatrix}d+1&d\\\ d&d+1\end{pmatrix}$

3: $R\leftarrow\begin{pmatrix}PC&PCP\\\ C&CP\end{pmatrix}+I_{4}$

4: if $R$ is not MDS then

5: continue

6: else

7: for $D\in\mathcal{D}_{4}$ do

8: Print matrix $D^{-1}RD$

9: end for

10: end if

11:end for

Now, we generate all $4\times 4$ involutory MDS matrices over
$\mathbb{F}_{2^{m}}$ for $m=3,4,5,6,7,8$ by searching through all involutory
MDS class representative matrices over these finite fields. The following
table gives the experimental findings.

Finite Field | Number of class representatives | Choices for $D$ | Total Count
---|---|---|---
$\mathbb{F}_{2^{3}}$ | $48$ | $7^{3}$ | $7^{3}\times 48$
$\mathbb{F}_{2^{4}}$ | $71856$ | $15^{3}$ | $15^{3}\times 71856$
$\mathbb{F}_{2^{5}}$ | $10188240$ | $31^{3}$ | $31^{3}\times 10188240$
$\mathbb{F}_{2^{6}}$ | $612203760$ | $63^{3}$ | $63^{3}\times 612203760$
$\mathbb{F}_{2^{7}}$ | $26149708368$ | $127^{3}$ | $127^{3}\times 26149708368$
$\mathbb{F}_{2^{8}}$ | $961006331376$ | $255^{3}$ | $255^{3}\times 961006331376$

Table 1: Count of all $4\times 4$ involutory MDS matrices over
$\mathbb{F}_{2^{m}}$

It is worth noting that in [9], the authors demonstrate a hybrid approach for
generating all MDS matrices, resulting in a search space of $2^{9m}$ to
generate all $4\times 4$ involutory MDS matrices. In another paper [23], the
search space for generating all $4\times 4$ involutory MDS matrices is
$2^{8m}$. Additionally, a recent study [21] proposes a hybrid method to
generate all $4\times 4$ involutory MDS matrices over $\mathbb{F}_{2^{m}}$,
reducing the search space to $(2^{m}-1)^{8}$ to identify all $4\times 4$
representative involutory MDS matrices. However, the enumeration of all
$4\times 4$ involutory MDS matrices in these papers is limited to finite
fields of small order; specifically, they can find the count
$\mathbb{F}_{2^{m}}$ for $m\leq 4$. In contrast, in our proposed method, the
search space for finding all the representative matrices is $(2^{m}-1)^{5}$.
Consequently, we can explicitly enumerate the total number of $4\times 4$
involutory MDS matrices over $\mathbb{F}_{2^{m}}$ for $m=3,4,\ldots,8$.

##  6 Conclusion

In this article, we have presented some characterizations of involutory MDS
matrices of even order. Our contributions include introducing a new matrix
form for obtaining all involutory MDS matrices of even order and comparing it
with other matrix forms available in the literature. Furthermore, we have
proposed a technique to systematically construct all $4\times 4$ involutory
MDS matrices over a finite field $\mathbb{F}_{2^{m}}$. This method
significantly reduces the search space by focusing on involutory MDS class
representative matrices. Specifically, our approach involves searching for
these representative matrices within a set of cardinality $(2^{m}-1)^{5}$.
Through this method, we provide an explicit enumeration of the total number of
$4\times 4$ involutory MDS matrices over $\mathbb{F}_{2^{m}}$ for
$m=3,4,\ldots,8$. This work may be expanded to suggest the generation of all
$n\times n$ involutory MDS matrices over $\mathbb{F}_{2^{m}}$, where $n>4$.

## References

  * [1] Joan Daemen and Vincent Rijmen.  The Design of Rijndael: AES - The Advanced Encryption Standard.  Information Security and Cryptography. Springer, 2002.
  * [2] Décio Luiz Gazzoni Filho, Paulo S.L.M. Barreto, and Vincent Rijmen.  The Maelstrom-0 hash function.  In Proceedings of the 6th Brazilian Symposium on Information and Computer Systems Security, pages 17–29, 2006.
  * [3] Praveen Gauravaram, Lars Knudsen, Krystian Matusiewicz, Florian Mendel, Christian Rechberger, Martin Schläffer, and Søren Thomsen.  Grøstl - a SHA-3 candidate.  Submission to NIST, 2008, Available at http://www.groestl.info/, 09 2008.
  * [4] Jian Guo, Thomas Peyrin, and Axel Poschmann.  The PHOTON Family of Lightweight Hash Functions.  In Phillip Rogaway, editor, Advances in Cryptology – CRYPTO 2011, pages 222–239, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg.
  * [5] Kishan Chand Gupta, Sumit Kumar Pandey, and Susanta Samanta.  On the Direct Construction of MDS and Near-MDS Matrices.  arXiv: 2306.12848, 2023.  https://arxiv.org/abs/2306.12848.
  * [6] Kishan Chand Gupta and Indranil Ghosh Ray.  On constructions of involutory MDS matrices.  In International Conference on Cryptology in Africa, pages 43–60. Springer, 2013.
  * [7] Kishan Chand Gupta and Indranil Ghosh Ray.  Cryptographically significant MDS matrices based on circulant and circulant-like matrices for lightweight applications.  Cryptography and Communications, 7:257–287, 2015.
  * [8] Gülsüm Gözde Güzel, Muharrem Tolga Sakalli, Sedat Akleylek, Vincent Rijmen, and Yasemin Çengellenmiş.  A new matrix form to generate all $3\times 3$ involutory MDS matrices over $\mathbb{F}_{2^{m}}$.  Information Processing Letters, 147:61–68, 2019.
  * [9] Yogesh Kumar, P. R. Mishra, Susanta Samanta, Kishan Chand Gupta, and Atul Gaur.  Construction of all MDS and involutory MDS matrices.  arXiv: 2403.10372, 2024.  https://doi.org/10.48550/arXiv.2403.10372.
  * [10] Jérôme Lacan and Jérôme Fimes.  A Construction of Matrices with No Singular Square Submatrices.  In Gary L. Mullen, Alain Poli, and Henning Stichtenoth, editors, Finite Fields and Applications, pages 145–147, Berlin, Heidelberg, 2004. Springer Berlin Heidelberg.
  * [11] Meicheng Liu and Siang Meng Sim.  Lightweight MDS Generalized Circulant Matrices.  In Thomas Peyrin, editor, Fast Software Encryption, pages 101–120, Berlin, Heidelberg, 2016. Springer Berlin Heidelberg.
  * [12] F.J. MacWilliams and N.J.A. Sloane.  The Theory of Error Correcting Codes.  North-Holland Publishing Co., Amsterdam-New York-Oxford, 1977.
  * [13] Meltem Kurt Pehlivanog̃lu, Muharrem Tolga Sakalli, Sedat Akleylek, Nevcihan Duru, and Vincent Rijmen.  Generalisation of Hadamard matrix to generate involutory MDS matrices for lightweight cryptography.  IET Information Security, 12(4):348–355, 2018.
  * [14] Ron M. Roth and Abraham Lempel.  On MDS codes via Cauchy matrices.  IEEE Transactions on Information Theory, 35(6):1314–1319, 1989\.
  * [15] Mahdi Sajadieh, Mohammad Dakhilalian, Hamid Mala, and Behnaz Omoomi.  On construction of Involutory MDS Matrices from Vandermonde Matrices in $GF(2^{q})$.  Designs, Codes and Cryptography, 64(3):287–308, sep 2012.
  * [16] Muharrem Tolga Sakalli, Sedat Akleylek, Kemal Akkanat, and Vincent Rijmen.  On the automorphisms and isomorphisms of MDS matrices and their efficient implementations.  Turkish Journal of Electrical Engineering and Computer Sciences, 28(1):275–287, 2020.
  * [17] Susanta Samanta.  On the Counting of Involutory MDS Matrices.  arXiv: 2310.00090, 2023.  https://doi.org/10.48550/arXiv.2310.00090.
  * [18] Sumanta Sarkar and Habeeb Syed.  Lightweight Diffusion Layer: Importance of Toeplitz Matrices.  IACR Transactions on Symmetric Cryptology, 2016(1):95–113, Dec. 2016.
  * [19] C. E. Shannon.  Communication Theory of Secrecy Systems.  The Bell System Technical Journal, 28(4):656–715, 1949.
  * [20] Siang Meng Sim, Khoongming Khoo, Frédérique Oggier, and Thomas Peyrin.  Lightweight MDS Involution Matrices.  In Gregor Leander, editor, Fast Software Encryption, pages 471–493, Berlin, Heidelberg, 2015. Springer Berlin Heidelberg.
  * [21] Gökhan Tuncay, Fatma Büyüksaraçoglu Sakalli, Meltem Kurt Pehlivanoglu, Gülsüm Gözde Yilmazgüç, Sedat Akleylek, and Muharrem Tolga Sakalli.  A new hybrid method combining search and direct based construction ideas to generate all $4\times 4$ involutory maximum distance separable (MDS) matrices over binary field extensions.  PeerJ Computer Science, 9:e1577, 2023.
  * [22] Dai Watanabe, Soichi Furuya, Hirotaka Yoshida, Kazuo Takaragi, and Bart Preneel.  A New Keystream Generator MUGI.  In Joan Daemen and Vincent Rijmen, editors, Fast Software Encryption, pages 179–194, Berlin, Heidelberg, 2002. Springer Berlin Heidelberg.
  * [23] Yumeng Yang, Xiangyong Zeng, and Shi Wang.  Construction of lightweight involutory MDS matrices.  Designs, Codes and Cryptography, 89(7):1453–1483, 2021.