$JetCurry$ I. Reconstructing Three-Dimensional Jet Geometry from Two-Dimensional Images

We present a three-dimensional (3-D) visualization of jet geometry using numerical methods based on a Markov Chain Monte Carlo (MCMC) and limited memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimized algorithm. Our aim is to visualize the 3-D geometry of an active galactic nucleus (AGN) jet using observations, which are inherently two-dimensional (2-D) images. Many AGN jets display complex structures that include hotspots and bends. The structure of these bends in the jet's frame may appear quite different than what we see in the sky frame, where it is transformed by our particular viewing geometry. The knowledge of the intrinsic structure will be helpful in understanding the appearance of the magnetic field and hence emission and particle acceleration processes over the length of the jet. We present the $JetCurry$ algorithm to visualize the jet's 3-D geometry from its 2-D image. We discuss the underlying geometrical framework and outline the method used to decompose the 2-D image. We report the results of our 3-D visualization of the jet of M87, using the test case of the knot D region. Our 3-D visualization is broadly consistent with the expected double helical magnetic field structure of the knot D region of the jet. We also discuss the next steps in the development of the $JetCurry$ algorithm.


Introduction
Relativistic jets transport mass and energy from subparsec central regions to Mpc-scale lobes, with a kinetic power comparable to that of their host galaxies and active galactic nuclei (AGNs). This profoundly influences the evolution of the hosts, nearby galaxies, and the surrounding interstellar and intracluster medium (Silk et al., 2012;Fabian, 2012). The generation of such flows is tied to the process of accretion onto (likely) rotating black holes, where the magneto-rotational instability can couple the black hole's spin and magnetic field to the disk or flow to produce highlatitude outflows at speeds close to the speed of light (Meier et al., 2001). While these jets have a dominant direction of motion (i.e., outward from the black hole), they often have bends and features (both stationary and moving) that are either perpendicular or aligned relative to the jet at some angle. Deciphering the true nature of these features and their geometry, relation, and dynamical meaning within the flow is a difficult problem, as any astronomical images we acquire are of necessity two-dimensional (2-D) views of three-dimensional (3-D) objects.
The problem of reconstructing 3-D information from 2-D images is common to many fields, but it is particularly critical in astronomy. In most other cases, e.g., medical imaging, one may take images of a source from multiple viewpoints to aid reconstruction. However, this is not possible in astronomy, so we must rely on other methods. For instance, Steffen et al. (2011), Wenger et al. (2012), Wenger et al. (2013), Cormier (2013), Sabatini et al. (2018), and Lagattuta et al. (2019) used symmetries inherent in, respectively, planetary nebulae and galaxies, plus 2-D images, to infer and reconstruct 3-D visualizations of these objects. This field is, in fact, rapidly growing in astronomy, as can be seen by the vast number of subjects explored on the 3DAstrophysics blog 2 .
In astrophysical jets, the problem is rather different. Unlike in galaxies or planetary nebulae, we cannot make assumptions such as spherical, elliptical or disk symmetry, or rotation. However, we can assume a dominant direction of propagation. As an example of the typical knotted structure of AGN jets, we show in Figure 1, a broad view of the M87 jet, one of the nearest of the class at 16.7 Mpc distance, taken from Meyer et al. (2013). In every single image, the M87 jet shows an amazing complexity of features, including knots, helical undulations, shocks, and a variety of other morphological structures, many of which are oriented at some odd angles with respect to the overall jet direction. As shown by Meyer et al. (2013), some of the features in the M87 jet seem to move with apparent velocities up to about 6 within the inner 12. ′′ 0 of the jet, with a general decline in apparent speed with increasing distance from the nucleus. However, there are some nearly stationary components that are largely located near the upstream ends of knots. Additionally, the polarimetric imaging of Avachat et al. (2016) shows apparent helical winding structures to the inferred magnetic field vectors in several knots. These features are clues to complex jet dynamics, but are difficult to interpret properly.
In this paper, we describe a geometrically based code that attempts to reconstruct the 3-D structure of jets, starting from 2-D images. This is an evolving project that in later stages will attempt to use kinematic information as well as incorporate special relativistic corrections so that foreshortening, Doppler boosting, and superluminal motion can be included. Our goal is to provide a firmer geometrical grounding to these modeling efforts by allowing reconstruction of a jet's structure in 3-D.

