Constraining the Gravitational Potential from the Projected Morphology of Extragalactic Tidal Streams

The positions and velocities of stellar streams have been used to constrain the mass and shape of the Milky Way's dark matter halo. Several extragalactic streams have already been detected, though it has remained unclear what can be inferred about the gravitational potential from only 2D photometric data of a stream. We present a fast method to infer halo shapes from the curvature of 2D projected stream tracks. We show that the stream curvature vector must point within 90 deg of the projected acceleration vector, in the absence of recent time-dependent perturbations. While insensitive to the total magnitude of the acceleration, and therefore the total mass, applying this constraint along a stream can determine halo shape parameters and place limits on disk-to-halo mass ratios. The most informative streams are those with sharp turns or flat segments, since these streams sample a wide range of curvature vectors over a small area (sharp turns) or have a vanishing projected acceleration component (flat segments). We apply our method to low surface brightness imaging of NGC 5907, and find that its dark matter halo is oblate. Our analytic approach is significantly faster than other stream modeling techniques, and indicates what parts of a stream contribute to constraints on the potential. The method enables a measurement of dark matter halo shapes for thousands of systems using stellar stream detections expected from upcoming facilities like Rubin and Roman.


INTRODUCTION
Low surface brightness imaging has revealed intricate tidal features in the halos of external galaxies (e.g., Shang et al. 1998;Ibata et al. 2001;Martínez-Delgado et al. 2009;Crnojević et al. 2016;Morales et al. 2018;Kado-Fong et al. 2018;Giri et al. 2023;Martínez-Delgado et al. 2023a). Similar tidal features have been precisely measured in the Milky Way, where the field of tidal debris can be studied on a star-by-star basis (e.g., Majewski et al. 2003;Odenkirchen et al. 2003;Belokurov et al. 2006;Koposov et al. 2010;Hasselquist et al. 2019;Malhan et al. 2022;Varghese et al. 2011;Koposov et al. 2022). The morphology of these systems is vast, ranging from streams, to shells, and continuous deformations of the two. These features originate from tidal stripping, when a lower mass satellite galaxy accretes onto the more massive host.

jnibauer@princeton.edu
The phase-space distribution (that is, the space of positions and velocities) of tidal debris is sensitive to the underlying gravitational potential of the host galaxy. Stellar streams, for instance, roughly trace orbits in the host gravitational potential, and provide a sensitive probe of both the baryonic and dark matter components of a galaxy. An extensive range of methods have been devised to measure the gravitational potential using tidal debris. These range from trial orbit integrations (e.g., Johnston et al. 1999;Fardal et al. 2006;Koposov et al. 2010;Price-Whelan et al. 2014) to fully generative models for stream formation (e.g., Küpper et al. 2012;Fardal et al. 2015;Bonaca et al. 2014;Gibbons et al. 2014). In the Milky Way some of these methods have been applied to real stellar streams, where kinematic information and distance tracks can be utilized to place tight constraints on the potential (e.g., Malhan & Ibata 2019;Koposov et al. 2022).
From an information theoretic perspective, Bonaca & Hogg (2018) found that kinematically cold streams provide a localized constraint on the galactic potential. In-deed, in Nibauer et al. (2022) it is shown that local galactic accelerations can be recovered from streams directly, without explicit reference to any potential model. The results from these analyses both imply that there must be some amount of information about the potential encoded in the geometry of tidal features (e.g., the light blue contours in Fig. 3 of Hogg 2018 andEq. 6 of Nibauer et al. 2022), though it has remained unclear what properties of the potential can actually be inferred from projected stream morphology. Using the generative "streakline" technique, Pearson et al. (2015) found that a spherical halo reproduces the curvature of Palomar 5 better than a triaxial halo. This method requires one to model the full time-evolution of the progenitor, introducing significant degeneracies with initial conditions, mass-loss rates, and integration times. If the object of interest is the shape of the gravitational potential, avoiding these degeneracies would be optimal. Especially in the context of external galaxies, attempting to model the full time-evolution of a progenitor from only the projected shape of a tidal feature introduces far more modeling degrees of freedom than the data can realistically constrain.
For the thousands of extragalactic tidal features expected to be discovered with upcoming surveys like The Rubin Observatory (Ivezić et al. 2019), Euclid (Laureijs et al. 2011), and The Nancy Grace Roman Space Telescope (Roman; Spergel et al. 2013), neither distance tracks nor velocity information will be available for most systems. In this limited data regime, modeling attempts will be hindered by degeneracies with projection effects and the unobserved phase-space dimensions. For instance, action-angle based methods also require measurements of the 6D phase-space distribution of the tidal feature (Sanders & Binney 2013a;Sanderson et al. 2015), severely limiting applicability to external galaxies. While a single extragalactic stream will not be measured to the same level of detail as in the Milky Way, measurements of the low surface brightness universe provide a unique opportunity to constrain dark matter halo properties at the population level. In this context, automated methods will be required to characterize the ensemble of tidal debris and translate observations to physical constraints on the dark matter distribution in external galaxies.
Recently, Pearson et al. (2022b) modeled the tidal feature surrounding the galaxy Centaurus A using the "particle spray" method (Fardal et al. 2015). This method provides a semi-analytic prescription for streamformation, and works by releasing test particles from the Lagrange points of a progenitor throughout a trial orbit integration. In the limited-data scenario of an external galaxy, Pearson et al. (2022b) demonstrates the significant degeneracies in modeling the projected morphology of an extragalactic stream without velocity information. Indeed, the shape of a stream is degenerate with mass enclosed, halo shape parameters, velocity, and projection effects. Pearson et al. (2022b) showed that if a single radial velocity is measured along the stream, some of these degeneracies (the halo mass, in particular) can be broken.
Other attempts to model the potential of external galaxies using extragalactic tidal features have relied on full N −body models for stream (or shell) formation (Fardal et al. 2013;Foster et al. 2014;Bílek et al. 2022;Martínez-Delgado et al. 2023b), additional applications of the "particle spray" technique (Amorisco et al. 2015;van Dokkum et al. 2019). In total, generative stream models have a large number of latent nuisance parameters characterizing the time-evolution of the progenitor. Without velocity information or a distance gradient along the tidal feature, 4 out of the 6 simulated phase-space dimensions cannot be compared to the measured data directly. Only the on-sky positions can be directly compared, requiring one to marginalize over the unobserved phase-space dimensions. This marginalization step introduces a range of additional degeneracies in attempting to reconstruct the potential.
Here, we aim to address what can possibly be learned about the host potential from only the projected track of a tidal stream. To this end, we connect the projected morphology (i.e., curvature) of a tidal feature to the gravitational acceleration field of its host galaxy. The resulting method for potential reconstruction is analytic, adaptable to any static potential model, and does not require trial orbit integrations. We demonstrate the performance of our method on N −body simulations of streams, and discuss what morphological properties of streams will be most useful for constraining the potential.
The paper is organized as follows. In §2 we introduce our modeling approach. In §3 we devise a likelihood. In §4 we apply the method to a series of N −body streams generated in potentials with different shape properties and viewing angles. In §5 we present preliminary results on the stream surrounding the galaxy NGC 5907. We discuss our results and future prospects in §6, and conclude in §7.

MODELING FRAMEWORK
We now develop a modeling framework that connects the projected track of a stream to the underlying gravitational potential of its host galaxy. We first consider the intuitive case of orbits, where the local curvature of Potential C Figure 1. A single orbit and acceleration vectors generated from different potential models. In each panel the same 2D projection of a 3D orbit is shown in black. The curvature vectors, which point orthogonal to the projected orbit, are the black arrows. The gray and red arrows depict acceleration vectors for three different potential models: A, B, and C, all projected onto the x − y plane. Which panel shows the correct acceleration vectors for the given orbit? The correct potential must produce acceleration vectors that fall within 90 deg of the curvature vector. Otherwise, the orbit would not curve in the observed direction. The correct potential is therefore Potential B, since potentials A and C do not satisfy this constraint over the full length of the orbit.
an orbit is related to the gravitational acceleration field ( §2.1). We then introduce a coherence condition which, if satisfied along the stream track, makes the same orbitbased framework applicable-even if the track of the stream does not trace a single stellar orbit ( §2.2). Our main assumptions are summarized below, and we refer the reader to the corresponding sections for further discussion.
• Stream Track and Coordinate System ( §2.2.1): The observed tidal feature has an elongated morphology in the plane of the sky which can be represented as a twice-differentiable track.
• Coherence Condition ( §2.2.2): The acceleration direction of a stream "fluid-element" in projection must not oppose the observed curvature direction of the projected stream-track. Otherwise, the orbits of stars populating the stream would completely decouple from the observed track. We derive an exact inequality that characterizes when this coherence condition is satisfied, and show that it does not require the stream to trace a single stellar orbit.
• Absence of Curvature ( §2.2.3): If the observed stream has a linear segment with zero projected curvature, we assume that the streamsegment is not accelerating perpendicular to its elongated axis in the plane of the sky.

Intuitive Picture: A Single Orbit
In each panel of Fig. 1 we plot the same 2D projection of an orbit generated in a 3D potential Φ(x, y, z). The red and gray vector fields represent accelerations under three logarithmic potential potential models (A, B, C) with different flattening parameters. The vector field for the correct potential is shown in the middle panel (Potential B). The black arrows indicate unit curvature vectors, which characterize the local concavity of the orbit. Unit curvature vectors are calculated locally, by differentiating the track of the orbit (black curve) in the x − y plane. For illustration purposes we assume a fixed trial distance z 0 when plotting all acceleration vectors, though our formal analysis considers a range of possible distances. Fig. 1 illustrates a connection between the curvature of the orbit in projection and the underlying acceleration field. Namely, under the correct potential model, acceleration vectors in the x − y plane evaluated at the location of the orbit always have a component along the local curvature vector. Otherwise, the orbit would not have curved in the observed direction. We mathematically derive this to be the case in the following section, though this simple principle forms the core of our analysis. In general, the projection of the acceleration vector along an orbit must always point within 90 deg of the projected curvature vector. If the orbit has no curvature in projection then we are viewing the orbital plane edge on, and all accelerations are along the line of sight.
If the 3D track of the orbit is known along with its speed, then the acceleration can be calculated directly (Nibauer et al. 2022). If the speed is not known and the orbit is assumed to be circular, then the correct acceleration vector must point in the same direction as the curvature vector. However, if only a random 2D-projection of the orbit is observed, then the topology of the orbit can be obscured by projection effects. For instance, circles can appear elliptical depending on their orientation with respect to the observer, and ellipses can appear more circular. Without knowing the speed, and without imposing additional priors (e.g., about the potential, circular or radial orbits, projection effects), we can only make the minimal assumption that the projected acceleration field along the orbit must point within 90 deg of the orbital curvature vector in the plane of the sky.
This principle can be used to rule out acceleration vectors that could not have possibly generated the observed curvature. For a given potential model, the direction of the acceleration field in a region is sensitive to shape parameters, e.g., how the mass distribution is oriented, flattened, or elongated. We devise a likelihood to perform this inference in §3. Note that the direction of acceleration vectors is not sensitive to the total mass of the system, since uniformly scaling the potential by a constant changes only the magnitude of accelerations. Because of this, orbital curvature alone cannot be used to determine an absolute mass estimate. However, varying the mass ratio between distinct potential components with different shapes (e.g, a disk and halo) does change the direction of accelerations, and therefore orbital curvature. We later show that limits on disk-tohalo mass ratios can be placed, assuming that the disk's potential is aligned with its stellar mass in shape and extent ( §4.2.3).
There are several caveats to the discussion above. For instance, streams are not orbits (Sanders & Binney 2013b). Instead, streams are populated by stars on an ensemble of orbits, with a typically small, albeit nonzero energy gradient along the stream (e.g., Johnston 1998;Johnston et al. 2001). Furthermore, the gravitational potential of the Milky Way and external galaxies is not expected to be static. This can affect the orbit of stars in a stream (e.g., Erkal et al. 2019;Shipp et al. 2021;Lilleengen et al. 2023). We next consider under what conditions the curvature of a stream can be linked to the acceleration field of its host. The resulting principle is the same as the orbit-based picture we have already discussed, though in §2.2 we show that this picture can be extended beyond the case of a single orbit.

