Edge percolation on $\mathbb{Z}^2$: probability that two neighbouring vertices are connected?

Indeed $P(1/2)=3/4$. With probability $1/2$, the edge $e$ joining the vertices is present, and they're connected. With probability $1/2$, the edge $e$ is not present, and then either the vertices are connected by a path not including $e$, or they're in separate clusters, and in the latter case, since there's no percolation in the critical phase, there's a closed path in the dual lattice that's composed of duals of absent edges, surrounds one of the vertices and includes the edge $e'$ dual to $e$. Removing $e'$ from such a path yields a path joining the endpoints of $e'$ that doesn't include $e'$. Thus these two situations (the endpoints of $e$ joined by a path and the endpoints of $e'$ joined by a path) are dual, and since the square lattice is self-dual, they must have the same probability. As exactly one of them obtains, that probability must be $1/2$. Thus $P(1/2)=1/2+1/2\cdot1/2=3/4$.

We can extend this to uncover a hidden symmetry in $P(p)$. Let $R(p)$ denote the probability that the vertices are connected by a path not including $e$; then $P(p)=p+(1-p)R(p)$. By the above duality argument, $R(p)=(P(p)-p)/(1-p)$ exhibits the symmetry $R(p)=1-R(1-p)$, which is antisymmetry about the point $R(1/2)=1/2$. Here's a plot of $R(p)$ that confirms this symmetry:

(Right-click on the image and open in a new tab/window to see the full resolution.)

The plot was obtained using this code, which associates each edge of a finite square lattice (with periodic boundary conditions) with a threshold uniformly distributed on $[0,1]$, then adds the edges in order of increasing thresholds and updates a histogram for $P'(p)$ when adjacent sites become connected. The red, green, blue and magenta lines show calculations on lattices of dimensions $5$, $10$, $20$ and $50$, respectively. The plot seems to be nearly converged up to its graphical resolution at dimension $50$. $10000$ runs were averaged for each lattice size, and the coincidence of the results for different sizes indicates that the statistical errors are at or below the plot's graphical resolution.

We can also obtain some analytical results using inclusion–exclusion. For small $p$, short paths between the vertices dominate, so we can derive a power expansion of $R(p)$ at $p=0$. We only have to consider self-avoiding paths, since other paths always contain a shorter path as a subpath. There are two paths of length $3$ connecting the vertices, which are disjoint from each other, so we arrive at $R(p)=2p^3+O(p^4)$ without having to deal with overlap between paths. There are six paths of length $5$, and merging them with one of the paths of length $3$ adds at least one more edge, so we can even get up to $5$th order without overlaps, $R(p)=2p^3+6p^5+O(p^6)$. There are six different unions of paths of length $3$ and $5$ that have $6$ edges, and one such union of the two paths of length $3$, so up to $6$th order inclusion–exclusion yields $R(p)=2p^3+6p^5-7p^6+O(p^7)$.

After this it becomes hard to avoid errors, so we'd better get some more help from our electronic friends. Here's code that finds all self-avoiding paths between the vertices and performs the complete inclusion–exclusion for all paths and unions containing up to $14$ edges. The result is

$$ R(p)=2p^3+6p^5-7p^6+28p^7-54p^8+140p^9-297p^{10}\\+690p^{11}-1534p^{12}+3710p^{13}-9285p^{14}+O(p^{15})\;, $$

and thus

$$ P(p)=p+2p^3-2p^4+6p^5-13p^6+35p^7-82p^8+194p^9\\-437p^{10}+987p^{11}-2224p^{12}+5244p^{13}-12995p^{14}+O(p^{15})\;. $$

Each of the self-avoiding paths corresponds to a self-avoiding loop obtained by adding the edge $e$. The number of self-avoiding loops of length $l$ is conjectured to be asymptotic to $B\mu^ll^{-\beta}$ with $\mu\approx2.63816$ (see Analytic Combinatorics, p. 365). The initial coefficients above suggest that inclusion–exclusion doesn't affect the exponential growth rate $\mu$, so we can expect the power expansion to converge up to $p=1/\mu\approx0.379052$. This is confirmed by a comparison of the numerical results with the expansions of $R(p)$ up to $5$th order and up to $13$th order:

The expansion up to $5$th order (in green) is already a very good approximation up to $p\approx1/\mu$, and around and beyond $1/\mu$ the expansion up to $13$th order (in blue) is actually worse. A detailed plot of values around $p=0.25$ indicates that the higher-order terms aren't simply wrong and do improve the approximation in the interior of the interval of convergence:

To summarize, for $p\lesssim0.4$ for practical purposes $R(p)$ is well approximated by

$$ R(p)\approx2p^3+6p^5\;, $$

and thus

$$ P(p)\approx p+(1-p)(2p^3+6p^5)\;. $$

(Note that this is derived manually above, so it doesn't rely on the correctness of the linked code.)