Start with the inclusion-exclusion argument. We compute the
cardinality of a set of arrangements where a given set of $q$ couples
or more sit together. There are $2^q$ possibilities of ordering the
$q$ couples. These are now fused and we have $2n-2q$ people left over
which are joined by $q$ couples in a circular permutation, so there
are $(2n-2q+q)!/(2n-2q+q) = (2n-q-1)!$ such arangements. This yields
by inclusion-exclusion the closed formula
$$\sum_{q=0}^n {n\choose q} (-1)^q 2^q (2n-q-1)!$$
which we divide by $(2n-1)!$ to get the probability. For the
asymptotics we use
$$(2n-q-1)! = \Gamma(2n-q) = \int_0^\infty x^{2n-q-1} e^{-x} dx.$$
This yields for the sum
$$\int_0^\infty x^{2n-1} e^{-x}
\sum_{q=0}^n {n\choose q} (-1)^q 2^q x^{-q}
dx
\\ = \int_0^\infty x^{2n-1} e^{-x}
\left(1-\frac{2}{x}\right)^n dx
= \int_0^\infty x^{n-1} (x-2)^n e^{-x} dx.$$
We will compute the asymptotic with Laplace's
method. To do this
put $x=nz$ to obtain
$$n^{n-1} \times n^n \times n
\times \int_0^\infty z^{n-1} (z-2/n)^n e^{-nz} dz
\\ = n^{2n}
\times \int_0^\infty z^{n-1} (z-2/n)^n e^{-nz} dz
.$$
The integrand is
$$\exp(n(\log(z)+\log(z-2/n)-\log(z)/n-z)) = \exp(n f(z)).$$
To find the saddle point we solve the saddle point equation $f'(z)=0$
which produces $$z_0 \approx 2 + \frac{1}{n^2}.$$
The asymptotic is then given by
$$\sqrt{\frac{2\pi}{n|f''(z_0)|}} e^{nf(z_0)}.$$
We have for $f(z_0)$
$$\log 2 + \log(1+1/n^2/2)
+ \log 2 + \log(1-1/n+1/n^2/2)
\\ - 1/n \log 2 - 1/n \log (1 + 1/n^2/2)
- 2 - 1/n^2.$$
The first two terns are
$$f(z_0) \approx 2\log 2 -2 - 1/n \log 2 - 1/n.$$
Furthermore we have
$$f''(z_0) \approx -\frac{1}{2}$$
so $f(z)$ does indeed have a maximum there.
Putting it all together we obtain
$$\sqrt{\frac{4\pi}{n}}
e^{2n\log 2 - 2n - \log 2 - 1}
= 2\sqrt{\frac{\pi}{n}}
\frac{2^{2n}}{2e^{2n+1}}
= \sqrt{\frac{\pi}{n}}
\frac{2^{2n}}{e^{2n+1}}.$$
This yields for our sum
$$n^{2n-1/2} \sqrt{\pi} \frac{2^{2n}}{e^{2n+1}}.$$
Now Stirling's formula yields
$$(2n-1)! \sim \sqrt{2\pi(2n-1)} (2n-1)^{2n-1} e^{-(2n-1)}.$$
This gives for the probability
$$\frac{1}{e^2} \frac{2\sqrt{n}}{\sqrt{2(2n-1)}}
\frac{2^{2n-1} n^{2n-1}}{(2n-1)^{2n-1}}
\sim \frac{1}{e^2}
\left(1+\frac{1}{2n-1}\right)^{2n-1}
\sim \frac{1}{e}$$
and the conjecture is proved.
Remark. These are the details of the inclusion-exclusion
argument. The nodes of the underlying poset represent arrangements
where some set of $q$ couples sit together, possibly more. The Mobius
function of the poset is $(-1)^q.$ A specific arrangement with exactly
$p$ couples sitting together is included in the ${p\choose q}$ nodes
where $q$ couples from the $p$ couples sit together, or more. It is
not included in any nodes representing more than $p$ couples sitting
together because this is impossible when the count is precisely
$p$. Therefore the total weight with the given Mobius function is
$$\sum_{q=0}^p {p\choose q} (-1)^q.$$
This is zero when $p\ge 1$ and one when $p = 0.$ These weights
correctly represent the requirement of arrangements with no couples
sitting together being counted once and all others being counted zero
times.
Addendum. An approach to this problem using the methods of
experimental mathematics might begin by computing the values for small
$n$ using total enumeration and hence verify the formula from
inclusion-exclusion, and then use that to compute the probabilities
for larger $n$, leading to a conjecture that we have $1/e.$ The
enumeration task is practicable to about $n=5$ and confirms the
correctness of the formula. This is the code, which consumes very
little memory. We get starting with one couple
$$0, 2, 32, 1488, 112512,\ldots$$
Observe that we got the correct value for $n=1$ unlike from the
formula, which does not work here because when there is only one
couple the two orderings of its constituents are the same, but the
factor $2^q$ produces a value of two instead of a value of one.
with(combinat);
X :=
proc(n)
option remember;
local perm, circ, count, pos;
count := 0;
perm := firstperm(2*n-1);
while type(perm, `list`) do
circ := [op(perm), 2*n, perm[1]];
for pos to 2*n do
if abs(circ[pos]-circ[pos+1]) = 1
and type(min(circ[pos], circ[pos+1]), `odd`)
then
break;
fi;
od;
if pos = 2*n + 1 then
count := count + 1;
fi;
perm := nextperm(perm);
od;
count;
end;
Q := n -> add(binomial(n,q)*(-1)^q*2^q*(2*n-q-1)!, q=0..n);