Why does subbing $s = j\omega$ into the Laplace transform of a cosine wave yield a purely imaginary result?
If I sub in $s=j\omega$, I get the Fourier transform of a cosine wave starting at $t=0$.


No, you don't. You can't just set $s=j\omega$ in an expression for the Laplace transform and expect the result to be the Fourier transform of the function. This only works if the imaginary axis is inside the region of convergence (ROC) of the Laplace transform. Since the given Laplace transform has poles on the imaginary axis, the imaginary axis is not part of the ROC (the ROC is $\textrm{Re}\{s\}&gt;0$).

The Fourier transform of the function

$$x(t)=\cos(\omega_0t)u(t)\tag{1}$$

exists in a distributional sense, i.e., if we allow Dirac impulses (and, possibly, its derivatives) in the expression of the Fourier transform.

Using the Fourier transform of the step function

$$\mathcal{F}\{u(t)\}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{2}$$

and the modulation property of the Fourier transform

$$\mathcal{F}\left\{e^{j\omega_0t}u(t)\right\}=\pi\delta(\omega-\omega_0)+\frac{1}{j(\omega-\omega_0)}=Y(\omega)\tag{3}$$

we can easily derive the Fourier transform of $(1)$ by noticing that $x(t)$ is the real part of the time-domain function on the left-hand side of  $(3)$, so its Fourier transform must be the even part of $Y(\omega)$:

$$\begin{align}X(\omega)&amp;=\frac12\left[Y(\omega)+Y^*(-\omega)\right]\\&amp;=\frac{\pi}{2}\left[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\right]+\frac{j\omega}{\omega_0^2-\omega^2}\tag{4}\end{align}$$