L Golinskii
email: golinskii@ilt.kharkov.ua

S Kupin
email: skupin@math.u-bordeaux1.fr

M Berezanskii

ON COMPLEX PERTURBATIONS OF INFINITE BAND SCHR ÖDINGER OPERATORS

Keywords: 2000 Mathematics Subject Classification. Primary: 35P20; Secondary: 47A75, 47A55 Schrödinger operators, finite-band spectrum, Lieb-Thirring type inequalities, relatively compact perturbations, resolvent identity

Let H0 = -d 2 dx 2 + V0 be an infinite band Schrödinger operator on L 2 (R) with a real-valued potential V0 ∈ L ∞ (R). We study its complex perturbation H = H0 +V , defined in the form sense, and obtain the Lieb-Thirring type inequalities for the rate of convergence of the discrete spectrum of H to the joint essential spectrum. The assumptions on V vary depending on the sign of Re V .

Introduction

Different characteristics of the distribution of the discrete spectrum for non-self-adjoint perturbations of model differential self-adjoint operators, e.g., a Laplacian on R d , a discrete Laplacian on Z d , etc., were studied in a number of papers (see Frank-Laptev-Lieb-Seiringer [START_REF] Frank | Lieb-Thirring inequalities for Schrödinger operators with complex-valued potentials[END_REF], Borichev-Golinskii-Kupin [START_REF] Borichev | A Blaschke-type condition and its application to complex Jacobi matrices[END_REF], Demuth-Hansmann-Katriel [START_REF] Demuth | On the discrete spectrum of non-selfadjoint operators[END_REF]). This paper focuses on complex perturbations of one dimensional Schrödinger operators with infinite band spectrum and certain behavior of the lengths of its gaps (the case of finite band Schrödinger operators was studied in [START_REF] Golinskii | On discrete spectrum of a complex perturbations of finite band Schrödinger operators[END_REF]).

So, consider a real-valued measurable function V 0 on R and denote by M V 0 a maximal multiplication operator by V 0 . The standing assumption is that the Schrödinger operator (0.1)

H 0 = -∆ + M V 0 , ∆ := d 2 dx 2
. is self-adjoint, H * 0 = H 0 , and its spectrum σ(H 0 ) is infinite band, i.e., (0.2)

σ(H 0 ) = σ ess (H 0 ) = I = ∞ k=1 [a k , b k ], a k → +∞.
We say that the gaps are relatively bounded if

(0.3) r = r(I) := sup k r k b k < ∞, r k := a k+1 -b k
is the length of k' gap in (0.2). A typical example here is the Hill operator with a periodic potential (see [START_REF] Reed | Methods of modern mathematical physics[END_REF]Section XIII.16]). It is well known (see [START_REF] Marchenko | A Characterization of the Spectrum of the Hill Operator[END_REF]) that r k → 0 as k → ∞ for potentials V 0 from L 2 on a period, so (0.3) obviously holds for such potentials. Furthermore, consider the form sum (0.4)

H = H 0 + M V ,
where V is a complex-valued potential. If M V is a relatively compact perturbation of H 0 , that is, dom(M V ) ⊃ dom(H 0 ), and M V (H 0 -z) -1 is a compact operator for z ∈ ρ(H 0 ), then, by the celebrated theorem of Weyl (see, e.g., [8, Section IV.5.6]), σ ess (H) = σ ess (H 0 ) and

σ(H) = I ∪ σ d (H) (disjoint union), the discrete spectrum σ d (H) of H, i.e.
, the set of isolated eigenvalues of finite algebraic multiplicity, can accumulate only on I. The main goal of the paper is to obtain certain quantitative bounds for the rate of this accumulation.

The assumption on the background V 0 looks as follows (0.5)

V 0 ≥ 0, V 0 ∈ L ∞ (R).
As for the perturbation V the conditions will vary depending on the sign of Re V . Precisely, for general V 's we assume that

(0.6) V ∈ L p (R), p ≥ 2,
and for accretive perturbations with Re V ≥ 0 we put

(0.7) V ∈ L p (R), p > 1.
Under assumptions (0.5)-(0.7) H is a well-defined, closed and m-sectorial operator, and there is a number ω 1 ≤ 0 such that

(0.8) σ(H) ⊂ N (H) ⊂ {z : Re z ≥ ω 1 },
where Theorem 0.1. Let H 0 be the Schrödinger operator (0.1) with V 0 satisfying (0.5). Assume that H 0 is the infinite band operator with the spectrum

N (H) = {(Hf, f ) : f ∈ dom(H), f 2 ≤ 1}
σ(H 0 ) = σ ess (H 0 ) = ∞ k=1 [a k , b k ], 0 ≤ a 1 < b 1 < a 2 < b 2 < . . . , a n → +∞,
and the lengths of gaps are relatively bounded (0.3). Then for the perturbation H (0.4) with V (0.6) and for each ω < ω 1 (0.8) the following Lieb-Thirring type inequality (0.9)

z∈σ d (H) d p (z, I) (|z -ω| + |ω|) 2p ≤ C(p, I) V p p (ω 1 -ω) p |ω| p-1/2 1 + V 0 ∞ a 1 + |ω| p ,
holds, where a positive constant C(p, I) depends on p and I = σ(H 0 ). Remark 0.2. If we take ω < ω 1 -1, bound (0.9) can be simplified. Indeed, now ω

1 -ω > 1, |z -ω| < |ω|(1 + |z|), 1 < a 1 + |ω| < |ω|(1 + a 1 ),
and so (0.10)

z∈σ d (H) d p (z, I) (1 + |z|) 2p ≤ C(p, I)|ω| p+1/2 (1 + V 0 ∞ ) p V p p .
There is an elementary way to specify ω and eliminate it from the final expression. The price we pay is an additional factor in the right hand side.

Theorem 0.3. Under assumptions (0.5), (0.6) (0.11)

z∈σ d (H) d p (z, I) (1 + |z|) 2p ≤ C(p, I)(1 + V 0 ∞ ) p (1 + V p ) p 2p+1 2p-1 V p p ,
where a positive constant C(p, I) depends on p and I = σ(H 0 ).

Denote D + = {|z| < 1}, D -= {|z| ≥ 1}.
Theorem 0.4. Let Re V ≥ 0. Under assumptions (0.5), (0.7) the following Lieb-Thirring type inequality holds for each 0 < ε < 1 (0.12)

z∈σ d (H)∩D + d p (z, I) |z| 1/2-ε + z∈σ d (H)∩D - d p (z, I) |z| 1/2+ε ≤ C(p, I, ε) V p p .
Remark 0.5. The only reason we restricted ourselves with the case of one dimensional Schrödinger operator H 0 as a background is that the class of multidimensional Schrödinger operators with spectra (0.2) is not well understood. Our technique works for any dimension d ≥ 1, and the corresponding problem will be elaborated on elsewhere.

Distortion for linear fractional transformations

The main analytic tool in the proof of Theorem 0.1 is the following distortion lemma for linear fractional transformations of the form

(1.1) λ ω (z) := 1 z -ω , ω ∈ R.
The argument here is quite elementary (though, rather lengthy).

Lemma 1.1. Let (1.2) I = I z = ∞ k=1 [a k , b k ], 0 ≤ a 1 < b 1 < a 2 < b 2 < . . . , a n → +∞,
and let λ ω (I) = I λ be its image under the linear fractional transformation (1.1)

λ ω (I) = I λ = ∞ k=1 [β k (ω), α k (ω)], β k (ω) = 1 b k -ω , α k (ω) = 1 a k -ω .
Then for ω < a 1 the following bounds hold: for Re z < a 1 or Re z ∈ I

