The Hilbert Transform idea would definitely work but may be overly complicated. An FFT would also provide phase and amplitude information for each bin, and the phase of the dominant bins where the tone resides can be compared to determine relative phase of each, assuming the data captures are synchronized.
In this case with a reasonably small phase another simple approach would be to multiply each waveform with a reference tone offset by a quarter cycle (90°) and filter. The waveforms which should be synchronously captured should be leveled to the same amplitude prior to taking the product since the phase result is both sensitive to amplitude and actual phase. Given two sinusoids, the average of the product is proportional to the cosine of the phase difference between the two sinusoids, by delaying one of the two by 90 degrees, the result will be proportional to the sine, and for small angles $\sin(\theta) \approx \theta$ when the angle is given in radians. For example an angle as large at 10 degrees which is 0.1745 radians has the result $\sin(0.1745) = 0.1736$ with an error that would most likely be below the noise of the result.
This is seen in the trigonometric relationship of the product of two sinusoids:
$$\cos(2\pi f t) \cos(2\pi f t + \theta) =\frac{1}{2} \cos(4\pi f t + \theta) + \frac{1}{2}\cos(\theta)$$
We see the product is the sum and the difference of the two frequencies and phase offsets. When the frequencies are the same we get a tone at double the frequency (which we filter out in the average, or low pass filter) and a tone at DC which is proportional to the cosine of the angle.
Building on this and the need to level each tone prior to taking the product (to remove amplitude sensitivity), a best approach would be to hard limit each received tone, ensuring the threshold is properly set for a 50% duty cycle in each resulting square wave. The low frequency result of this product (which is a simple XOR function that is averaged) will have a linear result versus phase, providing an unambiguous linear output versus phase over a $\pm 90$ degree phase range. Again, the reference tone should be delayed by 90° to center the result in the usable range.
The product and average solution has the benefit of minimizing noise in estimate versus any approach that would for example compare just one zero crossing, it is important however that the reference waveform be synchronized to the test waveform to eliminate inevitable drift errors between the two.
