Sampling with an alternating phase offset
Bottom Line
Given a sampled tone at normalized radian frequency $\omega_1$ and real amplitude $A$:
$$x[n] = A e^{j\omega_1n} \label{1} \tag{1}$$
If every other sample is shifted by a constant phase $\phi$, the result in the Fourier Transform will be two tones $x_1$ and $x_2$ (precisely so as the only non-zero values when the sampling rate is an integer multiple of $\omega_1$, otherwise $x_1$ and $x_2$ refer to the maximum bin in vicinity of the two peaks):
$$x_1 = |x_1|e^{j\phi_1}$$
$$x_2 = |x_2|e^{j\phi_2}$$
And $\phi$ is determined using the magnitude and phase of each of the two bins as follows:
$$\phi = 2\sin(\phi_1-\phi_2) \text{atan2}(|x_2|, |x_1|)$$

Further Details
The first tone will be at the bin associated with the same normalized radian frequency $\omega_1$ as the original but shifted in phase by $\phi/2$ and an additional tone will be given due to the phase modulation at the bin associated with the normalized radian frequency $\pi + \omega_1$ and initial phase such that the new tone will always be in quadrature phase with $x[n]$ as $-j (-1)^n x[n]$. The amplitude of the original tone will be reduced according to $A\cos(\phi/2)$ while the amplitude of the additional tone will be $A\sin(\phi/2)$. Given the modified original tone as $x_1[n]$, and the new tone as $x_2[n]$ this relationship to $x[n]$ ahd $\phi$ is:
$$x_1[n] = x[n]\cos(\phi/2)e^{j\phi/2}= A\cos(\phi/2) e^{j(\omega_1n + \phi/2)}\label{2} \tag{2}$$
$$x_2[n] = -j (-1)^nx[n]\sin(\phi/2)e^{j\phi/2} = jA\sin(\phi/2)e^{j((\pi+\omega_1)n+ \phi/2)}\label{3} \tag{3}$$
So given this and not knowing $\phi$, we can determine $\phi$ from the DFT results as follows: Given $x_1$ and $x_2$ as the complex values for two bins where $x_1$ is the frequency bin for the expected tone and $x_2$ is the frequency bin for the interference tone due to the introduced distortion of an unknown fixed phase shift every other sample, with:
$$x_1 = |x_1|e^{j\phi_1}$$
$$x_2 = |x_2|e^{j\phi_2}$$
$x_1$ and $x_2$ will always be in quadrature, and the phase between the $x_1$ and $x_1+x_2$ will always be $\phi/2$. Further, when $x[n]$ has no phase offset as we had introduced in \ref{1}, then the sum  of $x_1+x_2$ will also always be real and the phase can be determined from twice the phase angle of $x_1$ alone. However, what is developed here is the more general case applicable to any starting phase offset for $x[n]$. This general case can be solved for using the following:
$$\phi = 2 \times \text{atan2}(|x_2|\sin(\phi_1-\phi_2), |x_1|)$$
Where $\text{atan2}$ is the 2-element arctangent. The term $\sin(\phi_1-\phi_2)$ will always be either +1 or -1 depending if $\phi$ is positive or negative, since the relationship of the two tones will always be in positive or negative quadrature, so the above can also be written as:
$$\phi = 2\sin(\phi_1-\phi_2) \text{atan2}(|x_2|, |x_1|)$$
To extend this to a real waveform such as in the case of the OP, the original tone from the positive frequency axis would be used as $x_1$, and it's negative image would be used as $x_2$ as depicted in the graphic below:

Using a real waveform that is exactly at normalized frequency $\omega_1 = \pi/2$ would not be resolvable since the image from the added phase distortion and the primary negative frequency  component would both be on the same frequency location.

