Quit

<< NormalOrder1.m

a = AnnihilationOperator["c"];
b = AnnihilationOperator["d"];
A = CreationOperator["c"];
B = CreationOperator["d"];
vacR = Ket["0"];
vacL = Conjugate[vacR];
vacL ** vacR ^= 1;
AntiCommutator[A[p_], a[q_]] = 
        KroneckerDelta[p, q];
AntiCommutator[B[p_], b[q_]] = 
        KroneckerDelta[p, q];
AntiCommutator[a[_], b[_]] = 0;
AntiCommutator[a[_], B[_]] = 0;
AntiCommutator[A[_], b[_]] = 0;
AntiCommutator[A[_], B[_]] = 0;

aa[p_] = 
        Cos[θ[p]]*a[p] + 
          Sin[θ[p]]*B[-p];
BB[p_] = 
        Sin[θ[p]]*a[-p] - 
          Cos[θ[p]]*B[p];
AA[p_] = 
        Cos[θ[p]]*A[p] + 
          Sin[θ[p]]*b[-p];
bb[p_] = 
        Sin[θ[p]]*A[-p] - 
          Cos[θ[p]]*b[p];
θ[-p_] = θ[p];
Λ[-p_] = Λ[p];

H0[p_] = (AA[p] ** aa[p] + BB[p] ** bb[p])*E[p] + 
          AA[p] ** BB[-p]*Δ[p] + 
          bb[-p] ** aa[p]*Δ[
              p];
H[p1_, p2_, q1_, 
          q2_] = (AA[p1] ** aa[p1] + BB[p1] ** bb[p1])*E[p1] + 
          g*AA[p1] ** BB[q1] ** bb[q2] ** aa[p2]*Λ[
              p1]*Λ[q2]*Λ[p2]*Λ[
              q2]*KroneckerDelta[p1 + q1, 
              p2 + q2];
V[p1_, p2_, q1_, q2_] = H[p1, p2, q1, q2] - H0[p1];

NOSimp[x_, q_] := 
      FullSimplify[
        KronSimplify[
          NormalOrder[
            NormalOrder[
                x] /. (KroneckerDelta[-p_, p_] :> 
                  0)], 
q]]

NOSimp[
      vacL ** V[p1, p2, q1, q2] ** A[s] ** A[t] ** vacR, {}]

NOSimp[
      vacL ** V[p1, p2, q1, q2] ** A[s] ** A[t] ** 
        vacR, {}]
NOSimp[
      vacL ** V[p1, p2, q1, q2] ** B[s] ** B[t] ** 
        vacR, {}]
NOSimp[
      vacL ** V[p1, p2, q1, q2] ** A[s] ** A[t] ** A[u] ** A[v] ** 
        vacR, {}]
NOSimp[
      vacL ** V[p1, p2, q1, q2] ** A[s] ** A[t] ** A[u] ** B[v] ** 
        vacR, {}]
NOSimp[
      vacL ** V[p1, p2, q1, q2] ** A[s] ** B[t] ** B[u] ** B[v] ** 
        vacR, {}]
NOSimp[
      vacL ** V[p1, p2, q1, q2] ** B[s] ** B[t] ** B[u] ** B[v] ** 
        vacR, {}]

stateI1[p_, q_] := A[p] . B[q] . vacR;

NOSimp[
      vacL . V[p1, p2, q1, q2] . stateI1[p, q], {p1, p2, q1, q2}]

stateI2[p_, q_, r_, s_] := 
        A[p] . A[q] . B[r] . B[s] . vacR;

NOSimp[
      vacL . V[p1, p2, q1, q2] . A[s] . B[t] . vacR, {p1, p2, q1, q2, p, 
        t}]

$RecursionLimit = 1000

V2 = 
        FixedPoint[NormalOrder, 
          vacL . V . A[s] . A[t] . B[u] . B[v] . 
            vacR];

V3 = KronSimplify[V2, {p1, p2, q1, q2, p, q}];

KronSimplify[Expand[V3*V3], {s, t, u, v}]

f[p]*g[p]*h[q]*
        d[p, t] //. {Times[f__[p_], d[p_, t_], h___] :> 
          Times[f[t], h]}

Power[x, 2]