As listed your reasoning is indeed circular: you start by assuming $h \approx \sqrt{\epsilon}$, then use the results in point 3 to infer that we should take $h \approx \sqrt{\epsilon}$. So that's not convincing. So let's redo the analysis more carefully, then we'll look at the meaning of the constant.
A more careful analysis
As in your analysis, we'll assume $f'(x) \approx 1$ and $f''(x) \approx 1$.
Your analysis of truncation error is fine; the truncation error is indeed $\frac{h}{2} f''(x) + O(h^2)$, which given your assumptions is $\frac{h}{2} + O(h^2)$. In our analysis $h$ will be small enough that $O(h^2)$is negligible, so we'll call the truncation error $\frac{h}{2}$ for simplicity.
What about rounding error? The rounding error of representing $f(h)$ in floating point is $\epsilon$ in relative error, where $\epsilon$ is the machine epsilon; in absolute error, that is $\epsilon f(h)$. Since you assumed $f(h) \approx 1$, that is $\approx \epsilon$. Similarly, the rounding error of representing $f(x+h)$ is $\epsilon f(x+h)$ in absolute error; since $f(x+h) \approx f(x) \approx 1$, that is $\approx \epsilon$ as well. Now the absolute error of $f(x+h)-f(x)$ is at most about $2\epsilon + \epsilon (f(x+h)-f(x))$; the $2\epsilon$ comes from the error in each of the two subterms (and the triangle inequality), and the $\epsilon (f(x+h)-f(x))$ comes from roundoff error in representing the result. Since $f(x+h)-f(x) \approx h f'(x) \approx h$, that means the absolute error of $f(x+h)-f(x)$ is at most about $2\epsilon + \epsilon h$.
Then we divide by $h$. The absolute error becomes $2\epsilon/h + \epsilon$, plus a term due to representing the result in floating point; that extra term is $\epsilon ((f(x+h)-f(x)/h)$, and since $(f(x+h)-f(x)/h \approx f'(x) \approx 1$, we see that the extra term is $\approx \epsilon$. Therefore the absolute error in computing $(f(x+h)-f(x)/h$ is about
$2\epsilon/h + \epsilon + \epsilon$. So, the total roundoff error is about $2\epsilon/h + 2\epsilon$, in absolute error.
The total error is at most the sum of the truncation error and the roundoff error, or about
$$\frac{h}{2} + {2\epsilon \over h} + 2\epsilon.$$
We want to know what value of $h$ minimizes this expression. We can figure that out by letting $g(h) = \frac{h}{2} + 2\epsilon/h + 2\epsilon$ and minimizing $g(h)$; the minima must be at a point where the derivative is zero. So, differentiate:
$$g'(h) = \frac{1}{2} - {2\epsilon \over h^2}.$$
Setting $g'(h)=0$ and solving for $h$ yields $h^2 =4 \epsilon$, i.e., the function $g(h)$ is minimized at $h=2\sqrt{\epsilon}$; at that point, we have a total (absolute) error of at most $g(2\sqrt{\epsilon}) \approx 2 \sqrt{\epsilon} + 2 \epsilon \approx 2\sqrt{\epsilon}$.
The summary: assuming $f(x)\approx 1$ and $f'(x) \approx 1$, the total error is minimized when you choose $h$ to be about $h\approx 2\sqrt{\epsilon}$. If you choose $h$ too much larger than that, the truncation error increases and starts dominating. If you choose $h$ too much smaller than that, the roundoff error increases and starts dominating.
So, if you know $f(x)\approx 1$ and $f'(x) \approx 1$, then it's reasonable to choose $h\approx 2\sqrt{\epsilon}$.
The constant
Now, you might be wondering a bunch of things. Why the constant 2? Where did that come from? How critical is it for the constant to be 2? What if I use $h \approx \sqrt{\epsilon}$ instead of $h \approx \sqrt{2\epsilon}$?
Well, if you redo the calculation with a different value of $h$, you'll see that the total error is not too sensitive to the constant 2. For instance, if you choose $h=\epsilon$, then the total error will be at most about $g(\epsilon) \approx 2.5 \sqrt{\epsilon} + 2 \epsilon \approx 2.5 \sqrt{\epsilon}$, which is not that different. So don't worry too much about the exact value of the constant.
Where did the constant 2 came from? Well, partly it is coming from the fact that we are doing multiple operations, and you can get a little bit of roundoff error from each operation, so the error can accumulate. And, partly it is coming from the assumptions that $f'(x) \approx 1$ and $f''(x) \approx 1$. Of course, those assumptions are pretty sketchy. In practice, we don't know what $f'(x)$ is before doing the computation (that's the whole point of this computation; we are trying to estimate/compute $f'(x)$), and we don't know what $f''(x)$ is. So, that assumption also looks circular.
But here's why in this case it might be OK. You could redo the calculation with different assumptions, e.g., $f'(x) \approx 2$ and $f''(x) \approx 3$. If you do that, you'll find a slightly different optimal value of $h$; instead of $2 \sqrt{\epsilon}$, you'll have a slightly different constant replacing the 2. So we can see that the 2 comes from the assumptions about the magnitude of $f'(x)$ and $f''(x)$. But, crucially, if $f'(x)$ and $f''(x)$ aren't too large or too small, then the resulting constant in the optimal choice of $h$ won't be too different. It might be $5 \sqrt{\epsilon}$ or $0.5 \sqrt{\epsilon}$, but that's not so bad.
And, in practice, we're often dealing with smooth functions where the derivatives don't get too large, so while we don't know exactly what $f'(x)$ or $f''(x)$ are, we might be able to assume they won't be too large. As a result, taking $h \approx c \sqrt{\epsilon}$ for some small value of $c$ will be optimal. We don't know what the optimal value of $c$ is exactly, but we can just heuristically pick a small value for $c$, like $c=1$ or $c=2$, and even if that's not exactly optimal, we'll probably do reasonably well -- it won't be too much worse.
So, the final answer is: no, the 2 is not essential. It's perfectly reasonable to use $h = \sqrt{\epsilon}$ instead of $h = 2 \sqrt{\epsilon}$ -- but you probably want to choose $h$ somewhere in that neighborhood.
