As metioned in @assumednormal's comment, standard least squares regression requires linearity in the coefficients. In other words we need to be able to write the outcome as:
$Y = X\beta + \epsilon$
The matrix of independent variables $X$ however, can be non-linear. For example, the following belongs to this model
$Y_i =\alpha_y + f_1({X_1}_i)\beta_1 + f_2({X_2}_i)\beta_2 + {X_1}_i{X_2}_i\beta_3 + \epsilon_i$
where $f_1({X_1})$ and $f_2({X_2})$ are non-linear functions of $X_1$ and $X_2$ and clearly the interaction term $X_1X_2$ is non-linear.
Non-linear Regression on the otherhand is typically written as:
$y_i = x_i(\beta) + \epsilon_i$, where $x_i(\beta)$ is non-linear in the coefficients themselves. Your example $y_i = \beta_1 e^{(X_i\beta_2)} +\epsilon_i$ would be an example of this form. However under the assumption that $E[y_i|X_i] >0$ (which implies that $\beta_1 >0$) we could model $\log(E[y_i|X_i])$ as a linear model.
\begin{align}
E[y_i|X_i] &= \beta_1 e^{(X_i\beta_2)}\\
\log(E[y_i|X_i]) &= \log(\beta_1) + X_i\beta_2\\
\end{align}
Where $\log(\beta_1) \in \mathcal{R}$ is just some number just like a normal intercept term. So this is an example of non-linear least squares and could be treated and estimated as such, but under some conditions and goals can still be linearized and estimated. Chapter 6 of Econometric Theory and Methods (Davidson and McKinnon) discusses this and says more generally that many non-linear models can be reformulated into the form of a linear regression, but sometimes with non-linear restrictions on the coefficient themselves (if there are non-linear restrictions on $\beta$ we can't use the standard formula to estimate properly). In other words there can be slightly different definitions of what constitutes a linear model. Often implicitly people mean that with some transformation or reparametrization it can be estimated with ordinary least squares.
Which leads us to your leading example $Y = \beta_0 + X_1\beta_1^2 + X_2\beta_2 + \epsilon$, which is an interesting gray area in my mind. We can in fact reparametrize it to be a linear model, but we have to be careful with the parameter space and it cannot be estimated with ordinary least squares.
In this case, the only problem is that $\beta_1^2\geq 0$, which implies a restriction on the coefficient. We could reparametrize the model, with a new coefficient say $\beta_1^{\star} = \beta_1^2$ and write the model as:
\begin{align}
y = \beta_0 + \beta_1^{\star}1\{\beta_1^{\star} \geq 0\}X_1 + X_2\beta_2 +\epsilon
\end{align}
This is linear function in parameters over the parameter space $(\beta_0,\beta_1,\beta_2) \in (\mathcal{R},\mathcal{R}^{+},\mathcal{R})$. Ordinary least squares cannot guarantee a solution in this parameter space, but this is a special case of non-negative least squares, where we are solving the minimization problem
\begin{align}
\underset{(\beta_0,\beta_1,\beta_2):\beta_1^{\star} \geq 0}{\operatorname{argmin}} ||Y - \beta_0 - \beta_1^{\star}X_1 - \beta_2X_2||^2
\end{align}
This is a convex minimization problem and the solutions are well known (see for example this paper about it's application in high-dimensions https://arxiv.org/pdf/1205.0953.pdf)
So no, not technically OLS, but linear over a restricted parameter space and the restrictions are linear. But this is not usually what people mean when they say linearizable.
