A partial answer follows, contributions are welcome, and the question remains open.
Some of the questions that posed substantial roadblocks toward an understanding of what $R_{adj}^2$ is, are dealt with here in order to address the question "Would the real adjusted R-squared formula please step forward." What it is not is a proof of anything, it is merely an attempt to find what a solution to the adjustment problem is, and is not.
@EdM refers to this answer, which refers to Karch J (2019), which refers to the adjusted R-squared methods of Olkin and Pratt (1958), and in which it states that the method of Olkin and Pratt has not been implemented in any software because no one likes hypergeometric $\,_2F_1$ functions. Olkin and Pratt, as their last published function write that
$$R_{adj}^{2}=1-\frac{n-2 }{n-v}\left(1-R^2\right) \, _2F_1\left(1,1,\frac{1}{2} (n-v+2);1-R^2\right)\;,$$
but write, as identified by Karch, $\,_2F_1(\alpha,\beta,\gamma;x)$, as $F(\alpha,\beta;\gamma;x)$. Fortunately, Olkin and Pratt list the appropriate infinite sum definition for an $\,_2F_1$ function, $\sum _{k=0}^{\infty } \frac{x^k (\alpha )_k (\beta )_k}{k! (\gamma )_k}$. About this formula they write, "We cannot hope for a non-negative unbiased estimator, since there is no region in the sample space having zero probability for $\rho^2 = 0$ and positive probability for $\rho^2 > 0$. For the same reason there can be no positive unbiased estimator of $\rho$ either." What the Olkin and Pratt (O-P) solution does is find a minimum error estimator of $R_{\text{adj}}^2$. This results in a family of curves that looks quite similar to the real and imaginary plot to follow, below. However, the case in which $v=1$ has a region wherein the O-P correction factor is above the line of identity, which indeed agrees with O-P's caution to only use their formula for $v\geq2$.
The current most commonly used formula appears in the answer given by @EdM below as $$R^2_{adj} = 1 - (1 - R^2) \frac{n}{n-v}\;,$$ that is, if we count constant parameters as parameters (and not doing so would be indefensible). However, it has problems.
Suppose we examine what that formula means by taking its limits. First pro forma let us take the limit as $R\to 1$, that yields 1, and thus there is no correction, and that makes sense. However, let us now take the limit as $R\to 0$,
$$\underset{R\to 0}{\text{lim}}\left[R^2_{adj} =1-\frac{n \left(1-R^2\right)}{n-v}\right]\to R^2_{adj}=-\frac{v}{n-v}$$
That result is a negative number as $n>v$, and there are several problems with that.
(1) It means that a relationship to the square of Pearson correlation coefficient does not exist (see below), as the square of that coefficient is bounded on $(0,1)$, and the adjustment formula is bounded on $(-\frac{v}{n-v},1)$, which has no clear meaning.
(2) More generally, $R^2$ and the Pearson correlation coefficient squared should ideally have similar properties and indeed in the case of monovariate ordinary least squares linear regression (in x or y) R-squared is identical to r-squared, the square of Pearson correlation coefficient. This is not true if we assign a slope or intercept that is not the same as that implied by the Pearson correlation coefficient, or during bivariate regression.
(3) There have been attempts to justify the negative valued answers for both $R_{adj}^2$ and $R^2$ which can arise from its ANOVA definition, for example as in the following image of linear regression with intercept set at zero using Excel,
And we might reasonably ask how the $R^2$-value can be $-1.165$, and what that means. It could be argued that we are just subtracting a positive square, however, that is incorrect. Recalling the limit as $R\to 0$ above, $R^2_{adj}=-\frac{v}{n-v}$, it is clearly the square itself that is negative. The rules of mathematics thus allow one to take its square root which yields,
$$R_{adj}(R=0)=\pm\sqrt{\frac{v}{n-v}}\,i\;,$$
where $i=\sqrt{-1}$, an imaginary number. Let us plot the $R$-values, i.e., the square root of $R_{adj}^2$ from the formula above, which is then
$$R_{adj}=\pm\sqrt{1-\frac{n \left(1-R^2\right)}{n-v}}$$
to see just how bad it is as a family of curves for the positive root for $n=10$, and $v=1\text{ to }9$, as we cannot set $v=10$. The good news is that if $v=0$ the formula at least will return the identity. Plotting $R^2$ rather than $R$-values would change nothing, negative $R^2$-values are not real, that is, one cannot square a real number and obtain a negative value.
As we can see from the plot, imaginary number answers (Im) as dashed blue lines for $R_{adj}$ occur at smaller R-values when the degrees of freedom $(v)$ are fewest, and become progressively more problematic as their number increases. The real number answers (Re) are in solid blue lines. Trusting this type of adjustment appears to be questionable unless the $R_{adj}$-values are close to 1.
(4) It may be that what is being done during ANOVA is to ignore an "interaction term," a cross product, or "other" term in the ANOVA formula for $R^2$, For example, see this post Nonlinear regression SSE Loss.
Thus, we are stuck trying to squeeze blood from a turnip to make any real world sense of this. Can anything be done to breath life into this? Remember the $\pm$ sign when we took the square root? One could redefine adjusted R-values piecewise to be
$$\left\{
\begin{array}{@{}ll@{}}
\sqrt{1-\frac{n \left(1-R^2\right)}{n-v}}, & 1-\frac{n \left(1-R^2\right)}{n-v}>0 \\
-\sqrt{\frac{n \left(1-R^2\right)}{n-v}-1}, & 1-\frac{n \left(1-R^2\right)}{n-v}\leq 0 \\
\end{array}\right.
$$
This is not quite the same as adjusted R-values, because we would also have to redefine what it means to square this formula piecewise, but it would at least convert complex values to reals. In effect, this redefines when adjusted R-values are positive and when they are negative and is the result of sneaking in an absolute value to clear the complex numbers. It looks like this:
and now has the properties we might desire, except that the negative values are still inexplicable for R on $(0,1)$.
It turns out there is a way to define adjusted R-value, or at least a good start on one, although we do not follow through entirely on all its implications here, as one can only do so much at one sitting. That is, we can redefine $R_{adj}^2$ so as to remove the bias of correlation, which some authorities claim is not possible. That is, we can replace the sample means of $x$ and $y$ with their unbiased, population mean-values as follows; Adjusted r-squared is defined as the square of the population mean corrected Pearson correlation coefficient (r). That is,
$$r_{adj}^2=\left[\frac{\sum _{i=1}^n (x_i-w) (y_i-z)}{\sqrt{\sum _{i=1}^n (x_i-w)^2} \sqrt{\sum _{i=1}^n (y_i-z)^2}}\right]^2,$$
where $w$ and $z$ are the population mean values for $x$ and $y$ respectively. If we are not given the population mean values, we cannot produce a unique value for adjusted R-squared in a particular case, but we could produce a unique estimator under certain additional assumptions, for example, from simulation studies (see below), or if the population mean values are known. Using this new definition, in no case will adjusted R-squared be negative, and this is a distinguishing feature of this re-definition. For the example of three co-ordinate pairs (1,10}, {5,2} and {3,-5}, a population mean adjusted R-squared surface has the following appearance. (Note, this shows what the adjustment is for $w$ and $z$ but does not display how likely $(w,z)$ is as a bivariate model, see a bivariate Student's t-model as an example.)
The orange surface is population adjusted R-squared. The translucent blue flat surface is the ordinary Pearson correlation coefficient squared. Note that using a current formula the software I used listed the $𝑅^2$ value as $0.28$ and the adjusted $𝑅^2$ as $−0.43$, where the former is correctly the square of the Pearson correlation coefficient, and the latter the square of an imaginary number. Note that in a single individual case with only three realizations, that population mean adjusted r-squared will be larger than r-squared when the population mean values are far enough distant from the sample mean values (limited by 1 above), but that this is increasingly improbable when either the number of realizations, $(n)$, or for the mean of sample mean-values as the number of repetitions of a simulation increases for a fixed $n$ and given population mean-values. This may seem strange at first glance, but it has to be that way. It is well known that an outlier causes correlation magnitude to increase, and when an entire cloud of data are far from the population means in $x$ and $y$, those population means, which are known to be correct by assumption, then becomes an outlier with respect to the data cloud causing the adjusted R-squared value to increase. Such a situation becomes increasingly improbable when more data is collected. The formulas in current usage give no indication of such behaviour due to their one dimensionality. The problem type, however, only reduces to a one-dimensional one in the aggregate. We show this next in simulation studies. Doing simulations has the advantage of allowing us to assign population mean values a priori. Pearson correlation is:
$$
r=\frac{\sum _{i=1}^n \left(x_i-\bar{x}\right) (y_i-\bar{y})}{\sqrt{\sum _{i=1}^n \left(x_i-\bar{x}\right)^2} \sqrt{\sum _{i=1}^n (y_i-\bar{y})^2}}
$$
Notice that even if this arises from normal distributions that we do not know what the population means are, so we substituted the sample means $\bar{x}$ and $\bar{y}$. Wait one, that means that $r$ is biased and that is the bias we are trying to remove. To understand this (in case you do not) imagine you threw two darts at a dart board, in that case $(\bar{x},\bar{y})$ would be the midpoint of the line joining the two darts, that line is also the shortest distance between those points, and note that deviations are also distances. Thus, that is the absolute minimum root mean square value of all potential differences between the two dart positions and any population mean. To see that, throw a third dart and try to neatly split the positional difference between the first two darts. In fact, you would be lucky indeed to even hit a point within the minimum circle joining those first two darts. In effect, correlation using sample means typically (but not always, as above) overestimates correlation using population means, especially when the sample size is small, which is a good reason to ask how to adjust an R-value to reduce bias, and why the correction is a function of sample size. Now one way to address this problem might be to follow up with small number correction of standard deviation, see Why are we using a biased and misleading standard deviation formula for $\sigma$ of a normal distribution?, and work that through for what a correlation is in terms of small number formulas (a collection of gamma functions).
Finally the simulations. In this case, we simulate $y=2x$, where are generating data is from $x=0,1,2,\dots,n-1,n$ for $n=\{5,10,20,40,100\}$, and $v=\{1,2,3,4\}$. Then we inject noise into $y=2x$ in multiple ways using the independent standard normal distributions, $\mathcal{N}(0,1)$. Thus, we know that our population expected values are, for example, for $n=10$ with (when noiselessly progressing from 0 to 9) are $9/2$ for $x$, and $9$ for $y$. Such that, in that case, we would write an unbiased correlation as
$$
r_{adj}=\frac{\sum _{i=1}^n \left(x_i-\frac{9}{2}\right) (y_i-9)}{\sqrt{\sum _{i=1}^n \left(x_i-\frac{9}{2}\right)^2} \sqrt{\sum _{i=1}^n (y_i-9)^2}}
$$
That is, we modelled $r$ for a new $r_{adj}$ for up to 4 rv's in a linear equation and 2 rv's for $r$ itself using modelling for xy-data as follows, $\{x\text{-value}\to\text{ns}\,rv_1+x,y\text{-value}\to \text{ns}\, rv_2 +(\text{ns}\, rv_3+x) (\text{ns}\, rv_3+2)\}$, where all $rv$'s are independent standard normal random variates, ns is the "noise multiplier" and allowed to very from 1/15, to as much as needed to obtain r-values near zero (e.g., up to 13 to 1800), stepwise in a vaguely geometric progression, that defies easy description but was adjusted to make approximately similar spacing of mean r-values for up to about 30 r-values for each curve. To be clear, in case my notation is not self-evident, the random variables are either present in the noise equation above, or they are empty (0 would work as a marker for empty in this case, but actually the rv's were either used or not). So the $x$ value in the simulations sometimes were $x$ plus scaled (SND) noise, and sometimes not. Similarly, random variables were sometimes added to the intercept (of "0") the slope of 2, the x-value used for $y=f(x)$ (which is not the same as adding noise to the $x$-value independent variable, and conceptualizing that may seem paradoxical until one realizes that a injecting noise into $y=f(x)$ as $x$ plus noise, injects noise into $y$, not into $x$ itself, if we do not redefine what the $x$-values are explicitly and independently). Note well that although two points determine a line, we have assumed that noise can be independently present in each parameter or measurement of $y=m\, x+b$, for a total of 4 rv's. Now the perfectionist would point out that because we calculating a correlation that we cannot assume that the noise added in is uncorrelated between, for example $x$ and $y$. Although that is true, if we have correlated random variables for noise, each rv does not add a whole degree of freedom, thus the simulations to follow are a simplification that is not unlike the simplification of just adding up numbers of rv's to use the current formulas. Of these simulations let us examine one to see what this means.
This is a plot of r-values on the $x$-axis and $r-r_{adj}$ on the $y$-axis. The simulation dots in this case were made using $n=5$ and $v=2$ from 24 different noise multipliers each having 100000 trials and whose mean r-values ares shown as each blue dot. In the case shown, the blue line, which closely approximates the blue dots is from the formula
$r-r_{adj}\approx A0\left(z^{A1}-z^{A2}\right),$
where $z$ is a continuous version of discrete average $r$-values. Now notice how well the blue curve fits the blue dots. That is because we did curve fitting of the blue dots to obtain that fit. Note that the error from the fit equation (blue) is small. This appears to contradict O-P's statement that "We cannot hope for a non-negative unbiased estimator...."
Now a giant leap, the orange line is from estimating the parameters $A0,A1,A2$ from functions of $n$ and $v$, and that the error from the estimating equation (orange) is a maximum of about 10% of the magnitude of the $R-R_{adj}$. Explaining this procedure would take a while, because the fit function is actually scaled from an undocumented PDF that fits somewhat better than the beta distribution, and is algebraically of different form from the Kumaraswamy distribution. However, the simple fit parameter estimators were generated from regression of 24 simulation curves. This is heuristic, and more work is needed. Nevertheless, it is instructive that as an area under the curve (AUC) of r-value corrections from $R=0\text{ to }1$, the fit function, AUC*PDF, has a fairly high $r^2$ (0.9875, logarithmic fitting) to $\text{AUC}=\dfrac{0.12631}{n^{0.6351} (n-v)^{0.3582}}$, which strongly suggests that the correction $r-r_{adj}$ decreases for increasing $n$, and increases for increasing $v$. This is all over the simulation data, it is a real effect. The same implication arises from an examination of the current adjusted $R^2$ formulas, with the difference being that those current formulas are unrelated to density functions and do not span all r-values on $(0,1)$. More on this some other time, for now, the point is that one can estimate what $r-r_{adj}$ is from functions of $n$ and $v$ over the whole range of r-values at least on $r=(0,1)$, which is not what is being done now. So let us look at what is being done now to see how that stacks up.
Although the curve close to $R=1$, seems to fit, the formula overall seems unrelated to the problem, becomes complex numbered for smaller R-values, and has error so huge when it isn't complex-valued that it defies easy description as an error function. What it shows is that we would be better off redefining $R_{adj}^2$ to be $\left\{\Re\left[\sqrt{1-\frac{n \left(1-R^2\right)}{n-v}}\right]\right\}^2$ than the current formula.
There are umpteen other possible methods of creating r-value adjustments that do not have the problems listed above. The points here are that at the moment R-value formulas leave something to be desired. In sum, adjusted R-value needs work, the current formulas are range limited and have limited applicability.
A very large thank you to @whuber and @EdM for providing conceptual context for the work above.
References
Karch J (2019) Improving on adjusted R-squared. PsyArXiv, 16. Sept. 2019. (link)
Olkin, Ingram; Pratt, John W (1958) Unbiased Estimation of Certain Correlation Coefficients. Ann. Math. Statist. 29, 1, 201--211. doi:10.1214/aoms/1177706717. (link)
