How to justify the direction of the magnetic field formally, in symmetrical situations? I know the right hand rule and I know how to use it, what I am asking for is a formal argument for situations where symmetry is in our favour.

The right hand rule is a kind of mnemonic but it doesn't prove anything.

According to my teacher. To give a concrete case of discussion, I will throughout my question consider an infinite cylinder placed in a coordinate system such that the OZ axis passes through its geometric center.
For electrostatics situations, I found a little trick that works every time to justify the direction of the electric field, it is enough to show, by simple geometrical arguments, which are the invariances under rotation/translation and to arrive at the conclusion that $\bf E$ is only a function of $\rho$. For a long enough cylinder there is indeed invariance by translation along Z and by rotation around Z, thus $\bf E = \bf E(\rho)$ only. Then we can use $\bf E = -\nabla\phi$ where $\phi$ is the scalar potential. By writing the gradient in cylindrical coordinates,
$$
\nabla \phi := g^{ij}\partial_i \phi {\bf ê}_j := \frac{\partial\phi}{\partial\rho}{\bf ê_{\rho}} + \frac{\partial\phi}{\rho^2\partial\varphi} {\bf ê_{\varphi}} + \frac{\partial z}{\partial\rho}{\bf ê}_{z}
$$
and by the previous symmetry arguments, I know that $\bf E$ is only a function of $\rho$ thus the derivatives with respect to $z$ and $\varphi$ are therefore zero, so the gradient reduces to
$$
\nabla \phi = \frac{\partial\phi}{\partial\rho}{\bf ê_{\rho}}
$$
Since the electric field derives from the potential, I can easily conclude that $\bf E$ is directed along $\rho$ : ${\bf E} = E(\rho) {\bf ê}_{\rho}$. Note that we did not need to deduce the expression of the potential.
Unfortunately, my argument is no longer valid in magnetostatics. I looked for ways to do the same as in electrostatics, but I didn't find anything very convincing; on the Internet I always come across the right hand rule, and even if the result is the same (and of course faster to obtain), it is not the argument I am looking for, even if for a cylinder it seems clear that the field will be directed according to $\varphi$.
So my question is: how can I formally justify, as I did above, the direction of the magnetic field without invoking a ready-made argument? Simple discussion of symmetries are fine, I am not looking for formalism for formalism's sake, but something a bit more concrete than the rule quoted many times. If possible, it would be nice to make the justification applicable to other problems than a cylinder, e.g. my argument in electrostatics works equally well for spherical symmetry for example.
 A: Started as an edit, but it was too long. Your argument for the electric case is wrong. Take the simple example:
$$
\vec E = \vec e_z
$$
which is cylindrically symmetric and derives from the potential $\phi=-z$.
Your mistake is the most restrictive form of $\vec E$ under symmetry considerations is:
$$
\vec E = E_\rho(\rho)\vec e_\rho+E_\varphi(\rho)\vec e_\varphi+E_z(\rho)\vec e_z
$$
The irrotational condition gives:
$$
\vec E = E_\rho(\rho)\vec e_\rho+\frac{A}{\rho}\vec e_\varphi+B\vec e_z
$$
with $A,B$ any real constants. If you want it to be defined on the axis as well, or you add reflections whose invariant plane contain the $z$ axis (polar vector), then $A=0$. Finally, you can have $B=0$ if you add reflections whose invariant plane is perpendicular to the $z$ axis (polar vector). It is only after adding these extra conditions that you get the desired form:
$$
\vec E = E_\rho(\rho)\vec e_\rho
$$
Part of the argument can be salvaged for magnetostatics. You get once again:
$$\vec B = B_\rho(\rho)\vec e_\rho+B_\varphi(\rho)\vec e_\varphi+B_z(\rho)\vec e_z$$
This time you use the fact that it is divergence free. Using the expression in spherical coordinates, you get only the condition:
$$
(\rho B_\rho)’=0
$$
ie $B_\rho=\frac{A}{\rho}$ with $A$ a constant and that’s as far as you can get. First off, if you want it to be defined on the axis, you get $A=0$. If you add reflections whose invariant plane contain the $z$ axis you get $A=0$ and $B_z=0$, this time because the vector is axial. Finally, if you add reflections whose invariant plane are perpendicular to the $z$ axis, $A=0$ and $B_\varphi=0$. For example currents flowing along the $z$ axis fall in the first case while rotating currents around the $z$ axis fall in the second case (you can’t have both).
For usual groups of symmetry, similar arguments apply, once positioned in the appropriate coordinates. For spherical symmetry, you get:
$$
\vec E = E_r\vec e_r\\
\vec B = \frac{A}{r^2}\vec e_r
$$
and with reflections or if you want to extend it to the origin $A=0$.
The moral is to be careful with symmetry considerations when your objects don’t transform trivially under them. Also, don’t forget reflections whose effects will depend on whether your vector is axial or polar, hence the difference between electrostatics and magnetostatics. In fact thinking in terms of reflections allows you to directly identify the non zero components of the vectors and simplify the above arguments, especially for magnetostatics where it directly gives you the direction of the vector.
Hope this helps and tell me if something’s not clear.
