An obvious pattern to $i\uparrow\uparrow n$ that is eluding us all? Start with $i=\sqrt{-1}$.
This will be $a_1$.
$a_2$ will be $i^i$.
$a_3$ will be $i^{i^{i}}$.
$\vdots$
etc.
In Knuth up-arrow notation:
$$a_n=i\uparrow\uparrow n$$
And, amazingly, you can evaluate $\lim_{n\to\infty}a_n=\lim_{n\to\infty}i\uparrow\uparrow n=e^{-W(-\ln(i))}\approx0.4383+0.3606i$.
You can check this, it does indeed converge to this value.
In fact, I decided to make a graph of $a_n$ to show that it converges. (y axis is imaginary part, x axis is real part.)
And, to little astonishment, I quickly noticed that there is an apparent pattern to the graph.
Commonly, we define $x\uparrow\uparrow0=1$, which I have included in the graph.
So the pattern seems very obvious.  It follows a curved path that converges onto the point that was given above.
And, if you connect the dots, starting with the first point (given on the left as the first point) and trace a nice line to the second, third, and so fourth numbers, you will find an interesting spiral.  I thought that at first, this spiral was writable as an equation, but apparently, there are a few implications.
You will notice that the blue dots are way closer to the converging point and that the red and black dots are a little closer.  So whatever equation you can come up with should account that $a_{3n}$ is closest to the number you are trying to converge to.
I want (so desperately) to see if anyone can come up with an equation that allows the computation of $a_{0.5}$ that satisfies $$i^{a_{0.5}}=a_{1.5}$$a well known identity you can find on the Wikipedia.
At first glance of the graph I went on to think that perhaps, just perhaps, I (or you) could find a formula that allows us to define $i\uparrow\uparrow 0.5$.
If you are familiar with De'Moivres formula, it is a formula that allows us to perform compute
$$\sqrt{i}$$
with relative ease.  It was derived when De'Moivre noticed an interesting pattern to $(a+bi)^n$.  He proceeded to write his formula concerning the distance from zero and the angle from the positive real axis.
So I must tell you that I wish for the same to occur with $i\uparrow\uparrow n$.  Perhaps the answer lies in using a different coordinate system.  Perhaps the answer lies in calculating the distance one of the points on one of the lines (black, red, or blue) is from the converging spot and the adding in the angle at which the next point changes.
My progress on determining such a formula has gone nowhere.  The most I can say is that $a_n$ is probably not chaotic and does indeed converge in a way that is most certainly not random.
 A: Consider the functions $$f_n(z) = \frac{\sin \pi z}{z - n}$$ for $n\in\Bbb Z$. Note that $f_n(k) = 0$ if $n \ne k$ and $f_n(n) = 1$.
Define $$f(z) := \sum_{n=-\infty}^{\infty} a_{|n|}f_n(z) = \sin \pi z \sum_{n=-\infty}^{\infty} \frac{a_{|n|}}{z - n} = \sin \pi z\left(\frac{a_0} z + 2z\sum_{n=1}^{\infty} \frac{a_n}{z^2 - n^2}\right)$$
(for the appropriate definition of the double-infinite sum). Since the $a_n$ converge, they are bounded, so the last sum converges for all non-integer $z$, with simple poles that the integers, which are cancelled out by the multiplication of $\sin \pi z$, leaving $f(n) = a_{|n|}$ for all $n \in \Bbb Z$.
You can define $a_z = f(z)$ for arbitrary $z$.
However, this is just one possible function. For any entire function $h$, the function $g(z) = f(z) + h(z)\sin \pi z$ also satisfies $g(n) = a_{|n|}$ for $n\in \Bbb Z$, but generally has different values from $f$ off of $\Bbb Z$.
There are other choices as well. Maybe you want to define $a_n$ for negatives by the inverse recursion formula $a_n = -\frac {2\log a_{n+1}}{\pi}$. A similar construction can be done, though one has to take more care to ensure convergence (the reason I went with an even function for $f$).
So without some extra criteria, there is not a well-defined answer. Your formula $f(z + 1) = i^{f(z)}$ may resolve the issue, but I haven't worked that out yet.
A: The number $ ^\infty i \approx 0.4383+0.3606i$ where $ i^{^\infty i} = {^\infty i}$ is a hyperbolic fixed point. Let $\epsilon$ be a very small number such that $\epsilon^2 \approx 0$. Then 
$\large ^\infty i + \epsilon \rightarrow  i^{^\infty i+ \epsilon} = {^\infty i} \times i^\epsilon = {^\infty i} \times e^{Ln(i) \epsilon} = {^\infty i}(1 + Ln(i)  \epsilon) = {^\infty i} + Ln(^\infty i)  \epsilon$  
$\large ^\infty i + \epsilon \rightarrow  {^\infty i} + Ln(^\infty i)  \epsilon$ 
The number 
$ Ln(^\infty i) \approx Ln(0.4383+0.3606i) \approx -0.566386 +0.688444 i$ is called the multiplier in dynamics. 
In polar coordinates $ Ln(^\infty i),  r\approx 0.891487 ; \theta \approx 129.444°$. This describes the triangles in the plot and why they exist. The logarithm of the multiplier is called the Lyapunov exponent which is $-0.114865 +2.25923 i$ Because the Lyapunov exponent's x value is negative the system has a hyperbolic attractor at the fixed point.
A: You can find a non-trivial interpolation for the fractional iteration-height when you write down the consecutive iterates in log-polar-form (with center at the fixpoint). The nearer you come to the final fixpoint the log of the distance as well as the angle come nearer and nearer to a linear relation with the index and this suggests an obvious method of interpolation for fractional iteration-heights.             
I found it interesting that that type of interpolation agrees well with the solution, which you would find via the method invented by E. Schröder in the late 19 century. Although this log-polar/Schröder-interpolation gives a straightforward solution, there seems to be a better one (better in what sense? - too difficult to discuss it here) in the spirit of Kneser's analytic solution for the fractional interpolation of the $\exp()$-function (implementations available by tetrationforum). The latter can seemingly be approximated by a simple (but computationally much involved) procedure involving matrix-diagonalization and computing fractional powers of that matrix.              
You can find an introductory comparision of that mentioned methods (however for a different base for the exponentiation) in this small essay of mine

Here is an image for an interpolation to fractional heights starting at $z_0=1$ going to $z_1=î,z_2 \approx 0.2078,...,z_\infty \approx 0.438+0.361 î$ using the Schröder-mechanism. For instance for the half-iterate we find by this method $z_{0.5} \approx 1.1667+0.734 î$. The grey dotted line indicates the integer-iterates (should be the same as that of @GEdgar)                
 

I got Sh. Levenstein's Pari/GP-program "fatou.gp" (from the tetrationforum) for the (extended) Kneser-method working. Here is a comparision of the orbits produced by the two methods. For instance, the half-iterates differ even visually:       


