The term "hardest submodel" is appropriate because, firstly, any path is just a collection of probability distributions that are all in $\mathcal M$, though not all of them. Another term for "path" is therefore submodel, though "path" is more specific. The synonymous terms "hardest" or "least favorable" take a little more explaining. These have to do with the fact that estimating $\psi$ is "hardest" within the path that goes in the direction of the canonical gradient, relative to all other possible paths.

"Hardest" is defined in terms of the asymptotic variance of our best possible estimator. In other words, if we were to forget about the larger model $\mathcal M$ and only consider statistical estimation within the submodel $\{\tilde p_\epsilon: \epsilon \in[0,1]\}$, what could we say about the minimum attainable asymptotic variance if we were trying to estimate $\psi(\hat P)$?

As argued below, this is a parametric model with a single parameter $\epsilon$. For the moment, let's consider the problem of estimating the maximizer of the true log-likelihood under $\hat P$. In other words, we're looking for the distribution in the path that is most like $\hat P$. Obviously that will be $\hat P$ itself (i.e. $\epsilon=0$). But if we didn't know that a-priori and we had to estimate $\epsilon$ using the data, what could we say about the asymptotic variance of our estimate $\hat\epsilon$?

Since this is a parametric model, we can use the Cramer-Rao lower bound, which says that the minimum possible asymptotic variance that a RAL estimator of a univariate parameter $\epsilon$ can attain is $E\left[\frac{d^2}{d\epsilon^2} \log(\tilde p_\epsilon)|_{\epsilon = 0}\right]^{-1}$. This is the inverse Fisher Information for the parameter $\epsilon$ at 0 (i.e. the true value). Using the definition of the path you can see that this is the same as $E\left[h^2\right]^{-1}$ (for a path that goes in some direction $h$ / has score $h$).

The reason this is useful is because now we can consider the plug-in $\psi(\hat P_{\hat\epsilon})$. What is the asymptotic variance of this plug-in? First off, we have to realize that within the path, $\psi(\tilde P_\epsilon)$ is a univariate function of $\epsilon$ alone. The delta method tells us that the asymptotic variance of a function of an estimated parameter is given by $\frac{d}{d\epsilon} \psi(\tilde P_{\epsilon})|_{\epsilon=0}$ (i.e. evaluated at the truth) times the asymptotic variance of the estimate for the parameter. Using our fundamental identity relating the pathwise derivative and a covariance between score and gradient, we have that this term is the same as $E[h\hat\phi]$ (our path starts at $\hat P$).

Putting it together, our best-case scenario is that the asymptotic variance of $\psi(\hat P_{\hat\epsilon})$ is $E[h\hat\phi]/E[h^2]$, or, geometrically in $\mathcal L_2^0$, $\langle h,\hat\phi \rangle /\|h\|^2$. From this it is immediate (e.g. from Cauchy-Schwarz) that $h = \hat\phi$ gives the maximum value. Since this is the best-case asymptotic variance, and higher is worse in terms of estimation, we have that the path in the direction of $\hat\phi$ gives us the worst possible best-case asymptotic variance among all paths.

Intuitively, what makes certain model-parameter combinations in difficult to estimate in parametric models is that there are cases where you can change the distribution just a tiny bit, but the target parameter changes a lot. That's bad because if you try to fit the data and you get some estimated distribution, you still have a lot of uncertainty about what the parameter is: slightly different data would give you a slightly different estimated distribution, but possibly a big difference in the estimated parameter. So that is the "intuitive" meaning behind optimal asymptotic variance for a parameter in a parametric model and explains why "hardness" of the estimation problem is inherently linked to the minimum attainable asymptotic variance.