(1.3) d(λ ω (z), λ ω (I) d(z, I z ) > 1 3|z -ω|(|z -ω| + a 1 -ω) ; for b k < Re z < a k+1 , k = 1, 2, . . . (1.4) d(λ ω (z), λ ω (I) d(z, I z ) ≥ 1 2|z -ω| 2 1 + a k+1 -b k b k -ω -1
.

Moreover, if ω ≤ 0 and the gaps are relatively bounded (0.3), then the unique bound is valid

(1.5) d(λ ω (z), λ ω (I) d(z, I z ) ≥ 1 5(1 + r(I)) 1 |z -ω|(|z -ω| + a 1 -ω) , z ∈ C\I.
Proof. With no loss of generality we can assume that a 1 > 0.

We begin with the case ω = 0 and put

λ 0 = λ = z -1 . If z = x + iy and x = Re z ≤ 0, then Re λ = x|z| -2 ≤ 0 and so (1.6) d(λ, I λ ) d(z, I z ) = |λ| |z -a 1 | = 1 |z||z -a 1 | ≥ 1 |z|(|z| + a 1 )
.

Similarly, if x ∈ I z , then x ≥ a 1 and 0 < Re λ = x |z| 2 ≤ 1 x ≤ a -1 1 = α 1 , d(λ, [0, α 1 ]) = |Im λ| = |y| |z| 2 .
Since now d(z, I z ) = |y|, we have

(1.7) d(λ, I λ ) d(z, I z ) ≥ d(λ, [0, α 1 ]) d(z, I z ) = 1 |z| 2 > 1 |z|(|z| + a 1 )
. Consider now the case when

x = Re z / ∈ I z . Fix x in k's gap, (1.8) b k < x < a k+1 , k = k(x) = 0, 1, . . .
(we put b 0 = 0 and treat (b 0 , a 1 ) as a number zero gap). Then

d(z, I z ) = min(|z-b k |, |z-a k+1 |), k = 1, 2, . . . , d(z, I z ) = |z-a 1 |, k = 0.
Define two sets of positive numbers

u j = u j (x), v j = v j (x), j = k + 1, k + 2, . . . by equalities Re (λ(x + iu j )) = x x 2 + u 2 j = α j , Re (λ(x + iv j )) = x x 2 + v 2 j = β j ,
or, equivalently,

u j (x) = x(a j -x), v j (x) = x(b j -x).
We also put v k = 0, so

0 = v k < u k+1 < v k+1 < u k+2 < v k+2 < . . . , u n , v n → ∞, n → ∞.
While the point z traverses the line x + iy, y ∈ R, its image λ(z) describes a circle with diameter [0, 1/x]. We discern the following two cases.

Case 1. Assume that λ lies over the "gaps for λ". For each k = 0, 1, . . . there are two options for λ: interior gaps For gaps (1.9) we have

(1.9) Re λ ∈ (α j+1 , β j ) ⇐⇒ v j < |y| < u j+1 , j = k + 1, k + 2, . . .
(1.11) d(λ, I λ ) = min(|λ-α j+1 |, |λ-β j |) = 1 |z| min |z -a j+1 | a j+1 , |z -b j | b j .
Define an auxiliary function h on the right half-line

h(t) = h(t, z) := |z -t| t = x t -1 2 + y 2 , t > 0.
Clearly, h is monotone increasing on (x, +∞) and decreasing on (0, x) with the minimum h(x) = |y|. Hence (1.11) and (1.8) give

d(λ, I λ ) = min(h(a j+1 , z), h(b j , z)) |z| ≥ h(b j , z) |z| ≥ h(b k+1 , z) |z| ≥ h(a k+1 , z) |z| = |z -a k+1 | a k+1 |z| . Since by (1.8) d(z, I z ) ≤ |z -a k+1 |, we see that (1.12) d(λ, I λ ) d(z, I z ) ≥ 1 a k+1 |z| .
For gaps (1.10) let first k ≥ 1. Then as above in (1.11)

d(λ, I λ ) = 1 |z| min |z -a k+1 | a k+1 , |z -b k | b k , but it is not clear now which term prevails. If |z -a k+1 | ≤ |z -b k | then d(z, I z ) = |z -a k+1 | and d(λ, I λ ) d(z, I z ) = 1 |z| min 1 a k+1 , |z -b k | b k |z -a k+1 | = 1 a k+1 |z| . Otherwise |z -a k+1 | > |z -b k | implies d(λ, I λ ) d(z, I z ) = 1 |z| min 1 b k , |z -a k+1 | a k+1 |z -b k | ≥ 1 a k+1 |z| .
Next, for k = 0 one has 0 < x < a 1 , and in case (1.10)

d(λ, I λ ) = |λ -α 1 | = |z -a 1 | a 1 |z| , d(z, I z ) = |z -a 1 |,
and so

(1.13) d(λ, I λ ) d(z, I z ) = 1 a 1 |z| .
Finally, in the case of "gaps for λ" we come to the bound

(1.14) d(λ, I λ ) d(z, I z ) ≥ 1 a k+1 |z| , k = 0, 1, . . . .
A modified version of (1.14) will be convenient in the sequel.

For k ≥ 1 in view of |z| ≥ x > b k we have 1 a k+1 |z| ≥ b k a k+1 |z| 2
and so for k = 1, 2, . . .

(1.15) d(λ, I λ ) d(z, I z ) ≥ 1 |z| 2 1 + a k+1 -b k b k -1 = 1 |z| 2 1 + r k b k -1 , r k = a k+1 -b k is the length of k's gap. Similarly, for k = 0 one has from (1.13) (1.16) d(λ, I λ ) d(z, I z ) ≥ 1 |z|(|z| + a 1 )
.

Case 2. Assume that λ lies over the "bands for λ"

(1.17) Re λ ∈ [β j , α j ] ⇐⇒ u j ≤ |y| ≤ v j , j = k + 1, k + 2, . . . . Now d(λ, I λ ) = |Im λ| = |y| |z| 2 , d(z, I z ) ≤ |z -a k+1 | ≤ |y| + a k+1 -x = |y| + u 2 k+1 x ≤ |y| 1 + u k+1 x = |y| 1 + a k+1 -x x , so that (1.18) d(λ, I λ ) d(z, I z ) ≥ 1 |z| 2 1 + a k+1 x -1 -1
.

For k ≥ 1 (interior gap for z) inequality (1.18) can be simplified in view of

x > b k (1.19) d(λ, I λ ) d(z, I z ) ≥ 1 |z| 2 1 + r k b k -1
.

Let now k = 0, i.e., 0 < x = Re z < a 1 . In our case d(z,

I z ) = |z -a 1 | and |y| ≥ u 1 = x(a 1 -x).
If |y| ≥ 2x then |y| ≥ 2 3 |z| and so

(1.20) d(λ, I λ ) d(z, I z ) = |y| |z| 2 |z -a 1 | ≥ 2 3 1 |z|(|z| + a 1 )
.

Otherwise, |y| < 2x implies 2 √ x > √ a 1 -x, x > a 1 5 .
It follows now from (1.18) with k = 0 that

(1.21) d(λ, I λ ) d(z, I z ) ≥ 1 3|z| 2 > 1 3 1 |z|(|z| + a 1 )
.

We can summarize the results obtained above in the following two bounds from below. A combination of (1.6), (1.7), (1.16)

) d(z, I z ) ≥ 1 γ k |z| 2 , γ k = max 1 + r k b k , 1 + r k b k , b k < Re z < a k+1 , k = 1, 2, . . . . (1.23) Since γ k < 2(1 + r k /b k )
, the latter can be written as

(1.24) d(λ, I λ ) d(z, I z ) ≥ 1 2|z| 2 1 + r k b k -1 , b k < Re z < a k+1 , k = 1, 2, . . . .
To work out the general case ω = 0 and prove (1.3) and (1.4), it remains only to shift the variable and apply the results just obtained to the shifted sequence of bands

I z (ω) = k≥1 [a k -ω, b k -ω].
The final statement follows from a simple observation

r k b k -ω ≤ r k b k ≤ r.
The proof is complete.

Lieb-Thirring type inequalities

The key ingredient of the proof of our main statements is the following result of Hansmann [START_REF] Hansmann | Variation of discrete spectra for non-selfadjoint perturbations of selfadjoint operators[END_REF]Theorem 1]. Let A 0 = A * 0 be a bounded self-adjoint operator on the Hilbert space, and let A be a bounded operator with A-A 0 ∈ S p , p > 1. Then (2.1)

λ∈σ d (A) d p (λ, σ(A 0 )) ≤ K p A -A 0 p Sp ,
K p is an explicit (in a sense) constant, which depends only on p. We set

A 0 (ω) = R(ω, H 0 ) = (H 0 -ω) -1 , A(ω) = R(ω, H) = (H -ω) -1 ,
ω is defined above, and ω ∈ ρ(H 0 ) ∩ ρ(H) in view of (0.2) and (0.8). Let λ = λ ω (z) = (z -ω) -1 . The Spectral Mapping Theorem implies that

λ ∈ σ d (A(ω)) (λ ∈ σ(A 0 (ω))) ⇐⇒ z ∈ σ d (H) (z ∈ σ(H 0 )) .
Proof of Theorem 0.1. The second resolvent identity reads

R(z, H) -R(z, H 0 ) = -R(z, H)M V R(z, H 0 ), z ∈ ρ(H) ∩ ρ(H 0 ).
We wish to show that this difference belongs to S p and to obtain the bound for its S p -norm. First, we have

W = W (z) := M V R(z, H 0 ) = M V (-∆ -z) -1 (-∆ -z)(H 0 -z) -1 = M V (-∆ -z) -1 (1 -M V 0 (H 0 -z) -1 ), (2.2)
and so

W (z) Sp ≤ M V R(z, -∆) Sp I -M V 0 R(z, H 0 ) , z ∈ ρ(H) ∩ ρ(H 0 ).

It is clear that

I -M V 0 R(z, H 0 ) ≤ 1 + V 0 ∞ d(z, I) = 1 + V 0 ∞ |a 1 -z| , Re z < 0.
Next, write

M V (-∆ -z) -1 = V (x)g z (-i∇), g z (x) = (x 2 -z) -1 , x ∈ R.
By [START_REF] Simon | Trace ideals and their applications[END_REF]Theorem 4.1]

M V (-∆ -z) -1 Sp ≤ (2π) -1/p V p g z p , p ≥ 2.
Since 2|t -z| 2 ≥ (t + |z|) 2 for t ≥ 0 and Re z < 0, we have

g z p ≤ √ 2 g -|z| p
and so

M V (-∆ -z) -1 Sp ≤ C 1 |z| 1-1/2p V p , C 1 = C 1 (p) = √ 2 1 2π ∞ -∞ dx (x 2 + 1) p 1/p . (2.3) Thus, (2.4) W (z) Sp ≤ C 1 (p) V p |z| 1-1/2p 1 + V 0 ∞ |a 1 -z| , Re z < 0.
We put z = ω < ω 1 . Relation (0.8) implies in view of [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Theorem V.3.2] (2.5)

R(ω, H) ≤ 1 d(ω, N (H)) ≤ 1 ω 1 -ω ,
and the combination of (2.4) and (2.5) leads to the following bound for each

ω < ω 1 R(ω, H) -R(ω, H 0 ) p Sp ≤ R(ω, H) p W (ω) p Sp ≤ C 2 (p) V p p (ω 1 -ω) p |ω| p-1/2 1 + V 0 ∞ a 1 + |ω| p . (2.6)
We go back to (2.1) with

A 0 = A 0 (ω) = R(ω, H 0 ), A = A(ω) = R(ω, H),
so by the Spectral Mapping Theorem, in the notation of Lemma 1.1 we have (2.7)

