Second-order quadrupole shift for static crystal
The second-order quadrupole interaction is related to V(2,±1) and V(2,±2):

In static NMR experiment, V(2,k) is defined by
The second-order quadrupole interaction shifts the energy level |m> by an amount:

In the spectrum, the second-order quadrupole shift of the line position associated with the transition (m - 1,m) is

For the central-transition line with m = 1/2, the second-order quadrupole shift is

with

We also provide Mathematica notebook for calculating the three expressions 2V(2,1)V(2,-1), V(2,2)V(2,-2), and 2V(2,1)V(2,-1) + V(2,2)V(2,-2). Their analytical expressions can be determined as follows:
(1) Select and copy the following green lines; then paste them into a cell of Mathematica, a software for numerical and symbolic calculations.
(2) Press "Ctrl-A" for select all; then
press "Shift-enter" for evaluate cells.
(Or in the menu bar, select Kernel > Evaluation > Evaluate Cells)
Using Mathematica-5 running with a 3-GHz processor, these analytical expressions are obtained in 5 seconds.
(*V2pas is a row-matrix with 3 columns containing the 3 nonzero eigenvalues
of the EFG expressed as a 2-nd rang spherical tensor,in eq unit*)
V2pas={{eta/2,Sqrt[3/2],eta/2}};
(*D21 and D2m1 are reduced forms (3 rows x 1 column) of the 2-
nd order Wigner active rotation matrix*)
D21={{-(1/2)(1+Cos[beta1])Sin[beta1]*E^(-I(2alpha1+gamma1))},
{Sqrt[3/8]Sin[2beta1]*E^(-I gamma1)},
{(1/2)(1-Cos[beta1])Sin[beta1]*E^(I(2alpha1-gamma1))}};
D2m1={{-(1/2)(1-Cos[beta1])Sin[beta1]*E^(I(-2alpha1+gamma1))},
{-Sqrt[3/8]Sin[2beta1]*E^(I gamma1)},
{(1/2)(1+Cos[beta1])Sin[beta1]*E^(I(2alpha1+gamma1))}};
(*V21static and V2m1static are expressions in eq units*)
V21static=FullSimplify[V2pas.ComplexExpand[D21]];
V2m1static=FullSimplify[V2pas.ComplexExpand[D2m1]];
(*V1static is an expression in eq.eq units*)
V1static=Expand[2*V21static*V2m1static];
V1static=Expand[V1static /.{
Sin[beta1]^2 -> 1-Cos[beta1]^2,
Cos[alpha1]^2 -> (1+Cos[2 alpha1])/2,
Sin[alpha1]^2 -> (1-Cos[2 alpha1])/2}];
(*A1 is in eq.eq units*)
A1=(-3/2)*Collect[Expand[(-2/3)V1static],{Cos[beta1]^4,Cos[beta1]^2}];
(*************************************************)
(*D22 and D2m2 are reduced forms (3 rows x 1 column) of the 2-
nd order Wigner active rotation matrix*)
D22={{(1/4)(1+Cos[beta1])^2*E^(-2I(alpha1+gamma1))},
{Sqrt[3/8]Sin[beta1]^2*E^(-2I gamma1)},
{(1/4)(1-Cos[beta1])^2*E^(2I(alpha1-gamma1))}};
D2m2={{(1/4)(1-Cos[beta1])^2*E^(-2I(alpha1-gamma1))},
{Sqrt[3/8]Sin[beta1]^2*E^(2I gamma1)},
{(1/4)(1+Cos[beta1])^2*E^(2I(alpha1+gamma1))}};
(*V22static and V2m2static are expressions in eq units*)
V22static=FullSimplify[V2pas.ComplexExpand[D22]];
V2m2static=FullSimplify[V2pas.ComplexExpand[D2m2]];
(*V2static is an expression in eq.eq units*)
V2static=TrigExpand[V22static*V2m2static];
V2static=Expand[V2static /.{
Sin[beta1]^2 -> 1-Cos[beta1]^2,
Sin[beta1]^4 -> 1-2 Cos[beta1]^2+Cos[beta1]^4,
Sin[alpha1]^4 -> 1-2 Cos[alpha1]^2+Cos[alpha1]^4,
Sin[2alpha1]^2 -> 1-Cos[2alpha1]^2,
Sin[alpha1]^2 -> 1-Cos[alpha1]^2}];
V2static=Expand[V2static /.{
Cos[alpha1]^2 -> (1+Cos[2alpha1])/2,
Cos[alpha1]^4 -> (1+Cos[2alpha1])^2/4}];
(*A2 and B are in eq.eq units*)
A2=(3/2)*Collect[Expand[(2/3)V2static],{Cos[beta1]^4,Cos[beta1]^2}];
B=(-3/2)*Collect[Expand[(-2/3)(A1+A2)],{Cos[beta1]^4,Cos[beta1]^2}];
(***************************************************)
(*Suppression of the double curve brackets {{}} of A1, A2, and B*)
Print["2V21*V2m1 = ",(e^2)(q^2) A1[[1,1]]];
Print["V22*V2m2 = ",(e^2)(q^2) A2[[1,1]]];
Print["2V21*V2m1 + V22*V2m2 = ",(e^2)(q^2) B[[1,1]]];
Remove[V2pas,eta,D21,D2m1,D22,D2m2, alpha1,beta1,gamma1,V21static,V2m1static,
V1static,V22static,V2m2static, V2static, A1, A2, B];
