I found a better condition using Polya's characterization of non-zeroness of Fourier transform.
Theorem (Polya's, Thm. C of this paper). If $f(x)$ is a positive, continuous and decreasing 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). 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}$.
$\blacksquare$
Building on this, we have the following proposition that answers my question:
Proposition-2. Suppose that $X$ follows a unimodal distribution with a continuous density function $f_X$. Further suppose that $g$ is a continuous, increasing function with $g^\prime$ decreasing on $(-\infty,a)$, increasing on $(a,\infty)$, where $a$ is the mode of $X$. Then the random variable $Y:=g(X)$ is also unimodal and hence has a non-zero characteristic function.
Remark. An example for such $g$ could be the sigmoid function centered at $a$.
Proof. We have the density function $f_Y$ of $Y$ being
$$f_Y(y) = f_X(g^{-1}(y)) \cdot \frac{d}{dy} g^{-1}(y)\\ = f_X(g^{-1}(y)) \cdot \frac{1}{g^\prime(g^{-1}(y))}.$$
Since $g$ is an increasing function, so is $g^{-1}$. As a result, $f_Y(y)$ is increasing for $y\in (-\infty, g(a))$ and decreasing for $y\in (g(a),\infty)$.
$\blacksquare$
A cleaner condition is given by the following corollary:
Corollary-3. Suppose $X$ follows a unimodal distribution with mode zero and density function continuous. If $g$ is a S-shaped function satisfying $g(0)=0$. Then $Y:=g(X)$ has a non-zero characteristic function.
A function $g$ is called an S-shaped function if $g$ is continuous, increasing, and that $\exists a$ such that $g^\prime$ increases on $(-\infty,a)$ and decreases on $(a,\infty)$.
Remark. The strength of the above corollary is that it is closed under summation and function compostion. That is, for two S-shaped functions $g_1, g_2$, we have both $g_1(g_2(X))$ and $g_1(X)+g_2(X)$ have non-zero characteristic functions.
