© The Authors 2024

1. Introduction

The problem of galactic rotation is the empirical statement that rotational velocity around the galactic center seems to flatten out for a large fraction of the population of galaxies in which this has been measured at long-enough distances (Rubin et al. 1978). This creates a puzzle if the data is interpreted in terms of rotation external to a spherical source around which orbital equilibrium (Kepler’s third law written for the velocity, v) demands falling v for stellar objects or clouds of gas,

$$\begin{array}{c}\hfill \frac{v^2}{r}=\frac{\mathit{GM}}{r^2}\Rightarrow v=\sqrt{\frac{\mathit{GM}}{r}}·\end{array}$$(1)

The most accepted explanation is the existence of additional dark matter (DM) or a modification of the law of motion under gravity.

Because typical velocities in spiral galaxies are on the order of 200−300 km s⁻¹, v/c ∼ 10⁻³, relativity is a correction and Newtonian mechanics should get the bulk of the rotation right. Therefore, either a modification of mechanics, such as modified Newtonian dynamics (MOND; Milgrom 1983), or a modification of the gravity source, typically in the form of a spherical DM halo (Frenk et al. 1985), are invoked. MOND however runs into problems at larger, cosmological scales (Aguirre et al. 2001; Dodelson & Liguori 2006); and a spherical DM distribution has to be fine-tuned to have very nearly an isothermal ρ(r)∝1/r² profile to explain the flatness of the rotation curve.

If we inhabited a two-dimensional cosmos, however, the natural gravitational law would be $|F|\propto \frac{1}{r}$ instead of $\propto \frac{1}{r^2}$ and the observed rotational law would be v_2D ∝ constant, which is the law that the experimental data demands. We do not; but a cylindrical matter source achieves the same dimensional reduction by providing translational symmetry along the OZ symmetry axis of the cylinder (Slovick 2010; Llanes-Estrada 2021). If the linear density of the cylindrical DM source is λ, we can write

$$\begin{array}{c}\hfill \frac{v^2}{r}=\frac{2G\lambda }{r}\Rightarrow v=\sqrt{2G\lambda }.\end{array}$$(2)

That is, the constant velocity function, v(r), is natural for a filamentary source. Moreover, if the rotation curve is only measured to a finite r, it is obviously the case that the source does not need to be infinitely cylindrical: it suffices that it be prolate (elongated) instead of spherical, as is shown by detailed fits (Bariego-Quintana et al. 2023) to the SPARC database (Lelli et al. 2016) and as is consistent with simulations of DM halos (Allgood et al. 2006).

Observables in the galactic plane alone, such as detailed rotation curves, cannot distinguish between competing models such as spherical halos with nearly ${\rho }_{\mathrm{DM}}(r)\propto \frac{1}{r^2}$ profiles or elongated halos with arbitrary profile. To lift the degeneracy between shape and profile one needs to find adequate, simple observables from data outside the galactic plane.

For a while now, stellar streams (Ibata & Gibson 2007) in the Milky Way (MW) Galaxy have been a promising new source of information on the DM distribution (Nibauer et al. 2023), as they will eventually be for other galaxies (Pearson et al. 2022), or, more generically, for the effective gravitational potential (Sanders 2014). Streams have been proposed (Erkal et al. 2016) to detect DM clumps (for example, presumed subhalos) that may collide with a stream and leave a gap imprinted therein.

In the rest of this article, we develop what we think is a key observable to be measured on those streams that has a bearing on the question of the overall shape of the presumed halo. Section 2 is dedicated to reviewing the definition of torsion in differential curve geometry, because it is the parameter which locally quantifies the aplanarity of a curve in three dimensions, and to showing that, around a spherical halo, orbits as well as streams are torsionless. (The distinction is that streams can be thought of as a bundle of nearby orbits, one for each constituent body. If a whole stream is taken as a perfect orbit, biases in the full reconstruction of the gravitational potential can creep in Sanders & Binney 2013.)

Section 3 then shows how we expect tidal streams around elongated gravitational sources to show torsion if there is a component of the velocity parallel to the axis of elongation of the source. Section 4 makes a reasonable selection among the known stellar streams in our MW Galaxy and we plot the torsion calculated along each of them, showing that there seems to be a signal here. Section 5 then discusses how further studies can improve the conclusion.

2. Orbits and streams around central potentials are torsionless

2.1. Torsion quantifies separation from orbital planarity

Before explaining why we wish to propose torsion as a useful observable to probe the DM halo – to fix notation – let us recall a few concepts of differential geometry (do Carmo 2017).

In brief, the torsion of a curve measures how sharply it twists out of its osculating plane, instantaneously defined by the tangent and normal vectors (equivalently, by the velocity and normal acceleration).

To a curve, r(t), in three-dimensional space parameterized by an arbitrary variable, t, we can associate an arclength, $s(t)={{\displaystyle \int t}}^{}|\mathit{r}(t^{\prime })|\mathrm{d}{t}^{\prime }$, and the tangent vector, T = dr/ds. If at a certain point, P, the curvature is nonzero, then the normal vector at P is defined by N = dT/ds (its inverse modulus giving the radius of the circumference best approximating the curve at P), and the binormal vector (that completes the Frenet-Serret trihedron) by the vector product of both,

$$\begin{array}{c}\hfill \mathit{B}=\mathit{T}\times \mathit{N}.\end{array}$$(3)

If the curve is perfectly planar, the tangent and normal vectors will always lie in the same plane, and in such a case the binormal vector stays parallel to itself along the curve. Any natural definition of torsion will then yield zero.

But if the curve twists out of its osculating plane (like a uniformly advancing helix, which corresponds to constant torsion), the binormal vector will acquire a rotation. Torsion will then measure the speed of that rotation of the binormal – which is a locally defined vector at each point, P, along the curve, r(t) – as the scalar product of the intrinsic derivative of B and the normal vector (this discounts the change of the modulus of B and rather measures its twisting),

$$\begin{array}{c}\hfill \tau =-\frac{\mathrm{d}\mathit{B}}{\mathrm{d}s}·\mathit{N}.\end{array}$$(4)

If the arc length is not to hand and the arbitrary parameter, t, needs to be used, then a convenient formula (with the prime denoting d/dt) is

$$\begin{array}{c}\hfill \tau =\frac{({\mathit{r}}^{\prime }\times {\mathit{r}}^{″})·{\mathit{r}}^{″\prime }}{{|{\mathit{r}}^{\prime }\times {\mathit{r}}^{″}|}^2}·\end{array}$$(5)

Since up to three derivatives of the position along the curve need to be computed, several adjacent points of a discretized curve are needed to extract the torsion: but it is still quite a local observable that does not need long trajectory stretches.

We are going to demonstrate the use of this observable τ for stellar streams, particularly around the MW, to determine the shape of the gravitational potential of its DM Halo.

