  We present a comprehensive analysis of the shape of dark matter (DM) halos in
a sample of 25 Milky Way-like galaxies in TNG50 simulation. Using an Enclosed
Volume Iterative Method (EVIM), we infer an oblate-to-triaxial shape for the DM
halo with the median $T \simeq 0.24 $. We group DM halos in 3 different
categories. Simple halos (32% of population) establish principal axes whose
ordering in magnitude does not change with radius and whose orientations are
almost fixed throughout the halo. Twisted halos (32% of population), experience
levels of gradual rotations throughout their radial profiles. Finally,
stretched halos (36% of population) demonstrate a stretching in their principal
axes lengths where the ordering of different eigenvalues change with radius.
Subsequently, the halo experiences a "rotation" of $\sim$90 deg where the
stretching occurs. Visualizing the 3D ellipsoid of each halo, for the first
time, we report signs of re-orienting ellipsoid in twisted and stretched halos.
We examine the impact of baryonic physics on DM halo shape through a comparison
to dark matter only (DMO) simulations. This suggests a triaxial (prolate) halo.
We analyze the impact of substructure on DM halo shape in both hydro and DMO
simulations and confirm that their impacts are subdominant. We study the
distribution of satellites in our sample. In simple and twisted halos, the
angle of satellites' angular momentum with galaxy's angular momentum grows with
radius. However, stretched halos show a flat distribution of angles. Overlaying
our theoretical outcome on the observational results presented in the
literature establishes a fair agreement.
