Implemented calculations of decay-rate enhancement for point-dipole i…

…ncident beam. Total part is calculated with the help of a new function DecayCross() in crosssec.c/h (similar to Cext), non-radiative part - by normalizing Cabs. In case of surface, all parts of decay-rate enhancement (total, rad, nonrad) are normalized to that in free space, but additionally enhancement due to surface alone is shown.

Still, it is not clear whether the non-radiative part should be adjusted to account for absorption in the substrate (when it is absorbing).

Also, incident polarizations (Y and X) are not shown in the log for point-dipole incident beam. This is a bit problematic for standard LDR polarizability, since its factor S depends on both prop and incPol. But it probably makes little sense to use LDR polarizability for such incident field anyway (recommendation should be added to the manual).

One test was run for dipole near a sphere was run and matched exactly to first row of Table 1 in S. D'Agostino, F. Della Sala, and L.C. Andreani, “Dipole-excited surface plasmons in metallic nanoparticles: Engineering decay dynamics within the discrete-dipole approximation,” Phys. Rev. B 87, 205413 (2013). However, a few hacks were used to exactly match the simulation conditions (based on raw data provided by Stefania D'Agostino.

src/CalculateE.c

@@ -46,6 +46,8 @@ extern double * restrict muel_alpha;
 // defined and initialized in crosssec.c
 extern const Parms_1D phi_sg;
 extern const double ezLab[3],exSP[3];
+// defined and initialized in GenerateB.c
+extern const double C0dipole,C0dipole_refl;
 // defined and initialized in param.c
 extern const bool store_int_field,store_dip_pol,store_beam,store_scat_grid,calc_Cext,calc_Cabs,
 	calc_Csca,calc_vec,calc_asym,calc_mat_force,store_force,store_ampl;
@@ -716,6 +718,22 @@ static void CalcIntegralScatQuantities(const enum incpol which)
 			CCfile=FOpenErr(fname_cs,"w",ONE_POS);
 			if (calc_Cext) PrintBoth(CCfile,"Cext\t= "GFORM"\nQext\t= "GFORM"\n",Cext,Cext*inv_G);
 			if (calc_Cabs) PrintBoth(CCfile,"Cabs\t= "GFORM"\nQabs\t= "GFORM"\n",Cabs,Cabs*inv_G);
