SIMPLE 2 BIN SOLUTION FOR FREQUENCY AND PHASE - BOTTOM LINE
The following is the summary for the 2 bin solution to determine the exact frequency and phase given the largest DFT bin and the next largest adjacent bin for the case of a single tone. Details and resulting accuracy are further outlined below:
Step 1: Select the peak bin and next largest adjacent bin. Compute ratio of the magnitude of the two bins, the smaller bin in the numerator:
$|Y[k_{max}]|$ is the magnitude of the peak bin
$Y|[k_{max}+1]|$ is the magnitude of the next highest adjacent bin. here as abbreviated, the assumption is the bin to the right, but could also be $Y|[k_{max}-1]|$
$$r = \frac{Y[k_{max}+1]|}{|Y[k_{max}]|}$$
Step 2: Use the result $r$ with equation \ref{4} (further in the post and copied below) to compute the bin offset $k_\Delta$, which is the fractional frequency offset from bin center for the tone
$$k_\Delta \approx 0.63883277 - 0.6368504 e^{-1.475173r} $$
$k_\Delta$ will be a value between $0$ and $0.5$ corresponding to the maximum possible error between bins. If the next highest bin was to the right of the peak, $k_\Delta$ is to be interpreted as a positive quantity. However if the next highest bin was to the left of the peak, $k_\Delta$ as computed above should be multiplied by $-1$ to make it a negative quantity.
Step 3: Use the result $k_\Delta$ with equation \ref{2} (further in the post and copied below) to compute the phase correction due to the fractional frequency offset. The resulting $\phi$ is the correction to subtract from the phase for the peak bin to get the phase of the actual tone.
$$\phi = \pi k_\Delta \frac{N-1}{N}$$
DETAILS
Another solution is to simply zero pad the waveforms prior to taking the FFT. This will serve to interpolate more samples in the frequency domain from which the phase of the peak (or improved under noise conditions the average of the phase of the $M$ samples within the peak, where $M$ is a trade based on signal to noise of each of the DFT bins contributing to the answer). If we double the number of samples in the FFT by zero padding (using fft(x, M) where $M=2N$ and $N$ is the original length of the FFT), we will get one more sample in between each of the original samples. If we zero pad out to 5 times the original length, we will get 4 more samples in between each of the original samples, etc.
This isn't an ideal approach to compare the delay of the two signals. Alternatively consider a least squares based approach such as detailed here, which uses every sample optimally toward an estimate of the delay.
If it known that a single tone or narrow band signals with sufficient SNR will be used, delay can be estimated based on carrier phase and for that a simple phase detector approach can be used (multiply the two signals in the time domain and low pass filter). Further below I detail the development of an approach to determine the phase correction based on the DFT bin values.
Prediction of Phase Versus DFT Bin Offset For a Single Tone
The OP has asked for a mathematical relationship of the derived phase versus frequency offsets from bin center as is occurring. Assuming the SNR is high enough, the relationship between magnitude and phase for the samples in vicinity of a DFT bin for the DFT using a simple
Rectangular window (no additional windowing) would be dominated by the response given by the Dirichlet Kernel with frequency in units of DFT index as:
$$D(k) = e^{-j\pi (k-k_o)(N-1)/N}\frac{\sin(\pi (k-k_o))}{sin(\pi (k-k_o)/N)} \tag{1} \label{1}$$
Given an $N$ sample DFT, with $k$ as a continuous frequency variable ranging with $k \in [0, N)$ and $k_o$ as the discrete frequency index for a specific bin with $k_o \in 0 \ldots N-1$.
From this, and given a phase offset as the fractional portion of $k_o$, given as $k_\Delta = mod(k_0,1)$ we can get the relationship between bin offset and phase as:
$$\phi = \pi k_\Delta \frac{N-1}{N} \tag{2} \label{2}$$
The above provides the frequency response both in magnitude and phase for any DFT bin $k_o$ given an actual input frequency anywhere between $[0, N)$ on the continuous frequency domain. (With $N$ corresponds to the frequency of the sampling clock). Any additional windowing used would modify this response including spreading out the magnitude of the mainlobe and reducing the sidelobes.
For example the magnitude frequency response for the first bin ($k_o=0$ and sixth bin ($k_o=5$) for a 12 point DFT is plotted below. This shows the predicted magnitude for each of these bins for any single tone (a single tone is of the form $e^{j\omega t}$, not $\cos(\omega t)$!) at any frequency within the given range. Keep in mind that a sinusoid consists of two such tones ($2\cos(\omega t) = e^{j\omega t}+ e^{-j\omega t}$), so the cross interference of the response from each tone in the case of a sinusoid can have a significant effect on the computation for the case of lower and higher frequencies, or low number of total DFT bins.
The magnitude only is shown in the plot above, but significantly there is a considerable phase change as well versus variable $k$ as we move off of bin center toward the adjacent bins (where the magnitude goes to a null). Below is a plot for the result of bin 0 alone on a complex plane, where we can see the variation of both magnitude and phase concurrently:
Note in the unwindowed DFT without zero padding, the magnitude shown above for a single tone is going to the origin within the spacing of one DFT bin. The plot below provides further insight into the linear phase trajectory of the Dirichlet Kernel between bins and validates the formula given above for $D(k)$: The blue stem plot in magnitude and blue dots for phase are the actual results for a 12 point DFT with a single tone at $k=5.1$ as $y = e^{j2\pi(5.1)n/N}$ for $n = 0 \ldots 11$. Notably, the phase from bin center to the next bin is on a slope of $\pi N/(N-1)$, so goes nearly, but not quite $\pi$ radians in the spacing of one bin. With that, if we are able to predict (or estimate) the offset from bin center, we can then determine the phase correction.
Estimating $k_\Delta$
Under reasonable SNR conditions and for a single complex tone the offset in terms of a fraction of a bin can be estimated from the ratio of the two largest adjacent peaks, and then from that the phase correction in the linear range extending approximately from $-\pi/2$ to $+\pi/2$ can be determined. (A real tone can also be used with an insignificant error is not near DC or Nyquist. The real tone consists of two complex tones whose responses will overlap and result in an error more dominant for real tones close to DC or Nyquist. I haven't tested this case but suspect the result would approach that for $N=2$, which shown below also has a very low error.)
The ratio is chosen to be the smaller of the two peaks over the larger. Consider when we are on exact bin center; either adjacent sample will be at zero, so the ratio metric is better as the smaller value divided by the larger to avoid the divide by zero, with a ratio then going from 0 to 1 and mapping to $k_\Delta$. For example in the plot above, the next highest peak is at bin 6, and this magnitude would be divided by the magnitude of the peak at bin 5, and from that ratio (as detailed below) we can determine $k_\Delta$.
Below shows an example of such an approximator based on a simple curve fit that has been shown to provide reasonable accuracy for most applications. The vertical axis shows the actual fractional bin offset $k_\Delta$ from 0 to 0.5, and the horizontal axis is the inverse ratio of the magnitude of the peak bin to the magnitude of the next highest adjacent bin.
The exact solution is derived from equation \ref{1} as:
$$\frac{|Y[k_{max}+1]|}{|Y[k_{max}]|} = |D(k_\Delta-1)|/|D(k_\Delta)| \label{3} \tag{3}$$
Where
$|Y[k_{max}]|$ is the magnitude of the peak bin
$Y|[k_{max}+1]|$ is the magnitude of the next highest adjacent bin (here assuming the bin to the right resulting in a positive $k_\Delta$, or if the bin to the left as $Y[k_{max}-1]$ resulting in a negative $k_\Delta$.
Solving \ref{3} for $k_\Delta$ would be quite complicated, requiring an inverse Dirichlet Kernel. I curve fit the graph shown, which was specific to $N=12$, directly resulting in the following exponential fit, which was also included in the above plot:
$$k_\Delta \approx 0.63883277 - 0.6368504 e^{-1.475173r} \tag{4} \label{4}$$
Where
$r = \frac{|Y[k_{max}+1]|}{|Y[k_{max}]|}$
Thus the two bin solution is to determine $k_\Delta$ using equation \ref{4} and with that determine the phase correction using \ref{2}. To note, $\ref{4}$ is a fit specific to $N=12$, where $N$ is the size of the FFT. A new and more accurate fit could be determined for each $N$, however I tabulated the errors for all $N$ values up to $N=10,000$ using fit given in \ref{4} showing the error quickly going be $1E-5$ and converging to $3E-7$ for $N>6000$.
Zooming in on the range of $N<100$ shows a maximum error for $N=2$ to be $4E-2$.
Verification
I tested the algorithm using the waveform example plotted above, which was generated with a single tone with zero phase at $k=5.1$. The largest bin is at $k=5$ and the next largest is at $k=6$. The magnitude of each were as follows:
$|Y[5]| = 0.98374$
$|Y[6]| = 0.11031$
The angle of $Y[5]$ which is the error from true 0 phase is $0.28798$ radians.
The ratio $|Y[6]|/|Y[5]| = 0.11213$
The predicted $k_\Delta$ using \ref{4} is within 1.2% of the actual value of 0.1:
$k_\Delta =0.0988$
Using this estimate for $k_\Delta$ in equation {2} results in the predicted phase correction which is also within 1.2% of the true phase offset of $0.28798$ radians:
$\phi = \pi k_\Delta(N-1)/N = 0.28446$ radians.
The correction therefore is to subtract $0.28446$ radians from the angle for bin 5, resulting in a phase of:
$0.28798 - 0.28446 = 0.0035$ radians
Note that an rms angle of $0.0035$ is consistent with an SNR of $20Log_{10}(0.0035) = 49$ dB, suggesting the estimator as derived here to usable for a wide range of applications within that SNR limit. As noted earlier, the exponential curve fit given in equation \ref{4} was fitted for $N=12$. Other fits would be similar with simply a modification of the coefficients for equation $\ref{4}$ approaching an asymptotic limit as $N$ grows larger, but I haven't detailed that further.