Geometrical Framework
To visualize the geometry of the jet in 3-D, the key parameters to consider are the distance between any two features and the apparent angle between them with respect to the direction of the jet's axis. The distance between the core of the jet and the location of interest on the jet is defined as , and the angle of the location to the core is . Another parameter is the line of sight (LOS).
We assumed that the 2-D projection of the jet's axis lies along the -axis and measured the angles with respect to the positive direction. We considered and as our known parameters, which can be directly measured from the images, i.e., in the sky frame. A third parameter, assumed to be known (albeit from other information such as a vs. plot based on the observed superluminal motion, where is the space velocity and is the viewing angle with respect to the LOS) is the angle the jet's propagation axis makes with respect to the LOS.
Following Conway and Murphy (1993), Figure A-1 shows the relevant geometry for a single bend within a jet, and specifically how the 2-D sky frame can be related to the jet's frame, which is inherently 3-D. All primed points represent the observed, sky-frame projection we see, with the components lying at A ′ and B ′ in that frame but at points A and B in the jet's frame. The point D ′ is the projection of B ′ on the + axis, and the point D is projection of B on ( , ) plane. Angle is the angle between segment A ′ B ′ and the + axis. From point B draw a line BC perpendicular to the jet axis (OC) and set the ∠CAB as . Point C is then where line BC crosses line OC perpendicularly, so ∠BCA is 90 • . Segment BC makes an angle with the ( , ) plane (i.e., ∠BCD = ). This way, ΔABC is raised off the ( , ) plane through angle , while the segment AC still lies in the ( , ) plane. Segment AC makes an angle with the LOS, which is assumed to be along the + axis. The distance between points A and B in the jet's frame is , while is the projection of on the ( , ) plane, i.e., the distance between A ′ and B ′ . Angle is the apex angle of ΔBAD, and is the angle between triangles BAD and FAE. Finally, ΔAGH is the projection of ΔABD on the ( , ) plane.
To simplify the algorithm both computationally and physically, for now we assume that the jet is non-relativistic. The LOS effects in addition to various relativistic effects can enhance the intensity and shift the frequency observed, as well as change the comparison between geometry in the jet and observer frames (Böttcher, 2012). These effects will be included in the next version of JetCurry.

Non-linear Parametrized Equations
We used a set of non-linear parametrized equations containing the angles and distances described above. Assuming the non-relativistic jet flow and using the geometry in Figure A-1, we derived the following non-linear equations including three known parameters ( , , ), and five unknown parameters ( , , , , ).
If a local jet structure has a smaller bend; i.e., < ( 2 − ), the transformation is: If a local jet structure has a larger bend; i.e., ≥ ( 2 − ), the Equations (4) and (5) modify as: d cos cos + s cos tan = d sin cos sin (6) d 2 = s 2 sin sin 2 + s 2 cos cos 2 sin 2 ( − ) Our aim is to solve for the angles , , , , and the distance between the knots . The system is underdetermined, and so cannot be solved exactly. However, as we shall describe, by making use of the angle , distance and angle of LOS as known parameters, we can derive the solution space as well as the relative probability of various bend parameters.
Please note that Conway and Murphy (1993) did not give either these equations or a derivation of them, nor did they try to derive more detailed 3-D information about any individual jets. Their aim was to attempt to understand, in a geometrical sense, the misalignment of arcsecond-scale features in blazar jets with the features seen on milliarcsecond scales by VLBI arrays.