+			if (beamtype==B_DIPOLE) {
+				double self=1;
+				if (surface) self+=C0dipole_refl/C0dipole;
+				double tot=self+DecayCross()/C0dipole;
+				fprintf(CCfile,"\nDecay-rate enhancement\n\n");
+				printf("\nDecay-rate enhancement:\n");
+				PrintBoth(CCfile,"Total\t= "GFORM"\n",tot);
+				if (calc_Cabs) { // for simplicity we keep a single condition here
+					double nonrad=Cabs/C0dipole;
+					double rad=tot-nonrad;
+					// TODO: !!! how nonrad should account for absorption inside the surface???
+					PrintBoth(CCfile,"Radiat\t= "GFORM"\n",rad);
+					PrintBoth(CCfile,"Nonrad\t= "GFORM"\n",nonrad);
+				}
+				if (surface) PrintBoth(CCfile,"Surface\t= "GFORM"\n",self);
+			}
 			if (all_dir) fprintf(CCfile,"\nIntegration\n\n");
 			if (calc_Csca) {
 				Csca=ScaCross(f_suf);

src/GenerateB.c

@@ -46,6 +46,9 @@ extern const char *beam_fnameY;
 extern const char *beam_fnameX;
 extern const opt_index opt_beam;
 
+// used in CalculateE.c
+double C0dipole,C0dipole_refl; // inherent cross sections of exciting dipole (in free space and addition due to surface)
+
 // used in crosssec.c
 double beam_center_0[3]; // position of the beam center in laboratory reference frame
 /* complex wave amplitudes of secondary waves (with phase relative to particle center);
@@ -327,6 +330,31 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 					cvAdd(v1,b+j,b+j);
 				}
 			}
+			/* calculate dipole inherent cross sections (for decay rate enhancements)
+			 * General formula is C0=4pi*k*Im(p0*.G(r0,r0).p0) and it is used for reflection; but for direct
+			 * interaction it is expected to be singular, so we use an explicit formula for point dipole:
+			 * C0=(8pi/3)*k^4*|p0|^2. Thus we discard the choice of "-int ...", but still the used value should be
+			 * correct for most formulations, e.g. poi, fcd, fcd_st, igt_so. Moreover, it is also logical, since the
+			 * exciting dipole is really point one, in contrast to the dipoles composing the particle.
+			 */
+			C0dipole=2*FOUR_PI_OVER_THREE*WaveNum*WaveNum*WaveNum*WaveNum;
+			printf("C0=%g\n",C0dipole);
+			if (surface) {
+				r1[0]=r1[1]=0;
+				r1[2]=2*(beam_center[2]+hsub);
+				(*ReflTerm_real)(r1,gt);
+				double tmp;
+				/* the following expression uses that prop is real and a specific (anti-)symmetry of the gt
+				 * a general expression is commented out below
+				 */
+				tmp=prop[0]*prop[0]*cimag(gt[0])+2*prop[0]*prop[1]*cimag(gt[1])+prop[1]*prop[1]*cimag(gt[3])
+				   +prop[2]*prop[2]*cimag(gt[5]);
+//				v1[0]=prop[0]; v1[1]=prop[1]; v1[2]=prop[2];
+//				cReflMatrVec(gt,v1,v2);
+//				tmp=cDotProd_Im(v2,v1);
+				C0dipole_refl=FOUR_PI*WaveNum*tmp;
+				printf("C0refl=%g\n",C0dipole_refl);
+			}
 			return;
 		case B_LMINUS:
 		case B_DAVIS3:

src/crosssec.c

@@ -948,6 +948,29 @@ double AbsCross(void)
 
 //======================================================================================================================
 
+double DecayCross(void)
+// computes total cross section for the dipole incident field; similar to Cext
+// 4pi*k*Im[p0(*).Escat(r0)]
+{
+	double sum;
+	size_t i;
+
+	/* This is a correct expression only _if_ exciting p0 is real, then
+	 * (using G(r1,r2) = G(r2,r1)^T, valid also for surface)
+	 * p0(*).Escat_i(r0) = p0(*).G_0i.p_i = p_i.G_i0.p0(*) = p_i.G_i0.p0 = p_i.Einc_i
+	 * => Im(p0(*).Escat(r0)) = sum{Im(P.E_inc)}
+	 *
+	 * For complex p0 an efficient calculation strategy (not to waste evaluations of interaction) is to compute an array
+	 * of G_i0.p0(*) together with Einc and use it here afterwards.
+	 */
+	sum=0;
+	for (i=0;i<local_nvoid_Ndip;++i) sum+=cimag(cDotProd_conj(pvec+3*i,Einc+3*i)); // sum{Im(P.E_inc)}
+	MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm);
+	return FOUR_PI*WaveNum*sum;
+}
+
+//======================================================================================================================
+
 void SetScatPlane(const double ct,const double st,const double phi,double robs[static restrict 3],
 	double polPer[static restrict 3])
 /* Given theta (cos,sin) and phi, calculates scattering direction and unit vector perpendicular to the scattering plane.

src/crosssec.h

@@ -25,6 +25,7 @@ void CalcField(doublecomplex ebuff[static restrict 3],const double n[static rest
 void InitRotation(void);
 double ExtCross(const double * restrict incPol);
 double AbsCross(void);
+double DecayCross(void);
 double ScaCross(const char *f_suf);
 void ReadAlldirParms(const char * restrict fname);
 void ReadAvgParms(const char * restrict fname);

src/param.c

@@ -544,7 +544,7 @@ static struct opt_struct options[]={
 	{PAR(prognosis),"","Do not actually perform simulation (not even memory allocation) but only estimate the required "
 		"RAM. Implies '-test'.",0,NULL},
 	{PAR(prop),"<x> <y> <z>","Sets propagation direction of incident radiation, float. Normalization (to the unity "
-		"vector) is performed automatically.\n"
+		"vector) is performed automatically. For point-dipole incident beam this determines its direction.\n"
 		"Default: 0 0 1",3,NULL},
 	{PAR(recalc_resid),"","Recalculate residual at the end of iterative solver.",0,NULL},
 #ifndef SPARSE
@@ -2198,9 +2198,12 @@ void PrintInfo(void)
 		// log incident beam and polarization
 		fprintf(logfile,"\n---In laboratory reference frame:---\nIncident beam: %s\n",beam_descr);
 		fprintf(logfile,"Incident propagation vector: "GFORMDEF3V"\n",COMP3V(prop_0));
-		fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY_0));
-		fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n",COMP3V(incPolX_0));
-		if (surface) { // include surface-specific vectors
+		if (beamtype==B_DIPOLE) fprintf(logfile,"(dipole orientation)\n");
+		else { // polarizations are not shown for dipole incident field
+			fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY_0));
+			fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n",COMP3V(incPolX_0));
+		}
+		if (surface && beamtype==B_DIPOLE) { // include surface-specific vectors
 			fprintf(logfile,"Reflected propagation vector: "GFORMDEF3V"\n",COMP3V(prIncRefl));
 			fprintf(logfile,"Transmitted propagation vector: ");
 			if (msubInf) fprintf (logfile,"none\n");
@@ -2221,10 +2224,14 @@ void PrintInfo(void)
 				if (beam_asym) fprintf(logfile,"Incident Beam center position: "GFORMDEF3V"\n",
 					beam_center[0],beam_center[1],beam_center[2]);
 				fprintf(logfile,"Incident propagation vector: "GFORMDEF3V"\n",COMP3V(prop));
-				fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY));
-				fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n\n",COMP3V(incPolX));
+				if (beamtype==B_DIPOLE) fprintf(logfile,"(dipole orientation)\n");
+				else { // polarizations are not shown for dipole incident field
+					fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY));
+					fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n",COMP3V(incPolX));
+				}
 			}
-			else fprintf(logfile,"Particle orientation: default\n\n");
+			else fprintf(logfile,"Particle orientation: default\n");
+			fprintf(logfile,"\n");
 		}
 		// log type of scattering matrices
 		if (store_mueller) {