crystal symmetry and (χ^((2)))

demonstrate permutation symmetry

define chi tensor components and input fields

chiTi={{χi11,χi12,χi13},{χi21,χi22,χi23},{χi31,χi32,χi33}};
MatrixForm[chiTi]

e1={e1x,e1y,e1z};
e2={e2x,e2y,e2z};

calculate ith polarization component 
in intrinsic perm symm, can swap ω1 and ω2 along with second coordinate indices. This amounts to doing a transpose on the matrix:

Pi12=e1.chiTi.e2

Pi21=e2.Transpose[chiTi].e1

FullSimplify[Pi12-Pi21]

can we swap the fields without doing the transpose?

Pi21=e2.chiTi.e1

FullSimplify[Pi12-Pi21]

gather terms to see what we have
we can see this will be true only if the matrix is symmetric
Symmetry of the matrix is true in second harmonic generation (SHG) automatically. 
It is also true for nondispersive media (Kleinmann symmetry).

e1x ( e2y (χi12-χi21)  + e2z( χi13-χi31) ) +
e1y (e2x (-χi12+χi21) + e2z (χi23-χi32)) + 
e1z (-e2x( χi13-χi31)   -  e2y(χi23-χi32))

crystal symmetries: non-absorbing, non-dispersive

general case

chiTx={{χxxx,χxxy,χxxz},{χxxy,χxyy,χxyz},{χxzx,χxzy,χxzz}};
chiTy={{χyxx,χyxy,χyxz},{χyxy,χyyy,χyyz},{χyzx,χyzy,χyzz}};
chiTz={{χzxx,χzxy,χzxz},{χzxy,χzyy,χzyz},{χzzx,χzzy,χzzz}};
Print["χ's: ",MatrixForm[chiTx]," , ",MatrixForm[chiTy]," , ",MatrixForm[chiTz]];

class 222: many elements are zero, many of remaining are equal

chiTx=
{{0χxxx,0χxxy,0χxxz},
{0χxxy,0χxyy,1χxyz},
{0χxzx,1χxzy,0χxzz}};
chiTy={
{0χyxx,0χyxy,1χyxz},
{0χyxy,0χyyy,0χyyz},
{1χyzx,0χyzy,0χyzz}};
chiTz={
{0χzxx,1χzxy,0χzxz},
{1χzxy,0χzyy,0χzyz},
{0χzzx,0χzzy,0χzzz}};
Print["χ's: ",MatrixForm[chiTx]," , ",MatrixForm[chiTy]," , ",MatrixForm[chiTz]];

class 3m (BBO, LiNbO3)
note that if Kleinman's symmetry is not valid, there are two pairs of off-diagonal elements that aren't necessarily the same, meaning that you can't use the contracted notation.

chiTx=
{{0,-χyyy,χxxz},
{-χyyy,0,0},
{χxzx,0,0}};
chiTy={
{-χyyy,0,0},
{0,χyyy,χxxz},
{0,χxzx,0}};
chiTz={
{χzxx,0,0},
{0,χzxx,0},
{0,0,χzzz}};
Print["χ's: ",MatrixForm[chiTx]," , ",MatrixForm[chiTy]," , ",MatrixForm[chiTz]];

3m contracted notation (Kleinman valid)

d3m={
{0,0,0,0,d31,-d22},
{-d22,d22,0,d31,0,0},
{d31,d31,d33,0,0,0}};
Print["d matrix 3m: ",MatrixForm[d3m]];

in contracted notation, there is only one matrix, and the two electric fields are combined into a 6 element vector:
These combinations for sum frequency mixing are shown in eqn 1.5.27

E1E2Vector[E1_,E2_]:={(Subscript[E1, [[1]]]) (Subscript[E2, [[1]]]),(Subscript[E1, [[2]]]) (Subscript[E2, [[2]]]),(Subscript[E1, [[3]]]) (Subscript[E2, [[3]]]),
(Subscript[E1, [[2]]]) (Subscript[E2, [[3]]])+(Subscript[E1, [[3]]]) (Subscript[E2, [[2]]]),(Subscript[E1, [[1]]]) (Subscript[E2, [[3]]])+(Subscript[E1, [[3]]]) (Subscript[E2, [[1]]]),(Subscript[E1, [[1]]]) (Subscript[E2, [[2]]])+(Subscript[E1, [[2]]]) (Subscript[E2, [[1]]])};
E1={E1x,E1y,E1z};
E2={E2x,E2y,E2z};
E1E2General=E1E2Vector[E1,E2];
Print["general form for E1E2 vector: ",MatrixForm[E1E2General]];