JetCurry
JetCurry is a Python-based code to visualize the 3-D jet geometry of an AGN jet from its 2-D image. It incorporates the results of the Wavelet-based Image Segmentation and Evaluation (WISE, Mertens and Lobanov, 2014) method as a set of input parameters. These input parameters consist of the centroid locations of the identified components in the 2-D images. The code uses a nonlinear solving algorithm (emcee) to solve the non-linear parametrized equations for five unknown parameters: , , , , and (Equations 1-7). Then, it performs principal component analysis (PCA) and acquires the dominant direction of propagation of the knot D region from the M87's core with respect to our LOS. These tasks are discussed in detail in the following subsections. Our code is available for free as an interactive Jupyter Notebook on Github 3 .
We ran JetCurry on a radio image of knot D in the M87 jet, chosen as an example of a small region of a relatively complex jet. This is done because the knots have varying and complex flux structures that are separated by considerable distances (see Figure 2). Hence, the jet stream is better visualized by focusing on individual knot regions at a time instead of the entire jet. This restricts the projected jet track from wandering in regions in between knot complexes. 3 https://github.com/esperlman/JetCurry.

Image Decomposition and Component Identification
JetCurry employs the WISE method for multiscale structural decomposition and morphological identification in astronomical images (Mertens and Lobanov, 2014). The WISE algorithm implements segmented wavelet decomposition (SWD) to statistically extract significant structural patterns (SSP) at user-specified wavelet scales. It decomposes the image into a set of sub-bands (wavelet scales) and performs a wavelet transform to acquire significant wavelet coefficients against a noise threshold. Then, it uses these coefficients to extract the local maxima coordinates and applies watershed segmentation to retrieve the 2-D boundaries around the corresponding SSP features.
This methodology successfully identifies complex morphological components, including optically thin and partially overlapping structural features, that are otherwise undetected with standard object-recognition methods (e.g., Belongie et al., 2002;Lobanov et al., 2003;Bach et al., 2008). Figure 2 illustrates the centroid locations of the identified SSP features in the knot D region of the M87 jet at a wavelet scale of 50 milliarcseconds (mas).
The WISE method uses multiscale cross-correlation to detect, classify, and track different SSP features across a series of multi-epoch images. The results of different tests conducted by Lobanov (2014, 2016) demonstrate the robustness and reliability of the WISE method in identifying and tracking components of the M87 jet that undergo rotation, deformation, and segregation. For the next generation version of JetCurry, we will incorporate this capability to analyze the 3-D morphological and kinematic evolution of relativistic jets.

Assumptions, Constraints, and Non-linear Solvers
From an observed image in 2-D, we can only measure distances and angles , as shown in Figure A-1. We assumed that the LOS viewing angle for the M87 jet is 14 • with respect to the observer (Biretta et al., 1999;Perlman et al., 2011). It is important to that the results acquired from a parsec-scale jet kinematics study of AGNs as part of the 2 cm Very Long Baseline Array (VLBA) survey and Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) programs suggest that the values of between 70 • and 90 • are very unlikely (Lister et al., 2019). The bestfit Monte Carlo simulations of the 1.5 JyQC quasar sample indicate that the most jets are viewed at less than ∼ 10 • . This is also in agreement with the analysis carried out by Vermeulen and Cohen (1994), Lister and Marscher (1997), and Cohen et al. (2007).
The algorithm uses Equations 1-7 to solve for the most probable values of the unknown parameters, , , , , and using a non-linear solver, emcee, to explore the solution space. Emcee (Foreman-Mackey et al., 2013) is a highly efficient, open-source Python package that uses a Goodman & Weare affine-invarient MCMC Ensemble Sampler, aiming to find a global minimum solution. We used MCMC methods because of the underdetermined nature of the solution space, particularly because of the complex nature of our non-linear trigonometric equations. The number of times emcee is run depends on the user-preference (we are choosing three iterations). The initialization parameters for emcee are 5 dimensions, 1024 walkers, and 50 steps. For our current version of JetCurry, we are assuming that the observed jet structures have smaller bends. Therefore, the angles , , and are set to be explored from 0 to 1.57 radians, with a step size of 0.5 radians. Consequently, the angle to restricted to [20 • , ∕2− ]. The distance is set to be floor( ) to floor( )+ 81 parsec (pc). After this, the previous results are set as the initial guess for the next trial until a more defined range of possible solutions is found.
Once emcee has located the global minima/maxima regions, we used a nonlinear solver. To ensure the solutions of the variables are real, we took the logarithm of the equations. We found that the limited memory BFGS algorithm (Broyden, 1970;Fletcher, 1970;Goldfarb, 1970;Shanno, 1970) converges to the real solution faster and more accurately in our test cases compared to other nonlinear solvers. The algorithm is in the same class as Quasi-Newton methods and hill climbing with bounds on the solutions. It can be applied to a general non convex function that has continuous second derivatives. In JetCurry, the bounds in the nonlinear solver are set to be within 10 steps of the nonlinear solver by default. We present further details of our case study of knot D in § 5.

