Looking at subhalo descendents with merger trees

Nikita Agarwal
  • 3
  • 19 Nov '21

Hello,
To get the descendent subhalos for a particular subhalo for certain snapshots I used the merger tree.

tree = il.sublink.loadTree(basePath, 99,83280, fields=['SnapNum','SubfindID'],onlyMPB=True)
for i in (72,78,84,91,99):
    ind_in_tree = np.where(tree['SnapNum'] == i)
    print(tree['SubfindID'][ind_in_tree])

However, the central co-ordinate of the same subhalo for each snapshot is very different, i.e

Snapshot 99, (147.8469, 203.8065)
Snapshot 91, (70.7445, 90.4545)
Snapshot 84, (66.8501, 70.4983)
Snapshot 78, (45.0122, 72.8373)
Snapshot 72, (52.6200, 60.009)

This should not be the case, since I am looking at the same subhalo per snapshot. The coordinates for the subhalo center shouldn't vary so significantly.

Am I doing something wrong?

Regards,
Nikita

Dylan Nelson
  • 20 Nov '21

Hi Nikita,

I'm not sure what simulation you were looking at there, but I paste an example below. You can also load SubhaloPos from the tree for a positions. For the coordinates you posted, I am not sure why there are two numbers instead of three (x,y,z).

In [6]: basePath = 'sims.TNG/TNG50-1/output/'

