This is a theoretical question without much practical interest but still it can be nice to check the results and investigate if intuition still holds.
First of all, if the convolution is regarded an LTI operation between an input $x(t)=\cos(\omega_0 t)$ and a system $y(t)=h(t)=\sin(\omega_o t)$ then it's immediately obvious that since $h(t)$ is an unstable system, the output can be unbounded even for a bounded input signal.
Furthermore assuming that Fourier transforms of the input, the system and the output exists (the integrals converge!) then we can assert the property that: $$ x(t) \star y(t) \longleftrightarrow X(\omega)Y(\omega) $$
When the sinusodal signals $x(t)=\cos(\omega_0 t)$ and $y(t)=\sin(\omega_0 t)$ are considered, it's obvious that their Fourier integrals do not converge. Hence their formal Fourier transforms do not exist. The solution is the acceptance of the generalised impulse function to represent the Fourier transform of the sinusodial signals; i.e., $$ x(t)=\cos(\omega_0 t) \longleftrightarrow X(\omega) = \pi [\delta(\omega-\omega_0) + \delta(\omega+\omega_0)]$$ and $$ y(t)=\sin(\omega_0 t) \longleftrightarrow Y(\omega) = \frac{\pi}{j} [\delta(\omega-\omega_0) - \delta(\omega+\omega_0)]$$
Then we consider that the above theorem still holds for the impulse functions as well and apply it here:
$$ z(t)=\sin(\omega_0 t) \star \cos(\omega_0 t) \longleftrightarrow Z(\omega)=\frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)] \cdot \pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)] $$
Using the impulse sifting property;
$$f(x)\delta(x-a) = f(a)\delta(x-a)$$
(NOTE-WARNING: sifting property of $\delta(x-a)$ strictly requires that the function $f(x)=\delta(x)$ be sufficiently smooth around the discontinuity implied by the sifting function $\delta(x-a)$. Henceforth, in this application the property cannot be applied, as assuming that $\delta(x)$ is a sufficiently smooth function would then contradict the use of $\delta(x-a)$ as a sifting function in its own sifting property, and thus the rest of this manipulations is just an algebraically consistent illusion from mathematical point of view...)
we shall perform th multiplications as follows:
$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(\omega-\omega_0) \delta(\omega-\omega_0) + \delta(\omega-\omega_0) \delta(\omega+\omega_0) - \delta(\omega+\omega_0) \delta(\omega-\omega_0) - \delta(\omega+\omega_0) \delta(\omega+\omega_0)] $$
$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(0) \delta(\omega-\omega_0) + \delta(-2\omega_0) \delta(\omega+\omega_0) - \delta(2\omega_0) \delta(\omega-\omega_0) - \delta(0) \delta(\omega-\omega_0)] $$
Now noting that $\delta(2\omega_0)=\delta(-2\omega_0)= 0$ those terms cancel and $\delta(0)=\infty$ terms remain, hence
$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(0) \delta(\omega-\omega_0) - \delta(0) \delta(\omega-\omega_0)] $$
$$ Z(\omega)= \pi \delta(0) \frac{\pi}{j} [ \delta(\omega-\omega_0) - \delta(\omega-\omega_0)] $$
Which is recognized as the Fourier transform of the infinite amplitude sine wave as: $$ \boxed{ z(t) = \pi \delta(0) \sin(\omega_0 t) } $$
The conclusion is that the convolution between $x(t)=\sin(\omega_0 t)$ and $y(t)=\sin(\omega_0 t)$ produces an infinite amplitude sinudoidal wave of the same frequency $\omega_0$.
The time domain verification is as follows:
$$x(t) \star y(t) = \int_{-\infty}^{\infty} x(t-\tau)y(\tau) d\tau \leftrightarrow z(t) = \int_{-\infty}^{\infty} \cos(\omega_0(t-\tau)) \sin(\omega_0 \tau) d\tau$$
Using the trigonometric identity of $$\cos(x)\sin(y) = 0.5[\sin(y+x) + \sin(y-x)]$$ we can break the integral into two, where $x = \omega_0(t-\tau)$ and $y=\omega_0\tau$ , hence
$$z(t) = 0.5 \int_{-\infty}^{\infty} \sin(\omega_0(t-\tau)+\omega_0 \tau) + \sin(\omega_0 \tau - \omega_0(t-\tau) ) d\tau$$
$$z(t) = 0.5 \int_{-\infty}^{\infty} \sin(\omega_0 t)d\tau + 0.5 \int_{-\infty}^{\infty} \sin(2\omega_0 \tau - \omega_0 t) d\tau$$
Now the first integral becomes $$ 0.5 \sin(\omega_0 t) \int_{-\infty}^{\infty} d\tau$$ while the second integral can be shown to be zero after a suitable change of variables, assuming for a given (fixed) t, $\phi = 2\omega_0 \tau - \omega_0 t$ and $d\phi = 2\omega_0 d\tau$ then the second integral becomes $ \frac{1}{2\omega_0} \int_{-\infty}^{\infty} \sin(\phi) d\phi = 0$.
Theferore the result of the convolution is:
$$z(t) = \left( 0.5 \int_{-\infty}^{\infty} 1 d\tau \right)\sin(\omega_0 t) $$
Note that the integral that weights the sine wave has infinte value, moreover it can be shown that by the forward and inverse Fourier transform pairs one can recognize the integral as the forward Fouier transform $H(\omega)$ (evaluated at $w=0$, $H(0)$) of the signal $x(t)=1$ which is:
$$\mathcal{F} \{ 1\} = \int_{-\infty}^{\infty} 1 e^{-j \omega t} dt \equiv 2\pi \delta(\omega) $$ and therefore setting $w=0$ yields $$\int_{-\infty}^{\infty} 1 dt \equiv 2\pi \delta(0) $$
Finally plugging this into the result yields; $$\boxed{ z(t) = \pi \delta(0) \sin(\omega_0 t) }$$ which is the same as the result obtained from the Fourier method.