Derivation:
Relating \ref{1},  \ref{2} and \ref{3}:
$$x[n] = x_1[n] + x_2[n]$$
$$ = x[n]\cos(\phi/2)e^{j\phi/2} -j (-1)^nx[n]\sin(\phi/2)e^{j\phi/2}$$
$$ = x[n]e^{j\phi/2}(\cos(\phi/2) - j(-1)^n\sin(\phi/2))$$
$$ = x[n]e^{j\phi/2}e^{-j\phi/2(-1)^n}$$
$$ = x[n] \text{ for n even}$$
$$ = x[n]e^{j\phi} \text{ for n odd}$$
Intuition:
The above relationships are more intuitive to understand knowing that the Fourier Transform decomposes any single valued time domain analytic function into distinct phasors (as the basis functions) that each must rotate with time with a constant magnitude and at a constant rate. We can see the solution for how this would occur for the OP's problem by doing the decomposition geometrically.
First consider an arbitrary $x[n]$ as a rotating phasor on the complex plane in the time domain as shown on the left, where when sampled the normalized radian frequency corresponds to the phase shift from sample to sample. The Fourier Transform of this on the right is a single impulse in frequency, a positive frequency here given the counter-clockwise rotation of the phasor (positive increasing angle). The single impulse corresponds directly to the single spinning phasor.
x[n] with no distortion:

Now consider if every other sample was shifted forward in phase by a constant amount $\phi$:
x[n] with phase distortion:

To decompose this into phasors of constant magnitude and rate of rotation, the first phasor MUST be the first phasor rotated by $\phi/2$ since the error represents a modulation going back and forth by $\phi/2$ for every sample. And for this same reason the second phasor must be in alternating quadrature relationship such that it is modulating the phase back and forth from sample to sample (resulting in a normalize angular frequency of $\omega_1 - \pi$ when we consider it rotating opposite the original tone, or similarly within the range of $0$ to $2\pi$ at the normalized angular frequency of $\pi+\omega_1$.)
Fourier Decomposition of distorted x[n] into two tones:


And if we look at the relationship between any two samples at a given time instance in the decomposed waveform of two phasors, we see through the geometric relationship how the two must be always in quadrature with the cosine and sine relationship of magnitudes such that the normalized magnitude of the sum of the two phasors will be normalized (to one here) as in the original waveform. Consider the first two samples of the original waveform with the second sample offset by an additional $\phi$ as shown by the red phasors below, together with the first two samples of the decomposed tones showing how the quadrature relationship and relative magnitudes must be maintained:
First Sample of distorted x[n] and the Decomposed Tones:

Second Sample of distorted x[n] and the Decomposed Tones:

In summary, we shift the original tone by $\phi/2$ so that each sample it is midway between the sum of the two decomposed phasors. Then in order for the phase to jump back and forth by exactly $\phi/2$ from sample to sample, the second tone must be in alternate quadrature relationship from sample to sample (so will have rotated $-\pi + \omega_1$ since the original phasor modified with a phase shift of $\phi/2$ still rotates at the same rate as $+\omega_1$). Further for the sum of the two decomposed phasors to be on the unit circle (or whatever magnitude the original tone was), the magnitudes of each are related by the $\cos(\phi/2)$ and $\sin(\phi/2)$ since $\cos(\phi/2) + j\sin(\phi/2) = 1$.
This shows the relationship between the Fourier Transform components in the general case with a single complex tone, and taking the real of the above would result in two symmetric tones with the same relationships as is the case for the OP.  Further, to be precise, the descriptions above are specific to the Discrete Time Fourier Transform (DTFT) in which the time domain waveform is extended for all time and would be a continuous function in frequency with the exact solution shown as 2 weighted impulses in frequency corresponding to the two phasors in the time domain and zero everywhere else. To extend this derivation approach to the DFT, x[n] would be truncated in time (windowed with a rectangular window of duration $T$) which would result in convolving the impulses in frequency with a Dirichlet Kernel (aliased Sinc function) and then sampling the resulting frequency response at a spacing given by $1/T$ where $T$ is the window duration. When the original tone is on a bin center of the DFT (multiple of $1/T$), then the DFT would also result in just two tones for this case. For all other cases, multiple smaller tones would also appear (spectral leakage) but the same relationship would apply between the strongest tone closest to $\omega_1$ and the second tone at $\pi+\omega_1$, where the $N$ frequency bins given as integers $k$ from $0$ to $N-1$ are related to the normalized radian frequency $\omega$ as $\omega = 2 \pi k/ N$ .