2.2. Movement around a Newtonian spherical source

Newtonian gravity predicts, for motion around a spherical body,

$$\begin{array}{c}\hfill {\mathit{r}}^{″}=\frac{\mathit{F}}{m}=-GM\frac{\mathit{r}}{{|\mathit{r}|}^3},\end{array}$$(6)

with M the mass inside the sphere of radius |r|. The necessary third derivative can be computed in a straightforward manner, taking into account that $|\mathit{r}{|}^{\prime }=\widehat{\mathit{r}}·{\mathit{r}}^{\prime }$ is the modulus of the projection of the velocity along the visual from the origin,

$$\begin{array}{c}\hfill {\mathit{r}}^{″\prime }=\frac{-GM}{{|\mathit{r}|}^2}({\widehat{\mathit{r}}}^{\prime }-\frac{2}{|\mathit{r}|}(\widehat{\mathit{r}}·{\mathit{r}}^{\prime })\widehat{\mathit{r}}),\end{array}$$(7)

in terms of components along the velocity and along the position.

Because of Eq. (6),

$$\begin{array}{c}\hfill ({\mathit{r}}^{\prime }\times {\mathit{r}}^{″})\propto ({\mathit{r}}^{\prime }\times \mathit{r}),\end{array}$$(8)

and therefore, observing that both terms of Eq. (7) are proportional to either r′ or r, we see that (r′×r″)⊥r″′. Therefore, the scalar product in the numerator of Eq. (5) vanishes, and thus τ = 0 for motion around a spherical body.

The planarity of the orbit around a central potential is, of course, a textbook consequence (Goldstein et al. 2001) of the conservation of the direction of the angular momentum vector, $\widehat{\mathit{L}}$, which in this language is parallel to the binormal vector. And, additionally, the Newtonian gravity law is not strictly necessary: any central potential will yield the same result. This observation is of particular interest for the MOND explanation of the galactic rotation curves in Eq. (1) since, while the intensity of the acceleration induced by matter is different from Newtonian mechanics, the central direction of the force is respected: MOND likewise predicts no torsion.

2.3. Simulation of an N-point stellar stream around a spherical gravitational source

The discussion presented in Sect. 2.2 refers to the torsion of one test body under a central field. If the body lost dust grains, forming a kind of contrail, its shape through space would be a planar curve. But stellar streams are not quite of this nature, and rather the result of the tidal stretching of a globular cluster or dwarf galaxy (Noreña et al. 2019). Since each object in the cluster starts off at a different height, z, with respect to the galactic plane, its orbit around the center of force lies on a slightly different plane, the effect being that the cluster, in addition to stretching, contorts, with the upper particles passing under the center of mass (cm) to become the lower ones with each half orbit.

We here show that this effect is negligible and the torsion of a stream around a central potential can safely be neglected, as the cm of the stream follows one of the trajectories of Sect. 2.2, with τ = 0. Rather than entangling the discussion in detailed theory, a couple of simple simulations will illustrate the point.

We simulated a globular cluster of N (specifically 100, 200, 500, 700, and 1000 particles) point-like stars, randomly distributed at t = 0 over a sphere of radius R₀ (on the order of one or a few kpc, typically; in the following simulation, 2 kpc) at a distance, |r|, from the galactic center (of order 30 kpc in the following example), with a certain mass, $m_i^{*}$, and common initial velocity, v_* ∼ 220 kpc. An additional random velocity kick, Δv₀ = Gm_cluster/2r₀, in a random direction was given. We then let it evolve under the gravitational force of the central source with mass M ∈ (10⁹, 10¹²) M_⊙, standing for a galactic bulb or a spherical DM halo, and we allowed for a correction due to the inner binding forces of the cluster. This was small because the random masses were taken in the interval m_* ∈ (0, 20) M_⊙, and thus their mutual interactions were orders of magnitude smaller than those with the galactic center. The constant, GM, of the central source can conveniently be eliminated in terms of the typical velocity of circular orbits around the galactic center, from the orbital equilibrium, $v_{\text{rot}}^2/r=GM/r^2$. For the MW, this is typically 220 km s⁻¹.

We show the simulation in Fig. 1. The concentrated orange points in the stream visualization plots (three-dimensional views as well as Cartesian projections, as is marked on the axes) mark the initial cluster at t = 0. The evolved cluster at later times, the clouds of cyan, blue, and purple dots, is seen to stretch under tidal tensions.

Fig. 1. Motion of a star ensemble around around a spherical source. Top row: simulated movement of a 200-star cluster around a spherical gravitational source with parameters akin to those of the MW, in the three planar projections, extracted at times t = 0 Myr, t = 200 Myr, t = 400 Myr, t = 600 Myr, and t = 800 Myr. Lower left: three-dimensional view with the containing plane. The initial cluster has a radius of 2 kpc and is at a distance of 30 kpc from the galactic origin of coordinates. Member stars have slightly different velocities, v* = 220 kpc ± Δv, due to the additional random Δv. (All axes are in kpc). Lower right: extracted torsion as a function of time, from an analytical helicoidal fit to the resulting stream with χ2/Ngdl = 0.86, which yields a value of τ ≃ 3.5 × 10−4 kpc−1 (a one-star orbit would return exactly zero torsion).

The stream remains near the (slightly tilted) plane containing the initial velocity, without developing motion outside that plane, and thus with negligible torsion for each individual orbit. Nevertheless, a measurable τ ∼ O(10⁻⁴) is visible on the lower right plot, due to the need to fit a curve to the dispersed, crisscrossing points.

3. Orbits around elongated potentials and τ ≠ 0

3.1. Movement around a Newtonian cylindrical source

We now move on to quickly show how torsion is expected to look for an orbit around a perfectly cylindrical source of gravity, in a discussion paralleling that of Sect. 2.2. We naturally employ cylindrical coordinates (ρ, φ, z), so that, denoting the time derivative by a ′ sign as in Eq. (4),

$$\begin{array}{cc}\hfill \mathit{r}& =\rho \widehat{\mathit{\rho }}+z\widehat{\mathit{z}}\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\mathit{r}}^{\prime }& ={\rho }^{\prime }\widehat{\mathit{\rho }}+\rho {\phi }^{\prime }\widehat{\mathit{\phi }}+z^{\prime }\widehat{\mathit{z}}\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\mathit{r}}^{″}& =({\rho }^{″}-\rho {\phi }^{\prime 2})\widehat{\mathit{\rho }}+(2{\rho }^{\prime }{\phi }^{\prime }+\rho {\phi }^{″})\widehat{\mathit{\phi }}+z^{″}\widehat{\mathit{z}},\hfill \end{array}$$(11)

