Confirmation on how to calculate phase differences at every frequency point between two wideband signals?
What you seem to want is the phase response of the device that your optical system is passing through. Based on that, your high-level approach makes intuitive sense overall; I'm not sure why you chose a Hamming window and applied it three times. In any event, though, if you have a linear discrete-time system with (complex-valued) frequency response $H(e^{j\omega})$, then:

$$
Y(e^{j\omega}) = H(e^{j\omega})X(e^{j\omega})
$$

where $X(e^{j\omega})$ and $Y(e^{j\omega})$ are the discrete-time Fourier transforms (DTFTs) of the input and output signals, respectively. If you express the above in polar format (magnitude and phase), then the phase response of the system is just:

$$
\angle H(e^{j\omega}) = \angle Y(e^{j\omega}) - \angle X(e^{j\omega})
$$

that is, the phase response of the system just expresses the amount of phase shift that the system applies to the input signal as a function of (normalized) frequency $\omega$. Your approach takes advantage of the fact that for two complex numbers $X(e^{j\omega})$ and $Y(e^{j\omega})$:

$$
\angle Y(e^{j\omega}) - \angle X(e^{j\omega}) = \angle \left(X^*(e^{j\omega})Y(e^{j\omega})\right)
$$

So, given the DTFTs of the system's input and output signals, you can straightforwardly calculate the phase response that you seek. The devil is in the details of how one would obtain the DTFT representations of these signals. The nice part about the DTFT is that it is continuous in frequency; the bad part is that the DTFT sum is infinite in length. So, if you really want to evaluate the phase response at arbitrarily many frequencies, you better be patient, because you'll need to collect an infinite amount of data first.

The discrete Fourier transform (DFT) is often used as an approximation to the DTFT, which is valid under certain conditions. If the signals $x[n]$ and $y[n]$ are finite in length (which for a practical scenario like yours, they are), then the DFT can be thought of as evaluation of the DTFT on the finitely-spaced frequency grid $\frac{2\pi k}{N}, k = 0, 1, \ldots , N-1$ (where $N$ is the signal length). So, the larger you make $N$, the smaller your sample spacing across the frequency axis. 

Most physical devices will have a frequency response that is smooth to some degree (taken with the caveat that I'm not knowledgeable on optical device physics), so you should be able to find a DFT length that gives a decent approximation to the DTFTs for the input and output signals of your system. You can then apply your approach to estimate the phase response of the device.



You noted in your question that you see incorrect results for the estimated phase response for DFT bins in the output signal that didn't contain any content in the input signal; this is to be expected. Your assumption that you would just get zero phase for those bins is not likely to hold in a practical scenario; noise present in the sampled signals, spectral leakage, or other interference can cause you to observe some small amount of energy at frequencies that you wouldn't expect it. These extraneous sources can introduce distortion between the input and output that is not due to the device's response, and for complex numbers with small magnitude, small amounts of additive noise can drastically affect the number's resulting phase. 

Therefore, you can only consider the phase measurement valid at frequencies where there is significant spectral content in the input. Stated differently, if you want to make a precise phase measurement, you want high signal-to-noise ratio in the complex value that you're using to estimate the phase.