λ∈σ d (A(ω)) d p (λ, σ(A 0 (ω)) = z∈σ d (H) d p (λ ω (z), λ ω (I)) ≤ K p R(ω, H)-R(ω, H 0 ) p Sp .
We apply Lemma 1.1 in the form (1.5) to obtain

z∈σ d (H) d p (z, I) |z -ω| p (|z -ω| + a 1 + |ω|) p ≤ C 3 (p, I) V p p (ω 1 -ω) p |ω| p-1/2 1 + V 0 ∞ a 1 + |ω| p ,
and (0.9) follows. The proof is complete.

2

The proof of Theorem 0.1 shows that bound (0.9) essentially depends on the parameter ω. Roughly speaking, it comes from a bound from below of inf Re σ(H), and so it seems to be rather important to estimate this quantity in terms of V 0 and V only.

Proof of Theorem 0.3.

Put Ω = {Re z < 0} {|a 1 -z| > (1 + V 0 ∞ )} |z| 1-1/2p > 4C 1 (p)(1 + V p ) , C 1 is defined in (2.
3). We wish to show that Ω ⊂ ρ(H 0 ) ∩ ρ(H). Indeed, by

) W (z) ∞ ≤ W (z) Sp ≤ V p 2(1 + V p ) < 1 2 , z ∈ Ω, so I + W (z) is invertible and (I + W (z)) -1 < 2. An application of the identity H -z = (1 + W (z))(H 0 -z) completes the proof of our claim. (2.4), (2.8
Next, write the difference of the resolvents in another way

R(z, H) -R(z, H 0 ) = -R(z, H 0 )(1 + W (z)) -1 W (z) to obtain for z ∈ Ω R(z, H) -R(z, H 0 ) Sp ≤ R(z, H 0 ) (1 + W (z)) -1 W (z) Sp ≤ V p |a 1 -z|(1 + V p ) ≤ V p (1 + V 0 ∞ )(1 + V p ) .
It is clear by the definition of Ω that if t ∈ Ω, t < 0, then {Re z ≤ t} ⊂ Ω. Take z = ω < 0 so that (2.9)

|ω | 2 = a 1 2 + 1 + V 0 ∞ + (4C 1 (1 + V p )) 1/(1-1/2p) .
It is easy to check that {z : Re z < ω 2 } ⊂ Ω, so, in particular, ω ∈ Ω and hence

R(ω , H) -R(ω , H 0 ) Sp ≤ V p (1 + V 0 ∞ )(1 + V p ) .
Once again, (2.1) says (2.10)

λ∈σ d (A(ω )) d p (λ, σ(A 0 (ω ))) ≤ K p V p (1 + V 0 ∞ )(1 + V p )
p for p > 1 and, using Lemma 1.1, we come to

z∈σ d (H) d p (z, I) |z -ω | p (|z -ω | + a 1 + |ω |) p ≤ C 4 (p, I) V p (1 + V 0 ∞ )(1 + V p ) p .
By the choice of ω (2.9), we have Re z ≥ ω /2 for z ∈ σ d (H), and so

|z -ω | ≥ |ω | 2 > |a 1 | + |ω | 4 , and
|z -ω | + a 1 + |ω | < 5|z -ω |, |z -ω |(|z -ω | + a 1 + |ω |) < 5|z -ω | 2 .
Next As for the upper bound, we follow the line of reasoning from [6, Proof of Theorem 3.2], where such bound was proved in the case V 0 = 0. As a matter of fact, the argument goes through under assumption (0.5) as well. At any rate, we have and then integrate the latter inequality with respect to a ∈ (0, ∞). The proof is complete.

2

Figure 1 .

 1 Figure 1. Sets I = σ(H 0 ) and λ ω (I) with map λ ω (z) = 1 z-ω .

and the rightmost gap ( 1 . 10 )

 110 Re λ ∈ (α k+1 , 1/x) ⇐⇒ 0 < |y| < u k+1 .

2

 2 , |z -ω | ≤ (1 + |z|)(1 + |ω |) ≤ 2|ω |(1 + |z|) and hence z∈σ d (H) d p (z, I) (1 + |z|) 2p ≤ C 5 (p, I)|ω | 2p V p (1 + V 0 ∞ )(1 + V p ) p ≤ C 6 (p, I)(1 + V 0 ∞ ) p (1 + V p )Proof of Theorem 0.4. For accretive perturbations we have σ(H) ⊂ {z : Re z ≥ 0}, so one can take ω 1 = 0.The lower bound for the difference of resolvents is the same as above in Theorem 0.1. It is a consequence of the result of Hansmann and Lemma 1.1R(ω, H) -R(ω, H 0 ) p Sp ≥ z∈σ d (H) d p (z, I) (|z -ω| -ω) 2p , p > 1.

R

  (-a, H) -R(-a, H 0 ) p Sp ≤ C 5 (p) a 2p-1/2 V p p , a := -ω > 0.Since √ 2|z + a| ≥ |z| + a for Re z > 0 and a > 0, we come to the following inequality(2.11) z∈σ d (H) d p (z, I) (|z| + a) 2p ≤ C 6 (p) a 2p-1/2 V p p , a > 0.As in [6, Proof of Theorem 3.3], we multiply (2.11) through by (1 + a) 2ε , ε > 0, to obtainz∈σ d (H) d p (z, I) a 2p-3/2+ε (|z| + a) 2p (1 + a) 2ε ≤ C 6 (p) V p p a 1-ε (1 + a) 2ε , a > 0,

Acknowledgments. The authors are grateful to S. Denisov, F. Gesztesy and M. Hansmann for valuable discussions on the subject of the article. The paper was prepared during the visit of the first author, funded by the IdEx programm of University of Bordeaux, to the Institute of Mathematics of Bordeaux (IMB UMR5251). He would like to thank this institution for the hospitality.


Alban Derrien

Emmanuel Boutillon

Senior Member, IEEE Audrey Cerqueus
email: audrey.cerqueus@emse.fr

Additive, Structural and Multiplicative Transformations for the Construction of Quasi-Cyclic LDPC matrices

The construction of a Quasi-Cyclic Low Density Parity-Check (QC-LDPC) matrix is usually carried out in two steps. In the first step, a prototype matrix is defined according to certain criteria (size, girth, check and variable node degrees, etc.). The second step involves expansion of the prototype matrix. During this last phase, an integer value is assigned to each nonnull position in the prototype matrix corresponding to the rightrotation of the identity matrix. The problem of determining these integer values is complex. State-of-the-art solutions use either some mathematical constructions to guarantee a given girth of the final QC-LDPC code or a random search of values until the target girth is satisfied. In this paper, we propose an alternative/complementary method that reduces the search space by defining large equivalence classes of topologically identical matrices through row and column permutations using additive, structural and multiplicative transformations. Selecting only a single element per equivalence class can reduce the search space by a few orders of magnitude. Then, we use the formalism of constraint programming to list the exhaustive sets of solutions for a given girth and a given expansion factor. An example is presented in all sections of the paper to illustrate the methodology.

I. INTRODUCTION

Low Density Parity-Check (LDPC) codes are a family of error correction codes used in many communication standards (TV broadcasting [START_REF]DVB standard committee[END_REF], Wi-Fi [START_REF]IEEE standard committee[END_REF] and next generation cellular networks [START_REF]5G standard[END_REF] among others). LDPC codes were invented in the 1960s by Gallager [START_REF] Gallager | Low-density parity-check codes[END_REF]. Their success is due to the existence of a simple iterative decoding algorithm known as the Belief Propagation (BP) algorithm [START_REF] Mackay | Near shannon limit performance of low density parity check codes[END_REF]. The BP decoding algorithm (and its simplified versions [START_REF] Boutillon | Iterative Decoding of Concatenated Convolutional Codes: Implementation Issues[END_REF]) exchanges messages (belief on the value of a bit) in a bipartite graph composed of variable and parity-check nodes. LDPC codes were generalized by Davey and McKay in 1998 over a finite field arithmetic [START_REF] Davey | Low density parity check codes over GF(q)[END_REF]. These codes are called Non-Binary LDPC (NB-LDPC) codes. One of the main advantages of NB-LDPC codes is that good codes can be constructed with the degree of the variable nodes set to 2 (each variable is thus connected to exactly two parity-check nodes). NB-LDPC code has been adopted in a recent CCSDS standard [START_REF]Consultative Committee for Space Data Systems (CCSDS), Telecommand Sync and Channel Coding Specification using advanced Block Codes[END_REF].

To obtain a hardware-friendly error control code, some structures can be imposed on its structure. In particular, memory conflicts that can appear in a high speed decoder can be resolved using Quasi-Cyclic-LDPC (QC-LDPC) matrices. This type of matrices is used in many communication standards. There is a large literature on the construction of QC-LDPC matrices, some of them based on mathematical properties of particular sets, others on random generation, then selection, and others on exhaustive search [START_REF] Liu | Quasi-Cyclic LDPC codes: Construction and rank analysis of their parity-check matrices[END_REF], [START_REF] Divsalar | Capacityapproaching Protograph Codes[END_REF], [START_REF] Divsalar | Low-rate LDPC codes with simple protograph structure[END_REF], [START_REF] Thorpe | Methodologies for designing LDPC codes using protographs and circulants[END_REF], [START_REF] Ranganathan | On the girth of (3, L) Quasi-Cyclic LDPC Codes based on Complete Protographs[END_REF], [START_REF] Zhang | Construction of Girth-Eight QC-LDPC Codes from Greatest Common Divisor[END_REF], [START_REF] Wang | Construction of high-girth QC-LDPC codes[END_REF], [START_REF] O'sullivan | Algebraic construction of sparse matrices with large girth[END_REF], [START_REF] Zhang | Deterministic Girth-Eight QC-LDPC Codes with Large Column Weight[END_REF], [START_REF] Vafi | Combinatorial design based Quasi Cyclic LDPC codes with girth eight[END_REF], [START_REF] Tasdighi | Symmetrical Constructions for Regular Girth-8 QC-LDPC Codes[END_REF], [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF].

The main contribution of this paper is to propose a two step method to construct QC-LDPC matrices. The first step consists of partitioning the search space into equivalence classes of large size by means of additive, structural and multiplicative transformations. Thus, testing only one element per class allows us to reduce by a few orders of magnitude the number of solutions to be explored. The second step consists of using constraint programming tools to enumerate all solutions when they exist. This method can be applied to any type of QC-LDPC code. To explain and illustrate the method, we use an example of the construction of good QC-LDPC matrices from a given prototype matrix. Note that some of the ideas of the paper were independently developed by Tasdighi et al. and recently published in [START_REF] Tasdighi | Efficient Search of Girth-Optimal QC-LDPC Codes[END_REF]. More specifically, the idea of structural and multiplicative transformations was already mentioned in [START_REF] Tasdighi | Symmetrical Constructions for Regular Girth-8 QC-LDPC Codes[END_REF] and has been used to reduce the search space to find high-girth QC-LDPC matrices. But, as recognized by M. Tasdighi, the present approach is much more general than that previously proposed in [START_REF] Tasdighi | Symmetrical Constructions for Regular Girth-8 QC-LDPC Codes[END_REF].

The rest of the paper is organized as follows. Section II gives the background of QC-LDPC code construction. Section III presents an equivalence relation between QC-LDPC matrices. Section IV shows how to select a single element of each class. Section V presents a new equivalence relation based on multiplicative properties. Section VI gives the number of solutions of minimal girth for an L = 3 rows, J = 6 lines protograph with several expansion factors. Finally, section VII sets out the conclusions of this study.

II. CONSTRUCTION OF QC-LDPC MATRICES

In this section, we review the construction of a family of LDPC matrices well suited for hardware implementation called Quasi-Cyclic LDPC matrices. Then, we discuss the conditions that have to be satisfied to obtain QC-LDPC matrices with good topological properties. Finally, we provide some notation to describe the permutations.

A. Definition of a QC-LDPC matrix

An LDPC code can be represented by a bipartite graph that contains two types of nodes: the variable nodes and the check nodes. Variable nodes (respectively check nodes) are only connected to check nodes (respectively variable nodes). In the case of a regular LDPC code, the number d v of check nodes connected to a given variable node and the number d c of variable nodes connected to a given check node are constant (d v and d c are called the variable node degree and the check node degree, respectively).

LDPC codes were generalized to NB-LDPC codes in 1998 by Mackay and Neal [START_REF] Mackay | Near shannon limit performance of low density parity check codes[END_REF]. These authors show that good NB-LDPC codes can be constructed with a constant variable node degree d v = 2. The Belief Propagation (BP) algorithm is equivalent to the Maximum A Posteriori (MAP) decoder if the graph is cycle free [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF] and its performance is close to MAP decoder if the bipartite graph representing the code does not have small cycles. Thus, the main topological objectives, in the design of a bipartite graph, are, first, to maximize its girth g (i.e. the length of its minimal cycle), and second, to minimize its multiplicity (i.e. the number of cycles of length g) denoted by M(g)

The bipartite graph can be represented by a matrix H, called the parity-check matrix, where each variable of the code is associated with a column of H and each parity-check is associated with a row. If an edge exists between a check node i and a variable node j, then H(i, j) is not equal to zero, otherwise, H(i, j) = 0. In the case of a (d v , d c ) regular bipartite graph, each row of H contains d c non-zero values and each column of H contains d v non zero values. Note that permuting the rows and the columns of H does not affect the topology of the associated graph (see section III).

A QC-LDPC matrix is constructed in two steps. First, a protograph H of size J × L is constructed. Then, the protograph matrix is expanded, or lifted, by a factor N to obtain a (JN, LN ) matrix [START_REF] Thorpe | Low-density parity-check (LDPC) codes constructed from protographs[END_REF], [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF]. The simplest protograph matrix for

d v = 2, d c = 4 is the matrix H 2 of size J = 2, L = 4 defined as H 2 = 1 1 1 1 1 1 1 1 . (1)
Throughout this paper, we illustrate the proposed method using the example of QC-LDPC matrices with a rate 1/2 generated from the matrix H 3 of size J = 3, L = 6, defined as

H 3 =   1 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 1   . ( 2
)
This matrix is well suited for the NB-LDPC code, where d v = 2 gives a good code [START_REF] Mackay | Near shannon limit performance of low density parity check codes[END_REF]. Note that, for NB-LDPC code, we also need to optimize the edge values (i.e. value of non-null coefficients of H). This point is not discussed further here since the present study is focused on the topological properties of the graph associated with the expanded matrix. It can be noted that H 3 is unique up to rows/columns permutation (see section III).

The lifting process generates a matrix of size (JN, LN ) from H by replacing each 0 value in H by O N , the (N, N ) zero matrix, and each 1 value by an (N, N ) matrix I N a , where a is the shift value associated with the non-null position, see [START_REF] Gallager | Low-density parity-check codes[END_REF]. I N a is defined as

I N a = I N a (i, i + a mod N ) = 1, ∀i ∈ {0, . . . , N -1} 0, otherwise.
(3) Note that I N 0 represents the (N, N ) identity matrix. The matrix I N a , a ∈ {0, . . . , N -1} is called a circulant matrix [START_REF] Mitchell | Quasi-Cyclic LDPC Codes Based on Pre-Lifted Protographs[END_REF]. The index a is called the permutation shift. The QC-LDPC matrix H is thus fully described by the protograph matrix H, the factor of expansion N and by the shift values associated with each non-zero value of the circulant permutation matrices [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF].

For example, matrix H obtained from H 3 is defined as

H =   I N a1 I N b1 I N c1 I N d1 O N O N I N a2 I N b2 O N O N I N e2 I N f2 O N O N I N c3 I N d3 I N e3 I N f3   (4)
Matrix H defines a type-I QC-NB-LDPC code since non zero values of the prototype matrix are expanded to a single circulant matrix (a prototype matrix with at least one double circulant matrix is called a type-II protograph [START_REF] Smarandache | On regular quasicyclic LDPC codes from binomials[END_REF]). Since the protograph is of type-I, it is convenient to represent the expanded matrix H as

H =   a 1 b 1 c 1 d 1 -1 -1 a 2 b 2 -1 -1 e 2 f 2 -1 -1 c 3 d 3 e 3 f 3   , (5)
where

H(i, j) represents I N H(i,j) if H(i, j) is a positive integer, or O N if H(i, j) = -1.
In the sequel, -1 is replaced by a simple dot to simplify the notation.

B. Topological properties of QC-LDPC matrices

Due to the correspondence between the Tanner graph and the parity-check matrix, a cycle in the Tanner graph is equivalent to a cycle in the parity-check matrix. In the parity-check matrix, the cycle is composed of a circular consecutive sequence of horizontal and vertices vertices [START_REF] Gallager | Low-density parity-check codes[END_REF]. The construction objective in the lifting process is to maximize the girth of the Tanner graph.

Each cycle of length in the protograph lifts into one or several cycles of length in the expanded graph, with being a multiple of . The necessary and sufficient condition for a cycle of length to lift into cycles of length strictly greater than is given in [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF]. Let N be the expansion size and C a cycle of length in the protograph. Let (e 1 , e 2 ,. . . , e ) represent the sequence of edges of cycle C in H and let (a 1 , a 2 ,. . . , a ) be the corresponding permutation shifts. Then, it is possible to associate the cycle C with a value φ(C) defined as

φ(C) = i=1 (-1) i+1 a i mod N, (6)
where by convention, the first term of the cycle is on the upper left position of the cycle in the parity-check matrix (see Fig. 1 on page 5). Since the graph is bipartite, the length of a cycle is always even. Thus, the direction of the cycle does not impact the value of φ(C) since (-1) -i = (-1) i . Lemma 2.1 (Cycle breaking): Cycle C lifts into cycles of length > if and only if φ(C) = 0.

Proof: See [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF] The problem can now be modeled as finding the shift values of the circulant matrices to maximize the girth of the expanded matrix, or alternatively, finding the minimum expansion factor N required to obtain a given girth.

C. Equivalence relation

A permutation π of size n is a bijection function of set N = {0, 1, . . . , n-1} to itself: π(N ) = {π(0), π(1), . . . , π(n-1)}. The identity permutation is denoted Id. It is possible to represent the permutation π by a (n, n) matrix, called the permutation matrix P π , defined as P π (i, j) = 1 if j = π(i), 0 otherwise. The inverse permutation of π will be denoted π -1 and verifies π -1 (π) = π(π -1 ) = Id. Moreover, (P π ) -1 = P π -1 and P π -1 × P π = P π × P π -1 = I n 0 with I n 0 the identity matrix of size (n×n). Finally, (P π ) -1 = P t π , where P t means the transposed matrix of P .

Definition 2.2 (Equivalence of matrices): Two matrices A and B of size (n, m) are said to be "permutation equivalent" if and only if there exist two permutation matrices P r (index r for row) of size (n × n) and P c (index c for column) of size (m × m) such that A = P r × B × P c . The relation is denoted as A ≡ B.

Theorem 2.3: The permutation equivalent relation defines an equivalence relation.

Proof: The relation is reflexive

(A ≡ A), symmetrical (A = P r × B × P c implies that B = P -1 r × A × P -1 c ) and associative (A = P a r × B × P a c and B = P b r × C × P b c implies A = P a r P b r × C × P b c P a c
). Let A and B be two lifted matrices from the same prototype matrix H, then, if A ≡ B, then A and B are equivalent expanded matrices of H. This relation is denoted A ≡ H B, the subscript H being omitted when there is no ambiguity.

To summarize, if A ≡ H B, then A and B correspond to the same Tanner graph up to a renumbering of variable and check nodes, and thus, represent the same code. The aim of this paper is twofold: to define classes of equivalent solutions using permutations and present static symmetry breaking [START_REF] Chu | Dominance breaking constraints[END_REF] to identify a solution per class of equivalence.

III. ADDITIVE AND STRUCTURAL TRANSFORMATIONS

In the following, "equivalence class" should be understood in terms of the equivalence relation defined in Section II-C.

The main goals of this approach are (i) to determine the number of equivalence classes (possibly, how many elements there are in each equivalence class), and (ii) to obtain a unique representative for each equivalence class. To address the above issues, two transformations that preserve the equivalence relation are presented: the "additive transformations", already proposed in [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF] and the "structural transformation". In the sequel, the size N of the matrix is omitted when there is no ambiguity.

A. Additive transformation

In this section, we introduce the additive transformation, a state-of-the-art permutation [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF] which reduces the degree of freedom of the problem example from 12 (the number of variables in ( 5)) down to 4.

Definition 3.1 (Circulant permutation): Let π + a be the permutation over {0, 1, . . . n -1} defined as π + a (i) = i + a mod n. The permutation matrix P + a associated with π + a is also equal to the circulant matrix P + a = I n a . Property 3.2: The product of two circulant matrices P + a and P + b is equal to the circulant matrix P + a × P + b = P + a+b [START_REF] Mitchell | Quasi-Cyclic LDPC Codes Based on Pre-Lifted Protographs[END_REF]. Property 3.3 (Row rearrangement): Let M be an n × m matrix and π r a permutation of size n. Then, P πr × M is the (n, m) matrix obtained by permuting the rows of M according to π r .

Property 3.4 (Column rearrangement): Let M be a n × m matrix and π c a permutation of size m. Then M × P πc is the (n, m) matrix obtained by permuting the columns of M according to π -1 c . Proof: Since for any matrix A, A = (A t ) t , then M × P πc = (P t πc × M t ) t . The product P t πc × M t permutes the row of M t according to P t πc , i.e., according to permutation π -1 c . Thus, its transpose gives the permutation of the column of M according to π -1 c . Let us define C + r as the (JN, JN ) block diagonal permutation matrix defined as C + r = diag(P + r(1) , P + r(2) , . . . P + r(J) ), with r = (r(1), r(2), . . . r(J)) representing a vector of integers.

Property 3.5 (Row shift):

H r + = C + r ×
H c + = H × C + c is the expanded matrix defined by H c + with H c + (i, j) = H(i, j) + c(j) ∀i ∈ {1, 2, . . . , J}, j ∈ {1, 2, . . . , L}.
Let H N 3 denote the set of matrices obtained by an expansion of a factor N of the prototype matrix H 3 . Each matrix H ∈ H N 3 can be represented as [START_REF] Mackay | Near shannon limit performance of low density parity check codes[END_REF]. From r 1 = (-a 1 , -a 2 , 0) we can construct the (3N, 3N ) permutation matrix C + r1 = diag(P + -a1 , P + -a2 , P + 0 ) and H r + 1 = C + r1 × H. According to property 3.5, H r + 1 can be represented as

H r + 1 =   0 b 1 -a 1 c 1 -a 1 d 1 -a 1 . . 0 b 2 -a 2 . . e 2 -a 2 f 2 -a 2 . . c 3 d 3 e 3 f 3   (7)
Then, multiplying H r + 1 by the (6N, 6N ) permutation matrix

C + c with c = (0, a 1 -b 1 , a 1 -c 1 , a 1 -d 1 , a 2 -e 2 , a 2 -f 2 ), gives H r + 1 ,c + = H r + 1 × C + r . According to property 3.6, H r + 1 ,c + can be represented as H r + 1 ,c + =   0 0 0 0 . . 0 b2-a2 +a1-b1 . . 0 0 . . c3 +a1-c1 d3 +a1-d1 e3 +a2-e2 f3 +a2-f2   (8)
Finally, multiplying

H r + 1 ,c + on the left by C + r2 with r 2 = (0, 0, -c 3 -a 1 + c 1 ) gives H r + 1 ,c + ,r + 2 = C + r2 × H r + 1 ,c + . According to property 3.6, H r + 1 ,c + ,r + 2 can be represented as H r + 1 ,c + ,r + 2 , which is defined as H r + 1 ,c + ,r + 2 =   0 0 0 0 . . 0 b . . 0 0 . . 0 d e f   , (9)
where

         b = b 2 -a 2 + a 1 -b 1 , d = d 3 -d 1 -c 3 + c 1 , e = e 3 + a 2 -e 2 -c 3 -a 1 + c 1 , f = f 3 + a 2 -f 2 -c 3 -a 1 + c 1 . (10)
Property 3.7 (bdef-pattern): Any matrix H ∈ H N 3 is equivalent to a matrix of type [START_REF] Liu | Quasi-Cyclic LDPC codes: Construction and rank analysis of their parity-check matrices[END_REF].

Thus the study of the set H N 3 can be restricted to the study of matrices defined in [START_REF] Liu | Quasi-Cyclic LDPC codes: Construction and rank analysis of their parity-check matrices[END_REF]. To simplify notation, a matrix H ∈ H N 3 can be also represented by the 4-tuple H = b, d, e, f N . In this section, as mentioned before, we apply the results of [START_REF] Fossorier | Quasi-cyclic low-density parity-check codes from circulant permutation matrices[END_REF] to the case of the prototype matrix H 3 . This transformation is very general and can be applied to any QC-LDPC matrix. In fact, the additive transformation can be used to set the first non-negative coefficients of each column of H to zero and the first non-negative coefficient of each row of H to zero, thus reducing the dimension of the search space from Jd c (or Ld v , i.e., the number of non-negative values of H) down to Jd c -J -L + 1.

It is worth noticing that matrices of the form H = b, d, e, f do not give a unique representative of each equivalence class, i.e., two matrices H = b, d, e, f and H = b , d , e , f may be equivalent, even if (b, d, e, f ) = (b , d , e , f ). The next section presents the first contribution of this paper: the structural transformation that uses a different approach to identify equivalence between matrices.

B. Structural transformation

Let us consider the set P J,L of pair of permutations (π r , π c ) of size equal to L and J, respectively, that satisfies

H = P πr × H × P πc , (11)
where P πr and P πc are, respectively, the permutation matrices associated with π r and π c . For example, (π r , π c ) = ({2, 3, 1}, {4, 3, 6, 5, 2, 1}) belongs to P 3,6 since

P πr × H 3 =   0 1 0 0 0 1 1 0 0   ×   1 1 1 1 . . 1 1 . . 1 1 . . 1 1 1 1   =   1 1 . . 1 1 . . 1 1 1 1 1 1 1 1 . .   , then (P πr × H 3 ) × P πc gives   1 1 . . 1 1 . . 1 1 1 1 1 1 1 1 . .   ×         0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0         (12) =   1 1 1 1 . . 1 1 . . 1 1 . . 1 1 1 1   = H 3 . ( 13
)
Property 3.4 implies that the right product with P πc is equivalent to the permutation of the columns according to π -1 c . In the example given here, since in [START_REF] Ranganathan | On the girth of (3, L) Quasi-Cyclic LDPC Codes based on Complete Protographs[END_REF],

π c = {4, 3, 6, 5, 2, 1}, we have π -1 c = {6, 5, 2, 1, 4, 3}.
It is possible to interpret directly the effect of left multiplication of a matrix A by P πc : since π c (6) = 1, the last column of A would be the first column of A × P πc , and since π c (5) = 2, the fifth column of A would be second column of A × P πc and so on.

Let A ⊗ B represent the Kronecker product of A and B. Let P πr = P πr ⊗ I N 0 , P πc = P πc ⊗ I N 0 and H ∈ H N 3 , then H πr,πc = P πr × H × P πc is also an element of H N 3 satisfying H πr,πc ≡ H. Thus, structural transformations preserve the equivalence relation. In other words, several equivalence classes can be merged owing to the structural transformation, hence reducing the space search. Going back to the example, let H = b, d, e, f be an element of H 3 and (π r , π c ) the pair of permutations defined above, then H πr,πc can be represented as

H πr,πc =   0 0 b 0 . . f e . . d 0 . . 0 0 0 0   . (14)
Then, using [START_REF] Divsalar | Capacityapproaching Protograph Codes[END_REF], H πr,πc (and thus

H) is equivalent to e - f, b, b -d + f, b + f .
For each of the 6 permutations π r on the 3 lines of the prototype matrix H 3 , there are exactly 8 permutations π c on columns that give structural transformations. For example, if π r is equal to the identity, then permuting columns 1 and 2, columns 3 and 4 or columns 5 and 6 does not change the structure of H 3 [START_REF]IEEE standard committee[END_REF]. By combining these 3 permutations, a total of 8 structural permutations are obtained when π r is equal to the identity. The first column of Tables VI and VII gives an exhaustive enumeration of the 6 × 8 = 48 structural transformations. If b, d, f -e and N/2 are all distinct and non null, then each of the 48 permutations gives a distinct matrix after the transformation (see Tables VI andVII). If this condition is not fulfilled, then the number of distinct solutions is lower. For example, if f -e = 0, i.e., e = f , then permutations 1 and 2 of Table VI give the same matrix +b, +d, +e = +f, +f = +e . More generally, if e = f , then for any i ∈ {1, 2, . . . 24}, permutations presented in lines 2i -1 and 2i also give the same result: the number of distinct solutions is thus halved.

Tables VI and VII show the 48 equivalent 4-tuples of 1, 13, 15, 23 35 (subscript 35 indicates a factor of expansion N = 35) as well as the formal expression of the transformed 4-tuple.

IV. STRUCTURAL AND ADDITIVE DOMINANCE BREAKING

In the previous section, we show that, when applied to any H, the 48 permutation pairs of P 3,6 give a matrix H πr,πc in the same equivalence class. In this section, we present dominance constraints and use them to find one solution per class.

The constraints are based on the cycles in the protograph matrix. As shown in section II-B, a cycle can be defined by a sequence of consecutive moves (in a column, then in a row, then again in a column and so until going back to the starting point) between the non-null positions of the prototype matrix. In the example, the degree of the variable node is d v = 2, thus, there is no need to specify the row index to define cycles in H 3 . A column in H 3 is denoted by a letter, "A" being the first column and "F" the last column. By convention we denote a cycle by a sequence of columns considering that the first move is going down in the first column. The moves on subsequent columns (i.e., either down or up) are uniquely determined and do not need to be specified since the column weight is equal to 2.

Example C AB denotes the cycle going down on column A, then up in column B. Cycle C AB is thus going through nodes (a 1 → a 2 → b 2 → b 1 ) in Figure 1(a). Cycle C AEC denotes the cycle going down on column A then down on column E and finally up in column C (see Figure 1(b)). It is important to note that, for example, C ABAB (respectively C ABABAB ) denotes the cycle of length 8 (respectively 12) that makes 2 (respectively 3) round trips between columns A and B. According to lemma 2.1, the enumeration of all cycles of length that uses one or more columns more than once is also required.

  a 1 b 1 c 1 d 1 .
.

a 2 b 2 . . e 1 f 1 . . c 3 d 3 e 3 f 3   (a) Cycles C AB , C ABAB , . . .   a 1 b 1 c 1 d 1 . . a 2 b 2 . . e 2 f 2 . . c 3 d 3 e 3 f 3   (b) Cycles C AEC , C AECAEC , . . .
Φ(H) = φ(C AB ), φ(C CD ), φ(C EF ) N = b, d, f -e N ,
where the φ(C) is the value of the cycle, as defined in [START_REF] Boutillon | Iterative Decoding of Concatenated Convolutional Codes: Implementation Issues[END_REF].

A. Φ-order constraint: φ(C AB ) ≤ φ(C CD ) ≤ φ(C EF ).
Any matrix H = b, d, e, f of H 3 characterized by Φ(H) = b, d, f -e is equivalent to a matrix H 2 ∈ H 3 characterized by Φ(H 2 ) = b, f -e, d .
Proof: Same proof as Lemma 4.1 using the equivalence given in line 13 of Table VI.

Theorem 4.3 (Φ-Order): Any matrix H of H 3 is equivalent to a matrix H 0 = b, d, e, f , with Φ(H 0 ) = b, d, f -e satisfying b ≤ d ≤ f -e.
Proof: Lemmas 4.1 and 4.2 show that it is possible to generate equivalent solutions by swapping, respectively, the first and second terms of Φ(H), then the second and third terms of Φ(H). By combining these two swapping permutations, any order of Φ(H) can be obtained, in particular, the one where the values are in increasing order.

The Φ-order constraint defines a dominance relation over the solutions. Thus, one can restrict the search to solutions respecting this constraint, which reduces the size of the search space by almost a factor of 6.

C. Φ-order-N/2 constraint: φ(C AB ) ≤ φ(C CD ) ≤ φ(C EF ) ≤ N/2. Theorem 4.8 (Constraint Φ-order-N/2): Any matrix H of H 3 is equivalent to a matrix H 0 = b, d, e, f , with Φ(H 0 ) = b, d, f -e satisfying b ≤ d ≤ f -e ≤ N/2.
Proof: Using Theorem 4.7, H is equivalent to a matrix H 0 satisfying the constraint N/2. Then, the proof is completed by applying Theorem 4.3 on H 0 to sort the φ values of H 0 in increasing order.

Constraint φ(C AB ) ≤ φ(C CD ) ≤ φ(C EF ) ≤ N/2 defines a dominance relation, and the search of solutions can then be restricted to those respecting the constraint. This relation allows us to reduce by up to 48 the size of the space to explore. For example, in tables VI and VII, the only equivalent solution H 0 of H = 1, 7, 13, 17 [START_REF] Wang | Construction of high-girth QC-LDPC codes[END_REF] satisfying the Φ-order-N/2 constraint is given on line 13 (in a grey cell in the fifth column) by H 0 = 1, 4, 21, 28 [START_REF] Wang | Construction of high-girth QC-LDPC codes[END_REF] . The Φ value of H 0 is given by Φ(H 0 ) = 1, 4, 7 35 that satisfies

1 ≤ 4 ≤ 7 ≤ 35/2. D. Φ 0 -order-N/2 constraint: 0 < φ(C AB ) ≤ φ(C CD ) ≤ φ(C EF ) ≤ N/2.
Finally, we can also add an additional constraint to the Φorder-N/2 constraint by imposing 0 < b. In fact, b = 0 leads to weak expanded matrices with a girth of size 4 as explained in section VI-A. In the following, this new constraint is called the Φ 0 -order-N/2 constraint. We also define H N,0 is the subset of H N 3 with expanded matrices of girth g ≥ 6.

V. MULTIPLICATIVE TRANSFORMATION

In this section, we present a new transformation which relies on the multiplication of the values of the circulant matrices of H. The properties of this transformation are presented in Section V-A and used in Section V-B to identify new forms of equivalent matrices, thus reducing further the search space.

A. Properties Definition 5.1 (Product-matrix): Let k ∈ Z/N Z be an integer coprime with N (the greatest common divisor of k and N is equal to 1, i.e., gcd(k, N ) = 1), then P × k is the permutation matrix associated with the permutation

π × k (i) = k × i mod N .
Since k and N are coprime, then there exists an integer k -1 so that k × k -1 = 1 mod N and thus, (π

× k ) -1 = π × k -1 and (P × k ) -1 = P × k -1 . Theorem 5.2 (Multiplication over a circulant matrix): Let k ∈ Z/N Z coprime with N and I N a the a-circulant matrix, then P × k × I N a × P × k -1 is equal to the (ka mod N )-circulant matrix I N ak . Proof: P × k × I N a × P × k -1 is the permutation matrix associated with the permutation π = π × k (π + a (π × k -1 )). For all i ∈ Z/N Z, we have π × k -1 (i) = k -1 i, thus π + a (π × k -1 (i)) = k -1 i + a and therefor, π(i) = π × k (k -1 i + a) = k(k -1 i + a) = i + ka.
Thus, π is the ka-circular permutation, and therefore, its associated permutation matrix is I ka . Example 5.3: Consider N = 32, k = 7 and a = 3. Since gcd(7, 32) = 1, then 7 -1 exists over Z/32Z. Since 7 × 23 = 161 = 1 + 5 × 32 = 1 mod 32, then 7 -1 = 23 over Z/32Z. Thus, according to Theorem 5.2, P × 7 × I 32 3 × P × 23 = I 32 21 . Another way to reach this result has been presented in reference [START_REF] Tasdighi | Symmetrical Constructions for Regular Girth-8 QC-LDPC Codes[END_REF]. It can be extended to the whole prototype matrix to identify a new type of transformation (i.e., multiplicative transformation) that also preserves the equivalence relation.

Theorem 5.4 (Multiplicative rearrangement):

H × k = (I J 0 ⊗ P × k ) × H × (I L 0 ⊗ P × k -1
) is the expanded matrix defined as

H × k (i, j) = k × H(i, j).
Proof: The proof of this theorem comes directly from Theorem 5.2 since each (N, N ) expanded matrices of

H × k (see (4)) is equal either to P × k × 0 N × P × k -1 = 0 N or P × k × I N a × P × k -1 = I N ka .
Combined with additive and structural transformations, multiplicative transformation allows us to identify more matrices in the same equivalent class, thus reducing the search space since only one representative of each equivalence class needs to be tested.

B. Symmetry breaking

In this section, we use Theorem 5.4 to identify new subsets of equivalent matrices in H 3 . Note that this method can be applied to any prototype matrix. Proof: From Bezout's theorem, there exist two integers u, v so that ub + vN = gcd(b, N ), where gcd is the "greatest common divisor". For any integer ρ, let us define k ρ as k ρ = u + ρ × N/ gcd(b, N ). Thus, by construction, k ρ × b = gcd(b, N ) mod N . According to Diriclet's prime number theorem, since u and N/ gcd(b, N ) are coprime, the set {k ρ } ρ∈N contains an infinity of prime numbers, thus, there exists a k ρ coprime with N which complete the proof.

For example, let N = 30 and H = 12, 13, 4, 17 [START_REF] Marchand | Non-Binary WEB page of Lab-STICC[END_REF] . Since b = 12, we have gcd(b, N ) = 6 and N/ gcd(b, N ) = 5. Taking u = 3 and v = -1 gives 3 × 12 + (-1) × 30 = 6. Thus, k ρ = 3 + 5ρ takes respectively the values {3, 8, 13, 18, 23, 28} for ρ = 0, 1, . . . , 5. In this series, k 2 = 13 and k 4 = 23 are the only integers coprime with N = 30. Thus, using k 2 , H = 12, 13, 4, 17 30 is equivalent to 6, 19, 22, 11 30 , and using k 4 , H is also equivalent to 6, 29, 2, 1 30 .

Theorem 5.7 (Global Constraint): To conclude this section, we use multiplicative permutation to show that the number of equivalent classes (or equivalently, the space of search, since only one element per class needs to be tested) is upper bounded by ρ(N ), where ρ(N ) is the number of distinct tuples b, d, e, f respecting the Φ 0 -ordre-N/2 constraint and b ∈ N . The value ρ(N ) can be computed as

Any matrix H of H 0 3 is equivalent to a matrix H n = b n , d n , e n , f n , with Φ(H n ) = b n , d n , f n -e n verifying 0 < b n ≤ d n ≤ f n -e n ≤ N
ρ(N ) = b∈D(N )/N N/2 d=b d≤f -e≤N/2 1. ( 16
)
Let us focus on the last of the three sum operators of [START_REF] O'sullivan | Algebraic construction of sparse matrices with large girth[END_REF].

ρ(N ) = N b∈D(N )/N ( N/2 -b + 1)( N/2 -b + 2) 2 . ( 17
)
Note that if N is a prime number greater than 2, then D(N )/N = {1} and thus, ρ(N ) = (N 3 -N )/8. For example, for N = 5, the ρ(5

) = 15 solutions are 1, 1, 0, 1 , 1, 1, 0, 2 , 1, 1, 1, 2 , 1, 1, 1, 3 , 1, 1, 2, 3 , 1, 1, 2, 4 , 1, 1, 3, 4 , 1, 1, 3, 0 , 1, 1, 4, 0 , 1, 1, 4, 1 , 1, 2, 0, 2 , 1, 2, 1, 3 , 1, 2, 2, 4
, 1, 2, 3, 0 and 1, 2, 4, 0 .

C. Discussion

The links between structural additive equivalence and multiplicative transformation are still an open problem. In Tables VI and VII, the last column shows three equivalent matrices generated from the multiplicative transformation of 1, 7, 13, 17 [START_REF] Wang | Construction of high-girth QC-LDPC codes[END_REF] satisfying the Φ 0 -order-N/2 constraint, meaning that the three matrices (in grey in the last column) belong to the same equivalence class. So far, we have been unable to find an a priori method to determine a single element of this class. Nevertheless, once the matrices of interest are generated for a given girth, it is possible a posteriori to prune the space of solutions by seeking equivalent solutions using multiplicative transformation and keeping only one element per equivalent class. The present authors consider that, while the problem may be of interest to group theorist specialists, an exhaustive analysis of the group properties is beyond the scope of their paper. Finally, it is worth mentioning that the method used here can be generalized for any prototype matrix to reduce the search space. The difficulty is first to find the structural set P H of pair of permutations associated with the prototype matrix H, and then find a simple constraint to determine a single matrix per equivalence class. This is a challenging problem to solve in the general case! Finally, readers interested in the mathematical aspects of the present study are invited to consult the on-line note written of Xavier Giraud based on an early version of our paper [START_REF] Giraud | Comments on Additive, Structural and Multiplicative transformations for the construction of Quasi Cyclic LDPC matrices[END_REF]. M. Giraud also contributed to the paper by pointing out, in the global theorem, that b can be also be taken smaller or equal than gcd(d, N ) and than gcd(f -e, N ).

VI. CONSTRUCTION OF HIGH GIRTH MATRICES

In this section, we construct H matrices of minimum size from the H 3 prototype matrix with a girth ranging from g = 8 up to g = 14. Then, we give an explicit rule based on the multiplicity M(g) to select the matrix. Finally, to illustrate the efficiency of the proposed method, we present a very early promising result for another type of prototype matrix.

A. Optimal lifting of the H 3 protograph.

The problem of optimal lifting by a factor N of the H 3 protograph matrix can be formalized in two steps. The first step is to define the maximum achievable girth ḡN as ḡN (H) = max

H∈H N 3 g(H), (18)
where g(H) is the girth of matrix H. From ḡN (H), we can define H N,o 3 as the subset of H N 3 of matrices of maximum girth ḡN . Then the optimal lifted matrices H o is given by

H o = arg min H∈H N,o 3 M(ḡ N (H)), (19)
where M(ḡ N (H)) is the number of cycles of length ḡN (H) of the matrix H. Note that if several distinct matrices lead to the same M(ḡ N (H)) minimum value, then the value of M(ḡ N (H) + 2) can be used as a second criteria and so on.

In the following, we present some results on the inverse problem, i.e., to find the minimum value of N required to guaranty a given girth g 0 . To do so, we first express all the cycles of length = 4, 6, 8, 10, 12 and 14 of H 3 using the method presented in [START_REF] Asvadi | Design of Finite-Length Irregular Protograph Codes with Low Error Floors over the Binary-Input AWGN Channel Using Cyclic Liftings[END_REF]. Note that the software to perform this enumeration is available online [START_REF] Marchand | Non-Binary WEB page of Lab-STICC[END_REF]. The number N ( ) of cycles for each value of is given in Table I. Each cycle gives a new constraint on the lifted matrix according to equation ( 6). .

TABLE I NUMBER OF CONSTRAINTS TO BE TAKEN INTO ACCOUNT BY THE CONSTRAINT PROGRAMMING SOFTWARE

Then, using the Φ 0 -order-N/2 constraint, we obtained a well-defined constraint optimization problem to find the matrix b, d, e, f N giving a girth greater of equal than g 0 : according to lemma 2.1, any cycle C of the protograph matrix of length < g 0 should satisfy Φ(C) = 0. The total number of constraints Nl < g 0 is given in the last line of Table I.

It is noteworthy that some constraints are redundant. For example, the cycle C AECBEF B of length = 14 gives the same constraint as the cycle

C AEF EC of length = 10 since Φ(C AECBEF B ) = Φ(C AEF EC ) = f -2e. Nevertheless, introducing redundant constraints in a constraint programming tool is not a problem.
Table II presents the number of solutions for a targeted girth g 0 = 8 up to g 0 = 16, with respect to the expansion size N . In the first column, we give the highest value of N for which no solution can be found for the corresponding girth. We also give the number of solutions for the next 10 values of N to show the evolution of the number of solutions with increasing N . Line 2 gives the number of solutions n 0 when only the additive transformation is used to reduce the search space (matrices of type H = b, d, e, f ) with no particular constraint on b, d, e and f ), line 3 gives the number of solution n 1 when the Φ-order -N/2 is used, while line 4 shows the number of solutions n 2 when both Φ-order-N/2 constraint and multiplicative transformation are applied. The last line shows the final fraction of remaining search space R = n 2 /n 0 . For a given expansion factor N , if a solution with a girth g exist, this solution will be enumerated in the list of solutions obtained for a targeted girth g 0 = g -2. For example, in Table II, among the n 2 = 14 solutions obtained for a targeted girth g 0 = 8 with an expansion factor of N = 8, two have a girth equal to 10. When g 0 ≥ 10, any two cycles of length four have distinct values. The Φ 0 -order-N/2 constraint can thus be expressed with strict inequality as 0 < Φ(C AB ) < Φ(C CD ) < Φ(C EF ) < N/2 . Then, each of the 48 additive and structural transformations gives a distinct solution. In that case, n 1 is just equal to n 1 = n 0 /48. Moreover, due to the multiplicative transformation (Theorem 5.5), some of the leftover solutions can be proved equivalent, thus the reduction can be significantly greater than a factor 48. For example, in Table II when N equals 23, there are 23760 solutions when no equivalence is used. This number reduces to n 1 = 23760/48 = 495 when using Φ 0 -order-N/2 equivalence. Among those 495 solutions, 470 are proved to be equivalent due to multiplicative transformation. Thus, only a maximum of n 2 = 25 out of 23760 are distinct solutions (all solutions have a distinct multiplicity spectrum), which correspond to only R = 0.19% of the initial space. Note that the computation time is also significantly reduced: 3.9 seconds to determine n 2 and its related solutions as against 97.6 s to determine n 0 and its related solutions for N = 45 and g = 16.

