I found a better condition using Polya's characterization of non-zeroness of Fourier transform.

> **Theorem** (Polya's, Thm. C of [this paper](https://arxiv.org/abs/2110.01885)). If $f(x)$ is a positive, continuous and descreaing function on $x\in(0,\infty)$ such that the integral $\int_0^\infty f(x)dt$ exists, then
$$\int_0^\infty f(x)sin(kx)dt>0, \quad \forall k>0$$

We then have the following result:

**Proposition-1.** Any unimodal distribution with a continuous density function has its characteristic function (c.f.) non-zero.

**Proof.** Consider a random variable $X$ with density function $f$. We denote $F$ as the Fourier transform of $f$ (which is also the c.f. of $X$). We only show the case where $X$ has a unimodal distribution with mode zero, since translation of the density function does not affect the non-zeroness of $F$ (the Fourier transform of $f(x+a)$ is $e^{ika} F(k)$, if $F(k)$ is non-zero everywhere, then so is $e^{ika} F(k)$).

Due to the assumed unimodality, we have $f$ is increasing on $(-\infty,0)$ and decreasing on $(0,\infty)$. Hence, the function $h(x):=f(x)+f(-x)$ is decreasing on $(0,\infty)$. It is also continuous, positive and integrable. Using Polya's Theorem, we then have
$$\int_{-\infty}^\infty f(x) sin(kx)dx=\int_0^\infty f(x) sin(kx)dx + \int_{-\infty}^0 f(x) sin(kx)dx \\=\int_0^\infty f(x) sin(kx)dx + \int_0^\infty f(-x) sin(-kx)d(-x) \\=\int_0^\infty \{f(x)+f(-x)\} sin(kx)dx=\int_0^\infty h(x) sin(kx) dx>0,\quad \forall k>0$$

Till here, we have shown that the imaginary part of $F(k)$ is non-zero for all $k\in (0,\infty)$. We also have $F(0)=1$ (a property for all c.f.s, see Prop. 7 of [this note](https://www.stat.cmu.edu/~arinaldo/Teaching/36752/S18/Notes/lec_notes_11.pdf)). Moreover, we have from the symmetry of Fourier transform (for a real function $f$, its Fourier transform $F$ satisfies $F(k)=F^*(-k)$), the imaginary part of $F(k)$ is also non-zero for $k\in (-\infty,0)$. Together, these mean $F(k)\neq 0$ for all $k\in\mathbb{R}$.