Formal Derivation: Connecting Stream Curvature to the Gravitational Acceleration Field
In the remainder of this section we model streams in projection, and consider what can be learned from the S t r e a m T r a c k Figure 2. Time-evolution of a stream-element (gray) from the present day (t0) to a future time (t2). As the stream element evolves, its basis vectors are set by the minimum projected distance between the element's position and the stream track (black curve, x(γ)). Tangent (T ) and curvature (κ) unit vectors are set by the morphology of the stream track.
observed shape of a stream without making restrictive assumptions about its speed or shape along the line-ofsight. We show that the curvature-acceleration connection does not require the stream track to trace an orbit (i.e., an isoenergy curve).

Defining the Stream Track in On-Sky Coordinates
In this section we define the stream track in on-sky coordinates, and decompose the velocity of a stream "fluid element" into motion along the track and perpendicular to the track in the plane of the sky. We then write the acceleration of the fluid element in this frame, drawing a connection between projected stream-curvature and projected accelerations.
We consider streams that have a well-defined track, e.g., those resembling a series of loops rather than shells. We work in a coordinate system centered on the external galaxy where x − y is the plane of the sky and the z-axis is along the line-of-sight. The origin is at the center of the host galaxy.
Position along the projected stream-track is described by the vector-valued function x(γ) ∈ R 2 , where γ is a dimensionless phase-parameter encoding position along the stream. For simplicity, we work in physical coordinates so that x(γ) has units of length. This is not strictly necessary, since our acceleration analysis is scalefree. Instead, x(γ) could also represent pixel-based coordinates, or on-sky position angles (e.g., ϕ 1 , ϕ 2 ).
To characterize the morphology of a stream, we describe the stream track with a series of tangent vectors and orthogonal curvature vectors. The unit tangent vector is defined asT and the unit-curvature vector iŝ These vectors are orthogonal, satisfyingT ·κ = 0. An illustration of the coordinate system is shown in Fig. 2. The stream-track is depicted by the solid black curve, and the time-evolution of a stream element is shown in gray. Importantly,κ is invariant to the circulation direction along the track: it does not matter which tail is "leading" or "trailing" in our analysis. We now consider the motion of the gray stream element in Fig. 2, and under what conditions its acceleration vector is constrained by the curvature vectorκ. At the present epoch, t 0 , the velocity of the element in the plane of the sky is where v T and v κ are the tangent and perpendicular velocity components, respectively. For simplicity, we assume that v T > 0. 1 If the stream traces an orbit, then the stream element will evolve purely along the stream track with v κ = 0. More realistically, the stream element could have v κ ̸ = 0. As the stream element evolves in time, the vectorsT andκ are chosen by the minimum projected distance between the present-day stream track and the position of the stream element (dashed lines in Fig. 2). As a stream-element time-evolves, v planar and its components will change according to the planar acceleration where "over-dot" indicates a time-derivative, and a T and a ⊥ are the tangential and perpendicular accelerations in the plane of the sky, respectively. By construction, dT /dt will point alongκ, and dκ/dt will point along −T . The latter is sensitive to the circulation direction, so we do not use a T to constrain the potential. In §2.2.2, we discuss how and when Eq. 4 can be used to constrain properties of the potential using the morphology of a stream.
In subsequent sections we constrain a planar through a gravitational potential Φ m (x, y, z) with model parameters m. The coordinates (x, y) are on-sky position coordinates, and z is the unknown line-of-sight coordinate. The acceleration is related to the potential through its spatial gradient where a = a x , a y , a z . In these coordinates, the planar acceleration vector is

Connecting Curvature to Accelerations
We now use Eq. 4 to express a ⊥ in terms of the observable stream curvature,κ. The unit curvature vector is related to the time derivative of the unit tangent vector throughκ where γ is the dimensionless phase-parameter, parameterizing the track of the stream. Eq. 7 allows us to write an expression for a ⊥ (Eq. 4) in terms of the observed unit-curvature vector along the stream track, namely, By construction, a ⊥ points along ±κ: the sign is determined by the relative size between the two terms in parentheses. We use Eq. 8 to write a coherence condition, which specifies precisely when a ⊥ point along +κ. Namely, When Eq. 9 is satisfied, a ⊥ points along +κ and the observed stream-curvature provides a constraint on the acceleration direction. We now discuss whether we expect Eq. 9 to hold for streams. In static potentials, stream stars tend to sort out by their energies (see, e.g., Johnston et al. 2001), so that proximate regions of the stream are populated by stars on similar orbits without large oscillations away from the stream track. This behavior is seen in actionangle coordinates, where stars tend to sort out in angle by their frequency difference relative to the progenitor (Sanders & Binney 2013b). Then, at a given phase, the stream is populated by an ensemble of stars with locally similar frequencies. In this case, we do not expect stars to suffer large oscillations away from the stream track.
Then, Eq. 9 would be satisfied. This property of streams has been utilized to measure the solar velocity (Malhan et al. 2020), and to discover new streams in the Milky Way (Ibata et al. 2021).
Eq. 9 is not a particularly strict requirement: iḟ v κ << 0, the stream track would be totally decoupled from its trajectory, and the stream would begin the process of phase-mixing again from its new configuration. The mere existence of a stream implies a level of coherence, though time-dependent potentials can lead to perpendicular motion to the stream track (Erkal et al. 2019;Vasiliev et al. 2020). This does not necessarily violate Eq. 9, provided thatv κ is bounded below by −v T ∥dT /dt∥.
As long as Eq. 9 is satisfied, the stream-element's orbital trajectory has a perpendicular acceleration component alongκ. This means whereκ can be estimated from data directly (Eq. 7). If the stream track were to trace an orbit, then v κ = v κ = 0. However, the connection between perpendicular accelerations and stream-curvature (Eq. 10) holds even forv κ ̸ = 0. Thus, the track of the stream does not need to actually trace an orbit in order to connect streamcurvature to perpendicular accelerations.
We now cast Eq. 10 as a geometric constraint, which will be useful in interpreting our results. From Eq. 4, a ⊥ = (a planar ·κ)κ. The unit vector of a ⊥ is then where sgn(x) is the sign function, with sgn (x > 0) = +1, sgn (x = 0) = 0, and sgn (x < 0) = −1. The symbol θ(Φ) ∈ [−π, π] represents the angle between the planar acceleration, a planar , and the curvature vector to the stream track,κ. Planar accelerations depend on the underlying gravitational potential, Φ. The angle θ (Φ) is illustrated in Fig. 3, where the cyan arrows depict planar acceleration vectors along the stream track, and the green arrows are unit curvature vectors along the track (κ). When the coherence condition expressed in Eq. 9 is satisfied, a ⊥ points alongκ and sgn cos θ(Φ) = 1. This implies |θ (Φ) | < π/2, or that the angle between a planar andκ will fall within π/2 radians.T κ ∝ dT dγ a planar θ Figure 3. Illustration of a stream track (black curve), its tangent vectors (T ), and curvature vectors (κ). The planar accelerations along the track are the cyan vectors, a planar . The angle betweenκ and a planar is denoted with θ, which depends on the underlying gravitational potential, Φ.
Intuitively, Eq. 12 requires that a change in the unit tangent vector along the projected stream track must have been induced by an acceleration perpendicular to the track. The direction of this planar acceleration must fall within 90 deg of the curvature vector to produce the stream's morphology (Eq. 12). In the subsequent sections, we explore what limits this basic constraint, combined across the length of a stream, can place on the underlying gravitational potential.
The analysis presented here uses only the curvature direction and not curvature strength to constrain planar accelerations. While the latter does depend on the acceleration field, it also depends on the velocity of the stream (e.g., ∥dT /dt∥). In this work we aim to address what can be inferred from only the track of the stream without restrictive velocity assumptions. However, priors can be placed on the velocity of each stream segment to provide even tighter constraints. We leave this to future work.

Special Case: Zero Curvature
We have previously assumed that the projected stream track is non-linear. However, when ∥dT /dγ∥ ≈ 0, the stream-track is locally linear in the x − y plane. This could correspond to an inflection point where the stream changes concavity, or a radial orbit. We take a lack of curvature to indicate that the perpendicular acceleration component, a ⊥ , is equal to zero. We justify this below.
If we assume that the coherence condition in Eq. 9 is satisfied, then at an inflection point where the stream changes concavity it must be the case that a ⊥ = 0. This is required in order to satisfy Eq. 9 on either side of the inflection point, provided that a ⊥ is continuous along the stream. The result is that at an inflection point, a ⊥ = 0 and a planar points entirely along ±T for a planar ̸ = 0.
Alternatively, a linear stream segment in the plane of the sky could haven been observed at an epoch when the stream was undergoing a transition from having curvature in one direction, to zero curvature, and finally to curvature in the opposite direction (that is, from "ushaped", to linear, to "n-shaped"). We assume that such a scenario is unlikely to occur, and treat extended segments of a stream with zero curvature as having a ⊥ = 0.
In practice, due to statistical uncertainty in our estimation of the stream track, the derivative ∥dT /dγ∥ will not evaluate to exactly 0 even for streams with negligible curvature. Instead, the angular change in the tangent vectorT over an increment in angular arc-length, dℓ, is the dimensionless number: where ϕ is (e.g.) the angle between the tangent vector, T , and a fixed reference axis in the plane of the sky.
In Appendices A-B we discuss the stream track fitting routine and the threshold ϵ, which describes when we treat a stream segment as linear (|dϕ/dℓ| < ϵ).

Overview
In this section we develop a likelihood to apply the formalism developed in §2 to observational data. For our data inputs, we measure curvature unit vectors from a projected stream track that has been fit to the data (i.e., from imaging of the stream). We assume that N such vectors have been estimated from the data, and make up the set {κ i } N i=1 . The model component of the likelihood is the acceleration field evaluated along the stream, generated from a gravitational potential Φ m (x, y, z) with model parameters m (Eq. 5). The coordinates (x, y) are on-sky position coordinates, and z is the unknown line-of-sight coordinate. In our approach the acceleration field is allowed to vary along the line of sight, even in the absence of distance information along the stream. We handle the unknown line-of-sight coordinate by maximizing the likelihood over all distance tracks (discussed in §3.2).

Dealing with Unobserved Line-of-Sight Distances
In general, planar accelerations depend on both the on-sky position coordinates (x, y) and the line-of-sight coordinate, z. As a result, the angle between the stream curvature vectors and the planar component of the acceleration vector depends on the 3D location of the stream segment. We now discuss how our curvatureacceleration analysis treats the unobserved line-of-sight coordinate, by maximizing the likelihood over all possible continuous distance tracks in z. Then, the best possible potential model satisfies |θ i | < π/2 for all i while allowing for the possibility of a connected stream. This constraint is expressed mathematically in §3.3, though we provide a more qualitative description here.
In practice, each evaluation point along the stream track has some range of distances that yield |θ i | < π/2. For the i th evaluation point along the stream track, let D i represent the set of line-of-sight distances from the host galaxy in the interval [z min , z max ] for which |θ i | < π/2. Throughout this work we set z max = −z min ≈ 100 kpc. For some potential models we might find that the allowed distance ranges for subsequent stream segments (i.e., D i−1 and D i ) do not overlap. However, we expect distance to vary smoothly along the stream. This means that subsequent stream segments should have overlapping distance ranges within some reasonable threshold, ∆D. We impose this as a prior, by requiring that D i and D i−1 overlap within ∆D. We choose ∆D = 30 kpc, since this is a conservative upper-limit on the distance between subsequent evaluation points along a stream.
The distance sampling scheme has a few "edge cases," addressed in Appendix D. From the above discussion, it is evident that our analysis could, in principle, constrain the 3D shape of a stream track from only its projected morphology. Typically, we find that 3D tracks are poorly constrained, with the exception of a few special cases. We postpone exploring constraints on 3D stream tracks (i.e., measuring distance gradients) to a future work.
If we view the total mass distribution along a principal axis ( §4.2), then planar accelerations vary only in magnitude along the line of sight. In this case, all distance tracks are valid and requiring the stream to have a continuous distribution in z does not affect the likelihood.