Principal Component Analysis (PCA)
JetCurry applies PCA (see Ivezic et al., 2014, pg. 289-297, for derivation and algorithm) and extracts the underlying preferred direction of the knot D region with respect to the LOS. The PCA incorporates the M87 jet's core as the origin and finds a dominant direction of propagation in 3-D. This direction is mathematically equivalent to a regression line that minimizes the square of the perpendicular distances from the corresponding Cartesian coordinates (Ivezic et al., 2014, pg. 292). The solid red line in Figure 5 represents the preferred direction of propagation of the knot D region from the M87 jet's core with respect to the LOS ( = 14 • ).

Testing JetCurry
To test the accuracy of JetCurry, we simulated 100 bends using different combinations of , , , , and . For our simulated test cases, we assumed that bends greater than 90 • were most probably unrealistic. We chose = 50 pc and = 14 • . Additionally, we restricted the values of and to the range [0, ∕2] and [20 • , ∕2 − ], respectively. We then calculated a range of values of ( , ), which are measurable variables in our algorithm. These values of ( , ) were plugged back into our algorithm to test how well we could reproduce the given values of ( , ).
The resulting values of , , and are represented as absolute error distribution plots, which are shown in Figure 3.  The corresponding absolute errors are given by the color bars. The annotated circular markers are the acquired values of ( , ) for the identified components in the knot D region (see Figure 2). For the cases where the solutions for ( , ) converged successfully, JetCurry reproduced the expected bend parameters with the mean absolute errors of 18.16 • , 13.41 • , and 0.49 pc in , , and , respectively. However, it produced relatively large discrepancies for 0 • < < 20 • and 0 • < < 20 • . This is due to the nature of the equations, and particularly where the expressions in the denominators have singularities. This makes the probability distribution in these regions complex. As a result, MCMC methods do efficiently work backward to find the global minima. Figure 4 shows a corner plot of the marginalized total probability distribution of , , , , and for a simulated bend with = 44.5 pc, = 54.98 • , and = 14 • . The expected values for ( , , , , )

Results
We now present a case study of a small region of the M87 jet (i.e., knot D, shown in Figure 2), which displays a complex morphology. For this case study, we used the radio image from Avachat et al. (2016). It has a scale of 0. ′′ 025/pixel that translates to 2.02 pc/pixel at a distance of 16.7 Mpc. The knot has 3 sub-components, namely, D-east, D-middle and D-west, as well as multiple apparent bends. Knot D-east appears to be at an angle from the northern edge of the jet cross-section to the southern edge over a distance of about 25 pixels (∼ 0. ′′ 625). Near its downstream end, it appears to split into northern and southern branches, out of which the southern branch is brighter and has been identified as knot D-middle.
JetCurry uses the centroid coordinates of the identified components (see Figure 2) to calculate the values of and with respect to the jet's core ( core , core ) and the LOS angle (which we assume is 14 • ). It solves the non-linear Equations (1) through (7) and outputs the range of values for each unknown parameters, i.e., angles , , , and the distance . To plot the posterior probability distribution of the values of unknown parameters, we make corner plots for each bend in the knot (similar to Figure 4). Each corner plot is a multi-dimensional representation of projections of the posterior probability distribution of the parameter space (Foreman-Mackey, 2016). We interpreted the MAP values as the actual solutions for angles and distances, although this can be subject to irregularities near places where the function is discontinuous or nonlinear (see e.g., §5.2). The gray scale represents the output probability in parameter space, with higher probabilities corresponding to darker colors.