B. Elements of differentiation of the solutions

The above section presents the formulation of constraints forbidding cycles of a given size and shows the result of numerical experiments. From those experiments, we can highlight the reduction of the number of solutions found by using dominance breaking. Except for some rare cases, there exist several solutions for a given girth and expansion size and our method finds them all.

Due to the significant reduction in the number of solutions, it is now possible to analyze in more details the properties of each solution by computing the multiplicity M(g) of cycle of length g, and eventually, the multiplicity of greater length. Table III shows an example of the 5 solutions obtained for g = 14 and N = 30, giving the matrices and the values of M( 14) and M [START_REF] O'sullivan | Algebraic construction of sparse matrices with large girth[END_REF].

One can note that, in this example, all solutions are distinct since they have different cycle multiplicities and that solution 4 gives the code with the best topological properties.

One of the most surprising results on Table II is that the number n 2 of distinct solutions does not always increase when N increases. For g ≥ 8, the number n 2 of unique solutions for N = 5, 6, . . . 13 is respectively n 2 = 1, 7, 4, 17, 14, 37, 20, 99 and 35. Even more surprising, for g = 16, there exist n 2 = 4 solutions when N = 36, none when N = 37 and N = 38, then 2 for N = 39, 6 for N = 40 and one for N = 41! In other words, the girth is not necessarily an increasing function of the expansion factor of the lifting factor N , since g = 16 for N = 36, g = 14 for N = 37 and N = 38, g equals 16 again for N ≥ 39. This result is rather counter-intuitive so it has been carefully checked. It is worth mentioning that the method proposed here has been used to construct several good Non-Binary-LDPC matrices of several sizes and code rates available on line [START_REF] Marchand | Non-Binary WEB page of Lab-STICC[END_REF]).

