The hamster from $(i,n)$ goes to south to $(i,i)$ for sure, otherwise their paths intersect. This already has probability $(1-p)^{{n\choose 2}}$.

After that let the last point of the $i$-th hamster in the horizontal line $(\cdot,j)$ (where $j=1,2,\ldots,i$) have coordinates $(j,i+f(i,i+1-j))$. Then the function $f(x,y)$ increases both in $y$ (since they go only to south and to east) and in $x$ (the paths of the $i$-th and $(i-1)$-th hamsters do not intersect iff $i+f(i,i+1-j)>(i-1)+f(i-1,(i-1)+1-(j-1))$, equivalently, $f(i,i+1-j)\geqslant f(i-1,i+1-j)$, for all $j=2,\ldots,i$). And viceversa, any such increasing function $f$ corresponds to a collection of disjoint paths. For given $f$, the hamsters make $\sum_{i=1}^n f(i,i)$ moves to east and $1+2+\ldots+n=n+{n\choose 2}$ moves to south, thus the probability of this equals $p^{\sum f(i,i)}(1-p)^{n+{n\choose 2}}$. We want to prove that the sum over all $f$ equals $(1+p)^{-{n\choose 2}}$. This is equivalent to the formula $$ \sum_f p^{\sum f(i,i)}=(1-p)^{-n}(1-p^2)^{-{n\choose 2}}. \tag{$\clubsuit$} $$

The below proof of $(\clubsuit)$ is due to Fedor Ushakov.

Note that RHS of $(\clubsuit)$ is the generating function of the sum of elements of symmetric $n\times n$ matrices $A$ with non-negative integer elements: $$ (1-p)^{-n}(1-p^2)^{-{n\choose 2}}=\sum_{A=(a_{ij})}p^{\sum a_{ij}} $$ (each diagonal element corresponds to $(1-p)^{-1}$, each pair of symmetric off-diagonal elements to $(1-p^2)^{-1}$). The bijection from matrices $A$ to functions $f$ is based on RSK. The latter makes from $A$ a semistandard Young tableaux $\lambda$ with $\sum a_{ij}$ boxes and positive integer entries not exceeding $n$.

Now let $f(x,y)$ for $n\geqslant x\geqslant y\geqslant 1$be the number of elements not exceeding $n-x+y$ in the row number $n+1-x$ (where the first row is the longest, and the elements in every column strictly increase). This $f$ is non-strictly increasing both in $x$ and $y$, defines the semistandard Young tableaux uniquely, and the sum of $f(i,i)$ is the number of boxes.

The non-trivial part of the following argument is due to Fedor Ushakov (for that he gets A at my symmetric functions course).

The hamster from $(i,n)$ goes to south to $(i,i)$ for sure, otherwise their paths intersect. This already has probability $(1-p)^{{n\choose 2}}$.

After that we make a semistandard Young tableaux $T$ with at most $n$ rows and entries between 1 and $n$ (denote by $\lambda(T)$ the shape of $T$) from our collection of paths as follows: the hamster who is now at $(i,i)$ is responsible for the $(n+1-i)$-th row (which may be empty): when he makes a move to east being in the horizontal line $(\cdot,j)$, he adds a box with the number $n+1-j$ to this row, see an example.

The hamsters make $|\lambda(T)|$ moves to east and $1+2+\ldots+n=n+{n\choose 2}$ moves to south, thus the probability of this equals $p^{|\lambda(T)|}(1-p)^{n+{n\choose 2}}$. We want to prove that the sum over all Young tableaux $T$ with at most $n$ rows and entries not exceeding $n$ equals $(1+p)^{-{n\choose 2}}$. This is equivalent to the formula $$ \sum_T p^{|\lambda(T)|}=(1-p)^{-n}(1-p^2)^{-{n\choose 2}}. \tag{$\clubsuit$} $$

Note that RHS of $(\clubsuit)$ is the generating function of the sum of elements of symmetric $n\times n$ matrices $A$ with non-negative integer elements: $$ (1-p)^{-n}(1-p^2)^{-{n\choose 2}}=\sum_{A=(a_{ij})}p^{\sum a_{ij}} $$ (each diagonal element corresponds to $(1-p)^{-1}$, each pair of symmetric off-diagonal elements to $(1-p^2)^{-1}$). The bijection from matrices $A$ to semistandard Young tableaux $T$ is provided by RSK.

The hamster from $(i,n)$ goes to south for $(i,i)$ for sure, otherwise their paths intersect. This already has probability $(1-p)^{{n\choose 2}}$.