Solutions in 3-D Cartesian Space
We converted the data obtained from our main algorithm to 3-D Cartesian coordinates (i.e., , , and ). These coordinates make use of the most probable values taken from the corner plots obtained for each bend. We used the following transformation equations to calculate the required ( , , ) coordinates for our 3-D visualizations: = cos cos (8) = sin cos (9) = cos cos tan (10) where , , , and are the acquired angular bend parameters from JetCurry. Figure 5 shows 3-D visualizations for the knot D region of the M87 jet. These are snapshots of Plotly's interactive version from three different perspectives. Figure 6 illustrates the projection of these 3-D visualizations onto the ( , ) plane.

Visualizations
Each circular marker (dark red, blue, and purple) is a combination of the most probable values of , and for each bend. The dark red and blue markers belong to the northern and southern streams, respectively. However, based on the visual appearance of the knot D region (see Figure 2) and polarimetry data presented in Avachat et al. (2016), we chose a set of points to be in both streams, which are indicated as purple markers. All circular markers are annotated Figure 5: Snapshots of the 3-D visualizations for the knot D region of the M87 jet from three different perspectives. The circular markers represent the most probable solutions of the non-linear parametrized equations. The dark red and blue markers belong to the northern and southern streams, respectively. The purple markers belong to both streams. All circular markers are annotated with 1 , 2 , and 3 standard deviational ellipsoids, respectively. The dark red and blue dashed curves represent spline fits to the trajectories of the helical paths. The dark red, blue, and purple dotted lines act as a guide to the eye to follow the corresponding projections from 2-D to 3-D. The solid red line represents the preferred direction of knot D from the core ( = −84.23 pixels, = 20.69 pixels) of the M87 jet. The solid orange line corresponds to the LOS with respect to the z axis ( = 14 • ). The 3-D conical structure with a power-law index of 0.96 ± 0.1 illustrates the observed shape of the jet in this region (Asada and Nakamura, 2012). The scale on the 2-D radio image corresponds to 1 pixel = 0. ′′ 025.
with their corresponding 1 , 2 , and 3 standard deviational ellipsoids. These ellipsoids signify the spatial dispersion of the acquired bend parameters in the 3-D Cartesian space. The dark red and blue dashed curves are the spline fits to the northern and southern streams, respectively. We followed the regression analysis in Gilat and Subramaniam (2007, pg. 193-250); we found out that these curves are the best fits to the data since they yield the smallest total errors. The solid orange line represents the LOS with respect to the axis at an viewing angle of 14 • . The solid red line shows the visualized dominant direction of propagation of the knot D region from the M87 jet's core.
Based on the milliarcsecond-arcsecond structure of the M87 jet described in Asada and Nakamura (2012), we overlaid our visualizations with conical streamlines. Since the knot D region lies downstream of the HST-1 complex, we adopted a 3-D conical structure (shown in gray) with a power-law index of 0.96 ± 0.1 (Asada and Nakamura, 2012). This helped us further validate the 3-D bend parameters and their corresponding curve fits.
We also made a projection of the acquired 3-D Cartesian coordinates of the bend parameters onto the ( , ) plane. To acquire this projection, we used the same equations for and as in Equation (8), and assumed = 0. The dark red, blue, and purple dotted lines act as a guide to follow the corresponding projections from 2-D to 3-D. This 3-D interactive plot is available online.

