Hamiltonian with one constant of motion (besides the Hamiltonian itself) The background of my question is a well known fact: a Hamiltonian system with $n$ degrees of freedom with $n$ constants of motion is integrable. 
My question is about the case in which there are only two constants of motion, one being the Hamiltonian $H$ itself and the other a given function $J$:
$$\left\{H,J\right\}=0.$$
Books only make examples and do not provide theorems (e.g. Goldstein "Classical mechanics", third edition, sections 8.2 and 8.3, from page 343). One of the examples is when $J$ is a conserved momentum. In such examples, by means of a suitable canonical transformation, the Hamiltonian is taken to the form:
$$H(p_1, q_1, \dots, p_{n-1}, q_{n-1}, P)\tag{1}$$
in which $J=P$. Here $H$ still has $n$ degrees of freedom, i.e. the space is $2n$ dimensional, although $H$ does not depend on one of the $2n$ variables (let us call it $Q$).  The result is that the variable $Q$ is decoupled from the system, and another variable, $P$, is a constant.
The same can be expressed by saying that the system is rewritten in $n-1$ degrees of freedom, or $2(n-1)$ variables, with a reduced Hamiltonian 
$$H^r_J(p_1, q_1, \dots, p_{n-1}, q_{n-1})$$
Here the reduced Hamiltonian has $J$ as a parameter.
The books often lead the readers to think that this is a general fact.
Does anyone knows if this is a theorem? And where do I find the precise statement and proof? (edited: the answer should tackle the problem from the global point of view).
I think that the form of Eq. 1 is not general, i.e. I think that the presence of an additional constant of motion does not always lead to an expression like Eq. 1, unless some additional hypothesis is given. A counter-example is:
$$H(p_1, q_1, \dots, p_n, q_n) = H_a(p_1, q_1, \dots, p_m, q_m) + H_b(p_{m+1}, q_{m+1}, \dots, p_n, q_n)$$
where $H_a$ and $H_b$ are two Hamiltonians with only one constant of motion each (i.e. $H_a$ and $H_b$). Clearly, $H_a$ is a constant of motion for $H$, but $H$ cannot be taken to a form as in Eq.1.
 A: The specific theorem you are referring to in the background part of the question is the Liouville-Arnold theorem where an $n$-dimensional Hamiltonian system with $n$ involutory constants-of-motion is integrable. I've personally found it pretty tough to apply in practice because it usually just states a fact, but doesn't give you the recipe on how to find the explicit system. The Frobenius theorem is essentially a different version at the vector field level with generalizations to other topological spaces. There's one for Hamiltonian systems, or rather symplectic and Poisson manifolds.