where in the acceleration we recognize, from left to right, the radial, centrifugal, Coriolis, azimuthal, and vertical accelerations, respectively. The force law is the same as that for a line of charge in electromagnetism, except of course with the constant replaced, so that in terms of the linear mass density, λ,

$$\begin{array}{c}\hfill {\mathit{r}}^{″}=\frac{\mathit{F}}{m}=\frac{-2G\lambda }{\rho }\widehat{\mathit{\rho }}+0·\widehat{\mathit{\phi }}+0·\widehat{\mathit{z}}.\end{array}$$(12)

Comparing this with the general form in Eq. (11), we recover z″ = 0 ⟹ z′=constant (reflecting translational invariance along the OZ axis) and ρ²φ′=constant so that the third component, L_z/m, of angular momentum per unit mass is conserved, just as in the central force problem. However, now the direction of L is not conserved, so that the binormal vector changes and one expects a torsion. To obtain it, we started from Eq. (12) and took a further derivative to obtain

$$\begin{array}{c}\hfill {\mathit{r}}^{″\prime }=\frac{-2G\lambda }{\rho }(-\frac{{\rho }^{\prime }}{\rho }\widehat{\mathit{\rho }}+{\phi }^{\prime }\widehat{\mathit{\phi }})\end{array}$$(13)

(valid for the Newtonian force with cylindrical symmetry only). Calculating the cross product of Eqs. (10) and (12), while using the righthandedness of the trihedron $(\widehat{\mathit{\rho }},\widehat{\mathit{\phi }},\widehat{\mathit{z}})$ to evaluate each basis vector product, yields

$$\begin{array}{c}\hfill {\mathit{r}}^{\prime }\times {\mathit{r}}^{″}=(-2G\lambda )[-{\phi }^{\prime }\widehat{\mathit{z}}+\frac{z\prime }{\rho }\widehat{\mathit{\phi }}].\end{array}$$(14)

Next, we took the scalar product with r″′ and evaluated Eq. (5) to obtain the torsion, yielding

$$\begin{array}{c}\hfill \tau =\frac{z^{\prime }{\phi }^{\prime }}{{(\rho {\phi }^{\prime })}^2+z^{\prime 2}}=\frac{1}{\rho }\frac{v_z{v}_{\phi }}{v_{\phi }^2+v_z^2},\end{array}$$(15)

which we have cast in an easier-to-remember form in the second expression. Clearly, for there to be a torsion we need both azimuthal and vertical velocities (so that stellar streams in the galactic plane are not sensitive, as expected). Additionally, because $\left|v_z{v}_{\phi }\right|\le v_z^2+v_{\phi }^2$, torsion belongs to the interval τ ∈ [ − ρ⁻¹, ρ⁻¹], so its maximum magnitude is controlled by the distance from the stellar stream segment to the galactic axis.

3.2. Simulation of an N-point stellar stream around a cylindrical source

Next, we proceeded to repeat the exercise of Sect. 2.3 with the same starting data, but replacing the central spherical Newtonian source with a cylindrical source. The force, relegated to Eq. (B.4), now distinguishes the vertical acceleration from its two components parallel to the galactic plane.

Its first term is the acceleration caused by the cylindrical gravitational source (that along the OZ axis being zero by translational symmetry). Its linear mass density, λ = M/L, was obtained from the typical rotation curve around a galaxy, $v_{\mathrm{rot}}=\sqrt{2G\lambda }$ (Llanes-Estrada 2021). The second term of Eq. (B.4) is, again, the correction due to the tiny binding of the stellar streams stars among themselves, together with Eq. (B.4).

We show two typical computations in Figs. 2, a three-dimensional view, and 3, a more detailed analysis of a run.

Fig. 2. Spread of a 100-body stellar cluster after t = 0, 125, 250, 375, and 500 Myr. The resulting stream, with an initial vertical speed of vz = 50 km s−1, at different times clearly shows an upwardly advancing helix moving around the z axis.

Fig. 3. Motion of a star ensemble around around a perfectly cylindrical toy-model source. The initial velocity is vz = 90 km s−1, with the same conventions as in Fig. 1. Upper plots: a stream forms along a spiral trajectory. From a helicoidal fit to the cloud of points, as needs to be done on real data, with χ2/Ngdl = 1.16 we find a value for the torsion of τ ≃ 1.1 × 10−2 kpc−1 (lower plot). This is in agreement with the analytical estimate in the main text below Eq. (16).

If the starting velocity profile were perfectly set in the XY plane perpendicular to the cylinder, the torsion would still be zero, as per Eq. (15). We gave it a slight tilt and then the orbit started behaving as a helix. This can be appreciated by studying the originally compact cluster of stars at different times up to 1 Gyr, seeing it spread and ascend in a spiral, showing a torsion, which is a factor 30 times larger than that for the stream around a purely spherical source.

The value of the torsion extracted from the numerical simulation is τ ≃ 1.08 × 10⁻² kpc⁻¹. This is in excellent agreement with the analytical estimate that we now present. To obtain it, we first note that Eq. (15) for the torsion, with the parameterization of the stream track given in Eq. (A.2) and an almost circular planar projection (semiaxes satisfying a ≃ b), is

$$\begin{array}{c}\hfill {\tau |}_{\mathrm{helix}}\simeq \frac{c\omega }{a^2{\omega }^2+c^2}·\end{array}$$(16)

Here, c is the vertical “velocity”, v_z, in terms of the arbitrary curve parameter, t, and ω = v_φ/ρ is the rotational “angular velocity”, with ρ = a the radius of the orbit in the XY plane. With the setup data for the simulation in Fig. 3, namely c ∼ 80 kpc Gyr⁻¹, a ∼ 30 kpc (we note that the curvature scale is then a⁻¹ ∼ 0.03 kpc⁻¹) and ω ∼ 7 rad Gyr⁻¹, Eq. (16) yields τ|_helix ≃ 0.011 kpc⁻¹, in agreement with the numerical value, 0.0108, extracted from the simulation. The slight oscillation seen in Fig. 3 is due to the ellipticity of the helix (with the oscillation period of about half a Gyr corresponding to half of the actual period of the trajectory because of the reflection symmetry of an ellipse).

Unfortunately for observing this in a physical setup, at a fixed observation time the stream appears quite planar, and it is only in a computer in which we can simulate several time intervals at the same time that we see the helix. This is because of the translation symmetry along the cylinder, which does not allow for a vertical tidal stretching of the stream. On the upside, this makes the extracted torsion very insensitive to the number of bodies in the simulation, from 10² to 10³, since the additional bodies maintain reasonably close orbits in the z − v_z vertical phase space. Providing an initially larger speed dispersion along the OZ axis increases the length of the stream but also blurs it, making its reconstruction more difficult. Thus, we move onto the next case, in which both a sphere and a cylinder (a very idealized toy model of a galaxy and a prolate DM halo) respectively cause the vertical tidal stretching and the helicoidal motion.

