You seem want to introduce gauge invariance into a theory that doesn't appear to have the global symmetry need in the first place. One way to think of gauge invariance is that you 'gauge' the global symmetry, then you just change your derivative terms to covariant derivatives like you mentioned. In other words, we can only concern ourselves with the global symmetry for now and gauge it at the very end if we wish. Now, at the very least in the term you wrote down
$$\phi \bar{\psi} \psi$$
it's not clear how the indices are contracted. Do the $\phi$ and $\psi$ have indices? For example I could make the $\psi$ transform in under $SU(2)$ by introducing a second copy of the $\psi$ and sum over them:
$$ \phi \bar{\psi^a} \psi^a $$
but now I can't make the $\phi $ transform because there isn't anything left to contract the $\phi$ index with. That is,
$$ \phi^b \bar{\psi^a} \psi^a $$
isn't a singlet (a singlet doesn't transform under the symmetry ) so it doesn't make any sense as a term in your lagrangian. Or you could introduce a second spinor that is a singlet under the global symmetry so you have something to contract your scalar indices with:
$$ \phi^a \bar{\psi^a} \eta $$
Finally, if you want a spinor to transform as a triplet, or in the adjoint of $SU(2)$ you could introduce the generator of SU(2) and contract indices as follows:
$$ \phi^a (t^a)^{bc} \bar{\psi^b} \psi^c = \phi^a \bar{\psi} t^a \psi$$
where the $t^a$ is in the fundamental (or doublet) representation so that we can properly trade the adjoint index for two doublet indices and get an overall singlet for the lagrangian.
Now, as for the kinetic terms, like you mentioned, if you want to introduce gauge symmetry trade the regular derivative for a covariant derivative.
$$ \partial_\mu \phi^a \rightarrow D_\mu \phi^a =\partial_\mu \phi^a + i {A_\mu}^b (t^b)^{ac}\phi^c $$
where the $t^a$ is in whatever rep the $\phi^a$ transforms in. The same holds for the fermion fields.
There is one caveat to all this however - this prescription of naively gauging a global symmetry that I have outlined breaks down if the global symmetry is 'anomalous'. That is, quantum mechanical effects break the naive, classical global symmetry. I'm not going to get into what that is, but keep it in the back of your mind for the time being and read about it when you have a chance.
I have a feeling you might want more info than this, but I'll stop here for now and if you edit I will clarify/ add.
EDIT: In retrospect this seems to work more easily for $SU(2)$ reps easier than other groups since for $SU(2)$ the adjoint rep is the same as the triplet rep so I can trade triplet rep indices for doublet indices using the generators $(t^a)^{ij}$. I am not sure if you can do these sort of things for groups in general.
In general, the Wetterich equation reads
$$k \partial_k \Gamma_k = \frac{1}{2} {\rm Tr}\left[\left(\Gamma^{(2)}_k + R_k \right)^{-1} \cdot k \partial_k R_k \right],$$
where $\Gamma^{(2)}_k$ is a matrix of second-order derivatives of $\Gamma^{(2)}_k$ with respect to its arguments (here, $\phi$, $\psi$, and $\bar{\psi}$) and $R_k$ is the corresponding matrix of regulators. This is Eq. (4) in the paper you linked.
The original expression for the Wetterich equation missed an inverse on $\left(\Gamma^{(2)}_k + R_k \right)^{-1}$ before the edit. This inverse is important in general, because even if the regulator $R_k$ only couples to $\phi \phi$ and $\psi \bar{\psi}$ terms, the inverse involves all $9$ derivatives of $\Gamma_k$, not just $\delta^2 \Gamma_k/\delta \phi \delta \phi$ and $\delta^2 \Gamma_k/\delta \psi \delta \bar{\psi}$. This means that, in order to evaluate this inverse, in principle you will need to make your ansatz for $\Gamma_k$, evaluate all $9$ second derivatives, and then set $\phi$, $\psi$, and $\bar{\psi}$ to constant values, independent of space/momenta, so that you can actually evaluate the inverse. (The cited paper does not mention this point explicitly, but since they include a local potential $U_k(\phi^2/2)$ in their ansatz, they presumably take $\phi$ to be a scalar variable, and perhaps they set the $\psi$'s to be either zero or some other constant value). In this particular example, the inverse does simplify nicely when evaluated at zero fields (as shown for a simplified example below); the issue then comes from the fact that one needs to define the running fields in terms of derivatives of $\Gamma_k$ and differentiate the flow equation before evaluating at zero fields.
Note that the ansatz for $\Gamma_k$ that you have written didn't denote which quantities were assumde to run with $k$ (in the original post, before editing), so just to compare to the cited paper's ansatz (Eq. (5)),
$$\Gamma_k[\phi,\psi,\bar{\psi}] = \int dx~\left[ Z_k (\partial_\mu \phi)^2 + U_k(\phi^2/2) + Z_{\psi,k} \bar{\psi} i \gamma^\mu \partial_\mu \psi + i h_k \phi \bar{\psi} \psi \right].$$
It appears that you wish to explicitly truncate the authors' $U_k(\phi^2/2)$ at second order in $\phi^2$ (around $\phi=0$) and introduce an $i \bar{\psi} M_k \psi$. You will need to define how each coefficient is extracted from $\Gamma_k$ in order to compute your flow equations. For example, by defining $m^2_k = \lim_{p \rightarrow 0}\frac{\delta^2 \Gamma_k}{\delta \phi(p) \delta \phi(q)}\Big|_{\phi(p) \rightarrow \phi,\psi \rightarrow 0,\bar{\psi} \rightarrow 0}$, $M_k = \lim_{p \rightarrow 0}\frac{1}{i} \frac{\delta^2 \Gamma_k}{\delta \bar{\psi}(p) \delta \psi(q)}\Big|_{\phi(p) \rightarrow \phi,\psi \rightarrow 0,\bar{\psi} \rightarrow 0}$, $Z_k = \lim_{p \rightarrow 0}\partial_{p^2}\left[\frac{\delta^2 \Gamma_k}{\delta \phi(p) \delta \phi(-p)}\Big|_{\phi(p) \rightarrow \phi,\psi \rightarrow 0,\bar{\psi} \rightarrow 0}\right]$, etc. This will involve taking derivatives of the flow equation with respect to the fields in order to get $k \partial_k m^2_k$, $k \partial_k M_k$, etc. In doing so, it is again important to treat the inverse $\left(\Gamma^{(2)}_k + R_k \right)^{-1}$ carefully, since this will generate non-trivial terms in the flow equations for these quantities (using the identity $\frac{\delta}{\delta \phi} A^{-1} = -A^{-1} \frac{\delta A}{\delta \phi} A^{-1}$). See, for example, this paper for an example application to $\phi^4$ theory.
Edit: When multiple fields are involved, as in this case, or in a non-equilibrium problem, the inverse is a matrix inverse in the fields (when converted to momentum space). I'll sketch it for the cited paper's ansatz simplifed to $U_k(\phi^2/2) = m_k^2 \phi(x)^2/2$. In momentum space (assuming no mistakes...)
$$\begin{array}{c c}\Gamma_k &= \int dp_1 dp_2~\Big[(Z_k p^2 + m^2_k) \phi(p_1)\phi(p_2)\delta(p_1 + p_2) + Z_{\psi,k} i \bar{\psi}(p_1) \gamma^\mu (p_1)_\mu \psi(p_2) \delta(p_1+p_2)\Big] \\
& ~~~~~~~~ + \int dp_1 dp_2 dp_3~ i h_k \phi(p_1)\bar{\psi}(p_2) \psi(p_3) \delta(p_1 + p_2 + p_3)\end{array}.$$
(I am not very familiar with dealing with fermionic fields, so you'll want to double-check my math here).
Then, the matrix of derivatives is
$$\Gamma^{(2)}_k = \begin{bmatrix} \frac{\delta \Gamma_k}{\delta \phi(p) \delta \phi(q)} & \frac{\delta \Gamma_k}{\delta \phi(p) \delta \psi(q)} & \frac{\delta \Gamma_k}{\delta \phi(p) \delta \bar{\psi}(q)} \\ \frac{\delta \Gamma_k}{\delta \psi(p) \delta \phi(q)} & \frac{\delta \Gamma_k}{\delta \psi(p) \delta \psi(q)} & \frac{\delta \Gamma_k}{\delta \psi(p) \delta \bar{\psi}(q)} \\ \frac{\delta \Gamma_k}{\delta \bar{\psi}(p) \delta \phi(q)} & \frac{\delta \Gamma_k}{\delta \bar{\psi}(p) \delta \psi(q)} & \frac{\delta \Gamma_k}{\delta \bar{\psi}(p) \delta \bar{\psi}(q)} \end{bmatrix} = \begin{bmatrix} (Z_k p^2 + m_k^2)\delta(p+q) & \int dp_2~ i h_k \bar{\psi}(p_2) \delta(p+q+p_2) & \int dp_3~ i h_k \psi(p_3) \delta(p+q+p_3) \\ \int dp_2~ i h_k \bar{\psi}(p_2) \delta(p+q+p_2) & 0 & Z_{k,\psi} i \gamma^\mu p_\mu \delta(p+q) + \int dp_1~ i h_k \phi(p_1) \delta(p+q+p_1) \\ \int dp_3~ i h_k \psi(p_3) \delta(p+q+p_3) & Z_{k,\psi} i \gamma^\mu p_\mu \delta(p+q) + \int dp_1~ i h_k \phi(p_1) \delta(p+q+p_1) & 0 \end{bmatrix}$$
Then, one would add to this the matrix
$$R_k = \begin{bmatrix} R_{\phi,k} & 0 & 0 \\ 0 & 0 & R_{\psi,k} \\ 0 & R_{\psi,k} & 0 \end{bmatrix}\delta(p+q)$$
To calculate the inverse in practice, we need to choose functions $\phi(p)$, $\psi(p)$, and $\bar{\psi}(p)$ at which we wish to evaluate the inverse, since we cannot do it for arbitrary functions. Typically, one chooses momentum-independent values. If we had kept the full $\phi$ dependence of $U_k(\phi^2/2)$, we would typically set $\phi(p) = \phi$ (a scalar), but since we've truncated it, we could set all of $\phi(p) = \psi(p) = \bar{\psi}(p) = 0$, giving
$$\Gamma^{(2)}_k + R_k = \begin{bmatrix} (Z_k p^2 + m_k^2) + R_{\phi,k} & 0 & 0 \\ 0 & 0 & Z_{k,\psi} i \gamma^\mu p_\mu + R_{\psi,k} \\ 0 & Z_{k,\psi} i \gamma^\mu p_\mu + R_{\psi,k} & 0 \end{bmatrix}\delta(p+q),$$
which can now be evaluated as a matrix inverse (the inverse of the momentum-dependence here is just another $\delta(p+q)$). For this particular choice, the inverse does end up just being 1 over the non-zero elements. If we had chosen constant values $\phi = \phi_0$, $\psi = \psi_0$, and $\bar{\psi} = \bar{\psi}_0$, the inverse would not be as trivial.
Now, for the choice $\phi(p) = \psi(p) = \bar{\psi}(p) = 0$, you'll notice that the dependence on $h_k$ dropped out when we set the fields to $0$. That does not mean $h_k$ would not flow in this model because, as I mentioned above, one needs to define $h_k$ in terms of derivatives of $\Gamma_k$, e.g., $h_k = \frac{1}{i}\frac{\delta^3 \Gamma_k}{\delta \phi(p_1) \delta \bar{\psi}(p_2) \delta \psi(p_3)}\Big|_{\phi = \bar{\psi} =\psi = 0}$, and then differentiate the flow equation with respect to these fields before setting $\phi = \bar{\psi} =\psi = 0$. The result should be a non-trivial flow equation for $h_k$ (coupled to other running parameters).
Best Answer
1, 2: Yes. Putting aside issues of whether the Hilbert space of the interacting theory is well-defined, you can think of the Hilbert space here as the tensor product of the two Fock spaces. For example, the noninteracting vacuum is the state with no fermions, antifermions, or bosons.
3: No. Remember that $\phi$ has the capability of creating and annihilating bosonic particles. So the interaction term $\bar{\psi} \psi \phi$, when acting on a Fock state with some number of fermions/bosons, has the ability to implement one of the following processes: