My question can be summarized as:
Given a potential with a symmetry (e.g. $z\rightarrow-z$), should I expect orbits in that potential to exhibit the same symmetry? Below is the full motivation for this question.
A while back I came across an intriguing orbit (or rather class of orbits) that is not symmetric with respect to the $z=0$ plane in a potential with axial symmetry (axis along $z$) as well as symmetry with respect to reflection along the $z=0$ plane. Naively I would expect the symmetry of the potential to be reflected in the symmetry of the orbits, and this is the case in every textbook example I've seen, though admittedly I haven't done an exhaustive search through all classical mechanics texts.
The effective potential is:
$$\Phi_{\text{eff}} = \frac{1}{2}v_0^2\ln\left(R^2+\frac{z^2}{q^2}\right)+\frac{L^2}{2R^2}$$
It is drawn from Galactic Dynamics 2ed eq. 3.70 (see also fig. 3.3), and is a crude model of the potential of an oblate spheroidal galaxy with constant circular speed $v_0$. The parameter $q$ controls the axis ratio, basically how "squished" the spheroid is. For the examples I show below it's fixed to $q=0.7$.
The azimuthal motion of the particle is essentially ignored - since a test particle in this potential just goes around the $z$ axis with constant angular momentum, this is easy to do. The effective potential still accounts for the angular momentum due to this motion, of course. By "following the particle" in the azimuthal direction, the motion can be described as a trajectory in the $R-z$ plane. In general, the motion of a bound particle (i.e. negative total energy) is expected to be confined to some region in this plane, the shape of which depends on the $L_z$ and $E$ of the orbit. Typical regions look like the top two panels here:
The lower right panel is also a typical orbit, but it fills in its "box" very slowly, so it looks a bit peculiar here.
The bottom left panel is the strange orbit I'm interested in. Notice that it fills in a well-defined region that is not symmetric about $z=0$. What's going on here? Why isn't the region symmetric about $z=0$ as my intuition suggests it should be? The orbital parameters $(E,L_z)$ are given above the plot. An equivalent way of expressing the orbital parameters is in terms of initial conditions. Example initial conditions that give an orbit qualitatively similar to the one in the plot:
$$R_0 = 0.3; z_0=0.1; \dot{R}_0 = 0; \dot{z}_0 = 0$$
Footnote regarding the numerics: Orbits were computed using the RK4 algorithm with a timestep of $0.002$. I've played around with the parameters (timestep, initial conditions, different integrators, etc.) to be reasonably confident that this is not a numerics issue. Furthermore, one of my colleagues has independently computed the same type of orbit, the only thing I provided was the expression for the potential and the ICs. Of course we may both be doing something wrong, but the probability for that is much lower than for just me having screwed up...