3.3. Sphere and cylinder with an initial nonvanishing v_z, 0

To close this section, we combine both types of sources, a sphere (akin to a visible-matter galactic bulb) and a cylinder (mimicking an elongated DM halo). For the sphere, we took the typical mass of a galaxy, M_s ∈ (10⁹, 10¹²) M_⊙ (Busha et al. 2011; Binney & Tremaine 1987), and for the cylinder we used the expression for the linear mass density that we obtained from the asymptotic velocity at large r in the rotation curve, v_rot, of the MW, as was exposed in the previous Sect. 3.

The updated expression for the acceleration of the stars in the stream was obtained by combining Eqs. (B.3) and (B.4); that is,

$$\begin{array}{cc}& z^{″}=-\frac{v_0^2}{{|\mathit{r}|}^3}z\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\mathit{r}}_{\perp }^{″}=-\frac{{|\mathit{v}|}^2}{|{\mathit{r}}_{\perp }{|}^2}{\mathit{r}}_{\perp }-\frac{v_0^2}{{|\mathit{r}|}^3}{\mathit{r}}_{\perp },\hfill \end{array}$$(18)

where v is taken from the galactic rotation velocity when it has flattened out at large r, and v₀ estimated from the visible mass.

We added a small but appreciable (5 km s⁻¹) contribution to the average initial velocity in the z direction, so that $v_{\ast }^i=(220+\mathrm{\Delta }{v}_x^i,\mathrm{\Delta }{v}_{\mathit{y}}^i,5+\mathrm{\Delta }{v}_z^i)$ km s⁻¹, to induce sufficient vertical motion that would generate torsion as per Eq. (15). Additionally, there is a dispersion, Δvⁱ, which depends on the particular star and is randomized from a distribution between 1 and 1000 times the quantity, $\sqrt{G{\displaystyle \sum _j}{m}_{\ast j}/r_0}$, in terms of the mass of the small cluster and the characteristic length scale. There is no particular physical meaning to this initial dispersion except to be able to handle an ensemble of stars that are not all perfectly collinear to test the fitting procedure.

In Fig. 4, we clearly observe traits of the motion around cylinder+sphere sources described in Llanes-Estrada (2021). Along the symmetry OZ axis, a star will describe harmonic oscillations between the two hemispheres due to the Newtonian pull of the spherical part of the distribution acting toward the center (unless it is provided with escape velocity, in which case it will approach an asymptotic trajectory, a helix around the OZ axis). The orbit on the XY plane is not closed due to the additional 1/r force from the elongated DM halo. The net effect in three dimensions can be seen as a precession of the orbital plane around the OZ axis, with the trajectory creating complicated helicoidal patterns.

Fig. 4. Orbits around a source combining both a sphere and a cylinder. The first row represents the trajectory of a point mass with, alternatively, initial vertical velocities vz = 1 and 2 km s−1. The second row extends this to simultaneously show the motion of five stars.

The N-body simulation in Fig. 5 reflects this and clearly shows the appearance of nonvanishing torsion in the stellar stream. From a piecewise helicoidal fit with Eq. (A.2) to the stream at different times, t = 200 Myr, t = 400 Myr, t = 600 Myr, and t = 800 Myr, we obtained different values for the torsion above the value of τ > 10⁻⁴ kpc⁻¹ that we found on streams around a purely spherical source (basically, our irreducible background). The actual value happens to be apparently smaller than that of the purely cylindrical geometry presented in the previous section, but this is a contingent effect due to the initial conditions of the stream, and bears no significance. (The χ²/N_gdl values for each of the visible stream segments, one at each successive time, are 0.35, 0.32, 0.23, and 1.88, respectively). The three-dimensional rendering in the bottom left plot clearly shows the torsion causing a non-osculating plane movement.

Fig. 5. Stellar streams around a combined spherical+cylindrical source. This is a toy model representing the galaxy plus an elongated DM halo. Top row: formation of a stellar stream from an initial orbiting cluster, as in Fig. 4, with positions given at the times indicated in the legend. Lower left: three-dimensional view of the stream compared with a plane containing the origin and the initial stream’s position and velocity, showing that at subsequent times, the stream leaves that plane. Lower right: piecewise helicoidal fit. A nonvanishing torsion can be extracted for each of the short segments at different times, yielding a mean value of τ = 0.22 × 10−3 kpc−1.

3.4. Torsion in a galaxy with a spherical halo and a galactic plane

We wish to have a reference for a minimum torsion that we would consider “normal” in order of magnitude, so that if extensive studies of stellar streams show that their torsion exceeds that level, one could reject the hypothesis of a spherical halo. For this purpose, we propose here a toy model in which the halo is taken to be spherical, but we added a disk component. This adds a vertical (not radial) velocity outside of the galactic plane that points toward it.

The simplest (and coarsest) such model takes the galactic disk as being uniform and infinite. This is a reasonable approximation only for streams that do not elevate too much along the OZ axis; otherwise, it provides an upper bound to a more realistic torsion (since such additional vertical force will always be larger than that of a finite disk, whose effect will fall off with z). In that case, observed torsions above this bound would still entail an incompatibility with a spherical halo that should be studied further.

Therefore, in this minimum-torsion model we take the acceleration caused by the external force to be

$$\begin{array}{cc}& a_i=1,2=-GM\frac{x^i}{{(x^2+y^2+z^2)}^{3/2}}\hfill \\ \hfill & a_3=-GM\frac{z}{{(x^2+y^2+z^2)}^{3/2}}-2\pi G\sigma \mathrm{sign}(z).\hfill \end{array}$$(19)

In this equation, σ is a surface mass density of the disk, in the range of 50 − 100 M_⊙ parsec⁻², which is a usual estimate (Kuijken & Gilmore 1989, 1991) at 8 kpc from the galactic center.

Figure 6 shows the characteristic wobbling of movement near the galactic plane caused by the planar disk, which is qualitatively consistent with Cordoni et al. (2021).

Fig. 6. Simulation of a small tidal stream in the presence of a typical nonspherical galactic source. This is here taken as the galactic plane model of Eq. (19). From top to bottom: three-dimensional view, XZ, and YZ sections.

We can provide an analytical estimate of the torsion following the now-familiar reasoning. Since an instantaneous velocity that is parallel to any of the three coordinate vectors of the cylindrical base, $\{\widehat{\mathit{\rho }},\widehat{\mathit{\phi }},\widehat{\mathit{z}}\}$, will display zero torsion, we took a trajectory combining two of them,

$$\begin{array}{c}\hfill {\mathit{r}}^{\prime }=v_{\phi }\widehat{\mathit{\phi }}+v_z\widehat{\mathit{z}}.\end{array}$$(20)

