The numerical flux is chosen to ensure that information in the problem travels in the direction of the characteristic curves of the equation (upwinding). As mentioned in the comments, the numerical flux is necessary in order to couple the subproblems defined on each element.
One way to get an intuition for the role of the numerical flux is to consider the following simple example.
Consider the scalar advection equation (where for simplicity $a=1$) $$ \frac{\partial u}{\partial t} + \frac{\partial u}{\partial f} = 0 \qquad\text{on $\Omega$}, $$ where the domain is given by $\Omega = [0,1]$. Because this is a hyperbolic equation, and information is propagating from left to right, we need to enforce a boundary condition at $x = 0$ (but not at $x = 1$). For concreteness, suppose we enforce the Dirichlet condition $u(0,t) = g_D$ for some given $g_D$.
Suppose now we discretize this equation using the DG method, and we use two elements, $D_1 = [0,1/2]$ and $D_2 = [1/2,1]$. We could equally-well be discretizing the following set of two coupled PDEs, \begin{align*} \text{(PDE 1):}&& \quad v_t + v_x &= 0 \quad\text{on $D_1$},\\ \text{(PDE 2):}&& \quad w_t + w_x &= 0 \quad\text{on $D_2$}, \end{align*} where we will couple these equations to make them equivalent to the original equation.
To make the above equations well-posed, we need to enforce boundary conditions. As before, each equation is hyperbolic, and information is traveling from left to right. Therefore, we need to enforce a boundary condition for (PDE 1) on the left endpoint of $D_1$, and a boundary condition for (PDE 2) on the left endpoint of $D_2$.
The boundary condition on the left endpoint of $D_1$ must be chosen to be $v(0,t) = g_D$ in order to remain consistent with the original problem. We also look for smooth solutions, so the boundary condition on the left endpoint of $D_2$ must be chosen to enforce continuity. This condition reads $w(1/2,t) = v(1/2,t)$.
The DG method in this case chooses the numerical fluxes precisely to enforce the above boundary conditions. If we multiply by a test function $\psi$ and integrate by parts over each element $D_k$, we obtain boundary terms of the form \begin{align*} \int_{\partial D_1} \hat{n} \cdot v \psi \, dx &= \left[ v\psi \right]_0^{1/2}\\ \int_{\partial D_2} \hat{n} \cdot w \psi \, dx &= \left[ w\psi\right]_{1/2}^1 \\ \end{align*} In order to "weakly" enforce the boundary conditions, we replace $v$ and $w$ with the prescribed values at those points where boundary conditions are specified (i.e. the left endpoints of $D_1$ and $D_2$). This means we replace $v(0,t)$ by $g_D$ and $w(1/2,t)$ by $v(1/2,t)$ in the boundary integrals.
In other words, we define $u_h^* = g_D$ at $x = 0$, and $u_h^* = v(1/2,t)$ at $x=1/2$, and we recover exactly the standard upwind flux that is used in the DG method.
Looking at things this way, we can consider the numerical flux functions as weakly enforcing the boundary conditions on each element that are required to couple the equations in such a way that respects the characteristic structure of the equations.
For equations more complicated than constant-coefficient advection, information may not propagate always in the same direction, and so the numerical flux must be determined by solving (or approximating the solution to) a Riemann problem at the interface. This is discussed for linear problems in Section 2.4 of Hesthaven's book.