After that let the last point of the $i$-th hamster in the horizontal line $(\cdot,j)$ (where $j=1,2,\ldots,i$) have coordinates $(j,i+f(i,i+1-j))$. Then the function $f(x,y)$ increases both in $y$ (since they go only to south and to east) and in $x$ (the paths of the $i$-th and $(i-1)$-th hamsters do not intersect iff $i+f(i,i+1-j)>(i-1)+f(i-1,(i-1)+1-(j-1))$, equivalently, $f(i,i+1-j)\geqslant f(i-1,i+1-j)$, for all $j=2,\ldots,i$). And viceversa, any such increasing function $f$ corresponds to a collection of disjoint paths. For given $f$, the hamsters make $\sum_{i=1}^n f(i,i)$ moves to east and $1+2+\ldots+n=n+{n\choose 2}$ moves to south, thus the probability of this equals $p^{\sum f(i,i)}(1-p)^{n+{n\choose 2}}$. We want to prove that the sum over all $f$ equals $(1+p)^{-{n\choose 2}}$. This is equivalent to the formula $$ \sum_f p^{\sum f(i,i)}=(1-p)^{-n}(1-p^2)^{-{n\choose 2}}. \tag{$\clubsuit$} $$

The below proof of $(\clubsuit)$ is due to Fedor Ushakov.

Note that RHS of $(\clubsuit)$ is the generating function of the sum of elements of symmetric $n\times n$ matrices $A$ with non-negative integer elements: $$ (1-p)^{-n}(1-p^2)^{-{n\choose 2}}=\sum_{A=(a_{ij})}p^{\sum a_{ij}} $$ (each diagonal element corresponds to $(1-p)^{-1}$, each pair of symmetric off-diagonal elements to $(1-p^2)^{-1}$). The bijection from matrices $A$ to functions $f$ is based on RSK. The latter makes from $A$ a semistandard Young tableaux $\lambda$ with $\sum a_{ij}$ boxes and positive integer entries not exceeding $n$.

Now let $f(x,y)$ for $n\geqslant x\geqslant y\geqslant 1$be the number of elements not exceeding $n-x+y$ in the row number $n+1-x$ (where the first row is the longest, and the elements in every column strictly increase). This $f$ is non-strictly increasing both in $x$ and $y$, defines the semistandard Young tableaux uniquely, and the sum of $f(i,i)$ is the number of boxes.

The hamster from $(i,n)$ goes to south to $(i,i)$ for sure, otherwise their paths intersect. This already has probability $(1-p)^{{n\choose 2}}$.

After that let the last point of the $i$-th hamster in the horizontal line $(\cdot,j)$ (where $j=1,2,\ldots,i$) have coordinates $(j,i+f(i,i+1-j))$. Then the function $f(x,y)$ increases both in $y$ (since they go only to south and to east) and in $x$ (the paths of the $i$-th and $(i-1)$-th hamsters do not intersect iff $i+f(i,i+1-j)>(i-1)+f(i-1,(i-1)+1-(j-1))$, equivalently, $f(i,i+1-j)\geqslant f(i-1,i+1-j)$, for all $j=2,\ldots,i$). And viceversa, any such increasing function $f$ corresponds to a collection of disjoint paths. For given $f$, the hamsters make $\sum_{i=1}^n f(i,i)$ moves to east and $1+2+\ldots+n=n+{n\choose 2}$ moves to south, thus the probability of this equals $p^{\sum f(i,i)}(1-p)^{n+{n\choose 2}}$. We want to prove that the sum over all $f$ equals $(1+p)^{-{n\choose 2}}$. This is equivalent to the formula $$ \sum_f p^{\sum f(i,i)}=(1-p)^{-n}(1-p^2)^{-{n\choose 2}}. \tag{$\clubsuit$} $$

The below proof of $(\clubsuit)$ is due to Fedor Ushakov.

Note that RHS of $(\clubsuit)$ is the generating function of the sum of elements of symmetric $n\times n$ matrices $A$ with non-negative integer elements: $$ (1-p)^{-n}(1-p^2)^{-{n\choose 2}}=\sum_{A=(a_{ij})}p^{\sum a_{ij}} $$ (each diagonal element corresponds to $(1-p)^{-1}$, each pair of symmetric off-diagonal elements to $(1-p^2)^{-1}$). The bijection from matrices $A$ to functions $f$ is based on RSK. The latter makes from $A$ a semistandard Young tableaux $\lambda$ with $\sum a_{ij}$ boxes and positive integer entries not exceeding $n$.

Now let $f(x,y)$ for $n\geqslant x\geqslant y\geqslant 1$be the number of elements not exceeding $n-x+y$ in the row number $n+1-x$ (where the first row is the longest, and the elements in every column strictly increase). This $f$ is non-strictly increasing both in $x$ and $y$, defines the semistandard Young tableaux uniquely, and the sum of $f(i,i)$ is the number of boxes.