Primes as $X-Y$ with $\gcd(X,Y)=1$ and $\operatorname{rad}(XY)$ equal to the product of all primes below $p$ — is this known?

Title: Primes as $X-Y$ with $\gcd(X,Y)=1$ and $\operatorname{rad}(XY)$ equal to the product of all primes below $p$ — is this known?

While experimenting with multiplicative decompositions of primes, I arrived at the following question. I would like to know whether it is already in the literature, and whether it follows from known results on $S$-unit equations.

Setup. Let $p\ge 5$ be prime and let $N_p=\prod_{q<p,\ q\text{ prime}}q$ be the product of all primes strictly below $p$. Say that $p$ admits a partitioned representation if there exist positive integers $X,Y$ with $$X-Y=p,\qquad \operatorname{rad}(XY)=N_p,$$ where $\operatorname{rad}$ denotes the radical (product of distinct prime factors). Equivalently: the set of primes below $p$ is partitioned into two disjoint nonempty sets $A\sqcup B$ and $$p=\prod_{q\in A}q^{a_q}-\prod_{q\in B}q^{b_q},\qquad a_q,b_q\ge 1,$$ so every prime $q<p$ occurs in exactly one of the two products, none is omitted, and no other prime occurs. Note $\gcd(X,Y)=1$ is automatic: a common prime divisor would divide $X-Y=p$, but all prime factors of $XY$ are $<p$.

Examples. $$5=3^2-2^2,\quad 7=2\cdot 5-3,\quad 11=3\cdot 7-2\cdot 5,\quad 13=5\cdot 11-2\cdot 3\cdot 7,$$ $$17=2\cdot 7\cdot 13-3\cdot 5\cdot 11=182-165,\qquad 19=2^2\cdot3\cdot5\cdot17-7\cdot11\cdot13=1020-1001,$$ $$23=7\cdot11\cdot17^2-2\cdot3^2\cdot5\cdot13\cdot19=22253-22230,$$ $$29=3\cdot11\cdot13^2\cdot19\cdot23-2^{12}\cdot5\cdot7\cdot17=2437149-2437120.$$

Computations. Since every prime $q<p$ divides $X$ or $Y$, we have $XY\ge N_p$, and with $X=Y+p$ this forces $Y\gtrsim\sqrt{N_p}$; so no solution can be small. I enumerated all $p$-smooth integers up to a bound $B$ and looked for pairs $(Y,Y+p)$ with disjoint, complementary radicals. A second, independent method uses the fact that $X-Y=p$ and $XY=M$ force $p^2+4M$ to be a perfect square, and enumerates $M=N_p m$ over smooth $m$. Both agree.

The number of representations found is $$3,\;7,\;8,\;7,\;4,\;4,\;1,\;1 \quad\text{for } p=5,7,11,13,17,19,23,29,$$ and none at all for $p=31,37,41,43,47,53,59$ within the searched ranges. For $p=31$ ($N_{31}=6\,469\,693\,230$, so $Y\ge 80\,420$) I found no solution with $Y\le 10^{14}$, i.e. more than $10^{9}$ times the theoretical minimum $\sqrt{N_{31}}$.

Heuristic. Smooth numbers up to $T$ have polylogarithmic density, while $\sqrt{N_p}\asymp e^{p/2}$. The expected number of representations behaves roughly like $c(p)/\sqrt{N_p}$ with $c(p)$ subexponential, so it drops below $1$ around $p\approx 20$ and the sum over $p$ converges. This matches the observed decay to uniqueness at $p=23,29$ and then extinction — suggesting that $\{5,7,11,13,17,19,23,29\}$ may be the complete list.

What I would like to know.

1. Is this problem (or the sequence of representations) already known, perhaps under a different formulation?
2. Does the finiteness of the set of representable primes follow from a known theorem, or only heuristically? I am aware that for fixed $S$ the equation $X-Y=p$ with $X,Y$ $S$-units has finitely many solutions (Evertse), which gives finiteness for each individual $p$, but the union over $p$ is a different matter.
3. Most concretely: is $p=31$ decidable in practice? Writing $X=as^2$, $Y=bt^2$ with $a,b$ squarefree turns it into a family of generalized Pell equations $as^2-bt^2=31$ with $s,t$ smooth — e.g. the $p=29$ solution is $14421\cdot 13^2-595\cdot 64^2=29$. Would a Størmer–Lehmer style argument (primitive divisors in Lucas sequences, Bilu–Hanrot–Voutier) or an effective $S$-unit computation (Baker–Győry bounds plus LLL reduction à la de Weger) be feasible here, and is there existing software that handles the $11$-prime case?

Any pointers to the literature would be very welcome.