Finally, Table VIII in the Appendix gives the best expanded matrices found for values of N between 4 and 43. For each value of N , a NB-LDPC codes over GF(64) is constructed. The choice of the Galois Field coefficients is carried out following the rules given in [START_REF] Poulliat | Design of regular (2,dc)-LDPC codes over GF(q) using their binary images[END_REF]. Figure 2 compares the energy per symbol versus the Energy of the noise (E s /N 0 , in dB) required to obtain a Frame Error Rate (FER) of 10 -2 and 10 -3 , respectively. The decoding algorithm is the Extended Min Sum (EMS) algorithm of parameter n m = 20 [START_REF] Voicila | Lowcomplexity decoding for non-binary LDPC codes in high order fields[END_REF] with 10 decoding iterations. The modulation is the Binary Phase Shift Keying (BPSK) modulation over the Additive White Gaussian Noise (AWGN) channel. The channel outputs are quantified on 5 bits before entering the decoder.

C. Extension to the complete (3-L) QC-LDPC matrix

Before presenting the conclusion of this paper, we give an additional example to highlight the efficiency of the proposed method to construct QC-LDPC matrices. Let us consider the set of complete prototype matrices H (3,L) of size 3 × L, (H (3,L) is composed of 3 lines of L ones)). The problem of finding the minimum lifting factor N giving a QC-LDPC matrix of girth 8 and 10 is well studied in the literature. Table IV shows the evolution of the minimum value of N cited in the literature for L varying from 4 to 12 and girth g equal to 8.