Multiplying by the acceleration in Eq. (19) we obtained

$$\begin{array}{c}\hfill {\mathit{r}}^{\prime }\times {\mathit{r}}^{″}=\frac{GM\rho }{r^3}(v_{\phi }\widehat{\mathit{z}}-v_z\widehat{\mathit{\phi }})-G{v}_{\phi }(\frac{M}{r^3}z+(2\pi )\sigma \mathrm{sign}(z))\widehat{\mathit{\rho }}.\end{array}$$(21)

To construct the determinant ((r′×r″)⋅r″′) necessary for the torsion, we evaluated the third derivative outside the galactic plane (over which it is undefined),

$$\begin{array}{c}\hfill {\mathit{r}}^{″\prime }=-\frac{\mathit{GM}}{r^2}(\widehat{\mathit{r}}\prime -\frac{2}{r}(\widehat{\mathit{r}}·{\mathit{r}}^{\prime })\widehat{\mathit{r}}),\end{array}$$(22)

which is in the plane given by position and velocity, employing

$$\begin{array}{cc}\hfill \widehat{\mathit{r}}\prime & =\frac{1}{r}(\mathit{v}-(\widehat{\mathit{r}}·\mathit{v})\widehat{\mathit{r}})\hfill \\ \hfill & =\frac{1}{r}(v_{\phi }\widehat{\mathit{\phi }}+v_z\widehat{\mathit{z}}-\frac{z{v}_z}{r^2}(\rho \widehat{\mathit{\rho }}+z\widehat{\mathit{z}})).\hfill \end{array}$$(23)

A slightly tedious but straightforward calculation then yielded

$$\begin{array}{c}\hfill \tau =\frac{-6\pi (M\sigma )(|z|\rho )(v_z{v}_{\phi })}{v_{\phi }^2(M^{2}r+4\pi \sigma M|z|r^2+4{\pi }^2{\sigma }^2{r}^5)+v_z^2{M}^2{\rho }^2/r}·\end{array}$$(24)

The numerator has the mechanical dimensions of a squared momentum, and the denominator of squared momentum times length, yielding the correct 1/L dimensionality of the torsion. Moreover, the structure of the numerator shows that in the presence of a spherical source (M) alone, or a galactic plane (σ) alone, the torsion vanishes as it should. Likewise, both components of the velocity have to be nonvanishing, as in Eq. (15) for the cylindrical source; and the torsion is null both on the galactic plane (z = 0) and on its perpendicular axis through the center of the sphere (ρ = 0). We can then numerically evaluate Eq. (24) to obtain the floor value of the torsion that we should expect to be able to use in the galaxy. Taking into account that the galactic plane is not infinite, so that the elevation z will yield a diminishing multipolar field, it may be that galactic torsions from a spherical halo plus disk are even smaller; what we mean by this estimate is that those streams that may be found with larger values need to be further investigated as they may be teaching us something about the DM halo or about DM inhomogeneities.

Employing z ∼ 1 kpc, ρ ∼ r ∼ 10 kpc, v_z ∼ v_φ ∼ 220 km s⁻¹ (to take the most conservative floor to the torsion), M ∼ 10¹² M_⊙, and σ ∼ 10⁸ M_⊙ kpc² as was already discussed, the denominator of Eq. (24) is dominated by the M² terms, with the σr² correcting M only at the percent level. With these numbers we then find τ ∼ −9 × 10⁻⁴ kpc⁻¹.

We then conclude that torsions of stellar streams below 10⁻³ kpc⁻¹ in our MW Galaxy can be explained without resort to deformed DM halos or exotic phenomena. Of the few streams presently known at large radial distances, some present torsions below this level and are thus of no further interest to this application of the shape of the halos. It is those that reach τ at the percent level that deserve further scrutiny to bear on the halo shape, among the ones known and in future searches for streams.

4. Stellar streams in the Milky Way and their torsion

In this section we finally turn to some of the known stellar streams in the MW. We selected as relevant those found at distances d > 30 kpc from the galactic center, so that the internal structure of the galaxy, such as the disk and spiral arms, produces the minimum possible alteration in the stream. These streams (see Figs. 7 and 8) were extracted from Mateu (2023). The intention in this section is to extract the value of the torsion of the parameterized stream curves with Eq. (5) and check whether or not they vanished.

Fig. 7. Trajectories in galactocentric coordinates for some of the MW outer streams, extracted from Mateu (2023).

Fig. 8. Trajectories for some of the MW outer streams in Mollweide projection, extracted from Mateu (2023).

We took two of the streams out of further consideration; namely, those at Orphan-Chenab and Styx. The reason is that they may be influenced by gravity sources outside the MW. Due to the proximity of the Large Magellanic Cloud (LMC) to our Galaxy, the streams in its periphery in the angular direction of that cloud could suffer alterations due to this additional source of gravity (Conroy et al. 2021; Lilleengen et al. 2022). Moreover, time-dependent perturbations such as those caused by the Magellanic Cloud (Erkal et al. 2019; Vasiliev et al. 2021) make it difficult to interpret the stream track as that of an orbit, sowing further uncertainty (though they can be used in reverse to constrain the time-dependent potential of the perturbing body). This means that even those streams not removed from the analysis could be affected by time-dependent perturbations that we do, at present, not capture. Sometimes the effect can be very pronounced, such as a visible kink in the stream; at other times it can be more subtly manifested by, for example, a misalignment between the stream track and the velocity vectors (Vasiliev et al. 2021).

To obtain the torsion of the curves following the stream we employed both a smooth polynomial parameterization and a (piecewise) elliptical-helicoidal parameterization that we here report. Therefore, we first fit each of the galactocentric Cartesian coordinates tracking the individual streams in the data compilation of Mateu (2023) to order-four polynomial or elliptical-helicoidal parameteric curves. The parameter that describes each curve takes values in the interval t ∈ [0, 1]. The generic parameterizations are relegated to the appendix (see Eq. (A.1) and following). The parameters k or t therein are arbitrary coordinates that can be converted to arc length, which has a clearer geometric significance, by means of

$$\begin{array}{c}\hfill s={\displaystyle \int \sqrt{\sum _{i=1}^3{\left(\frac{\mathrm{d}{x}^i}{\mathrm{d}t}\right)}^2}}\mathrm{d}t;\end{array}$$(25)

the derivatives in Eq. (5) are to be taken with respect to the parameter, t, in general, or s, if a change in variables is effected (the outcome is the same, of course) to obtain the torsion.

To work with the streams in the database¹, we used the GALSTREAM library, and to perform the fits we used a least-squares method with a standard Python installation.