In [7]: tree = il.sublink.loadTree(basePath, 99, 83280, fields=['SnapNum','SubfindID','
   ...: SubhaloPos'], onlyMPB=True)

In [8]: tree
Out[8]:
{'count': 91,
 'SnapNum': array([99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88, 87, 86, 85, 84, 83,
        82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66,
        65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49,
        48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32,
        31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15,
        14, 13, 12, 11, 10,  8], dtype=int16),
 'SubfindID': array([  83280,   83620,   84004,   86622,   82201,   81549,   82277,
          83746,   81946,   82359,   81557,   84216,   76716,   69740,
          67530,  884748,  885558,  889999,   62147,  896395,  898557,
         894604,  894596,  890424,  887849,  883330,  875701,  865191,
         865550,  853731,  849347,  836392,  831886,  825239,  832945,
         820080,  829874,  828649,  831956,  828719,  828229,  831229,
         830840,  692189,  682215,  686668,  692903,  821559,  820682,
         815768,  820814,  820792,  823132,  822028,  825022,  831973,
        1576808, 1597414, 1647585, 1638560, 1666959, 1701017, 1680158,
        1671475, 1668466, 1645838, 1646854, 1615013, 1622570, 1575420,
        1658614, 1685624, 1727175, 1759638, 1780759, 2068963, 2243086,
        2345681, 2273652, 2443668, 2496936, 2633568, 2871972, 4351044,
        3932623, 4594749, 4760792, 4378543, 6826114, 4694826, 6399119],
       dtype=int32),
 'SubhaloPos': array([[23769.98046875, 15789.94726562,  2830.8605957 ],
        [23782.39648438, 15768.8671875 ,  2854.22998047],
        [23798.84960938, 15727.43457031,  2898.10400391],
        [23806.98046875, 15693.10253906,  2932.76367188],
        [23814.19726562, 15630.109375  ,  2993.17919922],
        [23813.4921875 , 15579.96777344,  3039.77001953],
        [23802.0859375 , 15491.34863281,  3117.8840332 ],
        [23784.65429688, 15421.38476562,  3174.55249023],
        [23736.63476562, 15298.10742188,  3260.17333984],
        [23686.93554688, 15205.98925781,  3307.51782227],
        [23595.046875  , 15072.78320312,  3346.32739258],
        [23502.92578125, 14959.45800781,  3352.59228516],
        [23445.99023438, 14894.25976562,  3343.13330078],
        [23366.13085938, 14809.12109375,  3314.76123047],
        [23316.97070312, 14759.03808594,  3289.02075195],
        [23249.32421875, 14693.76464844,  3242.57495117],
        [23188.390625  , 14637.68261719,  3188.96264648],
        [23151.45703125, 14605.1171875 ,  3149.96630859],
        [23102.109375  , 14563.16601562,  3087.32983398],
        [23059.06445312, 14528.55371094,  3020.70068359],
        [23034.36328125, 14509.26269531,  2974.62060547],
        [22992.04882812, 14479.3828125 ,  2880.21435547],
        [22974.70117188, 14468.49121094,  2831.30664062],
        [22952.984375  , 14457.02734375,  2758.06542969],
        [22936.84375   , 14450.4765625 ,  2683.98999023],
        [22928.44335938, 14449.01171875,  2634.08642578],
        [22919.38671875, 14449.36523438,  2558.72851562],
        [22914.34960938, 14454.40625   ,  2482.22680664],
        [22913.26757812, 14462.1015625 ,  2406.72460938],
        [22915.76757812, 14471.91015625,  2330.01342773],
        [22921.63085938, 14483.11328125,  2255.12768555],
        [22930.83789062, 14495.86425781,  2180.078125  ],
        [22943.81640625, 14508.61035156,  2106.13354492],
        [22961.30273438, 14521.49902344,  2033.88195801],
        [22982.20898438, 14534.31542969,  1961.88537598],
        [23015.72851562, 14552.3984375 ,  1869.21362305],
        [23043.85742188, 14566.73632812,  1801.11547852],
        [23074.69140625, 14581.98535156,  1734.04199219],
        [23106.81054688, 14597.59472656,  1668.12780762],
        [23150.40429688, 14619.15234375,  1582.68371582],
        [23184.64453125, 14635.79394531,  1519.72351074],
        [23230.56445312, 14658.86035156,  1438.82958984],
        [23264.53125   , 14676.86523438,  1379.35986328],
        [23310.36328125, 14702.23339844,  1302.63024902],
        [23346.66992188, 14721.72265625,  1246.61242676],
        [23392.671875  , 14749.390625  ,  1173.3996582 ],
        [23438.02929688, 14777.14550781,  1101.8515625 ],
        [23481.27734375, 14807.55273438,  1032.56469727],
        [23511.14453125, 14831.90625   ,   981.82641602],
        [23555.32617188, 14872.56054688,   901.33203125],
        [23590.19726562, 14901.57324219,   839.97583008],
        [23626.34960938, 14928.62304688,   781.06536865],
        [23662.32226562, 14954.8984375 ,   725.28881836],
        [23697.72851562, 14980.61132812,   671.31115723],
        [23742.5       , 15011.22167969,   606.89892578],
        [23777.6015625 , 15034.77734375,   557.66577148],
        [23820.95898438, 15064.33984375,   498.69729614],
        [23864.15625   , 15093.31933594,   442.6998291 ],
        [23906.82226562, 15120.8984375 ,   389.62753296],
        [23964.859375  , 15158.54492188,   318.61901855],
        [23988.64257812, 15174.47558594,   291.19625854],
        [24036.1953125 , 15205.17382812,   234.15309143],
        [24075.41601562, 15229.94824219,   189.74401855],
        [24119.45703125, 15258.17578125,   139.65345764],
        [24162.48632812, 15285.32617188,    91.94963074],
        [24204.5859375 , 15312.06152344,    44.97717667],
        [24251.22265625, 15342.54980469, 34991.828125  ],
        [24296.95898438, 15369.4765625 , 34942.7578125 ],
        [24340.79492188, 15394.80566406, 34896.37109375],
        [24382.81054688, 15420.22460938, 34850.171875  ],
        [24428.16601562, 15447.78222656, 34800.74609375],
        [24472.08984375, 15473.70214844, 34753.61328125],
        [24519.73242188, 15502.62597656, 34700.19140625],
        [24565.77929688, 15533.07617188, 34649.78515625],
        [24594.5546875 , 15551.22460938, 34618.140625  ],
        [24659.63671875, 15591.08398438, 34545.921875  ],
        [24702.421875  , 15617.61523438, 34496.93359375],
        [24742.80273438, 15643.13476562, 34455.8671875 ],
        [24790.56445312, 15676.76367188, 34401.66015625],
        [24813.6875    , 15693.28027344, 34375.7734375 ],
        [24847.8125    , 15714.11914062, 34341.63671875],
        [24876.2265625 , 15733.55175781, 34311.1796875 ],
        [24915.21875   , 15756.48828125, 34274.515625  ],
        [24939.47460938, 15770.890625  , 34250.32421875],
        [24968.13085938, 15788.0546875 , 34220.7890625 ],
        [24996.28320312, 15805.14648438, 34193.3984375 ],
        [25009.8984375 , 15813.30371094, 34180.5625    ],
        [25045.625     , 15835.35058594, 34146.2421875 ],
        [25077.70507812, 15855.21972656, 34113.51171875],
        [25091.25976562, 15863.09179688, 34100.20703125],
        [25123.984375  , 15882.2421875 , 34063.5234375 ]])}
Nikita Agarwal
  • 20 Nov '21

Hi Dylan,

Thanks, I am running TNG 100-1 simulation and the coordinates I had posted were just (x,y).

I want to create plots showing the entire subhalo, and require the 'Coordinates' field in the Snapshots catalog for that.

To get the correct coordinates, I am using:

    gas_pos = (gas['Coordinates'])*ckpc/h 
    sub_pos = (subhalo['SubhaloPos'])*ckpc/h

    x = (gas_pos[:,0] -  sub_pos[0])/kpc
    y = (gas_pos[:,1] -  sub_pos[1])/kpc
    z = (gas_pos[:,2] -  sub_pos[2])/kpc

Where x,y,z are the coordinates. However, on running this and getting the plots, the coordinates at the center of the same subhalo at different redshifts are:

Snapshot 99, (147.8469, 203.8065, 145)
Snapshot 91, (70.7445, 90.4545, 100)
Snapshot 84, (66.8501, 70.4983, 82)
Snapshot 78, (45.0122, 72.8373, 62)
Snapshot 72, (52.6200, 60.009, 65)

I fail to understand why this is the case, particularly for snapshot 99. Kindly request you to please help with this.

Regards,
Nikita

Dylan Nelson
  • 20 Nov '21

Hi Nikita,

Can you post a complete code example of the problem?

I am not sure what "the coordinates at the center of the same subhalo" means exactly. In TNG100-1 at snapshot 99, subhalo ID 83280 has SubhaloPos of (23867.8, 33076.1, 41411.2). If you load all the gas particles that belong to this subhalo, their positions will be at similar Coordinates.

Nikita Agarwal
  • 1
  • 21 Nov '21

Hi Dylan.

Below is the code:

   basePath = 'sims.TNG/TNG100-1/output'

   snap_id = [72,78,84,91,99]
   subhalo_id = [272444,293577,313393,337662,64]
   for (i,j) in zip(snap_id,subhalo_id):
    gas = il.snapshot.loadSubhalo(basePath, i, j, 'gas', fields=None) 
    subhalo = il.groupcat.loadSingle(basePath, i, subhaloID=j)

    if ('Coordinates' in gas.keys() ):
        header = il.groupcat.loadHeader(basePath,i)
    a = 1 / (1 + header['Redshift'])
        ckpc = kpc*a 

        gas_pos =(gas['Coordinates'])*ckpc/h 
        sub_pos =(subhalo['SubhaloPos'])*ckpc/h
        x = (gas_pos[:,0] -  sub_pos[0])/kpc
        y = (gas_pos[:,1] -  sub_pos[1])/kpc
        z = (gas_pos[:,2] -  sub_pos[2])/kpc

        xc = (np.max(x)-np.min(x))/2
        yc = (np.max(y)-np.min(y))/2
        zc = (np.max(z)-np.min(z))/2

        rho = (gas['Density'])*((1e10*Msun/h)/(ckpc/h)**3)

        B = gas['MagneticField']*(h/a**2)*np.sqrt(1e10*Msun/kpc)*1e5/kpc
        B_mean = np.mean(np.sqrt(B[:,0]**2+ B[:,1]**2+ B[:,2]**2))
        B_std = np.std(np.sqrt(B[:,0]**2+ B[:,1]**2+ B[:,2]**2))

    nPixels = [300,300]

        extent = [xc-100,xc+100,zc-100,zc+100]

        grid, _, _, _ = binned_statistic_2d(x,z, rho, 'sum', bins=nPixels)

        fig, ax = plt.subplots(figsize=[10,10],nrows=1, ncols=5)

        fig, ax = plt.subplots(figsize=[8,8])

        im = ax.imshow(np.log10(grid),extent=extent,cmap='viridis',interpolation='nearest', aspect =  nPixels[1]/nPixels[0],origin = 'lower', vmin = -26, vmax = np.max(np.log10(grid)))#,origin = 'lower')

        plt.xlim(xc-100,xc+100)
        plt.ylim(zc-100,zc+100)

        plt.title("Gas Density Plot with Magnetic Field Lines \n Redshift={r}\n Magnetic Field Mean (G) ={m}\n Magnetic Field Std(G) ={s}".format(r=header['Redshift'], m = B_mean, s = B_std))
        fig.colorbar(im,label = 'Gas Density (g/cm^3)')

The output for the code is as:

Screenshot from 2021-11-21 13-30-14.png

Screenshot from 2021-11-21 13-30-24.png

Screenshot from 2021-11-21 13-30-30.png

Screenshot from 2021-11-21 13-30-39.png

Screenshot from 2021-11-21 13-30-45.png

As you can see, the subhalo position significantly varies from z = 0.39 to z = 0

Regards,
Nikita

Dylan Nelson
  • 1
  • 21 Nov '21

I would suggest to simplify, to help understand if there is a problem. What about

basePath = 'sims.TNG/TNG100-1/output'
snap = 99
subhalo_id = 64

gas_pos = il.snapshot.loadSubhalo(basePath, snap, subhalo_id, 'gas', fields=['Coordinates'])
subhalo = il.groupcat.loadSingle(basePath, snap, subhaloID=subhalo_id)

sub_pos = subhalo['SubhaloPos']

# note: periodic boundaries of the simulation box need to be considered in general
dx = gas_pos[:,0] -  sub_pos[0]
dy = gas_pos[:,1] -  sub_pos[1]
dz = gas_pos[:,2] -  sub_pos[2]

fig, ax = plt.subplots(figsize=[8,8])
ax.set_xlabel('dx [ckpc/h]')
ax.set_ylabel('dz [ckpc/h]')

ax.plot(dx, dz, '.')
Nikita Agarwal
  • 22 Nov '21

Hi Dylan,
Thanks, I used your suggestion and checked for individual subhalos, but the problem still persists.
I shall try checking for other subhalo IDs as well.

Also, can you please clarify what you meant by the comment 'periodic boundaries of the simulation box need to be considered in general'

Regards,
Nikita

Dylan Nelson
  • 22 Nov '21

Dear Nikita,

Please see this previous thread about the periodic box.

Otherwise, I am not sure exactly what the problem is. The code above produces the output I would expect:

TNG100-1_99_64_gas_xz.png

Nikita Agarwal
  • 23 Nov '21

Hi Dylan,
Thanks, I am getting the same plot.

Your inputs helped solve my problem.

Regards,
Nikita

  • Page 1 of 1