Let's work out the Fourier transform. If I have a complex exponential in the time domain, as you noted, after this is
\begin{equation}x[n] = e^{j2\pi \frac{k_{0}n}{N}}\end{equation}
The Fourier transform of this complex exponential is
\begin{equation}X[k] = \sum_{n=0}^{N-1}e^{j2\pi \frac{k_{0}n}{N}}e^{-j2\pi \frac{kn}{N}} = \sum_{n=0}^{N-1}e^{j2\pi \frac{(k_{0}-k)n}{N}} \propto \delta(\frac{k_{0}-k}{N})\end{equation}
Now, let's apply a phase shift to the time domain signal. We get
\begin{equation}x'[n] = e^{j(2\pi \frac{k_{0}n}{N} + \phi)} = e^{j\phi}e^{j2\pi \frac{k_{0}n}{N}}\end{equation}
The DFT of this is
\begin{equation}X'[k] = \sum_{n=0}^{N-1}e^{j\phi}e^{j2\pi \frac{k_{0}n}{N}}e^{-j2\pi \frac{kn}{N}}\end{equation}
Because $e^{j\phi}$ is constant, we get
\begin{equation}X'[k] = e^{j\phi}\sum_{n=0}^{N-1}e^{j2\pi \frac{k_{0}n}{N}}e^{-j2\pi \frac{kn}{N}} \propto e^{j\phi}\delta(\frac{k_{0}-k}{N})\end{equation}
So, we get that the phase in the Fourier domain is the same as the phase offset in the time domain.
EDIT:
As pointed out in the comments, I wrote the DTFT on the left hand side and the DFT on the right. This has been fixed.
This is also important, as noted in the comments, for figuring out the phase in the Fourier domain. The phase of the DTFT, and the phase of the DFT in which the signal is periodic within the DFT window, is equal to the phase offset in the time domain. For a signal that is aperiodic within the DFT window, there will be an unknown phase offset due to the assumed frequency content by the aperiodicity. However, the linearity property still applies and will produce an equal phase offset. This can be demonstrated by the Matlab code below:
fs = 10;
f = 1;
N = fs/f;
T = 1/f;
t1 = (0:N-1)*T;
t2 = (0:N+1)*T;
c1 = exp(1i*2*pi*f/fs*t1);
c2 = exp(1i*2*pi*f/fs*t2);
c3 = exp(1i*(2*pi*f/fs*t1 + pi/2));
c4 = exp(1i*(2*pi*f/fs*t2 + pi/2));
Taking the angles of each of the 4 complex exponentials, one will see that when the signal is periodic within the DFT window (c1 and c3), the phase is exactly equal to the phase offset in the time domain (at the frequency of interest, excluding rounding errors in the other values). Taking the angle of c2 and c4, there will be an unknown phase offset applied to the FFT sample of interest, but the phase difference between the two is proportional to the phase difference between the signals in the time domain.
Hope this clarifies some things.
EDIT 2:
Adding the DTFT math per OP’s request.
If I have a complex exponential in the time domain, as you noted, it is defined as
\begin{equation}x[n] = e^{j\omega_{0}n}\end{equation}
The Fourier transform of this complex exponential is
\begin{equation}X(\omega)= \sum_{n=-\infty}^{\infty}e^{j\omega_{0}n}e^{-j\omega n} = \sum_{n=-\infty}^{\infty}e^{j(\omega_{0}-\omega)n} \propto \delta(\omega_{0}-\omega)\end{equation}
Now, let's apply a phase shift to the time domain signal. We get
\begin{equation}x'[n] = e^{j(\omega_{0}n+ \phi)} = e^{j\phi}e^{j\omega_{0}n}\end{equation}
The DFT of this is
\begin{equation}X'(\omega) = \sum_{n=-\infty}^{\infty}e^{j\phi}e^{j\omega_{0}n}e^{-j\omega n}\end{equation}
Because $e^{j\phi}$ is constant, we get
\begin{equation}X'(\omega) = e^{j\phi}\sum_{n=-\infty}^{\infty}e^{j\omega_{0}n}e^{-j\omega n} \propto e^{j\phi}\delta(\omega_{0}-\omega)\end{equation}
As you correctly noted, since it is a phase shifted dirac delta function, the phase is zero everywhere except the frequency of the complex exponential. This holds for the DTFT, and the DFT if the signal is periodic within the DFT window.