The definition of $f^*$ is given to me as below. But what is $f_*$? How can I justify $f_*, f^*$ is linear?
Definition. If $f: X \to Y$ is a smooth map and $\omega$ is a $p$-form on $Y$, define a $p$-form $f^*\omega$ on $X$ as follows: $$f^*\omega(x) = (df_x)^*\omega[f(x)].$$
The construction of lower star
For any tangent vector $v_a \in \mathbb{R}_a^n$, we can define a map $\tilde{v}_a: C^\infty(\mathbb{R}^n) \to \mathbb{R}$ by taking the directional derivative in the direction $v$ at $a$: $$\tilde{v}_af = D_vf(a)=\frac{d}{dt}\Bigg|_{t=0}f(a+tv).$$
The operation is linear and satisfies the product rule: $$\tilde{v}_a(fg) = f(a)\tilde{v}_a(g)+g(a)\tilde{v}_a(f).$$
A linear map $X: C^\infty(\mathbb{R}^n) \to \mathbb{R}$ is called a derivation at $a$ if it satisfies the following product rule: $$X(fg) = f(p)X_g + g(p)X_f$$ for all smooth function $f,g$ on $M$. If $N$ is smooth manifolds and $F: M \to N$ is a smooth map, for each $p \in M$ we define a map $F_*: T_pM \to T_{F(p)}N$, called the push-forward associated with $F$, by $$(F_*X)(f) = X(f \circ F).$$
The linearity of $f_*$
Consider linear maps $X, Y: C^\infty(\mathbb{R}^n) \to \mathbb{R}$, and scalar $\alpha \in \mathbb{R}$.
\begin{align} F_*(\alpha X + Y)(f) & = (\alpha X + Y)(f \circ F)\\ & = \alpha X (f \circ F) + Y (f \circ F)\\ & = \alpha (X (f \circ F)) + Y (f \circ F)\\ & = \alpha (F_*X)(f) + (F_*Y)(f)\\ \end{align}
The linearity of $f^*$
Suppose $A: V \to W$ is a linear map. Then the transpose map $A^*: W^* \to V^*$ extends to the exterior algebras, $A^*: \Lambda^p(W^*) \to \Lambda^p(V^*)$ for all $p>0$. If $T \in \Lambda^p(W^*)$, just define $A^* T \in \Lambda^p(V^*)$ by $$A^*T(v_1, \dots, v_p) = T(Av_1, \dots, Av_p)$$ for all vectors $v_1, \dots, v_p \in V$.
So following your definition, $$f^*\omega(x) = (df_x)^*\omega[f(x)]$$
Here we denote $x = (x_1, \cdots, x_p)$, where $x_i = \frac{\partial}{\partial x^i}$ is the $x^1$ coordinate derivation.
Hence $$(df_x)^*\omega[f(x)] = \omega(df_x x_1, \cdots, df_x x_p).$$
To show the linearity of $f^*$, we want to show $$f^*\omega(x + ay) = f^*\omega(x) + af^*\omega(y).$$
So we have $$f^*\omega(x + ay) = \omega(df_x (x_1 + ay_1), \cdots, df_x(x_p + ay_p)).$$ By linearity of $df_x$, we got
$$\omega(df_x (x_1 + ay_1), \cdots, df_x(x_p + ay_p)) = \omega(df_x x_1 + a (df_x)y_1, \cdots, df_xx_p + a(df_x)y_p)$$
$p$-forms are defined to be alternating $p$-tensors, which is multi-linear. Hence \begin{eqnarray*} &&\omega(df_x x_1 + a (df_x)y_1, \cdots, df_xx_p + a(df_x)y_p)\\ & = &\omega(df_x x_1, \cdots, df_x x_p) + a \omega((df_x)y_1, \cdots, (df_x)y_p)\\ & = & f^*\omega(x) + af^*\omega(y). \end{eqnarray*}