A few more interesting plots, including a spectacular one that I call The Eye of The Zeta Function (the last one on this post).
The series for $\phi_1$ and $\phi_2$ (see Appendix in my original question) converge very slowly and in a chaotic way depending on $\sigma$ and $t$. A better solution is to use the Riemann-Siegel formula.
Slow, chaotic convergence of the series
If for some large $n$ in the series, the quantity $t\log(n+1) - t\log n \approx t/n$ is close to an odd multiple of $\pi$, then around that $n$, a lot of terms in the series will have the same sign and similar value (as opposed to the regular alternating behavior), and if $n$ is not large enough, a sum that seems to have converged before $n$ will suddenly experience a huge shift. This is illustrated below.
In the above figure, the sum for $\phi_2$ with $\sigma=0.75$ and $t=18265.2$ after seemingly converging to around $-0.135$ using the first 5000 terms, experiences a massive drop between $n=5760$ and $5860$ before stabilizing at around $-0.292$ (the correct value confirmed by Mathematica). The next chart shows the cumulative sum of the term signs, for consecutive terms of the series up to $n=50000$. A huge drop or increase means many successive terms having the same sign around the $n$ in question (the X-axis represents $n$), explaining the mechanism at play in the first chart.
Note that when $n$ is large enough, the series (at least in this example, around $n=20,000$) eventually becomes almost perfectly alternating: the blue curve in the previous plot becomes an almost flat horizontal line.
The eye of the Zeta function
Below is the scaled orbit of the Zeta function for $\sigma=0.75$ and $0<t<3000$, featuring $300,000$ points $(X(t),Y(t))$, with $X(t)=\phi_1(\sigma,t)$, $Y(t)=\phi_2(\sigma,t)$, and using increments of $0.01$ for $t$. What you see below is just one slice corresponding to $\sigma=0.75$. The whole orbit (if you allow $\sigma$ to vary in $]\frac{1}{2},1[$) would densely cover a 3-D volume in space.
There is a hole (the "eye") around $(0,0)$ suggesting that $\zeta(\sigma+it)$ has no root if $\sigma=0.75$. It looks like there is an absolute positive constant $\epsilon_\sigma$ (not depending on $t$) such that $X^2(t)+Y^2(t) \geq \epsilon_\sigma$ for all $t$'s if $\sigma \neq \frac{1}{2}$. The million-dollar question is this: Do we have $\epsilon_\sigma \neq 0$? Of course, for the limiting case $\sigma=\frac{1}{2}$, we have $\epsilon_\sigma = 0$ since $\zeta$ has at least one root (indeed infinitely many) on the critical line.
Getting this plot right is not easy, since the series for $\phi_1,\phi_2$ converge in a chaotic way and require convergence boosting techniques (see here). Also using increments of $0.01$ is OK for small values of $t$, but as $t$ becomes larger, the oscillations in $\phi_1,\phi_2$ get exacerbated, requiring smaller and smaller increments to get a good coverage of the orbit, and to avoid missing a potential root.
My next step is to randomize $\phi_1,\phi_2$, possibly replacing $\cos(\gamma + t\log n)$ by $Z_n \cdot \cos(\gamma + t\log n)$, where $Z_1,Z_2\dots$ are i.i.d. Bernoulli-like random variables. I've already spent many hours on that, to no avail so far. Along similar lines, see this paper discussing the complete Riemann Zeta distribution.
There is a zero-free entire function bounded in every left half-plane, and such that $f-1$ is in $H^2$ in every left half-plane.
Let $\gamma$ be the boundary of the region $$D=\left\{ x+iy: |y|<2\pi/3, x>0\right\}
.$$ Consider the function
$$g(z)=\int_\gamma \frac{\exp e^\zeta}{\zeta-z}d\zeta,\quad z\in {\mathbf{C}}\backslash D.$$ The integral evidently converges and $g(z)=O(1/z)$ in ${\mathbf{C}}\backslash D.$ Now, $g$ has an analytic continuation to an entire function: deforming the contour to $\partial D_t$,
where $D_t=\{ x+iy:|y|<2\pi/3, x>t\},\; t>0$ does not change $g$ in $D$,
and shows that $g$ has an analytic continuation to ${\mathbf{C}}\backslash D_t$, and this is for every $t>0$, so $g$ is entire. Now $f(z)=e^{g(z)}$ is the desired function. If you want upper half-planes take $f(iz)$.
Remark. You can improve the estimate $g(z)=O(1/z)$.
Evidently, $g$ has infinitely many zeros $z_1,z_2,\ldots$. Then $g_k(z)=g(z)/((z-z_1)\ldots(z-z_k))$
satisfies $g_k(z)=O(z^{-k-1})$ as $z\to\infty$ outside $D_t$.
Remark 2. This construction is standard in the theory of entire functions, see, for example, Entire function bounded at every line
Sometimes this $g$ is called the Mittag-Leffler function.
Best Answer
Let $$f(z)=\frac{e^{iz}-1}{iz}.$$ This function is in the Hardy class for any upper half-plane, and has these properties: $f(0)=1,$ $f(2\pi n)=0$, $$|f(z)|\leq C\frac{e^y+1}{|z|+1},$$ (this evidently holds for large and small $|z|$, therefore there is a constant $C$ so that this holds everywhere).
Since the $L^2$ norm has the property $\| f(kx)\|^2=\| f(x)\|^2/k,$ and the norm is invariant under real translations, the following function $g$ belongs to the Hardy spaces for upper half-planes: $$g(z)=\sum_{n=0}^\infty f(n^2(z-a_n)),$$ where $a_n\in\{2\pi k: k\in {\mathbf{Z}}\}$. Moreover, it is entire, if the sequence $a_n$ grows fast enough. (It follows from the estimate for $f$ that the series is uniformly convergent on any compact subset of the plane). Now $g(a_n)=1$, so $h(z)=(g(z-i)-1)/(z-i)$, where has all properties that you required.
Some extra labor is needed to show that $h$ can be made of finite order.