Discussion
Based on the most probable values of our bend parameters, JetCurry has produced a very interesting 3-D visualization of the knot D region of the M87 jet. It is worth noting that our 3-D visualization is based purely on geometrical considerations. Thus it is different from physical modeling efforts, which while they start from well known physical principles, have the possible issue that they sometimes do not exactly represent the observed morphologies. A purely geometric visualization can be complementary to such modeling efforts, particularly once the appropriate effects produced by special relativity are included.
While not unique, this 3-D visualization allows us to make several comments regarding the nature of this knot complex. First of all, the acquired 3-D bend parameters are broadly consistent with the observed conical structure of the M87 jet downstream of the HST-1 complex (Asada and Nakamura, 2012). In fact, it appears likely that the bright regions of knot D lie on at least two filaments (i.e., the northern and southern streams) that are either wrapped around the jet or threading the jet cone, as the two streamlines ( Figures  5 and 6) broadly trace the outer edges of the cone. This supports the idea that the knot regions represented helical filaments was first proposed by Owen et al. (1989). Additionally, the radio and optical polarimetry images presented in Avachat et al. (2016) clearly show the apparent helical structure throughout the jet, particularly in the projected magnetic field vectors seen in knots HST-1, D, A and B. Although it is not the only possible way to visualize the jet, it has stood the test of time remarkably well, as discussed by Avachat et al. (2016).
The natural next step in improving JetCurry is to incorporate special relativity. If indeed the knot regions are helical filaments wrapped around or threading the jet, a number of special relativistic effects could be observed, such as brightening and apparent acceleration when components' direction of motion crosses the LOS, and dimming and apparent deceleration when components move away from us. This has the ability to break the limitations noted above, as different trajectories would make radically different predictions for the jet's brightness and apparent motion as a function of time. This could add context to the models that have been proposed (e.g., Hardee and Eilek, 2011) for certain features within knot D as structures downstream of the shock at the eastern end.
We also plan to use JetCurry to visualize the other regions of the M87 jet, especially knots A and B, where our polarimetric images point toward the presence of an intertwined double helix (Avachat et al., 2016). This kind of structure is also suggested by previous works. By adding the proper motion constraints (e.g., Meyer et al., 2013), the observed streamline structures (e.g., Asada and Nakamura, 2012), and the special relativistic considerations (e.g., LOS, Doppler boosting, flux variability, observed superluminal motions, and foreshortening), we can further restrict the ranges of our unknown parameters. These constraints will also help understand the differences in the morphology and hence the emission mechanisms in the inner and outer regions of the M87 jet.
We thank Dr. Jeremy Riousset 4 for useful discussions on PCA and 3-D curve fits. This research was supported by the National Science Foundation (NSF) grant AST-1716507.

Appendix: Derivation of Non-linear Parametrized Equations
We consider a single bend in the jet at a viewing angle as shown in Figure A-1. A and B are neighboring knots in the jet's frame, while A ′ and B ′ are the projections of knots A and B in the sky frame with the origin at O. The distance between A and B in the jet's frame is , and the apparent distance between A ′ and B ′ is . The projected bend angle relative to the -axis in sky frame is represented as . ∠BAC is set as .
When < ( 2 − ), the equations describing the bend geometry is derived below: Starting with the equation of Conway and Murphy (1993) and using the geometry shown in the top plot of where , and are labeled in Figure A-1, while is the angle of LOS to the observer. △AEF (same as △A ′ B ′ D ′ ) is the projection of △ABD onto the ( , ) plane (i.e., sky frame), and △ADE is the projection of △ABF onto the ( , ) plane. In the right △AEF: Because △AGH is the projection of △ABD onto the ( , ) plane, and △AHD is the projection of △ABG onto the ( , ) plane. In the right △ABG: = ∠BAG In the right △AGH: AG 2 = AH 2 + GH 2 (16) = = cos tan (17) = = ′ ′ = sin (18) 2 cos 2 = s 2 sin 2 + s 2 cos 2 tan 2 Equation (14) can be simplified by using Equation (12): cos 2 = cos 2 sin 2 + cos 2 cos 2 tan 2 (20) Equations (13) and (15) can be combined into: tan tan 2 = cos 2 In right triangle CID: = sin = d sin cos sin and AC cos = AH + CI therefore; cos cos = s cos tan + d sin cos sin Also, in right △ABC, AB 2 = AC 2 + BC 2 therefore; 2 = s 2 sin sin 2 + s 2 cos cos 2 sin 2 ( + ) From these, we get the following set of five equations, Equations (25) through (29), in five unknowns ( , , , , and ) along with the measurable parameters ( and ), and the assumed angle of LOS ( = 14 • ).
When the local jet structure has large , that is when ≥ ( 2 − ) (the geometry for this scenario is shown in the bottom plot of Figure The remaining three Equations (27) -(29), stay the same for both the cases.