(*This is a Wolfram Mathematica code to establish the relationships between various nonlinear susceptibility tensor components 
in zincblende crystals. Changing the generation matrices, as borrowed from the Ervin Hartman's pamphlet 'An Introduction to Crystal Physics,' University College Cardiff Press (2001), one can get similar relationships for other crystal classes. Nonlinear susceptibility orders 2, 3, 4, 5, and 6 are covered*) 

(*zincblende generators*)

M1={
{0,0,1},{1,0,0},{0,1,0}
};
M2={
{0,-1,0},{1,0,0},{0,0,-1}
};


(*χ2 constraints*)

(*generate empty 3^n array*)
χ=ConstantArray[0,{3,3,3}];

(*define generic tansor elements*)
For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,

χ[[i]][[j]][[k]]=ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]],

k++],
j++],
i++]

EQNS={};
Lin={};

(*Symmetries by generator 1*)

M=M1;

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,

S3=S3+M[[k]][[r]]χ[[p]][[q]][[r]],

r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

(*generate equations*)
AppendTo[EQNS,
χ[[i]][[j]][[k]]==
S1

];

(*generate variables*)
AppendTo[Lin,χ[[i]][[j]][[k]]],

k++],
j++],
i++]


(*Symmetries by generator 2*)

M=M2;
Lin={};

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,

S3=S3+M[[k]][[r]]χ[[p]][[q]][[r]],

r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];


AppendTo[EQNS,
χ[[i]][[j]][[k]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]]],
k++],
j++],
i++]

Print["Constraints on χ(2) elements"]
sol=Solve[EQNS,Lin]


(*χ4 constraints*)

χ=ConstantArray[0,{3,3,3,3,3}];

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,

χ[[i]][[j]][[k]][[l]][[m]]=ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]],

m++],
l++],
k++],
j++],
i++]

EQNS={};
Lin={};

(*Generator 1*)

M=M1;

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,
S5=0;
For[t=0,t<3,

S5=S5+M[[m]][[t]]χ[[p]][[q]][[r]][[s]][[t]],

t++];
S4=S4+M[[l]][[s]]S5,
s++];
S3=S3+M[[k]][[r]]S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]][[m]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]][[m]]],

m++],
l++],
k++],
j++],
i++];


(*Generator 2*)

Lin={};
M=M2;

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,
S5=0;
For[t=0,t<3,

S5=S5+M[[m]][[t]]χ[[p]][[q]][[r]][[s]][[t]],

t++];
S4=S4+M[[l]][[s]]S5,
s++];
S3=S3+M[[k]][[r]]S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]][[m]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]][[m]]],

m++],
l++],
k++],
j++],
i++];

Print["Constraints on χ(4) elements"]
sol=Solve[EQNS,Lin];

Print["Number of independent components"]
Dimensions[Lin][[1]]-Dimensions[sol][[2]]



(*χ6 constraints*)

χ=ConstantArray[0,{3,3,3,3,3,3,3}];

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,
For[o=0,o<3,

χ[[i]][[j]][[k]][[l]][[m]][[n]][[o]]=ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]<>ToString[n]<>ToString[o]],

o++],
n++],
m++],
l++],
k++],
j++],
i++];

EQNS={};
Lin={};

(*Generator 1*)

M=M1;

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,
For[o=0,o<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,
S5=0;
For[t=0,t<3,
S6=0;
For[u=0,u<3,
S7=0;
For[v=0,v<3,

S7=S7+M[[o]][[v]]χ[[p]][[q]][[r]][[s]][[t]][[u]][[v]],

v++];
S6=S6+M[[n]][[u]]S7,
u++];
S5=S5+M[[m]][[t]]S6,
t++];
S4=S4+M[[l]][[s]]S5,
s++];
S3=S3+M[[k]][[r]]S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]][[m]][[n]][[o]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]][[m]][[n]][[o]]],

o++],
n++],
m++],
l++],
k++],
j++],
i++];

M=M2;
Lin={};

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,
For[o=0,o<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,
S5=0;
For[t=0,t<3,
S6=0;
For[u=0,u<3,
S7=0;
For[v=0,v<3,

S7=S7+M[[o]][[v]]χ[[p]][[q]][[r]][[s]][[t]][[u]][[v]],

v++];
S6=S6+M[[n]][[u]]S7,
u++];
S5=S5+M[[m]][[t]]S6,
t++];
S4=S4+M[[l]][[s]]S5,
s++];
S3=S3+M[[k]][[r]]S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]][[m]][[n]][[o]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]][[m]][[n]][[o]]],

