The problem is that $k_x$ and $x$ do not commute.
The paper is using a hybrid representation for the atom's quantum state which is very useful and natural but which can also be disconcerting if you don't know what's going on. They are considering two degrees of freedom: one discrete (the three $m_F$ hyperfine levels) and one continuous (the atom's position). One way to represent such states is to use a column vector representation for the discrete d.o.f. and a wavefunction notation for the continuous one. Thus a quantum state $|\Psi\rangle$ is represented as
$$\begin{pmatrix}\psi_{-1}(x)\\\psi_0(x)\\\psi_{1}(x)\end{pmatrix},$$
where $\psi_m(x)=\langle m,x|\Psi\rangle$. Operators are similarly represented by three-by-three matrices whose entries are themselves operators on $L_2(\mathbb{R})$.
More to the point, they are transforming the hamiltonian
$$
H=\begin{pmatrix} \frac{k^2_x}{2m}+\epsilon_1 & \frac{\Omega}{2}e^{i2k_0x} & 0 \\ \frac{\Omega}{2}e^{-i2k_0x} & \frac{k^2_x}{2m} & \frac{\Omega}{2}e^{i2k_0x} \\0 & \frac{\Omega}{2}e^{-i2k_0x} & \frac{k^2_x}{2m}-\epsilon_2 \end{pmatrix}
$$
via the unitary transformation
$$
U=\left(\begin{array}{ccc}e^{-i2k_0 x} & 0 & 0 \\0 & 1 & 0 \\0 & 0 & e^{i2k_0x}\end{array}\right).
$$
What you get is something of a mess at first:
$$
H_{\text{eff}}=UHU^\dagger=\begin{pmatrix}
e^{-i2k_0 x}\frac{k^2_x}{2m}e^{i2k_0 x}+\epsilon_1 & \frac{\Omega}{2} & 0 \\
\frac{\Omega}{2} & \frac{k^2_x}{2m} & \frac{\Omega}{2} \\
0 & \frac{\Omega}{2} & e^{i2k_0 x}\frac{k^2_x}{2m}e^{-i2k_0 x}-\epsilon_2
\end{pmatrix}.
$$
As you can see, the $x$ dependence of the off-diagonal coupling terms has been "rotated" out, as it should be. The energies, however, need special handling, because $k_x$ does not commute with $x$. Instead, they obey the canonical commutation relation
$$[x,k_x]=i,$$
of course.
In this particular case, the operator $e^{-i2k_0 x}$ in the combination $e^{-i2k_0 x}\frac{k^2_x}{2m}e^{i2k_0 x}$ is acting as a displacement operator: multiplying by an exponential phase is equivalent to a translation in momentum space. (This should be obvious by how both operations relate to each other under a Fourier transform. For a wavefunction (instead of operator-viewpoint) look at this, see this answer.)
You can and should work the operator algebra to show that
$e^{-i2k_0 x}\frac{k^2_x}{2m}e^{i2k_0 x}=\frac{(k_x+2k_0)^2}{2m}$
as described, but the more general statement is that
$$e^{-i2k_0 x}f(k_x)e^{i2k_0 x}=f(k_x+2k_0).$$
This is enough to show that the transformed hamiltonian is
$$
H_{\text{eff}}=\left(\begin{array}{ccc} \frac{(k_x+2k_0)^2}{2m}+\epsilon_1 & \frac{\Omega}{2} & 0 \\ \frac{\Omega}{2} & \frac{k^2_x}{2m} & \frac{\Omega}{2} \\0 & \frac{\Omega}{2} & \frac{(k_x-2k_0)^2}{2m}-\epsilon_2 \end{array}\right).
$$