A word about the uncertainty in this extraction is warranted. The data points (x, y, z) of the extracted stream trajectories are quoted without (Δx, Δy, Δz) uncertainties in the original reference (Mateu 2023). Possible large systematic errors, particularly due to distance tracks, which are difficult to pinpoint with accuracy, can at present not be estimated. Thus, the understood uncertainty of our parameteric reconstruction stems in its entirety from our own fit interpreting the stream tracks in terms of Eq. (A.2). Until uncertainties in the track data on which we lean are compiled, the systematics will remain unclear. The reconstructed torsion will only become accurate when higher-precision estimates of the distance gradient along a stream with estimated uncertainties are provided.

After the reconstruction of the parameteric curves, we can obtain the expression for the torsion of each curve, r = (x, y, z), using once more Eq. (5) $\tau =\frac{\mathrm{Det}({\mathit{r}}^{\prime },{\mathit{r}}^{″},{\mathit{r}}^{″\prime })}{||{\mathit{r}}^{\prime }\times {{\mathit{r}}^{″}||}^2}$, where the derivative is taken respect to the t parameter, ′ = d/dt. We analytically express the derivatives in terms of the polynomial-elliptical parameterizations and then evaluate them as function of the t parameter that we used for the fit, taking values from 0 to 1. Because the torsion is parameterization-independent, the torsion can also be given as function of the arclength, τ(s), calculating the derivatives respect to either t or s.

The local torsion along many of the streams in Fig. 9 takes significant values, above the expected O(10⁻⁴) “noise” that we found in the simulations when the streams evolved around spherical sources, as well as the O(10⁻³) or less from simulations that included a galactic plane. Higher values for the torsion have been found for streams such as Cetus, Willka-Yaku, Cetus-Palca, and Sagittarius.

Fig. 9. Local torsion, |τ|, of the MW streams considered in this work, calculated using Eq. (5). The green and red lines and bands of the same shade represent the median value with the interquartile range and the mean value plus or minus one standard deviation, respectively. Stellar streams such as Eridanus and Jet show values of the torsion of the order τ ∼ O(10−3), so their torsion is to be taken as compatible with the influence of a source of the intensity of the galactic plane. The fit to the Sagittarius data was performed dividing the stream into five slices due to its irregular shape.

In fact, the two streams with the largest central value of the torsion (Willka-Yaku and Cetus) do not appear, when eyeballing the actual stream, to have such an important helicity. In the case of Cetus, the large fit uncertainty makes it irrelevant, and the case of Willka-Yaku, whose torsion is still near the distribution of the other streams, deserves being revisited in the future because its shape is very different from our fit functions.

In Table 1 we quantify the torsion in the stream sample: we quote the mean and median values of the torsion among the parameterized pieces of the stream, and also the standard deviation and the interquartile range of the torsion distribution. The high values, τ ∼ O(10⁻²), found for the torsion in several cases are compatible with a nonspherical gravitational source, and also perhaps with the interaction with other gravitational sources different from the overall galactic field, so the observable appears to merit future study.

Because torsion (as curvature) has dimensions of inverse length, we would expect stellar streams to perhaps show an inverse relation with respect to their distance from the galactic center, as is defined in Eq. (5). Irrespectively, in a galaxy such as the MW, the torsion of galactic streams should have a characteristic scale of (10 kpc)⁻¹. As per the discussion around Eq. (15), where we established $\tau \in [-\frac{1}{\rho },\frac{1}{\rho }]$, our selection of streams at 30 kpc or more means that we would consider torsion values of order 0.01 in units of inverse kiloparsecs to be sizeable and very different from zero, and Table 1 shows several such. Also, as per the discussion below Eq. (24), those with τ above 0.001 kpc⁻¹ could perhaps carry interesting information about the DM distribution.

Turning to the data, the torsions that we seem to observe in the MW streams show mean values above 0.001 that could give an insight into the DM distribution, but also information about the interaction of the MW with other gravitational sources. Several streams show a significant value for the torsion: only Eridanus and Jet seem to have torsions compatible with being seeded by the galactic plane, as small as τ ∼ O(10⁻³).

From the information we obtained from the streams in the database, we cannot explore specific details about the shape of the gravitational source, beyond establishing its probable nonsphericity.

5. Conclusions and outlook

The problem of galactic rotation curves suggests that galaxies are surrounded by significant amounts of DM, and the overall shape of these sources is yet to be ascertained. Whereas spherical DM distributions around galaxies have to be fine-tuned to explain the flatness of rotation curves, a cylindrical (or, generally, prolate) DM source can naturally explain that flattening. This avoids the fine-tuning of spherical DM halos to precisely follow the 1/r² fall-off for a large swath of r values. Observables inside the galactic plane cannot however distinguish between spherical (though fine-tuned) and cylindrical or elongated gravitational sources, but out-of-galactic-plane information could provide new strong discriminants.

The stellar streams around the MW have been extensively investigated for a while now, and are still nowadays a relevant subject of research. The track left by a stream, once extracted, can be used as the best proxy we have for a long-range trajectory of size commensurate to a galaxy, with all the attending dynamical uncertainties that this identification causes (Bovy 2014), and thus a tool to infer the geometry of the impacting gravitational sources. Stream tracks can be characterized by their torsion according to Eq. (5); around a central potential, orbits move in a plane and are expected to be torsionless (see Fig. 1). In addition, test masses around cylindrical sources are expected to follow helical orbits in which the torsion is nonzero (see Fig. 3) if given a vertical velocity. Another approach is to consider an ellipsoid-shaped halo, which is not perfectly cylindrical but rather elongated. The expected orbit of the stream components would arise from the combination of the orbits around central potentials and the helical orbits around cylinders, as is seen in Fig. 4. Torsion is then an observable directly tailored to assess the prolateness of a DM halo, as opposed to a multiparameter reconstruction that, though sensitive to the halo shape (Koposov et al. 2023), is affected by a degeneracy with many other features.

The streams of the MW have been a subject of research for a considerable time span, and many of the objects that constitute these streams have been catalogued. From the reconstruction of the stream tracks, we infer the torsion caused by the gravitational source in Fig. 9. In this work we only consider those streams that seem to be far enough away from the galactic center to (1) avoid large effects from the baryonic component of the galaxy and (2) have a bird’s eye view of the DM halo from outside a large fraction thereof. From the extraction of the torsion we see that it is non-negligible in some of the streams considered.

From our evaluation of the torsion we do not dare favor one or another interpretation of the DM halo shape in the view of current data; this article should be seen as a proposal for a new observable, τ, and a first exploratory study. We do find streams with significant torsion, which is encouraging and suggests that further scrutiny could eventually inform us about the shape of the halo. We have paid no attention to the sign of τ, since it only distinguishes left- from right-handed tracks. That opens another interesting study in itself as, with the predominant galactic rotation determining a sign, the preference of a given handedness would inform an asymmetry in the vertical motion. But this would require many more statistics and take us far afield at present.