These solutions are found by a heuristic search with no proof of optimality. Using the theory and the tools developed in the present study for this problem, we were rapidly able (in one day) to:

1) Enumerate all the cycles of length 4 and 6 of the H

Note that we do not have any proof that N = 40 is optimal (i.e. minimum expansion factor) for L = 11. Without a deeper study of this problem involving symmetry breaking, the running time required to demonstrate optimality is too long. Property 5.6 and Theorem 4.3 would help in finding new symmetries to reduce the size of the space to explore. We leave this as an open question.

For the girth g = 10, the minimum lifting factor N is greatly reduced compared to the state of the art. For example, for L =

VII. CONCLUSION

In this paper, we have defined the notion of equivalence classes between QC-LDPC matrices. We have proposed three transformations that preserve the equivalence between matrices: the additive transformation (already known), the structural transformation and the multiplicative transformation. As an example, we have applied these transformations to a particular 3-row, 6-column protomatrix with variable node degree of 2 and check node degree 4. For this problem, we proposed criteria to select a single element for each equivalence class of matrices, thus allowing a very fast exhaustive exploration. In fact, since an equivalent class can contain from a few tens up to a few thousand elements, the search space is reduced accordingly. Finally, we make use of a tool from the constraint programming literature to prove and find solutions to the construction problem.

We should note that the new equivalence relations between QC-LDPC matrices are very general and can be used for any construction of QC-LDPC matrices from the lifting of a prototype matrix. Finally, we conclude with an open question: is it possible to use additive, structural and multiplicative transformations on their own to generate the equivalent class of any QC matrix?