In case you're not familiar with it, maximum likelihood estimation is a general-purpose technique most commonly used in parametric model spaces. The goal is to find the distribution $\hat P$ in the model that maximizes the empirical log-likelihood of the observed data. We're trying to answer the question: "of all the possible probability distributions in the model, which is the most likely to have generated the data that I observe?". Or, formally, $\hat P = \argmax_{P\in \mathcal M} p(Z_{i, \dots n})$. "Likelihood" is just a term that we use to describe the density $p$ when we're talking about it as a quantity that itself isn't fixed (in the sense that we don't know what the true density is so we have to consider multiple options). If $Z_{i,\dots n}$ are assumed independent and identically distributed, the density factors over $Z_i$ and if we take a log (the argmax is the same either way) the product becomes a sum so we get $\hat P = \argmax_{P \in \mathcal M} \frac{1}{n}\sum_i^n p(Z_i)$. We can also use the more condensed notation $\hat P = \argmax_{P \in \mathcal M} \mathbb P_n \log p$.

This strategy is most often only usable for parametric models. If we have a parametric model, the maximum likelihood problem $\hat P = \argmax_{P \in \mathcal M} \mathbb P_n \log p$ can be equivalently stated as $\hat P = \argmax_{\theta \in \mathcal \Theta} \mathbb P_n \log p_\theta$. The advantage of this is that we've reduced the problem of looking at every possible distribution in $\mathcal M$ to a problem of looking at all possible values of the parameter $\Theta$, which is a finite-dimensional vector space like $\mathbb R^p$. If the log-likelihood function is differentiable in $\theta$, this is a classical optimization problem, the kind of which is solvable using tools from calculus. In particular, if $\mathbb P_n \log p_\theta$ as a function of $\theta$ is maximized at some value $\hat\theta$, then we know that the derivative at $\hat\theta$ must be 0. For simple models (e.g. linear regression) there is often an analytical solution which gives $0 = \mathbb P_n \frac{d}{d\theta}\log p_{\hat\theta}$, but generally you can use algorithms like gradient descent to find the maximum (or minimum if you put in a minus sign). The quantity $\frac{d}{d\theta}\log p_\theta$ is called the score with respect to $\theta$ and exactly corresponds to a path in the model space where we increment $\theta$ by some small amount. You can think of it as the sensitivity of the log-likelihood to small changes in the parameters (which is different at each point $p=p_\theta$). When $\theta$ has multiple components, the score is a vector of partial derivatives $[\frac{\delta}{\delta \theta_1}\log p_\theta, \dots \frac{\delta d}{\delta \theta_p}\log p_\theta]$, each of which may be called the score of $\theta_j$. We call $\mathbb P_n \frac{d}{d\theta}\log p_{\theta}$ a set of score equations (the free variable here is $\theta$) and we say that the maximum likelihood estimator $\hat\theta$ solves these score equations because at $\theta = \hat\theta$ we know we have to have $\mathbb P_n \frac{d}{d\theta}\log p_{\theta}$.

All of this falls apart very quickly if you work in nonparametric model spaces because we can't transform the maximum likelihood problem into a finite-dimensional optimization. How do you search over an infinite dimensional space? Even if you have a sensible way of writing the densities in your model as a function of an infinite dimensional parameter, (e.g. another function, say say $\theta(\xi)$), then for each "index" $\xi$ you have a score $\frac{\delta}{\delta \theta(\xi)}\log p_\theta$, which means you need to solve an infinite number of score equations.

Since we have freedom in choosing whatever form of path is most convenient for us, you might wonder if there is some automated way in which the estimation problem might be able to dictate what path is most convenient. This leads to the idea of a universally least-favorable submodel (or universally hardest path, etc.).

Up until now, the only constraints we've put on our paths have been local constraints in the sense that we only care about how the path behaves at its very beginning: it has to start at $\hat P$ and it has to start out in the direction of $\hat \phi$. After that we don't actually care what the path does. That means that if we want to pick among these many possible paths, we have to decide on what we want the path to look like in a universal sense (i.e. at all points along the path, not just at the beginning).

Well, why not just extend our local score condition to the rest of the path? We already have that at $\epsilon=0$ we want $\frac{d}{d\epsilon} \log \tilde p_\epsilon \big|_{\epsilon=0} = \hat\phi$, which we can write in a different way as $\frac{d}{d\epsilon} \log \tilde p_\epsilon \big|_{\epsilon=0} = \phi_{\tilde P_{\epsilon=0}}$. So why impose that condition universally instead of just locally? That would give us $\frac{d}{d\epsilon} \log \tilde p_\epsilon = \phi_{\tilde P_{\epsilon}}$. Any path starting at $\hat P$ that satisfies this condition is what we call the universally least-favorable submodel.

This path is universally least-favorable because it goes in the hardest direction at each point $\tilde P_\epsilon$. To be clear, though, what that direction is can and most often does change as we traverse the path. So it's not that the path is a "straight line" that keeps following the original gradient $\phi_{\hat P}$ no matter what. Instead, this is a path that reevaluates what the hardest direction is after each infinitesimal step and then course corrects to follow it at all times.

For a given estimation problem it might take a little math to be able to express what the universally least-favorable path actually is. But once we have it our lives do get easier in one sense: if we do an iteration of the targeting step along the universally least-favorable path we get to the final targeted estimate after just one step! Why? Because the universally least-favorable path that starts at the update $\tilde P_{\hat\epsilon}$ coincides with the continuation of the universally least-favorable path from $\hat P$. So if $\tilde P_{\hat\epsilon}$ was the maximizer of the log-likelihood along the latter path, it's still the maximizer along the path that starts from there because the two paths coincide. Thus subsequent updates will not change the distribution and this proves that the targeting step converges after just one iteration. For this reason, a TMLE that uses a universally least-favorable submodel is sometimes called a one-step TMLE.

Abbreviate $U = \left[ \frac{1_1(A)}{\hat\pi_1(X)} - \frac{1_0(A)}{\hat\pi_0(X)} \right]$ so that our definition of the path is $\text{logit}(\tilde p_\epsilon) = \text{logit}(\hat p) + \epsilon U$. Taking a derivative w.r.t. $\epsilon$ on both sides and applying the chain rule gives $\left(\frac{1}{\tilde p_\epsilon(1-\tilde p_\epsilon)} \right) \frac{d}{d\epsilon} \tilde p_\epsilon = U$, or, equivalently, $\frac{d}{d\epsilon} \tilde p_\epsilon = U\tilde p_\epsilon(1-\tilde p_\epsilon)$. Again by the chain rule we know $\frac{d}{d\epsilon} \log \tilde p_\epsilon \big|_{\epsilon=0} = \left( \frac{1}{\tilde p_\epsilon}\frac{d}{d\epsilon} \tilde p_\epsilon \right) \big|_{\epsilon=0} = \frac{1}{\hat p}\frac{d}{d\epsilon} \tilde p_\epsilon\big|_{\epsilon=0}$. Combining these two facts gives $\frac{d}{d\epsilon} \log \tilde p_\epsilon\big|_{\epsilon=0} = U(1-\hat p)$.

Our likelihood is $\hat p(Y|A,X) = \hat\mu_A(X)$ if $Y=1$ and $\hat p(Y|A,X) = 1-\hat\mu_A(X)$ if $Y=0$, so a little algebra shows that $1-\hat p$ can be equivalently written as $Y - \hat\mu_A(X)$. Combining this with the above gives the desired result.

We can write the conditional likelihood for a binary variable $Y$ where $P(Y=1|A,X) = \tilde\mu_A(X)$ like $\tilde p_\epsilon = \tilde\mu_A(X)^Y(1-\tilde\mu_A(X))^{(1-Y)}$. This shows that the log-likelihood is

Moreover we have from our construction of the path that $\text{logit}(\tilde p_\epsilon(Y|A,X)) = \text{logit}(\hat p(Y|A,X)) + \epsilon U$. Plugging in $Y=1$ and the definition of $\mu_A(X)$ (tilde or hat) immediately gives that

Compare this to the definition of logistic regression (e.g. eqs. 12.7 and 12.4 here). I've included underbraced quantities to the above so you can directly map the notation.

The only difference with the usual logistic regression is that in our case, we have a fixed intercept $\beta_0=0$ and a fixed coefficient $\beta_1=1$. This is known as an "offset" in the literature about logistic regression and options to use offsets exist in practically every software package that implements logistic regression. The coefficient $\beta_2 = \epsilon$ is what we fit using the "data" $(Y, U)_{1,\dots n}$, which is precisely the solution to $\hat\epsilon =\argmax_{\epsilon} \mathbb P_n \log \tilde p_\epsilon$, as desired.

We saw that the bias-corrected and estimating equations approaches give the same estimator for the ATE problem, which we called AIPW. The form of this estimator is a sample average over terms where one of the factors is $1/\hat\pi_a(X)$. If the estimated propensity score $\hat\pi_a(X)$ is close to 0 or 1 for some subjects in the study, that factor ends up being massive and it completely dominates the sample average. Of course, in an infinite sample, that couldn't happen because some other subject would have an equally large weight in the other direction or the large weight would be balanced with many subjects with smaller weights. But we don't have infinite data. We have finite data. That makes the AIPW estimator unstable in finite samples where there are strong predictors of treatment assignment. The only way to address this is to use ad-hoc methods like weight stabilization or trimming.

On the other hand, the TMLE estimator doesn't incorporate terms of the form $1/\hat\pi_a(X)$ directly into a sample average. Instead, it uses them as a covariate in a linear adjustment of the original regression model $\hat\mu_A(X)$. That makes the procedure more stable in the face of small or large propensity scores. This comes automatically. There are also additional variations on TMLE (e.g. collaborative TMLE) that are theoretically justified and which do even more to improve the finite-sample robustness of the resulting estimator.