Print["general form for output vector: ",MatrixForm[d3m.E1E2General]];

Example 1: only input along x
we get NL polarization along y and z

E1={A1x,0,0};
E2={A2x,0,0};
E1E2=E1E2Vector[E1,E2];
Print[MatrixForm[d3m],".",MatrixForm[E1E2]," = ",MatrixForm[d3m.E1E2]];

Example 2: only input along y
we get NL polarization along y and z

E1={0,A1y,0};
E2={0,A2y,0};
E1E2=E1E2Vector[E1,E2];
Print[MatrixForm[d3m],".",MatrixForm[E1E2]," = ",MatrixForm[d3m.E1E2]];

Example 3: if Ex1 and Ey2 are the only nonzero terms, we get output polarized along the x direction

E1={A1x,0,0};
E2={0,A2y,0};
E1E2=E1E2Vector[E1,E2];
Print[MatrixForm[d3m],".",MatrixForm[E1E2]," = ",MatrixForm[d3m.E1E2]];

Important comments: 
- just because there is a nonlinear polarization along a particular direction doesn't mean that we will see output there. Propagation has to be reasonably well phase matched to see output along that direction. 
- However, in surface harmonic- or sum-frequency generation, there is no phase matching, and you can use the harmonic output in reflection to test the symmetry of the material.
- coordinate system is for the crystal. Input polarization might be linear, but project along two crystal axes.

Effective d: type I phase matching
reference: Midwinter article

define the direction of the input beam wavevector 
the actual vector is in cartesian coordinates, but we specify the direction angles in spherical

kIn={Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]};

For type I phase matching, input waves are o-waves and parallel to each other:

The o-wave has no z-component (optic axis) of E-field and
E . (Subscript[k, in])=0
E . E=((Subscript[E, 0])^2)
We can use this to define the orientation of the E-field vector.

{Ex,Ey,0}.kIn

θ ≠ 0, or else the light that is generated would also be an o-wave. 
therefore,  (Subscript[E, x])Cos[ϕ] = -(Subscript[E, y])Sin[ϕ]

Ex1=E1 Sin[ϕ];
Ey1=-E1 Cos[ϕ];
Ez1=0;
E1Vect={Ex1,Ey1,Ez1};

Ex2=E2 Sin[ϕ];
Ey2=-E2 Cos[ϕ];
Ez2=0;
E2Vect={Ex2,Ey2,Ez2};

Make a vector of the E-field products that can be used with the d-matrix. These combinations for sum frequency mixing are shown in eqn 1.5.27

E1E2Lst:={Ex1 Ex2,Ey1 Ey2,Ez1 Ez2,
Ey1 Ez2+Ez1 Ey2,Ex1 Ez2+Ez1 Ex2,Ex1 Ey2+Ey1 Ex2};
MatrixForm[E1E2Lst]

calculate the induced polarization

Pinduced=FullSimplify[d3m.E1E2Lst];
MatrixForm[Pinduced]

Note that the nonlinear polarization, which is the source term for the sum frequency light, has components in each direction.
 
But not all of this polarization is in the right direction to propagate along k as it must. It is the e-wave component that will build up due to phase matching. 
Find the projection along this direction by dotting Pinduced with the direction of an e-wave. 
The e-wave is perpendicular to the o-wave and perpendicular to kIn

eWaveVect={eWavex,eWavey,eWavez};
eWaveVect.E1Vect

Setting this to 0, we can conclude that we can write the output wave in the form:

eWaveVect={axy Cos[ϕ],axy Sin[ϕ],eWavez};
eWaveVect.E1Vect

next, this field must be perpendicular to kIn

FullSimplify[eWaveVect.kIn]

again, setting this to 0, with the condition that this should be a unit vector, we have:

eWaveVect={-Cos[θ]Cos[ϕ],-Cos[θ]Sin[ϕ],Sin[θ]};

Finally, we dot this unit vector with the induced polarization to find the component of the NL polarization that drives the e-wave:
This gives a result that is consistent with eqn 1.5.30b

FullSimplify[Pinduced.eWaveVect]

This gives the effective d through the definition: 
P = deff E1 E2