Fig. 1 .

 1 Fig. 1. Graphical description of cyclesAs previously mentioned, any matrix in H 3 is equivalent to a matrix H = b, d, e, f . It can thus be characterized by the value of the 3 cycles of length four in matrix H 3 , i.e. C AB , C CD and C EF , i.e., by the tupleΦ(H) = φ(C AB ), φ(C CD ), φ(C EF ) N = b, d, f -e N ,where the φ(C) is the value of the cycle, as defined in[START_REF] Boutillon | Iterative Decoding of Concatenated Convolutional Codes: Implementation Issues[END_REF].

Lemma 4 . 1 (

 41 Swapping C AB and C CD ): Any matrix H = b, d, e, f of H 3 characterized by Φ(H) = b, d, f -e is equivalent to a matrix H 1 ∈ H 3 characterized by Φ(H 1 ) = d, b, f -e . Proof: According to line 26 of Table VII, H = b, d, e, f is equivalent to the matrix H 1 = d, b, -f, -e owing to the equivalence relation H 1 = P πr × H × P πc , with (π r , π -1 c ) ∈ P 3,6 defined as π r = {1, 3, 2} and π c = {3, 4, 1, 2, 5, 6}. Matrix H 1 is characterized by Φ(H 1 ) = d, b, f -e , which completes the proof. Lemma 4.2 (Swapping C CD and C EF ):

B. N/ 2 - 4 . 5 : 4 . 6 :

 24546 constraints: φ(C AB ) ≤ N/2, φ(C CD ) ≤ N/2 and φ(C EF ) ≤ N/2. Lemma 4.4: Any matrix H = b, d, e, f of H 3 characterized by Φ(H) = b, d, f-e is equivalent to a matrix H 3 ∈ H 3 characterized by Φ(H 3 ) = -b, d, f -e .Proof: Similarly to the proof of Lemma 4.1, this involves the transformation of line 5 of TableVI.Lemma Any matrix H = b, d, e, f of H 3 characterized by Φ(H) = b, d, f-e is equivalent to a matrix H 4 ∈ H 3 characterized by Φ(H 4 ) = b, -d, f -e .Proof: Similarly to the proof of Lemma 4.1, this involves the transformation of line 3 of TableVI.Lemma Any matrix H = b, d, e, f of H 3 characterized by Φ(H) = b, d, f-e is equivalent to a matrix H 5 ∈ H 3 characterized by Φ(H 5 ) = b, d, -(f -e) .Proof: Similarly to the proof of Lemma 4.1, this involves the transformation of line 2 of TableVI. Theorem 4.7 (constraint N/2): Any matrix H of H 3 is equivalent to a matrix H 0 = b, d, e, f , with Φ(H) = b, d, f -e satisfying b ≤ N/2, d ≤ N/2 and f -e ≤ N/2. Proof: Let us consider H = b, d, e, f where b > N/2, then using Lemma 4.4, H is equivalent to H 0 = -b, d, e, f . Since the operation is performed modulo N , b > N/2 ⇒ (-b mod N ) ≤ N/2. The proof is completed by Using lemma 4.5 and 4.6 in a similar way.

3 (or H 0 3 ,

 33 to simplify the notation) as the subset of all matrices of H N 3 equivalent to a matrix H = b, d, e, f respecting the Φ 0 -order-N/2 constraint, i.e. 0 < b ≤ d ≤ f -e ≤ N/2. Hence, H N,0 3

Property 5 . 5 (

 55 Coprime multiplication): Let k be a number coprime with N , then any matrix H = b, d, e, f N of H 3 is equivalent to the matrix kb, kd, ke, kf N Derived directly from Theorem 5.4.Once a solution is found, it is possible to build the class of all equivalent solutions for all k coprime with N . Property 5.6 (Divisors of N ): If b = 0, b can be restricted to values in the set D(N )/N of the divisors of N with N excluded.

  /2 and b n a divisor of N verifying b n ≤ gcd(d b , N ) and b n ≤ gcd(f n -e n , N ).Proof: The proof of this Theorem is performed by recursion, by showing the existence of a series of equivalent matricesH 0 = b 0 , d 0 , e 0 , f 0 , H 1 = b 1 , d 1 , e 1 , f 1 , . . ., H n = b n , d n ,e n , f n satisfying the Φ 0 -order-N/2 constraint, with decreasing value of b, i.e. (b 0 > b 1 > b 2 • • • > b n ≥ 1) until b n reaches the value of a divisor of N (eventually, b n = 1). Initial condition: From Theorem 4.8, any matrix H of H 0 3 is equivalent to a matrix H 0 = b 0 , d 0 , e 0 , f 0 satisfying the Φ 0 -order-N/2 constraint. Recursion: Let us assume that H is equivalent to a matrix H n = b n , d n , e n , f n satisfying the Φ 0 -order-N/2 constraint with b n ≥ 1. If b n is a divisor of N , the condition of the theorem is fulfilled. Otherwise, Theorem 5.4 shows that there exists k such that H n is equivalent to gcd(b n , N ), kd n , ke n , kf n . Note that gcd(b n , N ) < b n since b n is not a divisor of N . Due to the multiplication by k and the reduction modulo N , gcd(b n , N ), kd n , ke n , kf n may no longer respect the Φ 0 -order-N/2 constraint. Nevertheless, it is still possible to find an equivalent matrix H n+1 = b n+1 , d n+1 , e n+1 , f n+1 that respects the constraint. Since gcd(b n , N ) ≤ N/2, it will be only affected by swapping operations during the construction of an equivalent matrix respecting the Φ 0 -order-N/2 constraint (see proof of Theorem 4.8), and thus, 1 ≤ b n+1 ≤ gcd(b n , N ) < b n . Finally, if b n > gcd(d n , N ), then d n and b n can be swapped and the same process will lead to b n+1 = gcd(d n , N ) < b n . The same method can be applied if b n > gcd(f n -e n , N )

  For each value of e, there are N/2 -d + 1 possible values of f that fulfilled the constraint d ≤ (f -e) ≤ N/2, so the last summation term is equal to N ( N/2 -d + 1). According to the second summation term of (16), d varies from b to N/2 , thus ( N/2 -d + 1) takes all the values from 1 up to N/2 -b + 1, thus

2 )

 2 Prove the optimality of the existing solutions for L = 4 up to L = 10 1 . 3) Find a new solution for L = 11 with an expansion factor of N = 40. This new matrix H 40 (3,11) generated with the constraint programming tool is

COMPARISON

  OF EXECUTION TIME AND NUMBER OF SOLUTIONS WITH AND WITHOUT Φ-ORDER-N/2 CONSTRAINT FOR g 0 = 8, 10, 12, 14 AND 16.

				N	< 4	4		5	6	7	8	9	10	11	12	13
		g 0 = 8	n 0 n 1	0 0	2 2		16 2	92 14	288 18	702 55	1440 67	2648 150	4480 175	7130 335	10800 378
				n 2	0	1		1	7	4	17	14	37	20	99	35
				R	N/A	50 %	6.25 %	7.61 %	1.39 %	2.42 %	0.97 %	1.40 %	0.45 %	1.39 %	0.32%
				N	< 8	8		9	10	11	12	13	14	15	16	17
		g 0 = 10	n 0 n 1	0 0	48 1	336 7	480 10	1680 35	2160 45	5184 108	6096 127	12336 257	13968 291	24960 520
				n 2	0	1		2	3	4	15	10	24	40	45	35
				R	N/A	2.08 %	0.59 %	0.62 %	0.24 %	0.69 %	0.19 %	0.39 %	0.32 %	0.32 %	0.14 %
				N	< 14	14		15	16	17	18	19	20	21	22	23
		g 0 = 12	n 0 n 1	0 0	480 10	432 9	1056 22	1920 40	4032 84	5184 108	7296 152	12960 270	17520 365	23760 495
				n 2	0	3		2	4	3	18	7	25	30	41	25
				R	N/A	0.62 %	0.46 %	0.38 %	0.16 %	0.45 %	0.13 %	0.34 %	0.23 %	0.23 %	0.11 %
				N	< 28	28		29	30	31	32	33	34	35	36	37
		g 0 = 14	n 0 n 1	0 0	1440 30	1344 28	1632 34	5760 120	6528 136	12480 260	13824 288	30336 632	30048 626	60480 1260
				n 2	0	3		1	5	4	12	14	18	28	57	35
				R	N/A	0.21 %	0.07 %	0.30 %	0.069 %	0.18 %	0.11 %	0.13 %	0.09 %	0.19 %	0.058 %
				N	< 36	36		37	38	39	40	41	42	43	44	45
		g 0 = 16	n 0 n 1	0 0	1536 32		0 0	0 0	1728 36	2784 58	1920 40	5088 106	4032 84	15840 330	12096 252
				n 2	0	4		0	0	2	6	1	10	2	18	13
				R	N/A	0.26 %	N/A	N/A	0.12 %	0.21 %	0.052 %	0.20 %	0.050 %	0.11 %	0.11 %
										TABLE II
			#	b, d, e, f 30	M(14)	M(16)
			1	1, 3, 7, 19 30	180		825
			2 3	1, 5, 7, 16 30 1, 6, 25, 8 30	180 180		840 810	.
			4	1, 9, 2, 13 30	90		1080
			5	3, 5, 4, 13 30	240		750
					TABLE III
			THE FIVE SOLUTIONS FOR g = 14, N = 30
		3.4
		3.2							FER = 10 -3 FER = 10 -2
		3
		2.8
	0 (dB)	2.6
	E S /N	2.4
		2.2
		2
		1.8
		1.6
		0	100	200	300	400	500	600	700	800
					Information Blocksize
	Fig. 2. Required Es/N 0 to obtain a FER of 10 -2 and 10 -3 for NB-LDPC
	matrices over GF(64) constructed from matrices of Table VIII. The EMS
	decoding algorithm is used with 10 decoding iterations. Received samples are
	quantified on 5 bits.

RANDOMLY

  FOUND SOLUTION FOR A COMPLETE PROTOMATRIX, WITH J = 3 LINES AND L COLUMNS REQUIRED TO OBTAINED A GIRTH g = 10 FOR A COMPLETE (3,L) QC-LDPC MATRIX 11, N is reduced between 2008 and 2016 from N = 577 down to N = 560. In the present study, we propose a new value of N equal to 506! Smallest values of N are also obtained for L = 7, 8, 9 and 10 as shown in Table V. The corresponding matrices are

		L = 4	L = 5	L = 6	L = 7	L = 8	L = 9	L = 10	L = 11	L = 12
	2004	[20]	9	14	18	21	26	33	39	46	54
	2006	[16]	9	13	18	22	27	34	40	49	55
	2013	[15]	9	13	18	21	25	30	35	41	47
	2015	[33]	9	13	18	21	25	30	35	40	46
	2017	[19]	9	13	18	21	25	30	35	41	45
						TABLE IV
				L = 4	L = 5	L = 6	L = 7	L = 8	L = 9	L = 10	L = 11
	2013	Lower bound [34]	37	61	91	127	168	217	271	331
	2008	[35]		39	63	103	160	233	329	439	577
	2012	[36]		37	61	101	159	219	319	439	560	,
	2016	[33];[37]	37	61	91	155	227	323	429	571
	2017	[21]		37	61	91	-	-	-	-	-
	2018	Proposed	37	61	91	139	201	280	383	503
						TABLE V
	SMALLEST LIFTING FACTOR N

1

Taking the constraints H N (3,L) (2, 2) as a divisor of N thanks to multiplicative transformation and

H N (3,L) (2, j) ≤ H 40 (3,L)

