<filename>tex/discussion.tex
%!TEX root = ClementiBarba2020.tex

We have presented several results that replicate previous published findings in the general area of nanostructure responses to electromagnetic waves. 
Our field of interest is computational nanoplasmonics for applications in biosensors, and in a previous publication we developed the mathematical formulation and reported both solution verification activities, and an application demo with our software \pygbe, extended to treat complex dielectrics and imposed electric fields \cite{ClementiETal2019}.
The search for a physical context and published results that would allow us to undertake validation studies with \pygbe is what led us to this work. 
Even if we finally do have validation and replication cases, neatly presented here, the path to obtain these results was nonlinear, iterative, and arduous.

In the first case, we sought to replicate a result from Rockstuhl et al., 2005 \cite{rockstuhl2005}, where they computed the scattering cross-section as a function of wavelength for a silicon carbide rectangular nanostructure. 
They present their results as a plot (Figure 14 in their paper), and report the numeric value of the resonance wavelengths in the text. 
Lacking access to the secondary data behind the plots---computed from two-dimensional simulations with a boundary element solver---we were forced to manually digitize the values from the figure.
Our results are presented in Figure \ref{fig:rep_14}, together with the curve we obtain from digitizing the source image. 
We were successful at replicating the strong peaks reported  by Rockstuhl et al.\ at wavelengths
10.42 $\mu$m and 10.7 $\mu$ m,  when the electric field $E$ is parallel to the short side of the rectangle, and 10.42 $\mu$m and 10.82 $\mu$ m
when $E$ is parallel to the long side. Our results contain extra (small) peaks that are not present in the work of Rockstul et al.
The first one, located between the main two peaks, we attribute to the the effect of
sharp edges (see Figure \ref{fig:tri_reg_round_14}), as when we introduce a level of roundness, it diminishes. The second extra peak is the far right one, and we believe this
peak is a consequence of the 3D nature of our geometry; as observed in Figure \ref{fig:ext_y_14}, this peak intensity
decreases as the third dimension of the geometry lengthens.
The quantity of interest in these findings is the wavelength of the resonance peaks, and our results do indeed match the findings.

The second replication case is from a paper by Ellis and co-workers \cite{ellis2016} studying the effect of aspect ratio on the excitation of high-order modes in localized surface phonon-polariton nanostructures (Figure S4 of their supplementary material and Figure 2 of their paper).
Figure \ref{fig:rep_FS4_ellis} shows the results for this replication, where the relative errors between our computations and theirs is always smaller than 2$\%$. The smallest error is for the case with $AR=4$ and we proceeded to study the results on Figure 2 of their work, corresponding to this same aspect ratio.
Their results in Figure 2a include both experiments and simulations with the commercial software COMSOL, so we sought to both validate \pygbe using their experimental results, and replicate their computational findings. 
Again, we lack access to the data behind the plots, and we had to digitize the curves by hand. 
The quantity plotted in the original figures is reflectance as a function of wavenumber, whereas we compute the extinction---on the figures, they show inverted peaks, where we show positive peaks.
We can compare the results, nevertheless, because the quantity of interest is the wavenumber position of the peaks.
Figure \ref{fig:pygbe_vs_exp_2a} shows the results of our simulations using \pygbe on an isolated pillar,
compared with the experimental results of Ellis et al.\ on an array of pillars.
The results with \pygbe  do not account for the effect of coupling among the pillars, which explains the noticeable discrepancy on the wave numbers at which the peaks occur. 
We proposed a correction, based on the results reported on Figure S4 of Ellis et al.\ for the $AR=4$ case, where the shift on the wave number due to coupling is 12.17 cm$^{-1}$ (difference between black and red curves for $AR=4$ on figure S4 of the supplementary materials). 
We subtracted this value to that obtained with our simulations and we show the comparison of our corrected results against those of Ellis et al.: the experimental data on Figure \ref{fig:val_2a}, and the
computational data on Figure \ref{fig:rep_2a}. 
When comparing with their experiments, after the correction was applied, we observe a good match of the wavenumber for the lower (and stronger) mode, as well as a good match for the 
third and fourth peaks. The wave number of the second peak, related to a longitudinal excitation (mode $L_{000}$ in Ellis et al.), presents a discrepancy that we believe is related to the fact that our 
pillar does not have a substrate underneath. The remaining (fifth) peak, also presents a discrepancy, but in this case we can not identify the reason.
We did not analyze the peaks out of the range 864--961 cm$^{-1}$, since Ellis et al. \ describe these peaks to be associated with zone-folded LO (longitudinal) phonons of 4H-SiC,
and outside the scope of their study.
After considering all these details, we can say that we have validated our solver \pygbe against the experimental results of Figure 2a, as well as replicated their computational results.

%CHALLENGES 

Throughout all these activities aiming to replicate previous results and validate our computational model, we faced multiple challenges, starting with the complexity of the system under study. In both our source papers, we have systems that we were unable to fully 
model using \pygbe. For the case of Rockstuhl et al., even though they used a boundary element method, it was a 2D model instead of 3D like the one implemented in \pygbe. The computational work presented in Ellis et al.\ used a volumetric 
formulation (finite element method), whereas \pygbe implements a boundary integral one. Both works computed (or measured) quantities different from the extinction cross-section (computed in \pygbe), but given the relationship between their quantities (scattering cross-section and reflection)  
and ours (extinction cross-section), via the wavelength (wavenumber) at which resonance happens, this was not a problem. To produce the end results for comparing against these works, however, we went through an 
exhaustive process of modeling, making assumptions, and even corrections. We were lacking any information regarding the solvers, discretization, meshes, etc., in the original papers. In the case of Ellis et al., we benefited from multiple communications with the authors, 
who made available the geometry and mesh of the pillar, as well as the dielectric data used as input for the simulations. Even though we appreciate this helpful interaction, replication studies should ideally not depend on communications with the original authors. 
In both cases, the secondary data presented in their results were not publicly available and we relied on a manually digitized version, obtained from the plots. This work is tedious and introduces uncertainties that would be avoided by releasing 
the secondary data.

Validation has been a hard goal to achieve for our solver since the field does not have benchmarks that are meant to be used for this purpose. Multiple times,
we encountered experimental results that we aspired to use for validation, but due to the insufficient reporting of the experimental conditions, the setup and the geometries involved, 
we have been unable to validate until this moment. Even though we were not modeling the same experimental setup, the communications with the authors of Ellis et al., made 
the validation possible.

After all the challenges faced to achieve code validation and replication of results, we can say that the process was far from  linear, 
and the complexity involved increases if the work to be replicated or used for validation was not conducted and published using reproducible practices, and made open at the time 
of publication. We conclude that reproducible practices are needed not only for the work to be reproduced, but also 
replicated and even used for validation studies, if applicable.
