Sorry, I cant really offer you a solution in ArcGIS. If you think through the problem you will quickly realize that a one-to-one match is not possible in a flat attribute table. Each polygon will have it's own unique number of 1st, 2nd, ... order neighbors. You can do a many-to-one relate but need to generate the neighbor contingencies first.
Here is an example of creating Nth order Polygon contingencies using R and the spdep package.
library(spdep)
library(rgdal)
columbus <- readOGR(system.file("shapes/columbus.shp", package="spData")[1])
Using the columbus polygon data we create first-order neighbors (those touching a source polygon) using the spdep::poly2nb
function and second-order neighbors (those touching the polygons touching a source polygon) using the spdep::nblag
function. These neighbor objects are then coerced into sparse matrices.
nb.1st <- nb2listw(poly2nb(columbus), style = "B", zero.policy = TRUE)
W1 <- as(nb.1st, "CsparseMatrix")
nb.2nd = nblag(poly2nb(columbus), 2)
W2 <- as(nb2listw(nb.2nd[[2]], style = "B", zero.policy = TRUE), "CsparseMatrix")
Now we can index a polygon that we are interested (i) in and use the neighbor matrices to find the associated Nth order neighbors. The vectors created bellow are tracking the row.names of the original polygon object. As such, we are creating an index based on matching row.names between the polygon and matrix object given [W>0].
i = 5 # polygon index to check neighbors
nb.1st.idx <- colnames(W1)[which(W1[i,] > 0 )]
nb.1st.idx <- which(row.names(columbus) %in% nb.1st.idx)
nb.2nd.idx <- colnames(W2)[which(W2[i,] > 0 )]
nb.2nd.idx <- which(row.names(columbus) %in% nb.2nd.idx)
We now have the corresponding 1st and 2nd order polygons for polygon 5 and can plot the results.
plot(columbus, border="grey")
plot(columbus[nb.2nd.idx,], col="cyan", add=TRUE)
plot(columbus[nb.1st.idx,], col="blue", add=TRUE)
plot(columbus[i,], col="red", add=TRUE)
box()
title("First and Second Order Neighbors for polygon 5")
legend("topleft", legend=c("source polygon","1st order","2nd order"),
fill=c("red","blue","cyan"))
Keep in mind that this is only for one polygon. To do this for all polygons you would have to iterate through all polygons and store the Nth order polygon indexes in a list and then figure out how to relate them to the original data.