The most trustworthy streams should be those stemming from kinematically cold globular clusters (low-velocity dispersion), as those present the best-defined tracks. In future observational work it might be interesting to actively seek streams that show both vertical motion (along the axis perpendicular to the MW plane) and also azimuthal motion around that axis, as those with large v_z and v_φ will be most sensitive to the torsion. Should those streams show trajectories that are compatible with lying on a plane (zero torsion), a spherical halo would be preferred. Should they however appear helicoidal, with non-negligible torsion, they would be pointing to an elongated DM halo.

Streams that may be detected in nearby galaxies carry the same information about their respective halos.

However, for galactic-halo population studies with streams, which will require a 10-MPc reach, three-dimensional reconstruction of the streams appears more difficult; in that case, only the projected curvature is measurable (Nibauer et al. 2023), but not the torsion, and what can be learnt about the DM halo hangs on statistical analysis to reduce the impact of degeneracies due to the projection.

Finally, other observables can bear on the overall shape of the halo, and we are investigating the shape sensitivity of the gravitational lensing of both electromagnetic and gravitational radiation.

Acknowledgments

We thank our engineer David Fernández Sanz for maintaining a suitable computing environment at the theoretical physics departmental cluster, and Pablo Blanco-Mas for insightful conversations. Financially supported by spanish “Ministerio de Ciencia e Innovación: Programa Estatal para Impulsar la Investigación Científico-Técnica y su Transferencia”, Grants PID2021-124591NB-B-C41 and PID2022-137003NB-I00: funded by MCIN/AEI/10.13039/501100011033 and by the EU-ERDF, as well as U. Complutense de Madrid under research group 910309 and the IPARCOS institute.

References

Appendix A: Track fits

A simple parameterization (x(t),y(t),z(t)) that we employed is a piecewise fourth-order polynomial fit reading:

$$\begin{array}{c}\hfill x^i(t)=a_i{t}^4+b_i{t}^3+c_i{t}^2+d_{i}t+e_i\mathrm{for}\mathrm{i}=1,2,3.\end{array}$$(A.1)

It is futile trying to relate t to a Newtonian time since, not knowing a priori the dynamics of the system nor at which time a star was at what position along its trajectory: only the stream’ s instantaneous (present) geometry is known with certainty, and therefore an arbitrary parameter, t (or the arc length after computing it), should suffice. Another simple yet useful parameterization is that of an elliptical helix,

$$\begin{array}{cc}\hfill \mathbf{x}(t)& ={\mathbf{x}}_0+(acos(\omega t+\varphi ),bsin(\omega t+\varphi ),ct)\hfill \\ \hfill {\mathbf{x}}^{\prime }(t)& =(-a\omega sin(\omega t+\varphi ),b\omega cos(\omega t+\varphi ),c)\hfill \\ \hfill {\mathbf{x}}^{″}(t)& =(-a{\omega }^{2}cos(\omega t+\varphi ),-b{\omega }^{2}sin(\omega t+\varphi ),0)\hfill \\ \hfill {\mathbf{x}}^{″\prime }(t)& =(a{\omega }^{3}sin(\omega t+\varphi ),-b{\omega }^{3}cos(\omega t+\varphi ),0)\hfill \end{array}$$(A.2)

with (a, b, c, ω, ϕ, x₀) the parameter set, such that for t = 0 the curve cuts the XY plane with z = 0 with the initial phase, ϕ, in the x − y rotation, c controlling the advance of the helix, and a, b the elliptic projection on the XY plane. Should the rotation axis also need to be chosen to be different from OZ, an additional orthogonal matrix, O(φ, θ, ξ), parameterized by three Euler angles, should be introduced as x → Ox. We do not spell this out here for conciseness, and assume that OZ is an axis of symmetry for the DM halo, coinciding with the galactic rotation axis.

We then need to distinguish among two fit classes. First, when fitting simulated data, there is no obstacle to choose the coordinate axes so that OX lies along the initial visual from the galactic centre to the globular cluster that will become the stream.

In that case, the offset, x₀ in Eq. (A.2) is unnecessary. The elliptical parameterization is then, eliminating time,

$$\begin{array}{cc}& x(z)=Acos(\frac{\omega }{C}(z-z_0)+\varphi )+x_0\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& y(z)=Bsin(\frac{\omega }{C}(z-z_0)+\varphi )+y_0.\hfill \end{array}$$(A.4)

The length parameters, A and B, take values proportional to the distance of the stream to the galactic center in kpc. C is a vertical “velocity” and ω an “angular velocity” in the XY plane; but, because of the arbitrariness of the parameter t, they are only relevant in the combination, ω/C, that indicates the angle advanced per unit height in the corkscrew motion characterizing torsion.

(If one insists on thinking of the orbit of a single star instead of the track of a stream, then the angular velocity is of the order of magnitude $\omega \sim \frac{v_{\perp }}{r_{\perp }}\sim \frac{220\mathrm{km}/\mathrm{s}}{{10}^4\mathrm{Pc}}\times 3.15{10}^{16}\mathrm{s}/\mathrm{Gyr}\times \frac{1}{3.1{10}^{13}\mathrm{km}/\mathrm{Pc}}$, resulting in ω ∼ 22 rad/Gyr).

The second type of fit involves real data. Then, the osculating plane at the start of each piece of the orbit need not be aligned with the galactic Cartesian axes on account of initial conditions and precession in a noncentral field. We thus need a minimum of one additional angle, Θ, for each fit, to orient that starting plane,

$$\begin{array}{cc}& x\prime =x·cos\mathrm{\Theta }-y·sin\mathrm{\Theta }\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& y\prime =x·sin\mathrm{\Theta }+y·cos\mathrm{\Theta }.\hfill \end{array}$$(A.6)

Additionally, real streams may need an offset of x₀ = (x₀, y₀, z₀), bringing the total number of fit parameters to nine.

Once the parameterization has been chosen, we need to decide what is the optimal track through a cloud of them in a simulation. One strategy is to employ a squared-distance minimization. The function to be minimized is the following sum running over each of the stars in the cloud,

$$\begin{array}{c}\hfill D^2(a,b,c,t_0,\omega ,t_i)={\displaystyle \sum _{i=1}^N}{({\mathbf{x}}_i-\mathbf{x}(t_i))}^2,\end{array}$$(A.7)

where t_i is the point along the curve for which its distance to the ith point is minimized. Taking partial derivatives with respect to t_i or a_j results in the set of equations

