This question likely should be split into two (how to intersect overlapping polygons, and how to split results proportionally), but here is my method, inspired largely by Calculating zonal statistics of raster data in multiple overlapping zones and combining them into one table, which is more recent.
- Calculate population density for each census tract (people per square metre, people per square mile, etc.). Works best if density unit matches cell size used in next step.
- Convert census tract to raster, using density as cell value. Now, each cell represents the actual number of people contained within, on average.
- Run following Python code to calculate Zonal Statistics for each buffer, then merge all results together.
store_buffs = r"C:\junk\overlap\store_buffer.shp"
density_ras = r"C:\junk\overlap\density.tif"
table_list = 
with arcpy.da.SearchCursor(store_buffs, ["FID"]) as cursor:
for row in cursor:
exp = '"FID" = ' + str(row)
temp_table = r"in_memory\temp_table" + str(row)
temp_shp = r'in_memory\temp_shp'
arcpy.Select_analysis(store_buffs, temp_shp, exp)
arcpy.sa.ZonalStatisticsAsTable(temp_shp, 'FID', density_ras, temp_table, 'DATA', 'SUM')
final_table = r"C:\junk\overlap\pop.dbf"
- Join the output table to your store points to get total estimated population within store buffer