def1ch[1];

(* Note: we use normalization such that eta^2=1/2. *)
eta1u = (d[CR, UP] + d[AN, UP])/Sqrt[2];
eta2u = (d[CR, UP] - d[AN, UP])/(I Sqrt[2]);
eta1d = (d[CR, DO] + d[AN, DO])/Sqrt[2];
eta2d = (d[CR, DO] - d[AN, DO])/(I Sqrt[2]);
mu1u = (f[CR, 0, UP] + f[AN, 0, UP])/Sqrt[2];
mu2u = (f[CR, 0, UP] - f[AN, 0, UP])/(I Sqrt[2]);
mu1d = (f[CR, 0, DO] + f[AN, 0, DO])/Sqrt[2];
mu2d = (f[CR, 0, DO] - f[AN, 0, DO])/(I Sqrt[2]);

H1a = eps number[d[]] + U hubbard[d[]] + Bx spinx[d[]] + By spiny[d[]] + Bz spinz[d[]];
H1b = Vu I nc[eta1u, eta2u] + Vd I nc[eta1d, eta2d];

Hc1u = gammaPol1u I nc[eta1u, mu1u];
Hc2u = gammaPol2u I nc[eta2u, mu2u];
Hc1d = gammaPol1d I nc[eta1d, mu1d];
Hc2d = gammaPol2d I nc[eta2d, mu2d];

H = H0 + H1a + H1b + Hc1u + Hc2u + Hc1d + Hc2d;

params = Join[params, {
 gammaPol1u -> Sqrt[(1/Pi) thetaCh[1] (extraGamma1u) gammaA],
 gammaPol2u -> Sqrt[(1/Pi) thetaCh[1] (extraGamma2u) gammaA],
 gammaPol1d -> Sqrt[(1/Pi) thetaCh[1] (extraGamma1d) gammaA],
 gammaPol2d -> Sqrt[(1/Pi) thetaCh[1] (extraGamma2d) gammaA]
}];