$$\begin{array}{c}\hfill {\displaystyle \sum _{i=1}^N}({\mathbf{x}}_i-\mathbf{x}(t_i))·\mathbf{v}(t_i)=0={\displaystyle \sum _{i=1}^N}({\mathbf{x}}_i-\mathbf{x}(t_i))·\frac{\partial \mathbf{x}(t_i)}{\partial a_j}\end{array}$$(A.8)

which constrain the wanted parameters, a_j. The first of these equations is the condition of orthogonality of the curve’s tangent to the visual to the point. Small stretches of a stellar stream can be fit with the functional form of Eq. (A.2) for a straightforward extraction of the torsion. The torsion for this elliptical helix is then much simplified and becomes a t-independent constant for a circular helix with a = b, as in Eq. (16),

$$\begin{array}{c}\hfill \tau =\frac{abc\omega }{b^2{c}^2{sin}^2(\omega t+\varphi )+a^2{c}^2{cos}^2(\omega t+\varphi )+a^2{b}^2{\omega }^2}.\end{array}$$(A.9)

Appendix B: N-body equations of motion

Finally, we specify the (simplest available) numerical method employed to simulate an N-body stellar stream in section 2.3. The positions of the stellar objects composing it are updated in Cartesian coordinates. For this, we used Euler’s method, with the time step Δt = t_f/N_t, with the velocity updated via a once-improved Euler step,

$$\begin{array}{cc}& x_{j+1}^i=x_j^i+\mathrm{\Delta }t{v}_j^i\mathrm{for}\mathrm{i}=1,2,3\hfill \end{array}$$(B.1)

$$\begin{array}{cc}& v_{j+1}^i=v_j^i+\frac{1}{2}\mathrm{\Delta }t{f}^i({\mathbf{x}}_j+\frac{1}{2}h{\mathbf{v}}_j)\mathrm{for}\mathrm{i}=1,2,3,\hfill \end{array}$$(B.2)

where fⁱ is the function yielding each component’s acceleration.

The acceleration was calculated at each step from standard formulae such as those for a spherical central source,

$$\begin{array}{cc}\hfill a_{i=1,2,3}=& -GM\frac{x^i}{{(x^2+y^2+z^2)}^{3/2}}\hfill \\ \hfill & -{\displaystyle \sum _{j=1}^{N-1}}G\frac{m_j{x}^i}{{({(x-x_j)}^2+{(y-y_j)}^2+{(z-z_j)}^2)}^{3/2}}.\hfill \end{array}$$(B.3)

The first line of this expression is the acceleration under a central spherical source, and the second the force among the stellar-stream stars (which is too weak to avoid the tidal stretching). If we turn to a cylindrical source, the force external to the stream in Eq. (B.3) needs to be replaced, so that

$$\begin{array}{cc}\hfill a_{i=1,2}=& -2G\lambda \frac{x^i}{x^2+y^2}-{\displaystyle \sum _{j=1}^{N-1}}G\frac{m_j{x}^i}{{({(x-x_j)}^2+{(y-y_j)}^2+{(z-z_j)}^2)}^{3/2}}\hfill \\ \hfill a_{i=3}=& -{\displaystyle \sum _{j=1}^{N-1}}G\frac{m_j{x}^i}{{({(x-x_j)}^2+{(y-y_j)}^2+{(z-z_j)}^2)}^{3/2}}.\hfill \end{array}$$(B.4)

At last, we have simulated both cylinder+sphere and plane+sphere background fields, with respective forces given in Eq. (18) and (19).

All Tables

All Figures

Fig. 1. Motion of a star ensemble around around a spherical source. Top row: simulated movement of a 200-star cluster around a spherical gravitational source with parameters akin to those of the MW, in the three planar projections, extracted at times t = 0 Myr, t = 200 Myr, t = 400 Myr, t = 600 Myr, and t = 800 Myr. Lower left: three-dimensional view with the containing plane. The initial cluster has a radius of 2 kpc and is at a distance of 30 kpc from the galactic origin of coordinates. Member stars have slightly different velocities, v* = 220 kpc ± Δv, due to the additional random Δv. (All axes are in kpc). Lower right: extracted torsion as a function of time, from an analytical helicoidal fit to the resulting stream with χ2/Ngdl = 0.86, which yields a value of τ ≃ 3.5 × 10−4 kpc−1 (a one-star orbit would return exactly zero torsion).

In the text

	Fig. 2. Spread of a 100-body stellar cluster after t = 0, 125, 250, 375, and 500 Myr. The resulting stream, with an initial vertical speed of vz = 50 km s−1, at different times clearly shows an upwardly advancing helix moving around the z axis.
In the text

Fig. 3. Motion of a star ensemble around around a perfectly cylindrical toy-model source. The initial velocity is vz = 90 km s−1, with the same conventions as in Fig. 1. Upper plots: a stream forms along a spiral trajectory. From a helicoidal fit to the cloud of points, as needs to be done on real data, with χ2/Ngdl = 1.16 we find a value for the torsion of τ ≃ 1.1 × 10−2 kpc−1 (lower plot). This is in agreement with the analytical estimate in the main text below Eq. (16).

In the text

	Fig. 4. Orbits around a source combining both a sphere and a cylinder. The first row represents the trajectory of a point mass with, alternatively, initial vertical velocities vz = 1 and 2 km s−1. The second row extends this to simultaneously show the motion of five stars.
In the text

Fig. 5. Stellar streams around a combined spherical+cylindrical source. This is a toy model representing the galaxy plus an elongated DM halo. Top row: formation of a stellar stream from an initial orbiting cluster, as in Fig. 4, with positions given at the times indicated in the legend. Lower left: three-dimensional view of the stream compared with a plane containing the origin and the initial stream’s position and velocity, showing that at subsequent times, the stream leaves that plane. Lower right: piecewise helicoidal fit. A nonvanishing torsion can be extracted for each of the short segments at different times, yielding a mean value of τ = 0.22 × 10−3 kpc−1.

In the text

	Fig. 6. Simulation of a small tidal stream in the presence of a typical nonspherical galactic source. This is here taken as the galactic plane model of Eq. (19). From top to bottom: three-dimensional view, XZ, and YZ sections.
In the text

	Fig. 7. Trajectories in galactocentric coordinates for some of the MW outer streams, extracted from Mateu (2023).
In the text

	Fig. 8. Trajectories for some of the MW outer streams in Mollweide projection, extracted from Mateu (2023).
In the text

Fig. 9. Local torsion, |τ|, of the MW streams considered in this work, calculated using Eq. (5). The green and red lines and bands of the same shade represent the median value with the interquartile range and the mean value plus or minus one standard deviation, respectively. Stellar streams such as Eridanus and Jet show values of the torsion of the order τ ∼ O(10−3), so their torsion is to be taken as compatible with the influence of a source of the intensity of the galactic plane. The fit to the Sagittarius data was performed dividing the stream into five slices due to its irregular shape.

In the text