The theorem you are asking about in your question is instead a corollary of the theorem on the Marsden-Weinstein-Meyer reduction. Essentially, for a $2n$-dimensional Hamiltonian system, you can remove double the number of involutory constants-of-motion. If you have $m$ constants of motion, you reduce to $2n - 2m$ dimensions. So say you have $n$ constants-of-motion, then $2n - (2)(n) = 0$. A $0$-dimensional system has been integrated completely so that agrees with Liouville-Arnold.
Moving on to your specific example where $\{H,J\} = 0$, and then
$$
H_{\text{new}}(q_1, p_1, \cdots, q_{n-1}, p_{n-1}, J)
$$
You're absolutely correct that this isn't general. This can easily be seen by counting the number of variables. There's an odd number which means this system isn't Hamiltonian at all. However, it doesn't mean it's wrong. Since $\{J, H\} = 0$, then $\dot{J} = 0$ and therefore it shouldn't be included in the dynamical dimension of the system. It's an arbitrary parameter that affects how dynamics evolve, but is not a state itself.
Finally for the counter-example where $H = H_a(q_1, p_1, \cdots, q_m, p_m) + H_b(q_{m+1}, p_{m+1}, \cdots, q_n, p_n)$, you have
$$
0 = \{H_a, H_b\}
$$
This is really easy to show by just looking at what coordinates are involved in either Hamiltonian. The issue you're having is that the canonical transformation has to take place over the specific coordinates contained in the $H_i$ terms. This specific example, with $H_a = J_a$ and $H_b = J_b$, might look something like
$$
H_{\text{new}}(q_1, n_1, \cdots, q_{m-1}, p_{m-1}, J_a, q_m, p_m, \cdots, q_{n-1}, p_{n-1}, J_b)
$$
I personally prefer to write
$$
\begin{aligned}
H_{\text{new}} : \mathbb{R}^{2(n-c)} \times \mathbb{R}^{c} &\rightarrow \mathbb{R} \\
(Q_1, P_1, \cdots, Q_{n-c}, P_{n-c}, J_1, \cdots, J_c) &\mapsto (\text{the actual function})
\end{aligned}
$$
Where $c$ is the number of $J$'s you have. The first line in that equation tells you that, yes indeed we have an even dimensional system and therefore this is truly Hamiltonian since $2(n-c)$ is always even. It also says that we have $c$ arbitrary parameters. Your case has $c=2$, but this should hold for any number of $c$. The second line tells us that $q_i = Q_i$ and $p_i = P_i$ might not necessarily be true because our transformation might have been very complicated, which they almost always are.
Example with $c=1$
Apologies in advance for this not being a "pure" physics problem nor keeping everything in a Poisson language. I have very few real examples of symplectic reduction. Applying this procedure to the Brachistochrone problem, we start with the system
$$
\begin{aligned}
\min_{\theta} J &= \int_{0}^{t_f} 1 dt \\
\text{Subject to:} \; \dot{x} &= v\cos(\theta) \\
\dot{y} &= v\sin(\theta) \\
\dot{v} &= g\sin(\theta)
\end{aligned}
$$
plus some boundary conditions, but they are unimportant to symplectic reduction. Dualizing, or turning this into a Hamiltonian problem, we have
$$
\begin{aligned}
H &= \lambda_x v \cos(\theta) + \lambda_y v \sin(\theta) + \lambda_v g \sin(\theta) \\
\dot{\textbf{x}} &= \frac{\partial H}{\partial \textbf{p}}, \dot{\textbf{p}} = -\frac{\partial H}{\partial \textbf{x}}, \frac{\partial H}{\partial \theta} = 0
\end{aligned}
$$
where $\textbf{x} = [x,y,v]$ and $\textbf{p} = [\lambda_x, \lambda_y, \lambda_v]$. The total dimension of $\textbf{x}$ and $\textbf{p}$ here is 6. $\theta$ is what's called a control variable, and you can somewhat safely ignore this variable for the purpose of reduction. By crunching the dynamic equations, you can readily see that
$$
\begin{aligned}
\dot{\lambda}_x &= 0 \\
\dot{\lambda}_y &= 0
\end{aligned}
$$
So $c=2$, but let's just focus on $\int\ \dot{\lambda}_x = J_1$ for now. The prohibitive inversion in symplectic reduction is solving $J_1$ for some state, any state really. Here it is easy, we can solve $J_1$ for $\lambda_x$.
$$
\lambda_x = J_1
$$
That's step 1. Step 2 is to eliminate $x$ from the entire system. Because $x$ is also a symmetry, by definition it doesn't appear in the equations-of-motion and therefore we don't need to do any heavy lifting to remove it. Our final system becomes
$$
\begin{aligned}
H &= J_1 v \cos(\theta) + \lambda_y v \sin(\theta) + \lambda_v g \sin(\theta) \\
\dot{\textbf{x}} &= \frac{\partial H}{\partial \textbf{p}}, \dot{\textbf{p}} = -\frac{\partial H}{\partial \textbf{x}}, \frac{\partial H}{\partial \theta} = 0 \\
\dot{x} &= v \cos(\theta)
\end{aligned}
$$
where $\textbf{x} = [y, v]$ and $\textbf{p} = [\lambda_y, \lambda_v]$. The total dimension of $\textbf{x}$ and $\textbf{p}$ here is 4.
Wait a second, but $\dot{x}$ is still in the equations? Isn't this 5 dimensional? Yes and no. There is no closed form solution for $x$ because it is a function of other quantities that have no closed form solution. However, it is not a part of the dynamics, and can therefore be integrated using numerical quadrature rules like Simpson's rule, instead of propagators like RK45 or Hermite-Simpson collocation. This has a significant impact on computational resources since quadrature can be computationally parallelized whereas it is impossible to achieve complete parallelization on RK45 or collocation-based schemes. $J_1$ is a constant value that affects how the equations of motion evolve, but is also not an equation-of-motion, which is why I made the distinction above in $H_{\text{new}}$.
Slightly More Complicated $J$
In a different problem, say for instance, $J_1 = x^2 + y^2$. This represents the radius of a circle which is constant. Inverting $J_1$ to solve for any state, we get
$$
x = \pm \sqrt{J_1 - y^2}
$$
Clearly, you would be better off choosing radius and angle coordinates, however that requires prior knowledge about the structure of the solution.
Counterargument to your Counterargument
Consider two Hamiltonians, each representing a circle
$$
\begin{aligned}
H_1 &= x^2 + y^2 \\
H_2 &= u^2 + v^2
\end{aligned}
$$
Combining the system gives you
$$
H = x^2 + y^2 + u^2 + v^2 = H_1 + H_2
$$
where $\textbf{x} = [x, u]$ and $\textbf{p} = [y, v]$. In this case, you can choose action angle coordinates for both Hamiltonians, or even worse, you can do
$$
\begin{aligned}
x &= \pm\sqrt{J_1 - y^2} \\
u &= \pm\sqrt{J_2 - v^2}
\end{aligned}
$$
There is nothing preventing me from using the above transformation to bring the dimension down to a 0-dimensional, albeit horrible system. This process can be refined to use a purely Poissonian language, though I'm not that skilled.
A: Let there be given a $2n$-dimensional symplectic manifold $(M,\omega)$ with 2 globally defined functions $$H,J: M\to \mathbb{R},$$ such that
$$ \{H, J\}~=~0. $$
Then it is an easy exercise to prove the following two propositions.


*

*Proposition 1. 
Let $p_n\!\equiv\! J$. If there is given an open neighborhood $U\subseteq M$ and $2n-1$ functions 
$$q^1,\ldots, q^n, p_1,\ldots, p_{n-1}:U \to \mathbb{R}$$ 
such that
$$(q^1,\ldots, q^n, p_1,\ldots, p_{n-1}, p_n|_U),$$ 
is a local canonical/Darboux coordinate system, then
$$\frac{\partial H}{\partial q^n}~=~0. $$

*Proposition 2. Let $p_1\!\equiv\! H$ and $p_2\!\equiv\! J$. Given a point ${\rm pt}\in M$. If $$\mathrm{d}H({\rm pt})\wedge\mathrm{d}J({\rm pt})~\neq~ 0 ,$$ then (according to the Caratheodory–Jacobi–Lie theorem) there exists an open neighborhood $U\subseteq M$ of the point ${\rm pt}\in M$ and $2n-2$ functions 
$$q^1,\ldots, q^n, p_3,\ldots, p_n:U \to \mathbb{R}$$ 
such that
$$(q^1,\ldots, q^n, p_1|_U,p_2|_U,p_3, \ldots, p_n),$$ 
is a local canonical/Darboux coordinate system.