Likelihood
We now develop a likelihood that connects the model introduced in §2 to the data. Our strategy is to evaluate the likelihood of the observed stream curvature, given the planar acceleration field of a potential model with parameters m. We construct the likelihood to prefer curvature-acceleration angles within 90 deg (i.e., planar acceleration vectors pointing inside the stream's curve), and a vanishing perpendicular acceleration at linear segments of a stream.
Because we constrain the acceleration field along the stream track, we must also consider where in 3D space the potential model is evaluated. For the i th evaluation point along the stream track, the on-sky coordinates are (x(γ i ), y(γ i )). Let z(γ) represent the unknown distance track of the stream. The full 3D position vector is X(γ). For a trial distance track, we consider the three cases: • Case 1 (C 1 ): |θ i | < π/2 at X(γ i ).
• Case 3 (C 3 ): The stream segment has negligible projected curvature and θ i is undefined.
Each evaluation point at X(γ i ) corresponds to one of the three cases. This allows us to write the probability for the i th curvature vector, conditioned on the line-ofsight coordinate z(γ i ), model parameters m, and the specific case C j . For the first two cases, we write this conditional probability as can be replaced with either C 1 or C 2 . Note that "otherwise" in Eq. 14 is applicable if C {1,2} is not satisfied, or ifκ i is undefined in the absence of curvature. The proportionality sign is sufficient, since the normalizing constant is an overall scaling factor in the total likelihood. We treat Eq. 14 as a binary outcome, since the curvature either is or is not compatible with the given condition (C 1 or C 2 ).
For condition C 3 , the stream appears locally linear in projection andκ is undefined. In this case, rather than enforcing that a planar points entirely along ±T ( §2.2.3), we allow for deviations from this behavior by placing a Gaussian prior on the angle between a planar and the stream track. We adopt this approach since any given stream will have a non-negligible thickness, making it difficult to determine precisely what straight line path should be traced to construct the stream track. Placing a Gaussian prior on the angle between the planar acceleration vector and stream track in this case is more conservative, since we will not rule out potential models that produce the approximately correct behavior (i.e., planar accelerations nearly aligned with ±T ). The angle between the planar acceleration vector and the stream track is given by whereN ≡ −T y , T x is orthogonal toT = T x , T y in the plane of the sky. This definition ensures θ T,i ∈ [−π/2, π/2]. Note thatâ planar,i depends on z, so θ T,i also depends on the line-of-sight distance along the stream track. In this case, we write the conditional probability for a stream element with negligible curvature as where N is the Gaussian distribution with mean µ and standard deviation σ. Eq. 16 up-weights acceleration vectors that point tangent to the stream track (small θ T ), and down-weights acceleration vectors with large θ T . In practice, because θ T is set to a fixed interval N should really be a truncated Gaussian. However, because we never sample θ T outside of this range, the only difference between N and the properly truncated Gaussian is a normalization, which is why we have written a proportionality symbol instead of equality in Eq. 16. In this work we choose σ θ T = 10 deg, representing our prior uncertainty on the direction of the stream track. We discuss this choice of the hyperparameter further in Appendix B. By marginalizing over each case (C 1 , C 2 , C 3 ), the likelihood for a curvature vector,κ i , conditioned on the gravitational potential model parameters, m, and lineof-sight track, z(γ), is where P C j |m, z(γ) is the probability (weight) of each respective case with j P C j |m, z(γ) = 1. We assume that the probability of each respective case is constant along the stream track. This allows us to treat P C j |m, z(γ) for j = 1, 2, 3 as three parameters, independent of position along the stream. For the gravitational potential parameters m and distance track z(γ), the total likelihood taken over the length of the stream (i.e., as a product of Eq. 17 over all i) will be maximized when P C j |m, z(γ) is the fraction of times that each case, C j , occurs along the stream track. For instance, if C k is satisfied over the stream's length then P (C k |m, z(γ)) = 1 and P (C j̸ =k |m, z(γ)) = 0 will produce the most positive total likelihood. If case C k is satisfied for only a small subset of the evaluation points, then P (C k |m, z(γ)) will be the proportion of evaluation points for which C k is satisfied. Otherwise, giving higher weight to this case would then suppress the other two cases (since the sum of the weights is 1), despite the fact that they occur more frequently. The weights that maximize the likelihood are where the summation is over the N κ curvature vectors with non-zero curvature, and a planar depends on parameters m and the distance track z(γ). From Eq. 18, f 1 is the fraction of evaluation points with compatible curvature vectors and planar accelerations, f 2 is the fraction of evaluation points with incompatible curvature vectors and planar accelerations, and f 3 is the fraction of evaluation points with undefined curvature vectors (this is fixed for each stream track). The term in the summation of Eq. 18 is 1 if the curvature vector is compatible with the proposed planar acceleration, and 0 otherwise. With each term in Eq. 17 defined, we can now write an objective function-maximized over all possible distance tracks satisfying the constraints outlined in §3.2-as Eq. 19 is not a standard likelihood function due to the maximization. In the frequentist literature, a likelihood of this form is referred to as a "profile likelihood" (Sprott 2008), since it is maximized over nuisance parameters (the distance track, in this case). We remind the reader that we do not explicitly sample over all possible paths in the maximization step of Eq. 19. As discussed in §3.2, we compute the likelihood recursively along the stream track, so the range of distances that the i th evaluation point could occupy is near the range of distances compatible with the adjacent evaluation point indexed i − 1. The distance sampling procedure is discussed further in Appendix D.
Eq. 19 is degenerate with the f parameters. To break the degeneracies, we require This ensures that maximizing the likelihood (Eq. 19) entails finding the acceleration field that is most consistent with the measured curvature vector: i.e., the acceleration field for which Eq. 12 is most frequently satisfied across the stream.

TESTS ON SIMULATED STELLAR STREAMS
In this section we demonstrate the method on N −body streams generated in a ground truth background potential. We first discuss the N −body simulations in §4.1, then test our method using streams generated in an axisymmetric halo-only potential ( §4.2.1). In §4.2.2 we apply the method to streams generated in a more complicated disk + halo potential, and consider measuring the mass ratio of the two components in §4.2.3. In §4.2 we assume that the galaxy and its halo are viewed along a symmetry axis, and relax this assumption in §4.3.

Simulations
In this section we discuss the details of our simulations, including initial conditions ( §4.1.1), the potential models considered ( §4.1.2), viewing angles ( §4.1.3), and fitting a stream track to simulation snapshots ( §4.1.4).

Initial Conditions
We simulate stellar streams by following the tidal disruption of satellite dwarf galaxies with self-gravity. The satellites are initialized with a King distribution function (King 1966) with W 0 = 3 and r s ≈ 1 kpc using the software package AGAMA (Vasiliev 2019). We sample the distribution function to generate 5 × 10 6 particles which are then evolved forward using the N −body code PyFalcon, a python interface for gyrfalcON implemented in the NEMO framework (Dehnen 2002;Teuben 1995). The total mass of each cluster is 10 6 M Sun . For each simulation, the progenitor's location is initialized by randomly sampling a point within a galactocentric spherical radius r ∈ [20, 100] kpc. The initial velocity is randomly sampled in 0.3-1V c , where V c is the local circular velocity. This random sampling procedure gives rise to a diverse range of stream morphologies, providing a useful test-set for analyzing distinct tidal features generated in a known ground truth potential.

Potentials
We generate streams using N −body simulations with a static external potential. For a stream in the extended stellar halo of its host galaxy, we expect our constraints to be dominated by the dark matter component of the potential. Therefore, our simulations always include a halo potential, though in §4.2.2-4.2.3 we also incorporate a Miyamoto-Nagai disk (Miyamoto & Nagai 1975). We view the disk edge-on, so that its midplane is in the x−z plane, and its perpendicular axis is along the y direction. The disk's potential in our reference frame is The disk has scale-mass M Disk = 1.2 × 10 10 M Sun , scalelength a = 3 kpc, and scale-height b = 0.5 kpc. These parameter choices are similar to constraints derived from the Milky Way, though the scale mass is somewhat lower than recent estimates, typically based on Jeans modeling or rotation curve estimates (e.g., Sofue et al. 2009;Kafle et al. 2014;Bajkova & Bobylev 2017;Chrobáková et al. 2020). We have explored changing the disk parameters and find that our conclusions are robust to reasonable variations.
For the halo we use the logarithmic potential, which has a flat rotation curve, similar to the observations. The functional form is (22) We have introduced the primed coordinates x ′ , which are aligned with the principal axes of the triaxial halo potential. We relate the prime and unprimed frames in §4.1.3. When expressed in a common coordinate frame, the total potential consisting of a disk and halo is simply Φ = Φ Disk + Φ Halo . Note that for an axisymmetric logarithmic potential, the relationship between the potential flattening (q) and density flattening (q density ) is 3(1 − q) ≈ 1 − q density for moderate flattening (e.g., Binney & Tremaine 2008 pg. 77).
The parameter v c in Eq. 22 is fixed such that the circular velocity at 8 kpc is roughly 250 km/s, and r c = 16 kpc. These parameters are similar to Milky Way potential constraints derived from the GD-1 stream in Malhan & Ibata (2019). For most simulations we set q 1 = 1, in which case we refer to q 2 as q. We carry out multiple simulations with the flattening parameter, q 2 , varying in the range [0.75, 1.5]. In §4.3.2 we attempt to constrain the triaxial potential in Eq. 22 with q 1 ̸ = 1.

Observing our Simulations
In our experiments, the primed frame is aligned with the principal axes of the assumed halo potential, whereas the unprimed frame describes our reference frame. Both frames have the same origin, at the center of the host galaxy. In our reference frame x − y is the plane of the sky, and z is the line-of-sight coordinate.
We relate the primed and unprimed frames through the linear transformation where R is the rotation matrix x z x 0 y 0 z 0 y ↵ ∢ Figure 4. Illustration of the coordinate frames. The x − y plane is the plane of the sky, and the z-axis is along the lineof-sight. The primed frame (red) is rotated with respect to the unprimed frame (black), and is aligned with the principal axes of the dark matter halo potential. The angle β describes a rotation about the y-axis, while α describes a rotation about the z ′ axis.
The matrix R characterizes a counterclockwise rotation about the y-axis by an angle β, followed by a clockwise rotation about the (new) z ′ axis by an angle α (similar to Euler angles). The relationship between the prime and unprimed frames is illustrated in Fig. 4. We test the impact of projection effects on our inference by observing the same halo model at different viewing angles, defined by α and β ( §4.3). In §4.2.2 we attempt to constrain a disk-halo misalignment angle, which is the angle between the flattening axis of an axisymmetric halo potential (the y ′ axis in Eq. 22) and an axis orthogonal to the disk's midplane. For the simulations in §4.2.2, the disk-halo angle is simply α in Fig. 4, with β = 0. We use the same definition in §5.2, where we apply our method to the stream surrounding the edge-on galaxy NGC 5907.
In Table 1 we provide a summary of the simulations utilized in each section. Multiple simulations with different flattening parameters are indicated by the curly brackets {·}. For simulations with a disk, the disk column is labeled Y. Otherwise, the column is blank.

Fitting a Stream Track to Simulated Data
To fit a stream track to our simulation snapshots, we follow a similar approach to Nibauer et al. (2022). We first bin particle positions in the x−y plane over a coarse grid to prescribe an initial ordering to stars populating the stream. Bins that are populated with particles are connected to the nearest neighboring bin with a nonzero number of particles. Each bin is connected only once. This defines an initial jagged stream track. This procedure can fail if the stream intersects itself. In this Table 1. Parameters of the simulations analyzed in each section. The parameters q1 and q2 characterize the flattening of the halo potential, and the disk column indicates whether or not a disk is included (Y is yes, N indicates no disk). When q1 = 1, we simply refer to q2 as q in the text (since in this case the halo potential is axisymmetric). The rotation angles α and β specify the orientation of the halo potential with respect to our reference frame. When a disk is included, it is viewed edge-on.
case, we manually select neighboring bins to fit the elongated axis of the stream. Each particle is then assigned an ordering parameter γ ∈ [−1, 1] based on its distance to the initial track. The γ parameter is a monotonic phase, encoding progress along the stream track. A differentiable smoothing spline is then fit to the ordered stream, which is used to estimate unit curvature vectors to the track. Additional details of the fitting procedure are discussed in Appendix A.

Single-Component Model
We now apply our stream-curvature analysis to Nbody streams generated in a ground truth logarithmic potential (Eq. 22) with a single free parameter: the y ′axis flattening (labeled q), and no disk component. The x ′ flattening (q 1 in Eq. 22) is set to 1. We work in a coordinate system that is aligned with the x ′ − y ′ plane (i.e., α = β = 0), so that the halo is viewed along a principal axis. For a fixed point x 0 and y 0 in this frame, varying z changes only the magnitude of planar accelerations and not their direction. This can be seen by taking equipotential slices along z: for the viewing angle considered here, each equipotential slice describes an ellipse with the same axis ratios. Because planar accelerations are perpendicular to each equipotential ellipse, the direction of the planar acceleration vector at a reference point (x 0 , y 0 ) is independent of z for the configuration considered.
In the top row of Fig. 5, we show three distinct N −body streams generated in gravitational potentials with different flattening parameters, q = [0.8, 1.0, 1.5] (left to right). Snapshots are taken at t ≈ 10 Gyr. Each column corresponds to a different stream shown in the x − y plane, labeled Stream A, B, and C. Equipotential contours are illustrated as black dashed curves. The spline based stream track is overplotted as the colorful curve. The color-coding of the track corresponds to the γ-parameter, which encodes a monotonic phase along the stream. The true planar acceleration unit vectors along the track are shown as red arrows, and the curvature vectors are shown as the small green arrows. The latter are estimated from the data using derivatives of the stream track. For stream A, the tail of the stream near γ = −1 has negligible curvature, so curvature vectors are not shown in this region.
The middle row of Fig. 5 illustrates the central concept of our analysis: the correct potential model must produce accelerations that fall within 90 deg of the measured curvature vector (otherwise, the tidal feature would not curve in the observed direction). For a given flattening parameter, q, we plot the angle, θ, between the curvature vector and the proposed planar acceleration vector evaluated along the stream track. This yields a family of curves, which are color-coded by the flattening parameter q. Flattening parameters that produce |θ| > π/2 anywhere along the stream track are "ruled out", while those that produce stream-acceleration angles within the white space across the length of the stream are compatible. The θ versus γ profile for the true acceleration field is shown as the dashed red curve.
Streams A and C in Fig. 5, middle row, both have a true-potential profile (dashed red) that approaches the gray ruled out regions and abruptly changes. For Stream C this is due to an inflection point, where the curvature vector switches directions by 180 deg to demarcate a change in the stream's concavity. At an inflection point the planar acceleration vector should then point along the stream track, since a ⊥ = 0 ( §2.2.3). For Stream A, there is a similar inflection behavior in the θ versus γ profile towards the segment of the stream around γ = −1. However, this segment of the stream is nearly linear, and satisfies the tangent condition |dϕ/dℓ| < ϵ. Because this portion of the stream has negligible curvature, we treat its curvature vectors as undefined and up-weight accelerations withâ planar ∝ ±T in the likelihood (discussed in §2.2.3 and §3). Indeed, this can be seen in the top row of x input flattening parameter is shown as the red line. For the cyan likelihood curves the tangent condition is not implemented and stream segments with undefined curvature vectors are removed. For the black curves the tangent condition is implemented if |dϕ/dℓ| < ϵ anywhere along the stream track. When there is no inflection point or linear stream-segment (e.g. stream B), the two likelihoods are exactly the same. We expect that whether the tangent condition is included or excluded that the resulting likelihoods are consistent, though with the tangent condition providing additional constraining power. This is most clearly seen with Stream C: when the tangent condition is turned on, the likelihood becomes slightly more narrow and provides a more informative constraint on the shape of the potential. For Stream A, only a narrow range of potential models are compatible with the nearly linear segment of the stream towards γ ≈ −1 and the looping segment with γ > −0.5. In this case including the tangent condition makes a neg-ligible difference to the likelihood, since the full length of the stream already places informative constraints on the shape of the potential. Across all three cases, our curvature-acceleration analysis produces accurate constraints on the flattening of the potential with the input parameters supported by the highest likelihood regions. For Stream A in Fig. 5 the true flattening parameter is at the edge of the highest likelihood region. This is due to the linear stream segment in the top left panel around (x, y) ≈ (15, −15) kpc; potential models that are too oblate or prolate tilt planar accelerations away from the stream track in this region. This leads to a sharp peak in the likelihood function. Slightly different stream tracks fit using bootstrap resampling of stream stars can shift the highest likelihood region towards the true q, motivating a fully probabilistic treatment of the stream track rather than the fixed-track analysis we present here. We discuss the limitations of considering only a single best fit stream track in §6.1. Overall, the strength of the constraints in Fig. 5 varies depending on the particular stream or potential model considered, and the tangent condition can further limit the range of possible models that are compatible with the stream. Streams with nearly linear segments or inflection points are most useful, since for these systems the direction of the curvature vector changes by 180 deg over a small area. This discussion is extended in §6.2, where we highlight which stream morphologies will lead to the most informative constraints on the potential of external galaxies.

Two-Component Model
So far, we have only considered a single component potential model; however, we expect galaxies to have a stellar and a dark matter component. In this section we explore the possibility of detecting both components, and their relative orientation.
Cosmological simulations predict that dark matter halos can be aspherical and tilted with respect to the disk (Schneider et al. 2012;Chua et al. 2019;Emami et al. 2021;Petit et al. 2022). Indeed, in the Milky Way there is observational evidence for a disk-halo misalignment with an ellipsoidal dark matter halo that is tilted with respect to the galactic plane (Law et al. 2009;Han et al. 2022;Panithanpaisal et al. 2022). Motivated by these findings, we now treat the disk-halo misalignment angle as a free parameter, and consider constraining its value from projected stream morphology.
The simulations utilized in this section are outlined in §4.1. Importantly, we rotate the flattening direction of the axisymmetric halo by an angle α = 30 deg relative to a perpendicular axis passing through the midplane of the disk. We refer to this angle as the disk-halo angle. A 30 deg misalignment is similar to recent estimates of disk-halo misalignment in the Milky Way (Han et al. 2022). The disk-halo angle, α, in our frame is illustrated in the bottom left panel of Fig. 6. In this frame the disk is viewed edge-on, with the y-axis perpendicular to the midplane.
In the combined potential of the disk and rotated halo, we simulate two streams in exactly the same external potential. The initial conditions for these simulations were generated using the random sampling procedure described in §4.1. Snapshots of both N −body realizations are shown in the leftmost-column of Fig. 6. We purposely apply our analysis to two streams at very different stages in their evolution, in order to demonstrate how our constraints in this section depend on when the stream is observed throughout its evolution. The stream in the top row of Fig. 6 corresponds to a t ≈ 10.5 Gyr snapshot, and t ≈ 1.5 Gyr for the bottom row.
For both N −body systems in Fig. 6 the stream track is illustrated with the dashed red curve. For the stream in the top row, the tangent condition is satisfied along the bottom right segment of the track around (x, y) ≈ (20, −30) kpc. The shorter stream in the bottom row has no change in concavity nor a locally linear segment.
In the middle column of Fig. 6 the curvature vectors (black vectors) for each tidal feature are shown, along with trial accelerations in the x − y plane (blue and red vectors). The bottom middle row is zoomed-in around the smaller stream for clarity. The blue vectors correspond to the correct potential model, which has the halo component misaligned with the disk. The red vectors are generated from an incorrect potential with no diskhalo misalignment. Indeed, for both tidal streams there are several instances where the red vectors (incorrect model) subtend an angle greater than 90 deg from the estimated curvature vector. However, under the correct potential model with the halo oriented at an angle relative to the disk (blue vectors), the planar accelerations and estimated curvature vectors are compatible.
When constraining the shape of the potential, we set the mass ratio between the disk and halo component to the true value (where the enclosed mass within 40 kpc is M Disk ≈ 1.2 × 10 10 M Sun and M Halo ≈ 3.8 × 10 12 M Sun ). In §4.2.3, we discuss constraining the mass ratio as a free parameter. For the flattening parameter we sample q in the range [0.5, 2]. For the disk-halo angle, α, note that unit vector accelerations in the x − y plane are equivalent when the flattening axis is rotated by 90 deg with q − → 1/q. To reduce degeneracies when evaluating the likelihood, rotations are sampled in the interval of α ∈ [−45, 45] degrees. This approach still results in some degeneracies with rotation, as an oblate halo with α = 40 deg degrees is not significantly different from a prolate halo with α = −40 deg, even though the two models are not exactly the same. Constraints for the flattening parameter, q, and disk-halo angle, α, are shown in the rightmost column of Fig. 6. In the top right panel, the maximum likelihood estimate (MLE) is shown as the red point while the true parameters are the cyan point in both panels (top and bottom right).
For the case of the extended stream in the top row of Fig. 6, a spherical halo is ruled out at the ∼ 2σ level using a likelihood ratio test between the best fit model and a spherical halo model with q = 1. A non-zero diskhalo angle is slightly preferred, though the significance of the misalignment is below the 1σ level. However, the constraints imply that equipotential contours are flattened when measured along an axis with α = 30 deg, or stretched for α = −30 deg. As discussed above, these two solutions are complementary, since an oblate loga- rithmic potential can be approximated as a rotated prolate potential in Eq. 22. From the MLE for this stream (red point; Fig. 6), we find that there is a slight preference for the q < 1 mode over the q > 1 mode.
For the less extended stream in the bottom left, the likelihood surface is split into two regions: one where planar accelerations are compatible with the estimated curvature vectors (yellow; high likelihood) and one where they are not (black; low likelihood). Because the stream has a much smaller length compared to the extended stream in the top left, an incorrect potential model can produce incompatible planar accelerations for a large number of evaluation points along the stream (i.e., nearly 1/3 to 1/2 of the evaluation points have incompatible planar accelerations when sampling in the black regions of the likelihood surface). This leads to the large contrast between high and low (near zero) probability regions in the likelihood panel. Because there is no inflection point nor linear segment along the stream, the tangent condition ( §2.2.3) does not provide any additional improvement to the constraints.
While the constraints in the bottom right panel of Fig. 6 are far less informative than in the top right panel, we can still make interesting statements about the shape of the halo potential as a function of its relative orienta-tion to the disk. For instance, if we assume that the disk and halo are aligned (i.e., α = 0) then we find q ≲ 1.1. For α = −40 deg all flattening parameters in the sampled interval are compatible. In this case, incorporating more streams in the constraints can break degeneracies, though it is still possible to consider a meaningful range of halo geometries (i.e. q ≲ 1.1) as a function of diskhalo misalignment angles.
The results presented in this section demonstrate that while we can place tight constraints on the flattening of the dark matter halo (Fig. 6, top row), the disk-halo misalignment angle is more difficult to constrain from projected stream morphology. More extended streams with multiple loops or an inflection point represent cases where both parameters can be determined at a higher confidence level.

Limits on Mass Ratios
In the absence of dynamical information (e.g., velocity of the stream), the mass enclosed by the tidal feature is degenerate with viewing angles, distance gradients, and the 3D velocity distribution of stars populating the tidal feature (see, e.g., Sanderson & Helmi 2013;Pearson et al. 2022b). This is clear from Eq. 11: adjusting the total mass of the system by scaling the amplitude of the overall potential does not alter the output of the sgn function. This makes intuitive sense: unit vector accelerations do not change if the total mass of the system increases or decreases uniformly (that is, for a potential function of the form Φ(x) ≡ GM f (x) − → cGM f (x) for some scalar c, unit vector accelerations are invariant to this scaling:â − →â). However, for two-component potentials (e.g., a disk and halo) of the form the ratio m 2 /m 1 can be constrained, sinceâ depends on the ratio m 2 /m 1 . The amplitudes m 1 and m 2 are related to the characteristic mass of each potential component, with the dimensions of (e.g.) GM (so that f 1 and f 2 have dimensions of inverse length). Varying the mass ratio between potential components will change the direction of unit vector accelerations, as one component "pulls" unit accelerations away from the other. This is especially relevant for combined potentials with an anisotropic term (e.g., a disk or flattened halo), since aspherical mass distributions can significantly alter the direction of unit vector accelerations in the plane. Regardless of the geometry of the mass distribution, the correct mass ratio yields an overall potential model with planar accelerations that are consistent with the projected curvature vectors along the stream. This provides a means to estimate a mass ratio between a halo component and a disk component, even if the mass of both components are individually unknown. If the baryonic mass of the disk is determined independently, e.g., using mass-to-light ratios or other scaling relations, then the estimated ratio m 2 /m 1 can be used to solve for the mass of the dark matter component interior to the stream. If sufficiently precise, this can provide a "mass-ladder" for a single galaxy, where a stellar mass estimate can be used to calibrate the halo mass estimate derived from stream curvature.
We demonstrate the ability to place limits on mass ratios using the N −body stream in Fig. 6, top row. We fix q and the disk-halo angle to the best fit values from Fig. 6 (shown as the red scatter point). We constrain the mass ratio log 10 M Disk /M Halo for the two-component potential model used in §4.2.2. In this section M Disk /M Halo refers to the mass enclosed ratio within a radius of 40 kpc, since this is roughly the extent of the stream in the x − y plane. The scale length and height of the disk component are set to the true values used in our simulations, under the assumption that these can be roughly determined from the photometric properties of the disk. The core radius of the logarithmic  Fig. 6, top row, generated in an oblate halo potential. The stream is also reproduced in the top right panel. At the best fit flattening (q) and disk-halo angle (α), we evaluate the likelihood as a function of disk-to-halo mass ratios (defined as the enclosed mass ratio within a radius of 40 kpc). The true input mass ratio is shown as the red line. Bottom Panel: Same as the top row, but for a stream generated in a prolate halo potential with q = 1.25 and a slightly more massive edge-on disk. In both cases we infer the presence of a dark matter halo more massive than the disk (black dashed line).
potential component is also set to the true value used in simulations, though reasonable variations on the core radius do not alter our conclusions. At the best fit q and disk-halo angle (α), Fig. 7 (top panel) illustrates the model likelihood for a range of mass ratios. Using the true values for these parameters does not alter our discussion, since the MLE for q and α is so close to the input parameters (both in likelihood and distance to the truth). The true mass ratio is illustrated with the red line. Mass ratios to the right of the dashed black line have M Disk > M Halo . The likelihood sampled over mass ratios plateaus to the left of the black line, implying that the data (slightly) prefer a two-component mass model with a dominant halo component over the disk component. A likelihood ratio test shows that a two component model is preferred at the 1.5σ level, compared to a model with no dark matter halo or a dominant disk component. Note that we may only place an upper limit on M Disk /M halo from Fig. 7.
A limiting factor in constraining the disk-to-halo mass ratio is the degree of similarity between the acceleration field of the two mass components. For instance, the topology of the unit vector acceleration field pro-duced by an oblate halo potential is similar to that of a disk. In contrast, the unit vector acceleration field produced by a prolate halo potential is distinct from that of a disk. In the latter case, the direction of unit vector accelerations is more sensitive to the relative mass amplitude between the two components, so we expect our curvature-based analysis to provide more stringent limits. We demonstrate this point in the bottom panel of Fig. 7, where we apply the same mass ratio analysis to a stream generated in a prolate logarithmic potential (q = 1.25), with a major axis perpendicular to the disk's midplane. The true flattening parameter is used in our inference, and only the mass ratio is allowed to vary. A slightly larger disk mass is also used compared to the top panel, though our conclusions do not change for small variations of the disk mass. For the prolate halo potential, the likelihood function implies a strict upper bound of log 10 M Disk /M Halo ≈ −1.95, corresponding to the mass enclosed ratio within 40 kpc. This is a much more informative limit in comparison to the top panel of Fig. 7, where the degree of similarity between the oblate halo potential and the disk limits what can be inferred about the relative mass of the two components using only stream curvature. Still, in both cases the input mass ratio (red line) is recovered within the highest likelihood region.
The results presented in this section suggest that that extragalactic streams in halo potentials with significantly different shape parameters from the disk (i.e., prolate, or misaligned rather than strongly oblate) provide strong prospects for placing limits on disk-to-halo mass ratios. In future work we intend to explore a wider range of potential models and configurations, to determine the viability of constraining mass ratios from projected stream morphology.

General Viewing Angle
In the preceding sections we assumed that the observed galaxy is viewed along one of its principal axes (i.e., perpendicular to the flattening axis of the halo potential), so that the inferred flattening, disk-halo angle, and disk-to-halo mass ratio can be interpreted as global constraints. If this assumption is broken, then the bestfit flattening parameter, q, characterizes the axis ratios of the best fit projected equipotential slice taken over the length of the stream. This does not need to coincide with the global shape of the potential, but characterizes a cross section of the full equipotential surface. For example, slices of an oblate (prolate) ellipsoid are oblate (prolate), or circular at worst. We have verified that the best fit flattening, q, for an inclined oblate or prolate po-tential continues to favor oblate and prolate flattening parameters, respectively. Importantly, the best fit equipotential slice is specialized to the distance of the stream. For a highly elongated or flattened potential that is inclined with respect to the line-of-sight, distinct portions of the stream can prefer different flattening parameters, depending on the magnitude of the distance gradient along the stream. The likelihood devised in §3 accounts for this possibility by sampling acceleration vectors over the line of sight, allowing us to-in principle-constrain the 3D shape of the potential by stitching together the range of elliptical slices that each segment of the stream is compatible with. We now explore to what extent the 3D shape of the potential can be recovered from the projected track of a tidal stream. We consider an axisymmetric halo potential in §4.3.1, and a triaxial halo potential in §4.3.2.

Axisymmetric Potential
We assume the same logarithmic potential (Eq. 22) as in the previous sections with q 1 = 1 and q 2 free (simply referred to as q). We also attempt to constrain the rotation parameters (α, β) introduced in §4.1.3. These parameters characterize the 3D orientation of the dark matter halo relative to our reference frame (Eq. 23-24).
As a test case we use a snapshot (t ≈ 10 Gyr) of a stream generated in the axisymmetric logarithmic halo potential with flattening parameter q = 1.3. We rotate the system by the angles (α, β) = (150 • , 85 • ), so that equipotential ellipsoids are "tilted" significantly towards the line-of-sight, and rotated slightly away. This orientation represents a scenario where our previous assumption of viewing the halo along a principal axis is broken. The rotated system is illustrated in the left panel of Fig. 8, where we plot the 3D stream in black and an equipotential ellipsoid in blue. The plane of the sky is depicted as the gray slate. The "observed" stream projection is shown as the gray points. Our fit to the projected stream track is illustrated by the red curve, and unit curvature vectors to this track are plotted in black for a few locations. As before, the red track and its derivatives in the (x, y) plane are the only inputs to the analysis. When sampling accelerations over the line of sight, we consider a prior distance interval within ±100 kpc (see §3.2 for a discussion of distance gradients). In this section posterior samples are obtained using the dynamic nested sampling package, dynesty (Speagle 2020). Constraints on the 3D orientation and flattening of the projected stream are shown in the middle set of panels in Fig. 8. The true parameters are the black lines. Each contour (68 and 95% regions) is shaded according to the density of samples it contains. The islands of high probability within each 68% region comes from the nearly linear segment of the projected stream, where the tangent condition ( §2.2.3) is satisfied. There are several modes and high probability islands due to degeneracies in reconstructing the 3D shape of the potential from only the 2D morphology of the stream. For instance, a prolate (q > 1) potential can also be approximated by a rotated oblate (q < 1) potential. This is the source of the two modes in flattening: q < 1 and q > 1. Placing tighter priors on α and β can significantly alleviate degeneracies (e.g., from measurements of the disk inclination), though in the present work we sample over all possible orientations. Note however that a spherical halo with q = 1 is disfavored in this example. While the posterior distribution over these parameters is degenerate and multimodal, we can rule out a halo component with a flattening axis perpendicular to the line-of-sight (α = β = 0) at the ∼ 2σ level. Furthermore, the highest posterior density regions are indeed close to the ground truth logarithmic potential, both in its global flattening and 3D orientation.
In the right panel of Fig 8 we plot equipotential ellipsoids corresponding to two samples from the posterior distribution in (q, α, β) shown as the scatter points in the middle panel (green and pink points correspond to the green and pink ellipsoids, respectively). Posterior draws with q < 1-while less likely than those with q > 1-tend to be rotated to have prolate elliptical level sets when sliced along the line of sight. This is the expected behavior, since planar accelerations are orthogonal to line-of-sight slices of equipotential ellipsoids, and our analysis constrains a potential model through its planar acceleration field.
Still, there are significant uncertainties and degeneracies in reconstructing the 3D shape of the potential from only the projected track of a single stream. As we will discuss in §6.2, certain stream morphologies will be most useful for constraining the potential. Incorporating multiple streams in the 3D analysis presented here is also computationally feasible, and can break degeneracies to provide more stringent constraints. Perhaps surprisingly, our results show that the continuous nature of streams along the line of sight enables us to utilize projected stream curvature to estimate 3D properties of the observed galaxy's halo, even in the complete absence of line-of-sight information.

Triaxial Potential
We now explore constraining the axis ratios of a triaxial potential from the projected shape of tidal streams. The analysis works the same as before: by measuring the curvature to a stream, we can identity a range of potential models that are compatible with the observed curvature. We generate tidal features in a triaxial version of the logarithmic potential defined in Eq. 22. Mockobservations are produced by rotating the stream and potential by the angles α = 0 deg and β = 40 deg. The simulation procedure is further discussed in §4.1. A snapshot of the rotated system in the plane of the sky (x, y) is shown in Fig. 9 (left panel), corresponding to t = 8.85 Gyr. While we evolve a single satellite, several streams and shells form, each delineating its own family of orbits (this is expected for a triaxial model; see Yavetz et al. 2022). We use splines to fit tracks to two visually obvious, elongated streams (labeled Stream 1 and 2) after manually selecting a few representative points to construct the track. Unit curvature vectors to the two tracks are shown. For both streams there is a change in concavity.
We treat each stream independently in the likelihood, and constrain the two flattening parameters (q 1 , q 2 ) and the rotation angle β. When sampling acceleration vectors from the triaxial potential, we fix r c to the true value (16 kpc), and the choice of v c is arbitrary since we only constrain the direction of acceleration vectors and not their magnitude. Acceleration vectors are sampled along the line of sight within ±70 kpc of the host galaxy, which is roughly the extent of the tidal debris in the x−y plane. The posterior distribution over the model parameters are again drawn using dynesty (Speagle 2020).
We present constraints in Fig. 9 derived from each stream separately (middle column), and jointly (right set of panels). The combined constraints are produced by sampling from a joint posterior, which is constructed by taking the sum of the total log-likelihoods for each independent stream. The constraints in the middle column are marginalized over the rotation parameter β. Contours correspond to regions of 68 and 95% confidence. The true parameters are shown as the solid lines, and dashed lines correspond to q 1 = q 2 (for an axisymmetric potential). A key advantage of our method is that we are immediately able to identify which aspects of a stream inform constraints on the potential (e.g., §4.2.1 and the θ(Φ) v.s. γ plots). For Stream 1, most information comes from the nearly linear segment of the track around x = 34 kpc in the left panel of Fig. 9. For streams with nearly linear segments or inflection points, we expect planar accelerations to align with their track ( §2.2.3). For stream 2, the inflection point around x = 5 kpc again dominates the constraints, though the extended looping segment helps place limits on the orientation of the flattening axes. Perhaps surprisingly, Stream 1 appears to provide more information about the potential compared to Stream 2. This is because Stream 1 has a flat segment around x = 34 kpc, activating the tangent condition for several evaluation points. Stream 2 has an inflection point, but only satisfies the tangent condition at a single evaluation point (the inflection point). Streams with extended flat regions are especially useful for constraining the potential, since flat segments have a vanishing acceleration component. Furthermore, the change in concavity around the flat segment of Stream 1 also limits the range of possible orientations for the potential.
From the joint constraints, a spherical halo is ruled out (> 2σ) due to the presence of inflection points, and a triaxial halo with q 2 > q 1 is preferred. The correct input parameters are recovered within the 68% region. When combining information from both streams, we find only a marginal improvement in constraining power compared to the constraints from Stream 1 alone. However, there is a region at low q 1 in the joint constraints that receives less support due to the inclusion of Stream 2.
Unlike generative simulation based methods, including additional streams in our analysis is not computationally expensive, since we simply compare acceleration vectors to curvature vectors along each stream track. In principle, constraining the potential from, e.g., tens of streams simultaneously would be practical, provided that a differentiable track is fit to each stream. Future observations are expected to reveal a substantial population of streams in M31 (i.e., Pearson et al. 2022a): this work provides an efficient method to constrain the potential from an ensemble of tidal features.
We use this section to demonstrate that for galaxies with multiple streams or streams with constraining features (i.e., inflection points, flat segments, saddle points), it can be feasible to constrain shape parameters of the halo beyond axisymmetry. A more realistic potential model may also include a disk component, and parameters characterizing the relative orientation of the triaxial halo to the disk. In the absence of velocity information, constraints will typically be degenerate, though can still provide an indication of the range of halo geometries that are compatible with stream curvature.

MEASURING THE GRAVITATIONAL POTENTIAL OF NGC 5907
We now apply the method introduced in this work to low surface brightness imaging of the tidal feature surrounding the galaxy NGC 5907. NGC 5907 is an edge-on spiral galaxy at a distance of ∼17 Mpc (Tully et al. 2016), with stellar mass ∼ 8 × 10 10 M Sun (Laine et al. 2016). The detection of a stream associated with the galaxy is first discussed in Shang et al. (1998) and Zheng et al. (1999). Subsequent imaging of the galaxy and its tidal stream revealed a tentative second loop (Martínez-Delgado et al. 2008), which was reproduced in the context of N −body simulations. Recent imaging°1   Figure 11. Constraints on the gravitational potential of NGC 5907 for a single component model. Left: Using the curvature vectors to the stream track surrounding NGC 5907 (Fig. 10), we fit a single-component logarithmic halo potential with flattening parameter q. We allow both the flattening axis and parameter to vary. The disk-halo angle (α) is measured relative to an axis perpendicular to the disk. The likelihood surface for this two-parameter model is illustrated. The MLE is the red point. Middle: At the best fit disk-halo angle (α ≈ −7 • ), we plot the stream-acceleration angle θ versus the monotonic phase-parameter γ.
Each θ versus γ profile is color-coded by a trial flattening parameter, q. The θ profile for the best fit model is shown as the dashed red curve. Right: At the best fit parameters (with q ≈ 0.66), we plot equipotential contours. The stream track is also overplotted, color-coded by γ. The best fit equipotential contours are oblate, with a major axis that aligns closely with the disk. Fig. 10. The green segment represents a region where the estimated track satisfies the zero curvature condition, namely, |dϕ/dℓ| < ϵ ( §2.2.3 and §3). We treat the disk-halo angle, α, as a free parameter in §5.2. We present constraints on the full track (both the detection and tentative segment), and discuss where along the stream our constraints originate from.

Single Component Potential
We adopt the logarithmic halo potential in Eq. 22 with q 1 = 1 and q 2 as a free parameter (simply referred to as q). The core radius is fixed to r c = 26 kpc, the same as in van Dokkum et al. (2019). We also fit for the flattening direction, α (illustrated in Fig. 10, bottom panel). Because NGC 5907 is viewed nearly edge-on (i ≈ 87 • ; Sackett et al. 1994), we assume that the halo is viewed along a principal axis, with the flattening direction perpendicular to the line of sight, though possibly rotated with respect to the edge-on disk. We sample q ∈ [0.4, 2.5], and allow the flattening direction, α, to vary between −45 • and 45 • . While q < 1/ √ 2 ≈ 0.7 can produce negative mass densities near the origin for the logarithmic potential, the 3D geometry of the halo can bias the best fit q towards lower values if the halo is oblate and not viewed along a principal axis. Furthermore, q < 1/ √ 2 can still provide a useful description of dark matter halos far from the origin (∼ 10s of kpc), where a large portion of the stream resides.
The likelihood surface in the space of rotation angles, α, versus flattening, q, is shown in the left panel of Fig. 11. The best fit (maximum likelihood estimate; MLE) parameters occur at the red point. A rotation angle of 0 occurs when the halo flattening axis is per-pendicular to the disk. The best-fit rotation angle is α ∼ −7 deg, though there is a large range in α that may also produce the stream's curvature. If the curvature vectors along the stream are mostly due to the disk's potential, then this is the expected behavior for a single component potential model (i.e., zero misalignment).
At the best-fit rotation angle, we plot the streamacceleration angle, θ, in the middle panel of Fig. 11 as a function of position along the stream track, γ. The θ versus γ curves are color-coded by trial flattening parameters, q. For reference, the stream track is color-coded by γ in the right panel. From the middle panel, we can determine where along the stream the constraints originate from. Namely, the regions of the stream closest to the disk (in projection) and the intermediate regions with γ ≈ 0 rule out q ≲ 0.5. The extended regions of the stream towards γ ≈ 0.25 rule out q ≳ 1.7. Additional constraints come from the segment of the stream close to the disk (in projection), where we find that the tangent condition |dϕ/dℓ| < ϵ is satisfied. This indicates that planar accelerations pointing along the stream track should be preferred in this region ( §2.2.3). The best fit potential (dashed red curve) is shown to simultaneously satisfy these constraints across the stream's length.
We overplot equipotential contours of the best-fit model on the image of NGC 5907 in Fig. 11, right panel. We emphasize that our analysis knows nothing about the existence of the disk in photometry: the only input is the track of the stream, and the orientation of the potential model is allowed to vary along with its flattening. At the best fit parameters, equipotential contours align closely with the projected orientation of the disk. We find that our fits to a single-component potential model prefer a strongly oblate potential (q ≈ 0.6 − 0.8), indicating a possibly substantial influence of the disk's potential on the orbit of the tidal stream. This is further supported by the small relative angle between the equipotential contours and the disk.
There is an additional class of high-likelihood models in the left panel of Fig. 11 at more prolate q ≳ 1.2 with a rotation angle α ≈ 45. We have verified that when sampling α > 45 deg that this mode continues until equipotential contours align with the disk. We take this to indicate a degeneracy in the likelihood surface, where an oblate potential can be expressed as a rotated prolate potential. We find similar behavior in §4.2.2 on simulated data. Either mode suggests that equipotential contours are quite flattened when measured along an axis perpendicular to the disk. In the following section we explore a two-component potential model, with separate disk and halo components.

Two Component Potential
In this section we follow the same approach presented in §4.2.2 on simulated data, and attempt to constrain the parameters of a two-component potential model using the stream track of NGC 5907. We work in the same parameter space as in §4.2.2. Following Martínez-Delgado et al. (2008) and van Dokkum et al. (2019), we represent the disk with a Miyamoto-Nagai potential (Miyamoto & Nagai 1975) with scale-length 6.24 kpc and scale-height 0.26 kpc. The disk potential is rotated so that it is viewed edge-on and closely aligned with the major axis of the disk in observations, traced by its stellar mass (Fig. 10). For the halo, we use the logarithmic potential (Eq. 22) with scale-radius r c = 26 kpc (following van Dokkum et al. 2019), q 1 = 1, and q 2 = q (a free parameter). The halo potential is therefore axisymmetric. In what follows, we constrain the 2D orientation of the flattening axis (the disk-halo angle; α), the flattening (q), and the mass ratio log 10 (M Disk /M Halo ), where M Disk /M Halo is the mass enclosed ratio within a 65 kpc radius, since this is roughly the extent of the stream in the x−y plane. We again assume that the halo is viewed along a principal axis, so that the flattening axis is perpendicular to the line-of-sight. It is possible to relax these assumptions, but we leave this to future work.
Posterior samples are obtained using the likelihood discussed in §3 and the dynamic nested sampling package, dynesty (Speagle 2020). Constraints for the twocomponent potential are illustrated in Fig. 12, where contours correspond to regions of 68 and 95% confidence. The mock stream in the top right panel is discussed in , and mass ratio between the two potential components are all free parameters. Dark and light blue contours correspond to regions of 68 and 95% confidence. We again find a preference for an oblate halo component, with q ≈ 0.8. At lower disk-to-halo masses, the flattening parameter prefers slightly higher values. In the top right panel, we generate a mock stream in the two-component potential with parameters set to the trial point (red). The morphology of the mock stream generated at the red point is well-matched to the stream surrounding NGC 5907. Both tentative and detection segments are captured by the mock stream with an oblate halo potential. §5.3. When incorporating the disk we find a preference for a slightly higher q = 0.76 +0.10 −0.13 at 68% confidence, when conditioning on log 10 M Disk /M Halo < −1 (reasonably, we expect the halo to have a dominant mass over the disk by at least an order of magnitude). Our fits are consistent with q ≈ 0.85 at the 1σ level, corresponding to q density ≈ 0.55.
The mass ratio constraints neither confirm nor ruleout a substantial contribution from the dark matter halo surrounding NGC 5907. This is, in part, due to degeneracies induced by sampling over a wide range in stream line-of-sight distances (within ±150 kpc in this case). At large distances from the disk, the scale mass of the halo does not need to be dominant to reproduce the curvature of the stream, since the contribution of the disk to the acceleration field diminishes with distance.
Still, it is interesting to consider the shape of the dark matter halo as a function of realistic disk-to-halo mass ratios. For instance, the top-left panel of Fig. 12 implies an oblate halo component with q ≈ 0.85 for log 10 M disk /M Halo ≈ −1. The bottom right panel implies that there is not a significant misalignment between the halo and disk component. Both constraints indicate a preference for an oblate halo component, with a flattening direction within ± ∼ 20 • of an axis orthogonal to the disk's midplane. These results do not necessarily imply that we are detecting the dark matter halo of NGC 5907. Instead, it could also be the case that the stream's curvature cannot rule out an oblate halo potential. Our analysis is consistent with both explanations.
On the other hand, our constraints disfavor q > 1 at the 2σ level within the disk-halo angle interval considered. This indicates that we do not expect the halo component of the potential to be significantly different from the disk, at least interior to the stream. This is similar behavior to Milky Way potential constraints presented in Koposov et al. (2010), derived from the 6D phase-space distribution of the GD-1 stellar stream. In their work, the flattening of the halo model is attributed mostly to the influence of the disk. As we will discuss in §6.4, we anticipate that kinematic data (e.g., radial velocity measurements) from the steam surrounding NGC 5907 will help disentangle the influence of the disk from that of the halo, and break degeneracies in our interpretation of Fig. 12.

Validation with a Forward Model
We now use a forward model to generate a simulated stream in a potential consistent with our constraints in Fig. 12. We compare our results to van Dokkum et al. (2019), who found that a simulated stream generated in a mildly prolate halo (q = 1.1) reproduces the most prominent parts of the stream in observations. However, their model fails to reproduce the full extent of the stream (γ ≳ 0.25). In contrast, our analysis, which fits the entire observable length of the stream, prefers an oblate halo.
To validate that an oblate halo can produce the tidal feature surrounding NGC 5907, we select a trial point contained within the 68% credible interval (red point in Fig. 12) and generate a mock stream at the trial point. The trial point is q, log 10 (M Disk /M Halo ), α ≈ (0.77, −0.6, 8.5 • ). In order to generate mock streams, we calibrate the mass of the disk to yield a combined (disk plus halo) maximum circular velocity consistent with that used in van Dokkum et al. (2019), derived from Posti et al. (2019). For the disk-halo mass ratio considered here, this corresponds to M Disk ≈ 1.5 × 10 11 M Sun . The resulting halo mass at the trial point is then M Halo ≈ 10 0.6 M Disk = 6 × 10 11 M Sun . For comparison, this corresponds to M 200 ≈ 3 × 10 12 M Sun for the total potential, where M 200 is defined as the mass within a radius enclosing 200 times the critical density of the universe. This is similar, by construction, to rotation curve measurements from Posti et al. (2019), who found M 200 ≈ 10 12 M Sun .
We do not claim that these are in fact the absolute masses of the disk and halo components. Instead, we aim to show that the stream curvature is consistent with such absolute masses, based on our constraint of the mass ratio M Disk /M Halo . A larger mass difference is also supported by our constraints, e.g., with log 10 M Disk /M Halo = −1. The same M 200 can be achieved with this mass ratio, by calibrating the absolute masses of each component (e.g., with the Posti et al. 2019 measurements).
Using the trial point and the absolute masses for the two-component potential previously discussed, we generate a mock stream using the particle-spray method implemented in Gala (Price-Whelan 2017). We sample over a coarse grid in initial phase-space conditions, until a qualitative fit between the observed stream and the generated stream is achieved. The integration time is ∼ 6 Gyr and the progenitor mass is taken to be ∼ 5 × 10 8 M Sun .
A mock stream generated using this approach is shown in the top right panel in Fig. 12. We find a strong qualitative agreement between the morphology of the stream imaged with the Dragonfly telephoto array and the mock stream generated in this section. The mock stream appears to have a similar tentative loop to the one shown in Fig. 10. This tentative loop has not been previously reproduced in the literature, though an oblate halo gives rise to this feature.

DISCUSSION
We have shown that the curvature of extragalactic tidal streams can be used to constrain properties of the host potential. In this section we discuss our results, starting with the limitations of our analysis in §6.1. We then consider which geometrical aspects of extragalactic tidal features will be most useful for constraining the potential of external galaxies in §6.2. We draw comparisons to forward modeling techniques in §6.3, and discuss future directions in §6.4.

Caveats
We first consider the stream track fitting routine discussed in §4.1 and Appendix A. The fitting routine currently implemented is not unique, nor optimal. For a given stream, there could be several tracks with slightly different curvature vectors that are compatible with photometry. One solution is to sample a posterior distribution over tracks, and propagate errors through to the model parameters of the potential. In Nibauer et al. (2022) an unsupervised stream track fitting algorithm is developed, though the algorithm utilizes velocity information that we are unlikely to measure in external galaxies. Starkman et al. (2022) developed another approach for fitting stream tracks, utilizing self-organizing maps and Kalman Filters. Both methods operate in the context of the Milky Way, and we leave the task of extending them to photometry of external galaxies to future work.
We now discuss whether a reliable distance to the host galaxy is required to carry out our analysis. The curvature-acceleration connection discussed in §2 is independent of distance to the host galaxy. However, a potential model is adopted to interpret the curvature of a stream. As long as the total potential model depends on distance to the host galaxy only through a scalar multiple (e.g., Φ − → aΦ when scaling distance), then unit vector accelerations in the plane of the sky are scalefree. This is the case for the Miyamoto-Nagai potential and logarithmic halo potential considered in this work. For a composite potential consisting of both a disk and halo component, the relative scale between characteristic lengths in each potential model can still influence the direction of unit vector accelerations in the plane of the sky (e.g., the scale length of the disk versus the core radius of the halo). The possibility of constraining the ratios of these parameters is left to future work.
Throughout this work we make the assumption that a coherence condition is satisfied along the stream track ( §2.2.2), allowing us to connect stream curvature to planar accelerations. In the Milky Way the gravitational potential is evolving, and this can influence the kinematics and morphology of tidal debris (e.g., Erkal et al. 2019;Shipp et al. 2021;Vasiliev et al. 2021;Lilleengen et al. 2023). It is possible that the coherence condition can be satisfied even in time-dependent potentials, since we do not assume that proper motions along the stream are aligned with its track. Future work will be needed to determine the effect of a time-dependent potential on projected stream morphology.
Throughout this work we have focused on stellar streams, though stellar shells have also been observed around external galaxies (e.g., Schweizer & Seitzer 1988;Martínez-Delgado et al. 2008;Duc et al. 2015). Our method for constraining the potential of external galaxies can be extended to the case of shells, where the curvature of the projected shell boundary places limits on the direction of planar accelerations using similar arguments to those developed in §2. We defer an exploration of this extension to future work.

Sensitivity of Analysis to Stream Morphology:
What features are most informative?
In this section we discuss what morphological properties of tidal streams will be most useful for constraining the potential of external galaxies, based on our analytic model in §2 and the tests presented throughout the paper. For instance, in Fig. 5, stream B seems to suggest that our method has poor constraining power for spherical potentials since the likelihood function supports a wide range in flattening parameters. Alternatively, it could also be the case that the constraining power of our analysis is tied to the specific geometrical properties of the tidal feature considered, and therefore when it is observed throughout its evolution. We now carry out an exercise to test how constraints on the potential vary as a function of stream morphology, by analyzing the same stream at two different stages in its evolution.
In Fig. 13 we test how the observed morphology of a stream affects our ability to place informative constraints on the underlying potential. In each column we plot the same stream, though "observed" at two different epochs indicated by times t A and t B . The particular snapshots are chosen to demonstrate how the likelihood function can change significantly as a function of different stream lengths and geometries. The background potential varies across streams D-F, though for a given column the potential model is exactly the same.
For each snapshot, the spline-based stream-track is shown by the red dashed-curve. We plot the corresponding likelihood functions for each system in the bottom row. We utilize the full likelihood introduced in §3, with the tangent condition implemented when the stream passes through a region without curvature. The blue likelihood curve corresponds to snapshot t A , while the dashed red curve corresponds to snapshot t B . The true flattening parameter for each case is shown as the black vertical line.
For the same stream observed at two different epochs, the width of the likelihood function can vary significantly. This is most evident for Streams D and E in Fig. 13, which display a significant change in morphology between snapshots t A and t B . While all snapshots yield constraints that are consistent with the true halo shape, streams with segments of negligible curvature, inflection points, or abrupt apocentric passages tend to provide much stronger constraints on halo shape parameters. This aligns with the analytic treatment introduced in §2: streams with tangent vectors that change direction quickly along the track sample the largest range of curvature vectors in the plane. The range of compatible accelerations for a given potential will be severely limited, since the constraint |θ| < π/2 becomes increasingly restrictive if the stream samples several directions over a small area. This is typically the case for streams with sharp apocenters, as for streams D and F at their respective t B . Our analysis is also powerful when a segment of the stream has negligible curvature (e.g., the tidal feature surrounding Centaurus A forward modeled in Pearson et al. 2022b). In this case, we employ the tangent condition discussed in §2.2.3, which constrains the planar acceleration vector to point along the local, curvaturefree segment of the tidal feature. Our analysis is most informative when a stream has segments with negligible curvature, and with significant curvature. This occurs for stream D at time t B in Fig. 13. Towards the left end of the track, the curvature becomes negligible (i.e., |dϕ/dℓ| < ϵ) and the tangent condition is applied. The intermediate region of the track has a sharp turn, so the tangent vector changes rapidly over a small spatial scale. The result is a tight constraint on the shape of the gravitational potential, since collectively, there are only a small range of continuous acceleration fields that can satisfy these constraints simultaneously across a single stream. Alternatively, multiple streams surrounding one galaxy with different constraining features can be used to fit the potential through a joint likelihood ( §4.3.2).
Finally, we highlight a third regime that is illustrated in Fig. 13, stream E at t A . The snapshot is taken at a moment when the stream passes near the origin at an angle inclined from the x-axis. In this snapshot, curvature vectors are pointing up and to the right. For a potential that is strongly oblate, the bottom right segment of the stream curves in the correct direction (towards the central mass density), while the top left segment curves away from the central mass density. In this case, the top left segment of the stream is incompatible with an oblate potential. For a strongly prolate potential, the converse is true: the top left segment of the stream curves towards the central mass density, while the bottom right segment curves away. Therefore, the bottom right segment of the stream is incompatible with a strongly prolate potential. From the likelihood plot in the bottom row, this configuration is enough to rule out either extreme (i.e., an oblate or prolate potential), and is best satisfied by q ≈ 1; in between prolate and oblate. Thus, extended streams passing near the central mass density of the host (in projection) provide another case where the potential is well constrained.
Collectively, the above scenarios can occur over a wide range of potentials. The constraining power of our analysis is tied to the presence of morphological features along a stream, which is predominantly determined by the stream's length and our viewing angles, rather than the halo shape. As long as the constraining morphologies discussed above are present, we can, in principle, constrain even complicated potential models with many shape parameters.

Comparison to Forward Modeling
In this section we discuss how our method compares to forward modeling techniques for stream formation in a specified potential.
A common method for forward modeling streams is the "particle spray" technique (Fardal et al. 2015), where test particles are released from a parent orbit near the Lagrange points of a progenitor with a specified mass. Other forward models adopt similar techniques, though with different stream formation prescriptions (e.g., Johnston et al. 1999;Fardal et al. 2006;Koposov et al. 2010;Price-Whelan et al. 2014;Küpper et al. 2012;Bonaca et al. 2014;Gibbons et al. 2014).
Forward models simulate the full time-evolution of the released particles, and can reproduce a wide range of stream morphologies. Additionally, these techniques can utilize the observed width of the stream to constrain the potential, along with its surface brightness profile (and therefore number density). The effects of a time-dependent potential can also be incorporated, along with perturbations from substructure (e.g., dark matter subhalos; Bonaca et al. 2019).
On the other hand, forward models have a large number of nuisance parameters that specify the full 6D phase-space initial conditions of the progenitor, progenitor mass, mass loss rates, integration times, the orientation of the system with respect to the line of sight, and the parameters of the adopted potential model. If the object of interest is the gravitational potential of the host galaxy, attempting to infer all of these quantities simultaneously from only projected stream morphology introduces strong degeneracies with the parameters that we wish to constrain. These degeneracies are only lifted with the inclusion of kinematic data. In the absence of this information, it is likely that the number of modeling degrees of freedom are higher than what the data can realistically constrain (e.g., the same stream morphol-ogy can be achieved for different halo masses; Pearson et al. 2022b). The method also scales poorly with the number of streams, and does not indicate which parts of the stream(s) considered are driving constraints on the potential.
In contrast, the method presented in this work does not model initial conditions, nor the unobserved velocity components along the stream. The only required input is the shape of the current stream track, and an estimate of the track's curvature along the stream. We do not rely on orbit integrations, stream formation prescriptions, nor prior knowledge of the progenitor's mass or phase-space position. Instead, the computational cost is front loaded to the track-fitting routine. Once a track is fit to the stream(s) of interest, varying the parameters of a potential model and the orientation of the potential is inexpensive, since we only compare the direction of acceleration vectors to the observed stream curvature. Furthermore, our approach also indicates where along the stream constraints on the potential originate from. Our inference can be performed simultaneously for multiple streams around a single galaxy (e.g., Crnojević et al. 2016, and some of the candidate mergers in Giri et al. 2023), or for many galaxies with streams (e.g., Shang et al. 1998;Ibata et al. 2001;Martínez-Delgado et al. 2009;Crnojević et al. 2016;Morales et al. 2018;Kado-Fong et al. 2018;Sola et al. 2022;Giri et al. 2023). This gives rise to the possibility of constraining dark matter halos at the population level, as we will discuss in §6.4.
A distinct advantage of our approach over forward modeling is that our model in §2 is very general, though still provides meaningful insights on what can actually be inferred about the potential from extragalactic stream morphology. For instance, the curvatureacceleration connection derived in §2.2.2 tells us that from stream curvature, we can only constrain properties of the potential that influence the direction of the acceleration field when projected onto the plane of the sky. In this case, it is clear that the absolute mass of a galaxy and its halo cannot be determined (without imposing priors on the progenitor's initial conditions), but a mass ratio between the two components can be constrained. Shape properties such as the flattening direction, and even 3D orientation of the halo or disk can also alter the direction of planar acceleration vectors making these good parameters to target. We can also predict which geometrical aspect of tidal features will be most useful for constraining halo shape parameters, such as inflection points or segments of a stream that appear linear in projection.
There are limitations to our method as well. For example, we do not constrain the initial conditions of the progenitor, nor its present-day location. However, if the only available measurements are of the stream's projected morphology, we do not expect that initial conditions can be determined with high confidence. Our approach does not impose priors on the velocity of the progenitor, or stars released from the progenitor. This means we are more agnostic to the progenitor mass, though stream-width could encode information about this quantity (e.g., Erkal et al. 2016). We have also neglected modeling the stream's surface brightness profile, since we have only considered stream curvature. Incorporating surface brightness information is possible in our analysis, though we defer this to future work. Finally, we have considered the realistic scenario of having no kinematic information of the stream but only its projected morphology. For a limited number of extragalactic streams, radial velocities or distance information has been measured (e.g., from the progenitor globular cluster; Veljanoski et al. 2014; surface brightness fluctuations; Toloba et al. 2016; and tip of the red-giant branch measurements; Crnojević et al. 2019). When available, these measurements can be incorporated in our constraints as discussed in §6.4. However, the present work considers the more ubiquitous scenario of having only measured the stream morphology from imaging.

Future Directions
The fundamental properties of dark matter halos are their mass and shape; while there are many techniques to measure halo masses (e.g., Rix et al. 1997;Napolitano et al. 2009;Watkins et al. 2010;Posti et al. 2019), measurements of halo shapes beyond spherical models have been more limited. In the extragalactic context, weak lensing has been the primary method for constraining halo shapes by stacking galaxy images and measuring the resulting shear field (e.g., Natarajan & Refregier 2000;Mandelbaum et al. 2006). Extragalactic tidal features provide a new method to measure halo shape parameters at the level of individual galaxies, and across an ensemble. Pearson et al. (2022b) showed that with only a single radial velocity measured from an extragalactic stream (in addition to stream morphology), the halo mass is much better constrained. While the analysis presented in this work is entirely concerned with stream curvature and its connection to halo shape parameters, radial velocity information can also be included in our approach. By calibrating the rotation curve of our morphology-based model against a radial velocity measurement from the stream, there is a clear path for-ward to self-consistently incorporate kinematic information in our constraints. Indeed, we have already carried out a similar analysis in §5, where we use previous measurements of the rotation curve from NGC 5907 (Posti et al. 2019) to calibrate the inferred mass ratio M Disk /M Halo . If a radial velocity is supplied to our model at the stream's radius, we anticipate that this information could also be utilized to better inform halo shape parameters, and help disentangle the effect of the halo from the disk. We leave this to future work.
Perhaps the most exciting future direction to consider is a population level description of dark matter halo shapes (e.g., axis ratios and orientations), probed by the morphology of extragalactic streams. From the 1000s of extragalactic tidal features expected to be observed over the next decade with The Rubin Observatory (Ivezić et al. 2019), Euclid (Racca et al. 2016), and Roman (Spergel et al. 2013), constraints on dark matter halo shape parameters can be combined at the population level to obtain a global view of dark matter halo properties. While stream morphology based constraints on the dark matter halos of individual galaxies can be weak or degenerate with viewing angle, a population level (e.g., hierarchical) model that combines information across many galaxies and streams will be better constrained. Our analysis is suited for this purpose, since all that is required is the track of a stream and its curvature vectors. From a population level description of dark matter halo morphologies, comparisons can be drawn to cosmological simulations which make important predictions for dark matter halos shapes in ΛCDM (e.g., Allgood et al. 2006;Schneider et al. 2012;Chua et al. 2019;Petit et al. 2022).
We have worked under the assumption of cold dark matter, though our analysis could be sensitive to other descriptions that alter the acceleration profile of a galaxy (e.g., Modified Newtonian Dynamics; Milgrom 1983). We leave an exploration of extragalactic streams and alternative descriptions of the gravitational acceleration field to future work.

SUMMARY AND CONCLUSION
We have drawn an analytic connection between the curvature of tidal streams and the underlying acceleration field of the host galaxy. The curvature-acceleration connection can be used to identify gravitational potential models that are consistent with the projected 2D morphology of the tidal feature considered. We have shown that in order to connect stream curvature to accelerations-and therefore the potential-the stream does not need to trace a single stellar orbit. In §3 we devised a likelihood to identify the range of potential models that are compatible with the local curvature direction(s) of a stream. The only required input to our analysis is a differentiable stream track, from which curvature vectors can be estimated. The method is computationally inexpensive, and does not rely on forward modeling techniques. Because of this, our analysis scales well with the number of streams considered.
We have demonstrated the method on stellar streams generated in N −body simulations with a static background potential. From 2D spatial snapshots of these simulations, we recover flattening parameters of the host halo, disk-halo misalignment angles, and in some cases mass ratios between different potential components. The best constrained parameter is the halo flattening, while constraints on misalignment angles and mass ratios are weaker. We have also shown that in some cases the 3D orientation of the halo can be constrained, and that limits can also be placed on the axis ratios of a triaxial halo model from only projected stream morphology. These constraints have significant degeneracies, though still allow us to rule out a range of possible geometries. Modeling multiple streams with constraining features can help place more stringent constraints on the 3D potential.
Our analytic treatment implies that streams with a change in concavity or a linear segment are most informative, since in these regions planar accelerations align with the stream track. Streams with tightly wound loops are also informative, since their curvature vectors sample many directions over a small area in the sky.
We apply our analysis to the stellar stream surrounding the edge-on galaxy NGC 5907, and find a preference for an oblate halo potential (q = 0.76 +0.10 −0.13 when conditioning on M Disk /M Halo < 0.1). Mock streams generated in potentials consistent with our constraints successfully recover the observed stream morphology in full.
There are already several observations of extragalactic streams (e.g., Shang et al. 1998;Ibata et al. 2001;Martínez-Delgado et al. 2009;Crnojević et al. 2016;Morales et al. 2018;Kado-Fong et al. 2018;Giri et al. 2023), with 1000s more expected from The Rubin Observatory (Ivezić et al. 2019), Roman (Spergel et al. 2013), and Euclid (Laureijs et al. 2011). Traditional methods for constraining the potential using extragalactic streams will be hampered by the lack of kinematic information and distance tracks, since we will often only have measured the projected stream morphology from imaging. The approach developed in this work is well-suited to this limited data scenario, and is computationally inexpensive. The method can be applied to a large number of galaxies with tidal streams, providing the means to measure the 3D shapes of dark matter halos at the population level. This enables future tests of ΛCDM using the morphology of extragalactic streams, by comparing the inferred distribution of dark matter halo geometries to theoretical predictions from cosmological simulations. The modeling approach introduced in §2 relies on estimates of unit curvature vectors along the projected stream track. The unit curvature vector is related to the second derivative of position along the projected track, as expressed in Eq. 7. We estimate these derivatives by fitting a smooth, differentiable spline to a series of representative points along the N −body stream in the plane of the sky. Our fitting routine works in two steps. For a given snapshot, we first bin particle positions in the x − y plane and take the average (x, y) coordinate in each bin. For the binned coordinates, we select nearby neighbors to construct a fiducial stream track, or automatically join neighboring bins based on the discussion in §4.1.4. This fiducial track is used to assign an ordering to the stars in the stream, by measuring the minimum distance between stream particles and the fiducial track. Each star is assigned an ordering parameter γ ∈ [−1, 1] based on where along the fiducial track they are closest to. The difference in γ-values for subsequent stars is set to a constant, such that the interval γ ∈ [−1, 1] is partitioned into N − 1 equal segments (where N is the number of stars belonging to the tidal feature). For instance, stars towards one end of the track will have γ closer to −1, while stars located at the opposite end will have γ close to 1. The final track used in the analysis is determined by fitting a spline to the now ordered stream. We utilize the Scipy interpolation library to fit a smooth spline to the projected stream, and differentiate the track using a representation of the same spline in PyTorch (Paszke et al. 2019) with automatic differentiation enabled. We fit a smoothing spline of degree k ∈ [3, 4, 5], with a large (e.g., 10 5 − 10 6 ) smoothing factor to ensure that the projected track characterizes the average shape of the stream rather than passing exactly through the data points.
This approach works for observations that resolve individual stars belonging to a tidal feature, however this is generally not the case. For unresolved streams, binned (i.e., average) photometric measurements can be utilized to fit a track to the tidal feature. As long as photometric measurements can capture the average shape of the stream, we expect our analysis to be applicable. This is the approach we take in §5, where we analyze a real tidal feature surrounding the galaxy NGC 5907. In a future work, we will consider more rigorous and automated techniques for fitting flexible and differentiable curves to the observed population of tidal debris (e.g., Starkman et al. 2022).

B. HYPERPARAMETERS
We now discuss the various hyperparameters used throughout this work, and their values. Most hyperparameters relate to the stream track and our fitting routine. We could consider a Bayesian determination of the stream track, allowing us to treat the hyperparameters in a probabilistic framework. However, the present work is aimed at introducing and validating the basic principles discussed in §2, so we postpone a more thorough exploration of how to effectively fit the stream track to future work.
In §3 we introduced the hyperparameters ∆D, ϵ, and σ θ T . These correspond to the maximum change in line-of-sight distance between adjacent evaluation points along the stream track, the zero curvature threshold (|dϕ/dℓ| < ϵ), and our prior uncertainty on the direction of the stream track when the tangent condition is satisfied, respectively. For ∆D, we adopt 30 kpc, since this is a reasonable upper bound for the line-of-sight distance between adjacent stream segments in our N −body simulations. For a stream at ≈ 17 Mpc (i.e, NGC 5907 considered in §5), a 0.035 deg change in the direction of the stream track per kpc corresponds to an ϵ ∼ 10. We adopt this value for ϵ, and σ θ T = 10 deg. In our tests on simulated data, we find this value of ϵ successfully identifies regions of the stream track with negligible curvature. Our usage of a relatively large σ θ T ensures that the tangent condition will not strongly disfavor planar acceleration vectors with angular separations of ∼ 1 − 4 deg from the local stream-track.

C. EXPANDED LIKELIHOOD
In this appendix we write out the full likelihood utilized in the analysis. The likelihood is already derived in §3, here we simply expand Eq. 19 to demonstrate a relationship with the trinomial distribution. Each of the conditional probabilities in Eq. 14 and Eq. 16 are zero if the condition (C j ) is not satisfied. That is, for each evaluation point along the stream (eachκ i ) only one term in the summation of Eq. 17 survives. In this way these probabilities have similar behavior to indicator variables, which are either zero or one. Eq. 18 is the fraction of times each conditional probability is non-zero. This means the full likelihood (Eq. 19) can be expanded as where N is the total number of evaluation points along the stream track and the product is over all N × f 3 evaluation points with undefined curvature vectors. The parameter f 1 is the fraction of evaluation points along the stream for which the curvature and planar acceleration vectors are consistent, f 2 is the fraction of evaluation points with inconsistent curvature and planar acceleration vectors, and f 3 is the fraction of evaluation points with no curvature. Note that f 3 is independent of the model parameters m, since the number of evaluation points without curvature does not change with the choice of a potential model. The terms with f 1 , f 2 , f 3 in Eq. C1 collectively form the trinomial distribution, since the analysis has three possible binary outcomes: either a trial acceleration is consistent with the curvature vector, not consistent, or the curvature vector is undefined. The last case introduces an additional probability component that is continuous, since in this scenario accelerations should align close to the stream track.

D. SPECIFYING EDGE CASES IN DISTANCE SAMPLING
In §3 we discuss the likelihood used in this work, which performs a maximization over the line-of-sight coordinate, z. In §3.2 we discuss how the line-of-sight coordinate is sampled to ensure the possibility of a continuous distance distribution along a stream. In general, there will be some distance interval, D i , over which the stream segment indexed i has |θ i | < π/2. This interval for the preceding segment, indexed i − 1, must overlap within some threshold ∆D of D i in order to ensure the possibility of a continuous distance track with a well-behaved gradient. Therefore, the range of allowable distances for the i th evaluation point is set by the preceding evaluation point with index i − 1. We refer to the range of acceptable distances for the i th evaluation point with the ∆D "padding" as the acceptable distance interval, represented by A i−1 (since the interval depends on the preceding segment).
We now discuss a few edge cases in this sampling scheme, and how they are addressed: • The first evaluation point (i = 0): For the first evaluation point the acceptable distance interval that we sample over is the full distance interval z ∈ [z min , z max ]. We usually take this interval as the maximum extent of the stream in the plane of the sky. For NGC 5907 in §5, we use a more conservative interval with z ∈ [−150, 150] kpc.
• Extended Stream Segment with Negligible Curvature: If both evaluation points i and i − 1 have negligible curvature, we definê which is the most probable line-of-sight distance for the i th evaluation point. The acceptable distance interval A i−1 is taken to be within 100 kpc of the preceding evaluation point with distanceẑ i−1 . The choice of 100 kpc is conservative, and ensures that we sample over a range of possible distances to identify a potential model that is compatible with the stream track.