Replicating Simulation Study of Gary King in R In this article, the authors describes the following simulation study of Propensity Score Matching,

For each of two covariates, we randomly and independently draw 100 control units from a Uniform(0,5) and 100 treated units from Uniform(2,6). This leaves the overlapping [1, 5] x [1, 5] square as a completely randomized experiment and observations falling outside adding imbalance. We generate the outcome as $Y = 2T_i + X_{i1} + X_{i2} + \varepsilon$, where $\varepsilon \sim N(0,1)$

However, there is no clear description of how T is defined here. $X_{i,k} ~ \forall ~k$ is, as far as the paper goes, a vector of length 100, while T is also a (naturally) vector.

My first interpretation of this, was of the following,
# First Attempt

dgp <- function(n = 100) 
{
  
  
  # Generate Covariates
  x1 <- runif(n, min = 0, max = 5)
  x2 <- runif(n, min = 1, max = 6)
  
  # Generate data with treatment indicators
  # and outcomes
  
  tibble(
    x1 = x1,
    x2 = x2
  ) %>% mutate(
    treatment = as.numeric(x1 >1 | x2 > 1),
    y         = 2*treatment + x1 + x2 + rnorm(n, 0, 1)
  )
}

However, this attempt assigns some of the obs in the untreated runif(n, 0,5) to the treatment group, which is not what the authors intends.
My second attempt, includes cbind(x1 = runif(100, 0,5), treatment = 0) for each covariate - however, then the length of the treatment-variable exceeds that of the covariates, and therefore breaks down.
How do I interpret the simulation study mentioned in the paper, and replicate it in R?
Edit: Quoted exactly from the article after @whuber's comment.

I attempted to replicate this study two ways using the data-generating process described in the text and what I hoped to demonstrate was an equivalent data-generating process. My replication is here.
Here is a snippet of the code I used to generate the data, where A is the treatment:
n <- 200
A <- rep(0:1, each = n/2)
X1 <- ifelse(A==1, runif(n, 1, 6), runif(n, 0, 5))
X2 <- ifelse(A==1, runif(n, 1, 6), runif(n, 0, 5))

Y <- 2*A + X1 + X2 + rnorm(n)

Note that in the paper, they generate the treatment groups first, and then assign covariate values to the assigned units. That is, treatment doesn't depend on the covariates. Instead, the treatment is randomized and the covariates are actually mediators.
A colleague thought this undermined their entire paper, and I disagreed, which is why I set out to formulate an equivalent data-generating process where the treatment was created based on the covariates but the treatment and covariates had the same joint distribution as in the original formulation. Here is code to do that:
G1 = sample(LETTERS[1:5], n, replace = TRUE, 
            prob = c(5/50, 4/50, 32/50, 4/50, 5/50))
X1 = runif(n, 
           sapply(G1, switch, A = 1, B = 0, C = 1, D = 5, E = 0),
           sapply(G1, switch, A = 6, B = 1, C = 5, D = 6, E = 5))
X2 = runif(n, 
           sapply(G1, switch, A = 5, B = 1, C = 1, D = 1, E = 0),
           sapply(G1, switch, A = 6, B = 5, C = 5, D = 5, E = 1))

#P is propensity score.
P = 1*(X1 > 1)*(X2 > 1) - .5*(X1 > 1)*(X1 < 5)*(X2 > 1)*(X2 < 5)

A <- rep(0, n)
A[P == 1] <- 1
A[P == .5] <- sample(c(rep(1, n/2 - sum(P==1)), 
                       rep(0, n/2 - sum(P==0))), 
                     sum(P == .5), replace = FALSE)

Basically, you create an exogenous auxiliary grouping variable G1, then generate X1 and X2 based on G1. The propensity score is defined in an unusual way; values can be either 0, .5, or 1. Normally, treatment would be generated using rbinom() for those with propensity scores other than 0 or 1, but to ensure the number of treated and control units is equal, we use sample() instead. This ensures both the propensity score is accurate and the group sizes are balanced. Those with a propensity score of .5 correspond to the "completely randomized experiment" and those with a propensity score of 0 or 1 correspond to those "adding imbalance".
You can verify that the joint distributions of the variables are the same regardless of the method, indicating that the messages of the paper apply regardless of which data-generating process is used as long as the joint distribution of the variables is the same. This is because the paper makes a statistical argument, not a causal one, and the statistics are agnostic to the causal assumptions one would typically invoke (i.e., about whether the treatment causes the covariates or the covariates cause the treatment).
This has to do with the statistical estimand, which for matching is $E[E[Y|X,A=1]] - E[E[Y|X,A=0]]$. Whether this corresponds to a causal estimand is another question. Under a confounding data-generating process (the second one I describe), the causal estimand is the de-confounded average treatment effect $E[Y^1] - E[Y^0]$. Under the mediation data-generating process (the one in the paper), the causal estimand is the direct effect of treatment not through the covariates. The authors don't clarify this because it is tangential to their argument, but for those investigating deeply, it is worth knowing what is actually going on.
