Models for testing non-linear dynamics (e.g., threshold effects), where threshold is not necessarily known Any recommendations for models (or applications in papers) to test theoretical predictions about non-continuous effects, e.g., based on a threshold, rather than marginal effects? Some examples of the sort of threshold theoretical predictions I'm interested in testing:
The simplest example would be where the threshold is known. e.g.,
\begin{equation}
Y_i=1 \ \text{if} \  X_i>0.5 \ \text{and}\  Y_i=0 \ \text{otherwise} \ \ \ \ \ \ \ \text{(1)}
\end{equation}
A more useful example would be where we theorize a threshold, but we don't know its value, or even assume that it's constant across individuals.
\begin{equation}
Y_i=1 \ \text{if} \  X_i>T_i \ \text{and}\  Y_i=0 \ \text{otherwise} \ \ \ \ \ \ \ \text{(2)}
\end{equation}
Where $T_i$ is a threshold that is allowed to differ across units indexed in $i$.
I don't just want to test whether $X$ has a positive effect on $Y$ (e.g., $Y$~$X$), but rather I want to test for evidence for the threshold prediction.

I think you're describing change point detection, and you can find an enormous amount of information on that topic once armed with the right term to search for. By the way, this kind of thresholding does not constitute "non-linear dynamics", since there is no non-linearity, and no dynamics over time.
There are lots of existing tools for doing change point detection, but here's a simple approach using a noisy version of your data:
library(tidyverse)
x = 1:100
changepoint = 50
y_mean = ifelse(x < changepoint, 0, 1)
y = rnorm(length(x), y_mean, .5)
plot(x, y)
abline(v=changepoint)


#' Log-likelihood of x under the best-fitting parameters
gaussian_loglik = function(x){
  m = mean(x)
  s = sd(x)
  dnorm(x, m, s, log = T) %>% sum()
}

evaluate_changepoint = function(proposed_changepoint){
  y_pre = y[x < proposed_changepoint]
  y_post = y[x >= proposed_changepoint]
  gaussian_loglik(y_pre) + gaussian_loglik(y_post)
}

changepoint_liks = map_dbl(x, evaluate_changepoint)
estimated_changepoint = x[which.max(changepoint_liks)]
plot(x, changepoint_liks, 'l',
     xlab = 'Value',
     ylab = 'Log lik. of changepoint')
abline(v=estimated_changepoint)