(2, j + 1), j = 1, . . . , L -1 thanks to column permutations of the matrix.

  + f, -b, +e, -d + e N 4, 34, 13, 6 35 1, 26, 12, 19 35 18 {2, 3, 1} {3, 4, 5, 6, 2, 1} +e -f, -b, +f, -d + f N 31, 34, 17, 10 35 1, 9, 22, 15 35 19 {2, 3, 1} {3, 4, 6, 5, 1, 2} -e + f, -b, -d + e, +e N 4, 34, 6, 13 35 1, 26, 19, 12 35 20 {2, 3, 1} {3, 4, 6, 5, 2, 1} +e -f, -b, -d + f, +f N 31, 34, 10, 17 35 1, 9, 15, 22 35 21 {2, 3, 1} {4, 3, 5, 6, 1, 2} -e + f, +b, +b + e, +b -d + e N 4, 1, 14, 7 35 1, 9, 21, 28 35 22 {2, 3, 1} {4, 3, 5, 6, 2, 1} +e -f, +b, +b + f, +b -d + f N 31, 1, 18, 11 35 1, 26, 13, 6 35 23 {2, 3, 1} {4, 3, 6, 5, 1, 2} -e + f, +b, +b -d + e, +b + e N 4, 1, 7, 14 35 1, 9, 28, 21 35 24 {2, 3, 1} {4, 3, 6, 5, 2, 1} +e -f, +b, +b -d + f, +b + f N 31, 1, 11, 18 35 1, 26, 6, 13 35

	#	πr	πc	b, d, e, f N ≡	1, 7, 13, 17 N	b divisor of N
	1	{1, 2, 3}	{1, 2, 3, 4, 5, 6}	+b, +d, +e, +f N	1, 7, 13, 17 35	1, 7, 13, 17 35
	2	{1, 2, 3}	{1, 2, 3, 4, 6, 5}	+b, +d, +f, +e N	1, 7, 17, 13 35	1, 7, 17, 13 35
	3	{1, 2, 3}	{1, 2, 4, 3, 5, 6}	+b, -d, -d + e, -d + f N	1, 28, 6, 10 35	1, 28, 6, 10 35
	4	{1, 2, 3}	{1, 2, 4, 3, 6, 5}	+b, -d, -d + f, -d + e N	1, 28, 10, 6 35	1, 28, 10, 6 35
	5	{1, 2, 3}	{2, 1, 3, 4, 5, 6}	-b, +d, +b + e, +b + f N	34, 7, 14, 18 35	1, 28, 21, 17 35
	6	{1, 2, 3}	{2, 1, 3, 4, 6, 5}	-b, +d, +b + f, +b + e N	34, 7, 18, 14 35	1, 28, 17, 21 35
	7	{1, 2, 3}	{2, 1, 4, 3, 5, 6}	-b, -d, +b -d + e, +b -d + f N	34, 28, 7, 11 35	1, 7, 28, 24 35
	8	{1, 2, 3}	{2, 1, 4, 3, 6, 5}	-b, -d, +b -d + f, +b -d + e N	34, 28, 11, 7 35	1, 7, 24, 28 35
	9	{2, 1, 3}	{1, 2, 5, 6, 3, 4}	-b, -e + f, -e, +d -e N	34, 4, 22, 29 35	1, 31, 13, 6 35
	10	{2, 1, 3}	{1, 2, 5, 6, 4, 3}	-b, +e -f, -f, +d -f N	34, 31, 18, 25 35	1, 4, 17, 10 35
	11	{2, 1, 3}	{1, 2, 6, 5, 3, 4}	-b, -e + f, +d -e, -e N	34, 4, 29, 22 35	1, 31, 6, 13 35
	12	{2, 1, 3}	{1, 2, 6, 5, 4, 3}	-b, +e -f, +d -f, -f N	34, 31, 25, 18 35	1, 4, 10, 17 35
	13	{2, 1, 3}	{2, 1, 5, 6, 3, 4}	+b, -e + f, -b -e, -b + d -e N	1, 4, 21, 28 35	1, 4, 21, 28 35
	14	{2, 1, 3}	{2, 1, 5, 6, 4, 3}	+b, +e -f, -b -f, -b + d -f N	1, 31, 17, 24 35	1, 31, 17, 24 35
	15	{2, 1, 3}	{2, 1, 6, 5, 3, 4}	+b, -e + f, -b + d -e, -b -e N	1, 4, 28, 21 35	1, 4, 28, 21 35
	16	{2, 1, 3}	{2, 1, 6, 5, 4, 3}	+b, +e -f, -b + d -f, -b -f N	1, 31, 24, 17 35	1, 31, 24, 17 35
	17	{2, 3, 1}	{3, 4, 5, 6, 1, 2}	-e

  TABLE VI EQUIVALENT MATRICES THROUGH STRUCTURAL AND ADDITIVE TRANSFORMATIONS. PART I: LINES 1 TO 24. +b, +d -e, +d -f N 28, 1, 29, 25 35 7, 34, 6, 10 35 28 {1, 3, 2} {3, 4, 2, 1, 6, 5} -d, +b, +d -f, +d -e N 28, 1, 25, 29 35 7, 34, 10, 6 35 29 {1, 3, 2} {4, 3, 1, 2, 5, 6} +d, -b, -b -e, -b -f N 7, 34, 21, 17 35 7, 34, 21, 17 35 30 {1, 3, 2} {4, 3, 1, 2, 6, 5} +d, -b, -b -f, -b -e N 7, 34, 17, 21 35 7, 34, 17, 21 35 31 {1, 3, 2} {4, 3, 2, 1, 5, 6} -d, -b, -b + d -e, -b + d -f N 28, 34, 28, 24 35 7, 1, 7, 11 35 32 {1, 3, 2} {4,3, 2, 1, 6, 5} -d, -b,-b + d -f, -b + d -e N 28,34, 24, 28 35 7, 1, 11, 7 35 +d, -b + d -e, +d -e N 31, 7, 28, 29 35 1, 7, 28, 19 35 48 {3, 2, 1} {6, 5, 4, 3, 2, 1} -e + f, +d, -b + d -f, +d -f N 4, 7, 24, 25 35 1, 28, 6, 15 35

	#	πr	πc	b, d, e, f N ≡	1, 7, 13, 17 N	b divisor of N
	25	{1, 3, 2}	{3, 4, 1, 2, 5, 6}	+d, +b, -e, -f N	7, 1, 22, 18 35	7, 1, 22, 18 35
	26	{1, 3, 2}	{3, 4, 1, 2, 6, 5}	+d, +b, -f, -e N	7, 1, 18, 22 35	7, 1, 18, 22 35
	27 -d, 33 {1, 3, 2} {3, 4, 2, 1, 5, 6} {3, 1, 2} {5, 6, 1, 2, 3, 4} -d, +e -f, +e, +b + e N	28, 31, 13, 14 35	7, 4, 22, 21 35
	34	{3, 1, 2}	{5, 6, 1, 2, 4, 3}	-d, -e + f, +f, +b + f N	28, 4, 17, 18 35	7, 31, 18, 17 35
	35	{3, 1, 2}	{5, 6, 2, 1, 3, 4}	+d, +e -f, -d + e, +b -d + e N	7, 31, 6, 7 35	7, 31, 6, 7 35
	36	{3, 1, 2}	{5, 6, 2, 1, 4, 3}	+d, -e + f, -d + f, +b -d + f N	7, 4, 10, 11 35	7, 4, 10, 11 35
	37	{3, 1, 2}	{6, 5, 1, 2, 3, 4}	-d, +e -f, +b + e, +e N	28, 31, 14, 13 35	7, 4, 21, 22 35
	38	{3, 1, 2}	{6, 5, 1, 2, 4, 3}	-d, -e + f, +b + f, +f N	28, 4, 18, 17 35	7, 31, 17, 18 35
	39	{3, 1, 2}	{6, 5, 2, 1, 3, 4}	+d, +e -f, +b -d + e, -d + e N	7, 31, 7, 6 35	7, 31, 7, 6 35
	40	{3, 1, 2}	{6, 5, 2, 1, 4, 3}	+d, -e + f, +b -d + f, -d + f N	7, 4, 11, 10 35	7, 4, 11, 10 35
	41	{3, 2, 1}	{5, 6, 3, 4, 1, 2}	+e -f, -d, -e, -b -e N	31, 28, 22, 21 35	1, 28, 12, 21 35
	42	{3, 2, 1}	{5, 6, 3, 4, 2, 1}	-e + f, -d, -f, -b -f N	4, 28, 18, 17 35	1, 7, 22, 13 35
	43	{3, 2, 1}	{5, 6, 4, 3, 1, 2}	+e -f, +d, +d -e, -b + d -e N	31, 7, 29, 28 35	1, 7, 19, 28 35
	44	{3, 2, 1}	{5, 6, 4, 3, 2, 1}	-e + f, +d, +d -f, -b + d -f N	4, 7, 25, 24 35	1, 28, 15, 6 35
	45	{3, 2, 1}	{6, 5, 3, 4, 1, 2}	+e -f, -d, -b -e, -e N	31, 28, 21, 22 35	1, 28, 21, 12 35
	46	{3, 2, 1}	{6, 5, 3, 4, 2, 1}	-e + f, -d, -b -f, -f N	4, 28, 17, 18 35	1, 7, 13, 22 35
	47	{3, 2, 1}	{6, 5, 4, 3, 1, 2}	+e -f,

  TABLE VII EQUIVALENT MATRICES THROUGH STRUCTURAL AND ADDITIVE TRANSFORMATIONS. PART II: LINES 25 TO 48).

			0	0	50	125	850	975	3000
	12	1, 2, 6, 15 26	0	0	39	182	676	1274	2704
	12	1, 2, 6, 18 27	0	0	27	135	891	1062	2916
	14	1, 6, 23, 8 28	0	0	0	168	882	1288	2478
	14	1, 4, 6, 15 29	0	0	0	174	812	1421	2349
	14	1, 9, 2, 13 30	0	0	0	90	1080	1050	3195
	14	1, 3, 9, 14 31	0	0	0	124	992	1147	2697
	14	1, 8, 11, 26 32	0	0	0	64	1176	832	3264
	14	1, 10, 2, 14 33	0	0	0	99	1023	1155	2871
	14	1, 3, 8, 23 34	0	0	0	102	918	1326	2601
	14	1, 4, 12, 18 35	0	0	0	70	1015	1085	3150
	16	1, 11, 6, 19 36	0	0	0	0	1296	648	3834
	14	1, 3, 12, 17 37	0	0	0	37	1184	851	3219
	14	1, 3, 10, 15 38	0	0	0	114	988	1140	2565
	16	1, 12, 7, 21 39	0	0	0	0	1209	897	3393
	16	4, 10, 7, 19 40	0	0	0	0	1330	0	5760
	16	1, 3, 14, 19 41	0	0	0	0	1312	574	3567
	16	1, 3, 14, 19 42	0	0	0	0	1260	728	3528
	16	1, 3, 14, 19 43	0	0	0	0	1204	774	3483
					TABLE VIII
		BEST MATRICES FOUNDED FOR SEVERAL EXPANSION FACTORS N .

ACKNOWLEDGMENT

The authors would like to thank several people that have also help for the realization of this paper: Ahmed Abdmouleh for his proposal of the prototype matrix, Marc Sevaux for his reading of a very early version of the paper, Hassan Harb, Cédric Marchand, Titouan Gendron and Franklin Cochachin to have developed some useful tools and performed simulations. A special thanks to Xavier Giraud, Valentin Savin, Alireza Tasdighi and the anonymous reviewers for their reading and their suggestions to improve the paper, Laura Conde-Canencia for her reading of the paper and Michael Carpenter for copyediting of the final version. Part of the work uses ressources funded by the Région Bretagne and the ANR through the CPER Sophie. APPENDIX