The problem likely does not have a closed form solution, but asymptotically, the probability is $Cn^{-3/2} e^{-an}(1±o(1))$ where $C$ is a constant. ($o,O,Θ,Ω,ω$ are per Big O Notation.)
$a = 4(\operatorname{ln} 7 - h)$, where 7 is the number of move choices (away from the border), 4 corresponds to $≈4n$ total moves, and $h$ is the (base $e$) entropy of the maximum entropy distribution for mean $v=3/4$ (speed), with 3 choices with value -1 (wrong direction), one choice with value 0 (stand still), and 3 choices with value 1.
This also applies for stopping time $cn$ and $k$ dimensions (assuming constant $c>k/2$) by modifying $a$ and $C$ accordingly and changing $n^{-3/2}$ to $n^{-k/2}$.
Approximate analysis
To meet each other, the flies must choose an appropriate direction at least a net $≈v$ fraction of the time. Thus, the collision probability is exponentially small. Each wasted move reduces the probability in $1+Ω(1)$ times. For a typical collision path (for the pair), the border may only be hit during the first $O(1)$ steps, and except at $O(1)$ steps before the time limit, all coordinates of the first fly will be smaller than the corresponding coordinates of the second fly.
The above maximum entropy condition leads to random paths with the desired expected displacement. The distribution (and thus $h$) can be calculated using $3p_{-1}+p_0+3p_1=1$ (probability normalization), $-3p_{-1}+3p_1 = v$ (mean), and $p_{-1} p_1 = p_0^2$ (maximum entropy; $p_d$ increases exponentially in $d$).
Random fluctuations cause a $Θ(\sqrt n)$ uncertainty of the actual endpoint, and thus we have the $n^{-3/2}$ factor to match the position in three dimensions.
Computing C
A careful analysis yields probability $C n^{-3/2} e^{-an} (1±o(1))$ for a constant $C$. Briefly, $C$ can be estimated with arbitrary precision as the product of three components:
Borderless calculation of collision probability at (exactly) the stopping time, assuming $n→∞$. Consider the first fly starting at $(0,0,0)$. Under the above maximum entropy distribution with $+x$, $+y$, $+z$ being desired directions (and each move independent of the previous ones), we get an approximately Gaussian distribution of the end points, with $Θ(n^{0.5})$ variances. The second fly is symmetric, so we can consider it to be stationary and give twice as much time to the first fly. The desired collision probability is the product of the probability density at the desired end point and $e^{an}$. (The above maximum entropy distribution overestimates the probability in $e^{an}$ times here (and more generally the difference in log probability is linear in $(x,y,z,t)$.)
Note: The starting point being $(1,1,1)$ rather than $(0,0,0)$ can be handled by multiplying the probability by $e^{3c_1}$ (for the $c_1$ below).
Correction for the border. For this, consider a random path starting at $(0,0,0,0)$ (initial location of the first fly (after translating the corner to the origin) and time 0) and ending when we are far enough from the border, or the path becomes sufficiently irrelevant. Let $(x,y,z,t)$ be the stopping coordinates. The existence of the border multiples the probability by $E(e^{c_1(x+y+z)-c_2t})$ where $c_1=-dh/dv$ and $c_2 = vc_1 - (\operatorname{ln} 7 - h)$. For the second fly, the border correction factor is the same (also multiplicative).
Correction for the possibility to collide strictly before the stopping time. For this, for each time $t$ close enough to the end, we can compute the probability of collision at time $t$ with no later collisions, and add up across $t$. If given a collision (and assuming there are no borders) the probability of no further collisions for $m$ steps is $f(m)$, this factor ends up being (multiplicative) $\sum_{m=0}^∞ e^{-c_3 m} f(m)$ for $c_3=2c_2$ (as $c_2$ (above) corresponds to saving one time step for just one fly).
Other start locations and borders: Many other starting locations for the flies can be handled similarly but without the symmetry between $x$, $y$, and $z$. However, because in your description, the probability of not moving is $1/(k+1)$ rather than $1-k/7$, where $k$ is the number of neighbors, it is sometimes faster to move along the borders. For example, if the flies start at $(εn,εn,εn)$ and $(n,n,εn)$, then conditional on meeting in time $2n$ (or $10n$), typical journeys spend most of the time with $z=Θ(1)$ (and correspondingly for the probability, we get a $n^{-1}$ instead of $n^{-3/2}$ factor, and a complicated calculation for the exponent base). This effect dissipates as $v→0$; a typical path near the border would spend only $Θ(v^2)$ time fraction exactly at the border (assuming $d^{-1/4}≪v≪1$ where $d$ is the source-destination distance).
Low speed: If $\sqrt{t} ≪ d ≪ t$, the collision probability is $d^{-O(1)} e^{-\frac{7}{8} \frac{d^2}{t}(1±o(1))}$, where $d$ is the Euclidean distance between the flies and $t$ is the time limit. (Here, $7/8 = 1/(2σ^2)$ where $σ^2=4/7$ is the approximate unit-time variance of the Euclidean distance for faraway flies.)
Unlimited Time
A natural variation of the problem is to allow unlimited time, and ask for the expected time of the first collision. Let $k$ be the number of dimensions. The expected collision time is:
- $(1±o(1)) \, C_1 n^2$ for $k=1$
- $(1±o(1)) \, C_2 n^2 \operatorname{ln} n $ for $k=2$
- $(1±o(1)) \, C_k n^k$ for $k≥3$
For $k≥3$, this holds (with unchanged $C_k$) for arbitrary starting points of the flies at initial distance $ω(1)$, and similarly for $k=2$ and initial distance $n^{1-o(1)}$. For $k≥2$ and $t=ω(n^2)$, the probability distribution of the time $t$ of the first collision is approximately exponential.
The mixing time is $Θ(n^2)$, so at larger time scales the expected number of collisions per unit time is $≈n^{-k}$ (independently of the other, each fly is at an approximately random location). To approximate the expected time to the first collision ($k>1$), one needs to multiply $n^k$ by the expected number of collisions until the distributions become well-mixed again. For $k≥3$, $C_k=1/(1-p_k)$, where $p_k$ is the probability that the random walk of the displacement between the flies (ignoring the border) will return to the starting point.
For $k=2$, each time step of that walk has variance 4/5 in both $x$ and $y$ directions, with 0 covariance. A random collision after well-mixing of locations is likely far enough from the border that we can ignore the border for $m≪n^2$ steps. Now, after $m≪n^2$ steps past a collision, we can approximate the displacement as a 2-dimensional Gaussian, leading to density at origin (collision probability) $(2π \frac{4}{5} m)^{-1}$, which integrates (up to $m=Θ(n^2)$) to $\frac{5}{4π} \operatorname{ln} n ± O(1)$, so $C_2 = \frac{5}{4π}$.
