In section 2.2 of the paper https://arxiv.org/abs/1807.07112, there appears a Fourier transformation named $F_k^n$ that comes out of a matrix called $F_2$,
$$ F_2 = \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1/\sqrt{2} & 1/\sqrt{2} & 0\\ 0 & 1/\sqrt{2} & -1/\sqrt{2} & 0\\ 0&0&0&-1 \end{pmatrix} $$
This matrix seems to be explained in the appendix of https://doi.org/10.1103/PhysRevA.79.032316, but I don't understand where the $c^\prime s$ appeared and how it is related to the Fourier transformation.
It's the beamsplitter unitary (a.k.a. the QFT in two dimensions, a.k.a. the Hadamard gate), represented via its action on two-mode Fermion states.
With two modes, there are four possible fermionic states: $|11\rangle\equiv c_1^\dagger c_2^\dagger |\text{vac}\rangle$ (one fermion per mode), $|01\rangle\equiv c_2^\dagger |\text{vac}\rangle$ (a fermion in the second mode), $|10\rangle\equiv c_1^\dagger |\text{vac}\rangle$ (a fermion in the first mode), and $|00\rangle\equiv|\text{vac}\rangle$ (no fermion at all).
A beamsplitter will act on these states as follows: \begin{align} |00\rangle &\to |00\rangle, \\ |10\rangle &\to \frac{1}{\sqrt2}(|10\rangle+ |01\rangle), \\ |01\rangle &\to \frac{1}{\sqrt2}(|10\rangle- |01\rangle), \\ |11\rangle &\to -|11\rangle. \end{align} To see this you just need to consider that the beamsplitter acts on the fermionic modes as $$c_1^\dagger\to\frac{1}{\sqrt2}(c_1^\dagger+c_2^\dagger), \qquad c_2^\dagger\to\frac{1}{\sqrt2}(c_1^\dagger-c_2^\dagger),$$ so that for example $$|11\rangle\equiv c_1^\dagger c_2^\dagger |\text{vac}\rangle \to \frac12(c_1^\dagger+c_2^\dagger)(c_1^\dagger-c_2^\dagger)|\text{vac}\rangle \to -c_1^\dagger c_2^\dagger|\text{vac}\rangle.$$ The rest of the rules is similarly derived. The matrix reported in the paper is simply the matrix representation of these rules.
This is a special case of a more general problem: given a unitary $U$, how does it act on many-fermion states? The general result is that the scattering amplitude between an $n$-fermion, $m$-mode input $$|r_1,...,r_m\rangle\equiv c_1^{r_1\dagger}\cdots c_n^{r_n \dagger}|\text{vac}\rangle$$ and an output $|s_1,...,s_m\rangle$ is given by the determinant of the matrix obtained from $U$ by taking its first column $r_1$ times, its second column $r_2$ times, etc., and similarly taking rows according to the occupation numbers in $|s_1,...,s_m\rangle$. For example, applying this to the above case with $U=H$ and $|r_1,r_2\rangle=|11\rangle$, we get that the probability amplitude of $|11\rangle$ evolving into $|11\rangle$ is the determinant of $H$ itself, which is $-1$, consistently with what we found by direct analysis before.