o++],
n++],
m++],
l++],
k++],
j++],
i++];

Print["Constraints on χ(6) elements"]
sol=Solve[EQNS,Lin]

Print["Number of independent components"]
Dimensions[Lin][[1]]-Dimensions[sol][[2]]

(*Calculations of nonlinear polarizations generated by an external field as a function of the crystal lattice orientaion (see Methods in the paper for the description). All nonzero components are assumed to be = 1 *)

conds={};

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,
For[o=0,o<3,

AppendTo[conds,Solve[ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]<>ToString[n]<>ToString[o]]==1,ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]<>ToString[n]<>ToString[o]]][[1]][[1]]],

o++],
n++],
m++],
l++],
k++],
j++],
i++];

Clear[θ];

R=(Cos[(π θ)/180]	-(Sin[(π θ)/180]/Sqrt[2])	-(Sin[(π θ)/180]/Sqrt[2])
Sin[(π θ)/180]/Sqrt[2]	1/2 (1-Cos[(π θ)/180])+Cos[(π θ)/180]	1/2 (-1+Cos[(π θ)/180])
Sin[(π θ)/180]/Sqrt[2]	1/2 (-1+Cos[(π θ)/180])	1/2 (1-Cos[(π θ)/180])+Cos[(π θ)/180]

);

Ein=R.{0,1,0}
IH6=Norm[N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein.Ein/.conds][[1]]]^2;
IH6proj=Norm[N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein.Ein/.conds][[1]]];

H6x=N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein.Ein/.conds][[1]].{1,0,0};
H6y=N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein.Ein/.conds][[1]].{0,1,0};
H6z=N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein.Ein/.conds][[1]].{0,0,1};

Ex=H6x Cos[θ π/180]+H6y Cos[θ π/180]+H6z Sin[θ π/180];
LogPlot[IH6,{θ,0,90},
Frame->True,FrameLabel->{"E_loc to [001] axis tilt (degrees)","H6 polarization (arb. un.)"},FrameStyle->Black,
LabelStyle->Directive[FontFamily->"Helvetica",FontSize->16,FontColor->Black],PlotStyle->{Directive[Black]},AspectRatio->1/2.5,PlotLegends->Placed[{"P^2","Subscript[P, x]^2,Subscript[P, y]^2","Subscript[P, z]^2"},{Right,Top}],FrameTicks->{Automatic,{{0,30,60,90},None}},PlotRange->{0.01,200}]

LogPlot[
{IH6, H6x^2,H6z^2},{θ,0,90},
Frame->True,FrameLabel->{"[001] axis tilt (degrees)","H6 polarization (arb. un.)"},FrameStyle->Black,
LabelStyle->Directive[FontFamily->"Helvetica",FontSize->16,FontColor->Black],PlotStyle->{Directive[Red],Directive[Blue,Dashed],Directive[Black,Dotted]},AspectRatio->1/1.4,PlotLegends->Placed[{"P^2","Subscript[P, x]^2,Subscript[P, y]^2","Subscript[P, z]^2"},{Right,Bottom}],FrameTicks->{Automatic,{{0,30,60,90},None}},PlotRange->{0.01,200}]



(*χ3 constraints*)

χ=ConstantArray[0,{3,3,3,3}];
For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0, l<3,
χ[[i]][[j]][[k]][[l]]=ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]],
l++],
k++],
j++],
i++]

EQNS={};
Lin={};

(*Generator 1*)

M=M1;

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,

S4=S4+M[[l]][[s]]χ[[p]][[q]][[r]][[s]],

s++];
S3=S3+M[[k]][[r]] S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]]],

l++],
k++],
j++],
i++]


(*Generator 2*)

M=M2;
Lin={};

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,

S4=S4+M[[l]][[s]]χ[[p]][[q]][[r]][[s]],

s++];
S3=S3+M[[k]][[r]] S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]]],

l++],
k++],
j++],
i++]


Print["Constraints on χ(3) elements"]
sol=Solve[EQNS,Lin]

Print["Number of independent components"]
Dimensions[Lin][[1]]-Dimensions[sol][[2]]



(*χ5 constraints*)

χ=ConstantArray[0,{3,3,3,3,3,3}];

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,

χ[[i]][[j]][[k]][[l]][[m]][[n]]=ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]<>ToString[n]],

n++],
m++],
l++],
k++],
j++],
i++];

EQNS={};
Lin={};

(*Generator 1*)

M=M1;

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,




S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,
S5=0;
For[t=0,t<3,
S6=0;
For[u=0,u<3,

S6=S6+M[[n]][[u]]χ[[p]][[q]][[r]][[s]][[t]][[u]],

u++];
S5=S5+M[[m]][[t]]S6,
t++];
S4=S4+M[[l]][[s]]S5,
s++];
S3=S3+M[[k]][[r]]S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]][[m]][[n]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]][[m]][[n]]],

n++],
m++],
l++],
k++],
j++],
i++];

M=M2;
Lin={};


For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,




S1=0;
For[p=0,p<3,
S2=0;
For[q=0,q<3,
S3=0;
For[r=0,r<3,
S4=0;
For[s=0,s<3,
S5=0;
For[t=0,t<3,
S6=0;
For[u=0,u<3,

S6=S6+M[[n]][[u]]χ[[p]][[q]][[r]][[s]][[t]][[u]],

u++];
S5=S5+M[[m]][[t]]S6,
t++];
S4=S4+M[[l]][[s]]S5,
s++];
S3=S3+M[[k]][[r]]S4,
r++];
S2=S2+M[[j]][[q]]S3,
q++];
S1=S1+M[[i]][[p]]S2,
p++];

AppendTo[EQNS,
χ[[i]][[j]][[k]][[l]][[m]][[n]]==
S1

];

AppendTo[Lin,χ[[i]][[j]][[k]][[l]][[m]][[n]]],

n++],
m++],
l++],
k++],
j++],
i++];

Print["Constraints on χ(5) elements"]
sol=Solve[EQNS,Lin]

Print["Number of independent components"]
Dimensions[Lin][[1]]-Dimensions[sol][[2]]


conds={};

For[i=0,i<3,
For[j=0,j<3,
For[k=0,k<3,
For[l=0,l<3,
For[m=0,m<3,
For[n=0,n<3,

AppendTo[conds,Solve[ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]<>ToString[n]]==1,ToExpression["f"<>ToString[i]<>ToString[j]<>ToString[k]<>ToString[l]<>ToString[m]<>ToString[n]]][[1]][[1]]],

n++],
m++],
l++],
k++],
j++],
i++];

Clear[θ];
(*θ=π/12;*)
(*ϕ=π/4;*)

R=(Cos[(π θ)/180]	-(Sin[(π θ)/180]/Sqrt[2])	-(Sin[(π θ)/180]/Sqrt[2])
Sin[(π θ)/180]/Sqrt[2]	1/2 (1-Cos[(π θ)/180])+Cos[(π θ)/180]	1/2 (-1+Cos[(π θ)/180])
Sin[(π θ)/180]/Sqrt[2]	1/2 (-1+Cos[(π θ)/180])	1/2 (1-Cos[(π θ)/180])+Cos[(π θ)/180]

);

Ein=R.{0,1,0};


IH5=Norm[N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein/.conds][[1]]]^2;
IH5proj=Norm[N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein/.conds][[1]]];

H5x=N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein/.conds][[1]].{1,0,0};
H5y=N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein/.conds][[1]].{0,1,0};
H5z=N[Dot[χ/.sol, Ein].Ein.Ein.Ein.Ein/.conds][[1]].{0,0,1};

LogPlot[
{IH5},{θ,0,90},
Frame->True,FrameLabel->{"[001] axis tilt (degrees)","H5, 0th order emission (arb. un.)"},FrameStyle->Black,
LabelStyle->Directive[FontFamily->"Helvetica",FontSize->16,FontColor->Black],PlotStyle->{Directive[Black]},AspectRatio->1/1.4,PlotLegends->Placed[{"P^2","Subscript[P, x]^2,Subscript[P, y]^2","Subscript[P, z]^2"},{Right,Top}],FrameTicks->{Automatic,{{0,30,60,90},None}}]

LogPlot[
{IH5, H5x^2,H5y^2},{θ,0,90},
Frame->True,FrameLabel->{"[001] axis tilt (degrees)","H5 polarization (arb. un.)"},FrameStyle->Black,
LabelStyle->Directive[FontFamily->"Helvetica",FontSize->16,FontColor->Black],PlotStyle->{Directive[Red],Directive[Blue,Dashed],Directive[Black,Dotted]},AspectRatio->1/1.4,PlotLegends->Placed[{"|P|^2","Subscript[P, x]^2,Subscript[P, z]^2","Subscript[P, y]^2"},{Right,Bottom}],FrameTicks->{Automatic,{{0,30,60,90},None}},PlotRange->{0.1,50}]