A powerful method in theoretical physics are ladder operators. They are used in QM to solve problems like the harmonic oscillator and the hydrogen atom. The idea is to solve with their help the groundstate problem in order to get the full spectrum and eigenfunctions afterwards by successively applying them to the preceding states. Now the thing is, that most oftens these problems were solved somehow earlier so that somebody came up with a good way to choose the operators so that it works out. I have never heard of an analytical way to choose them.

My question is: Do we know how to construct them for periodic potentials?

So, if I have a Hamiltonian $H \psi = (- \frac{d^2}{dx^2} + V)\psi $, where $V$ is a smooth($\in C^{\infty}$) $2 \pi $ periodic function. If it helps, we could assume that it has a finite Fourier series expansion. In that case spectral theory (see for example Simon/Reed Analysis of Operators) tells us that $H$ has a completely discrete spectrum. Hence, in principle the chances should not be that bad that such operators exist. Furthermore, we are able to choose eigenvectors $(\psi_n)$ so that they are analytic in their argument for $x \in (0,\pi) \cup (\pi,2\pi)$ and continuous at $\{0,\pi,2\pi\}$ with $\psi_n(0) = \psi_n(2\pi)$.

Now the question is: Can we find operators $A$ and an adjoint version $A^*$ such that $A^* \psi_n = \lambda_n \psi_{n+1}$ and $A \psi_{n+1} = \mu_n \psi_n$ in this general setting without solving the problem completely?

EDIT:

I just want to include one example(where this method would fail I think) in response to the answer below, as it is not that clear to me, when this method now actually works:

Let's take the Hamiltonian on $[0,2\pi]$. $H = -d^2_x - a \cos(x) - a^2 \cos(x)$, where $a$ is some real number and $d_x:=\frac{d}{dx}$.

You find that $W(x)=a \sin(x)$ and $E_0 = -a^2$. (Maybe I should add that this superpotential is just a particular solution and the other ones are constructable from the Riccati equation, though this does not seem to be possible explicitely).

You construct $A^+ = d_x + a\sin(x)$ and $A^-=-d_x+a\sin(x)$ according to your answer.

Hence, you get something in $[A^+,A^-]$ that is dependent on $x$, so you need the generalized method that is described in your paper.

I guess that in their paper they want to choose: $V_-(x) = W^2(x)-W'(x)+E_0 = -a \cos(x) - a^2 \cos^2(x)$ and $V_+(x) = W^2(x)+W'(x)+E_0 = a \cos(x) - a^2 \cos^2(x)$, correct?

Now they define a function (see equation 2.7) $R(x,a_1):= V_-(x,a_0)-V_+(x,a_1) = -(a_0+a_1)\cos(x) +(a_1^2-a_0^2) \cos^2(x)$. Hence this one is independent of $x$ iff we have $a_0 = -a_1$. But in that case $R(x,.)=0$ and you cannot get a new energy eigenvalue accoding to (2.19), cause you stay exactly at the energy level where you already are. What does this tell me now? Did I apply this method correctly or is this method just not applicable in that case?