With using ?NumericQ
we have
\[Gamma]21 = 0.1; \[Gamma]31 = 0.1; \[Gamma]41 = 1; w = 1; kp = 1; \
\[Beta] = 1; \[CapitalDelta]2 = 0; \[CapitalDelta]3 = 0; \
\[CapitalDelta]1 = 1;
\[Phi] = Pi/6; \[CapitalDelta]33 = \[CapitalDelta]3 + kp*v;
\[CapitalDelta]11 = \[CapitalDelta]1 +
kp*v; \[CapitalDelta]22 = \[CapitalDelta]2 +
kp*v; d = 1; M1 = 4; N1 = 4; \[Lambda] = 0.25;
\[Phi] = Pi/6; L = 6;
m1 = 0; n1 = 0;
\[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1; r = 1; l = 3;
\[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)]*(r/w)^l;
\[CapitalOmega]3 = \[CapitalOmega]f1*
Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22; \[CapitalLambda] = 1;
\[CapitalDelta]13 = \[CapitalDelta]11 - \[CapitalDelta]33;
\[CapitalOmega]2 = \[CapitalOmega]20*(Sin[Pi*(x/\[CapitalLambda])] +
Sin[Pi*(y/\[CapitalLambda])]);
\[Chi]p = (I*(\[Gamma]21 - I*\[CapitalDelta]12)*(\[Gamma]31 -
I*\[CapitalDelta]13))/((\[Gamma]41 -
I*\[CapitalDelta]11)*(\[Gamma]21 -
I*\[CapitalDelta]12)*(\[Gamma]31 -
I*\[CapitalDelta]13) + (\[Gamma]31 -
I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21 -
I*\[CapitalDelta]12)*\[CapitalOmega]3^2);
a1 = ComplexExpand[Im[\[Chi]p]];
aa[v_] := (1/(d*Sqrt[2*Pi]))*Exp[-v^2/(2*d^2)];
aaa[x1_?NumericQ, y1_?NumericQ] :=
NIntegrate[a1*aa[v] /. {x -> x1, y -> y1}, {v, -10, 10}];
b1 = ComplexExpand[Re[\[Chi]p]];
cc1[x1_?NumericQ, y1_?NumericQ] :=
NIntegrate[b1*aa[v] /. {x -> x1, y -> y1}, {v, -10, 10}];
tt[x_, y_] := Exp[(-aaa[x, y])*L];
TT1[x_, y_] := Exp[(-aaa[x, y])*L + I*cc1[x, y]];
q[\[Theta]1_, x_, \[Theta]2_, y_] :=
Exp[(-aaa[x, y])*L + I*cc1[x, y]]*
Exp[(-2*Pi*I*\[CapitalLambda]*x*Sin[\[Theta]1])/\[Lambda]]*
Exp[(-2*Pi*I*\[CapitalLambda]*y*Sin[\[Theta]2])/\[Lambda]];
dd[\[Theta]1_?NumericQ, \[Theta]2_?NumericQ] :=
NIntegrate[
q[\[Theta]1, x, \[Theta]2, y], {x, 0, 0.001}, {y, 0, 0.001}]
y[\[Theta]1_, \[Theta]2_] := Abs[dd[\[Theta]1, \[Theta]2]]^2
s[\[Theta]1_, \[Theta]2_] := (Sin[(M1*Pi*\[CapitalLambda]*
Sin[\[Theta]1])/\[Lambda]]^2/(M1^2*
Sin[(Pi*\[CapitalLambda]*
Sin[\[Theta]1])/\[Lambda]]^2))*(Sin[(N1*
Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2/
(N1^2*Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2))
z[\[Theta]1_, \[Theta]2_] :=
y[\[Theta]1, \[Theta]2]*s[\[Theta]1, \[Theta]2]
DensityPlot[
Evaluate[z[\[Theta]1, \[Theta]2]], {\[Theta]1, -0.6,
0.6}, {\[Theta]2, -0.6, 0.6}, PlotRange -> All, PlotPoints -> 50,
ColorFunction -> "Rainbow",
PlotLegends ->
Placed[BarLegend[Automatic, LegendMarkerSize -> 280], Right],
ImageSize -> 300, Background -> Transparent, FrameStyle -> Black,
LabelStyle -> {Black, FontSize -> 18}, RotateLabel -> True]

Update 1. In a case of integration limits {v, -Infinity, Infinity}
and {x, 0, 1}, {y, 0, 1}
we modified code as follows (we also update function s
to exclude singularities)
\[Gamma]21 = 0.1; \[Gamma]31 = 0.1; \[Gamma]41 = 1; w = 1; kp = 1; \
\[Beta] = 1; \[CapitalDelta]2 = 0; \[CapitalDelta]3 = 0; \
\[CapitalDelta]1 = 1;
\[Phi] = Pi/6; \[CapitalDelta]33 = \[CapitalDelta]3 + kp*v;
\[CapitalDelta]11 = \[CapitalDelta]1 +
kp*v; \[CapitalDelta]22 = \[CapitalDelta]2 +
kp*v; d = 1; M1 = 4; N1 = 4; \[Lambda] = 0.25;
\[Phi] = Pi/6; L = 6;
m1 = 0; n1 = 0;
\[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1; r = 1; l = 3;
\[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)]*(r/w)^l;
\[CapitalOmega]3 = \[CapitalOmega]f1*
Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22; \[CapitalLambda] = 1;
\[CapitalDelta]13 = \[CapitalDelta]11 - \[CapitalDelta]33;
\[CapitalOmega]2 = \[CapitalOmega]20*(Sin[Pi*(x/\[CapitalLambda])] +
Sin[Pi*(y/\[CapitalLambda])]);
\[Chi]p = (I*(\[Gamma]21 - I*\[CapitalDelta]12)*(\[Gamma]31 -
I*\[CapitalDelta]13))/((\[Gamma]41 -
I*\[CapitalDelta]11)*(\[Gamma]21 -
I*\[CapitalDelta]12)*(\[Gamma]31 -
I*\[CapitalDelta]13) + (\[Gamma]31 -
I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21 -
I*\[CapitalDelta]12)*\[CapitalOmega]3^2);
a1 = ComplexExpand[Im[\[Chi]p]];
aa[v_] := (1/(d*Sqrt[2*Pi]))*Exp[-v^2/(2*d^2)];
aaa[x1_?NumericQ, y1_?NumericQ] :=
NIntegrate[
a1*aa[v] /. {x -> x1, y -> y1}, {v, -Infinity, Infinity}];
b1 = ComplexExpand[Re[\[Chi]p]];
cc1[x1_?NumericQ, y1_?NumericQ] :=
NIntegrate[
b1*aa[v] /. {x -> x1, y -> y1}, {v, -Infinity, Infinity}];
tt[x_, y_] := Exp[(-aaa[x, y])*L];
TT1[x_, y_] := Exp[(-aaa[x, y])*L + I*cc1[x, y]];
q[\[Theta]1_, x_, \[Theta]2_, y_] :=
Exp[(-aaa[x, y])*L + I*cc1[x, y]]*
Exp[(-2*Pi*I*\[CapitalLambda]*x*Sin[\[Theta]1])/\[Lambda]]*
Exp[(-2*Pi*I*\[CapitalLambda]*y*Sin[\[Theta]2])/\[Lambda]];
dd[\[Theta]1_?NumericQ, \[Theta]2_?NumericQ] :=
NIntegrate[q[\[Theta]1, x, \[Theta]2, y], {x, 0, 1}, {y, 0, 1},
PrecisionGoal -> 2, AccuracyGoal -> 2]
y[\[Theta]1_, \[Theta]2_] := Abs[dd[\[Theta]1, \[Theta]2]]^2
s[\[Theta]1_, \[Theta]2_] :=
If[Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]1])/\[Lambda]]^2 <= 10^-10,
1, (Sin[(M1*Pi*\[CapitalLambda]*
Sin[\[Theta]1])/\[Lambda]]^2/(M1^2*
Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]1])/\[Lambda]]^2))]*
If[Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2 <= 10^-10,
1, (Sin[(N1*Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2/
(N1^2*
Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2))];
z[\[Theta]1_, \[Theta]2_] :=
y[\[Theta]1, \[Theta]2]*s[\[Theta]1, \[Theta]2]
To plot z
we first prepare interpolation function
lst = Table[{\[Theta]1, \[Theta]2, {y[\[Theta]1, \[Theta]2]}}, {\
\[Theta]1, -0.6, 0.6, .1}, {\[Theta]2, -0.6, 0.6, .1}]
Y = Interpolation[Flatten[lst, 1]];
Finally we plot functions s, Y, z
{DensityPlot[
Evaluate[s[\[Theta]1, \[Theta]2]], {\[Theta]1, -0.6,
0.6}, {\[Theta]2, -0.6, 0.6}, PlotRange -> All, PlotPoints -> 100,
ColorFunction -> "Rainbow",
PlotLegends -> Automatic],
DensityPlot[Y[t1, t2], {t1, -.6, .6}, {t2, -.6, .6},
PlotRange -> All, ColorFunction -> "Rainbow",
PlotLegends -> Automatic, PlotPoints -> 100],
DensityPlot[
s[\[Theta]1, \[Theta]2] Y[\[Theta]1, \[Theta]2], {\[Theta]1, -0.6,
0.6}, {\[Theta]2, -0.6, 0.6}, PlotRange -> All, PlotPoints -> 100,
ColorFunction -> "TemperatureMap",
PlotLegends ->
Placed[BarLegend[Automatic, LegendMarkerSize -> 280], Right],
ImageSize -> 300, Background -> Transparent, FrameStyle -> Black,
LabelStyle -> {Black, FontSize -> 18}, RotateLabel -> True]}

Update 2. This code is answer to the question on this page
Clear["Global`*"]
plotset = {FrameStyle -> Directive[Thickness[0.004]],
TicksStyle -> Directive[Black, 18]};
plotset2 =
FrameTicksStyle -> {{Directive[Black, 18],
Directive[FontOpacity -> 0, FontSize -> 0]}, {Directive[Black,
18], Directive[FontOpacity -> 0, FontSize -> 0]}};
\[Gamma]21 = 0.1; \[Gamma]31 = 0.1; \[Gamma]41 = 1; w = 1; kp = 1; \
\[Beta] = 1; \[CapitalDelta]2 = 0; \[CapitalDelta]3 = 0; \
\[CapitalDelta]1 = 1;
\[CapitalDelta]33 = \[CapitalDelta]3 + kp*v;
\[CapitalDelta]11 = \[CapitalDelta]1 +
kp*v; \[CapitalDelta]22 = \[CapitalDelta]2 +
kp*v; d = 1; M1 = 4; N1 = 4; \[Lambda] = 0.25;
L = 6;
m1 = 0; n1 = 0;
\[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1; r = 1; l = 1;
\[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)]*(r/w)^l;
\[CapitalOmega]3 = \[CapitalOmega]f1*
Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22; \[CapitalLambda] = 1;
\[CapitalDelta]13 = \[CapitalDelta]11 - \[CapitalDelta]33;
\[CapitalOmega]2 = \[CapitalOmega]20*(Sin[Pi*(x/\[CapitalLambda])] +
Sin[Pi*(y/\[CapitalLambda])]);
\[Chi]p = (I*(\[Gamma]21 - I*\[CapitalDelta]12)*(\[Gamma]31 -
I*\[CapitalDelta]13))/((\[Gamma]41 -
I*\[CapitalDelta]1)*(\[Gamma]21 -
I*\[CapitalDelta]12)*(\[Gamma]31 -
I*\[CapitalDelta]13) + (\[Gamma]31 -
I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21 -
I*\[CapitalDelta]12)*\[CapitalOmega]3^2);
a1 = ComplexExpand[Im[\[Chi]p]];
aa[v_] := (1/(d*Sqrt[2*Pi]))*Exp[-v^2/(2*d^2)];
aaa[x1_?NumericQ, y1_?NumericQ, phi_?NumericQ] :=
NIntegrate[
a1*aa[v] /. {x -> x1, y -> y1, \[Phi] -> phi}, {v, -Infinity,
Infinity}];
b1 = ComplexExpand[Re[\[Chi]p]];
cc1[x1_?NumericQ, y1_?NumericQ, phi_?NumericQ] :=
NIntegrate[
b1*aa[v] /. {x -> x1, y -> y1, \[Phi] -> phi}, {v, -Infinity,
Infinity}];
tt[x_, y_, \[Phi]_] := Exp[(-aaa[x, y, \[Phi]])*L];
TT1[x_, y_, \[Phi]_] :=
Exp[(-aaa[x, y, \[Phi]])*L + I*cc1[x, y, \[Phi]]];
q[x_, y_, \[Phi]_] :=
Exp[-aaa[x, y, \[Phi]]*L + I*cc1[x, y, \[Phi]]]*Exp[-2*Pi*I*m1*x]*
Exp[-2*Pi*I*n1*y];
pp1[\[Phi]_?NumericQ] :=
NIntegrate[q[x, y, \[Phi]], {x, 0, 1}, {y, 0, 1},
PrecisionGoal -> 2, AccuracyGoal -> 2];
yy[\[Phi]_] := (Abs[pp1[\[Phi]]])^2;
Y = Interpolation[Table[{x, yy[x]}, {x, 0, 4, .1}]]
Visualization
p1 = Plot[Y[\[Phi]], {\[Phi], 0, 4}, PlotRange -> {All, All},
PlotStyle -> {Blue, Thickness[0.007]}, GridLines -> Automatic,
Frame -> True,
FrameLabel -> {Style["", 18, Bold], Style["", 18, Bold]},
PlotLegends ->
Placed[LineLegend[{"\!\(\*SubscriptBox[\(I\), \
\(p\)]\)(\!\(\*SubsuperscriptBox[\(\[Theta]\), \(x\), \
\(0\)]\),\!\(\*SubsuperscriptBox[\(\[Theta]\), \(y\), \(0\)]\))"},
LegendMarkers -> Automatic,
LegendMarkerSize -> {{30, 25}}], {After, Top}],
Evaluate@plotset, Evaluate@plotset2, Axes -> True] //
AbsoluteTiming

In[*]
entries. Also provide the definitions forx
andv
. Thanks. $\endgroup$