I think you want to do something like Generate Near Table or just Near which will allow you to join attributes between the points for each census tract and the near table to find the distance to the city centres.
If you use closest only on the near table you will only get the closest feature, Near will always return just the nearest feature and distance from one feature class to another, in addition it should add Near_FID to the input points so you can identify which one by joining Near_FID to the FID on the city centeres.
Addendum
To find the distance between census tracts and the matching code number:
Generate near table using tracts as input and city centres as near features with all distances and a very large or not entered tolerance; the documentation says If no value is specified, all near features will be candidates which should mean that every census tract to every city will be represented in the table should no tolerance be entred, if it is not the case then use a very large tolerance.
From here there are many ways to trim that down, but this is the method I use as it's less confusing (for me anyway):
Add two fields to the near table, let's call them TractCode and CityCode. Load up the table in ArcMap (you will need to use add data as tables don't want to drag 'n drop), add the census tracts and city centres as well.
Open the attribute table for the near table and Add a join from the near table to the city centres using NEAR_FID to FID/OBJECTID (depending on data type). Using field calculator calculate the value of the CityCode in the near table from the code field in the city centres data. Now remove that join and add another one to the census tracts and calculate the TractCode from the census tracts and then remove that join too.
Now your near table has the codes from the census tracts and city centres. Select by attributes on the near table where CityCode = TractCode which will select the rows where the code is the same on city centre and census tract, export that to a new table for simplicity.
This newly exported table will contain the distances between census tracts and city centres where the code matches.
Best Answer
If you need to use the road distance as the cost of travel, then you would indeed need to use Network Analyst if you are using ArcGIS and want to take advantage of the restrictions imposed on the roads in the network dataset (one-way, turns, etc.).
The solver you would need for this is OD Cost Matrix. You basically need to generate a matrix where each point is routed to all other points. The OD matrix is usually used for modelling business cases where you have a couple of origins (warehouses) and many consumers (retail network points) and you want to report the travel costs between the origins and the destinations.
In your case, however, because everything is stored in the same feature class, you need to load your feature class points both into origins and destinations (i.e., they will be identical). This is a valid approach, and I am using the OD cost matrix this way quite often myself, too.
The number of lines generated for a feature class with
Npoints will be(N*N) - N. This is because after the solve, you would exclude those pairs where the source and the origin are the same. You can filter them out easily by using the SQL query on the output cost (travel time or distance) field (traveltimefield > 0).Performance wise, it would be wise to choose None for the Output shape type option. I don't know how many points you have in your feature classes, but ~1500 points both in origins and destinations (resulting in ~2.4 mln lines) never was a problem for me to run in just a couple of minutes (I have a decent laptop with i7 CPU though). If you have more points in your feature class, you can run out of memory in ArcMap (it's a 32bit application that cannot take more than just a couple of GBs of RAM).
If you will hit a RAM limit of ArcMap, you might consider running this solve from Python 64bit process. This is what I usually do when processing large road networks (the RAM consumption for the running Python process was about 12GB at peak, no problem).
There are other options of generating the OD matrix outside of ArcGIS if you are OK with just one-way restrictions and no access to other ArcGIS network dataset based restrictions (vehicle height, turns, etc). networkx - a Python package for working with graphs is something I use often for network analysis, but you would need to program Python then.