Let $\{a_{i,j}\}_{i\ge0,\;j\ge1}$ be defined by
$$
a_{0,j}=\frac1j,\qquad
a_{i,j}=j\bigl(a_{i-1,j}-a_{i-1,j+1}\bigr)\qquad(i\ge1,\ j\ge1).
$$
It is known that $a_{i,1}=B_i$ (Bernoulli numbers, with the usual convention), so the “trivial zeros’’ are $a_{2m+1,1}=0$ for $m\ge1$. The goal is to understand zeros with $j>1$.
I start by deriving a convenient generating function in the column index. Define
$$
P_i(x):=\sum_{j\ge1} a_{i,j}\,x^j\qquad(|x|<1).
$$
From the recursion,
$$
\begin{aligned}
P_i(x)
&=\sum_{j\ge1} j\bigl(a_{i-1,j}-a_{i-1,j+1}\bigr)x^j\\
&=\underbrace{\sum_{j\ge1} j\,a_{i-1,j}\,x^j}_{=\,xP'_{i-1}(x)}
-\underbrace{\sum_{j\ge1} j\,a_{i-1,j+1}\,x^j}_{=:S(x)}.
\end{aligned}
$$
For the second sum, write $m=j+1$. Then
$$
S(x)=\sum_{m\ge2} (m-1)a_{i-1,m}x^{m-1}
=\frac1x\sum_{m\ge2}(m-1)a_{i-1,m}x^{m}
=\frac{1}{x}\bigl(xP'_{i-1}(x)-P_{i-1}(x)\bigr).
$$
Hence
$$
P_i(x)=(x-1)P'_{i-1}(x)+\frac{P_{i-1}(x)}{x}
\qquad\text{with}\qquad
P_0(x)=\sum_{j\ge1}\frac{x^j}{j}=-\log(1-x).
\tag{1}
$$
A more structural representation comes from the observation that
$$
\frac1j=\int_0^1 x^{\,j-1}\,dx.
$$
Using this and the recursion once, and then integrating by parts, one proves by induction on $i$ that there exist polynomials $Q_i(x)$ such that
$$
a_{i,j}=(-1)^i\int_0^1 x^{\,j-1}\,Q_i(x)\,dx,
\tag{2}
$$
where $Q_0(x)=1$ and
$$
Q_i(x)=x\,\frac{d}{dx}\Bigl[(1-x)\,Q_{i-1}(x)\Bigr]\qquad(i\ge1).
\tag{3}
$$
For the convenience of later algebra, it is best to write $Q_i$ explicitly in terms of Stirling numbers of the second kind $S(i,\ell)$. Using the standard recurrence $S(i,\ell)=\ell\,S(i-1,\ell)+S(i-1,\ell-1)$ and checking the base case $i=1$, induction on $i$ shows
$$
Q_i(x)=\sum_{\ell=1}^{i}(-1)^{\ell}\,\ell!\,S(i,\ell)\,x^{\ell}.
\tag{4}
$$
Plugging (4) into (2) and integrating term‑by‑term yields the exact partial‑fraction formula
$$
\boxed{\quad
a_{i,j}=(-1)^i\sum_{\ell=1}^{i}(-1)^{\ell}\,\ell!\,S(i,\ell)\,\frac{1}{\,j+\ell\,}
\quad}
\qquad(i\ge1,\ j\ge1).
\tag{5}
$$
The sum in (5) always collapses to a single rational function in $j$. If we put $(j)_{i+1}:=j(j+1)\cdots(j+i)$, then there is a unique monic polynomial $p_i(j)$ of degree $i$ with integer coefficients such that
$$
\boxed{\quad
a_{i,j}=\frac{p_i(j)}{(j)_{i+1}}\quad}
\qquad(i\ge0,\ j\ge1).
\tag{6}
$$
Clearing denominators in (5) gives an explicit numerator
$$
p_i(j)=(-1)^i\sum_{\ell=1}^i (-1)^\ell\,\ell!\,S(i,\ell)\,
\prod_{\substack{0\le m\le i\\ m\ne \ell}}(j+m).
\tag{7}
$$
There is also a simple recurrence for these numerators. Substituting (6) into the defining recursion,
$$
\begin{aligned}
a_{i,j}
&=j\Bigl(\frac{p_{i-1}(j)}{(j)_i}-\frac{p_{i-1}(j+1)}{(j+1)_i}\Bigr)\\
&=\frac{j(j+i)\,p_{i-1}(j)-j^2\,p_{i-1}(j+1)}{(j)_{i+1}},
\end{aligned}
$$
so
$$
\boxed{\quad
p_0(j)=1,\qquad
p_i(j)=j\,(j+i)\,p_{i-1}(j)-j^2\,p_{i-1}(j+1)\quad(i\ge1). \quad}
\tag{8}
$$
Several consequences follow quickly.
First, for all $i\ge2$ the factor $j^2$ divides $p_i(j)$. This is immediate by induction from (8): the first term has a factor $j^3$ if $j^2\mid p_{i-1}$, and the second term is visibly divisible by $j^2$.
Second, evaluating (6) at $j=1$ gives
$$
a_{i,1}=\frac{p_i(1)}{(i+1)!}=B_i,
$$
so when $i$ is odd and at least $3$, one gets $p_i(1)=0$ and hence $(j-1)$ divides $p_i(j)$. These are the trivial zeros.
Third, if we define $q_i(j):=p_i(j)/j^2$ for $i\ge2$, then dividing (8) by $j^2$ and letting $j\to0$ gives
$$
q_i(0)=-\,p_{i-1}(1)=-\,i!\,B_{i-1}.
\tag{9}
$$
In particular, when $i$ is even and at least $4$ one has $B_{i-1}=0$, so $q_i(0)=0$ and $j^3$ divides $p_i(j)$.
Fourth, evaluating the explicit sum (7) at negative integers yields exact values at the “Lagrange nodes’’:
$$
p_i(-\ell)=(-1)^i\,(\ell!)^2\,(i-\ell)!\,S(i,\ell)\qquad(1\le\ell\le i).
\tag{10}
$$
The special zero in the row $i=4$ is now easy to see directly from (8). Starting from $p_0(j)=1$, the first few steps are
$$
p_1(j)=j,\qquad
p_2(j)=j^2,\qquad
p_3(j)=j^2(j-1),\qquad
p_4(j)=j^3(j-5),
$$
hence
$$
a_{4,j}=\frac{j^3(j-5)}{(j)_5}
=\frac{j^2\,(j-5)}{(j+1)(j+2)(j+3)(j+4)}.
$$
Therefore $a_{4,5}=0$, and no other positive integer $j$ is a zero in that row.
To understand the overall sign pattern, expand the summands in (5) for large $j$:
$$
\frac{1}{j+\ell}=\frac1j-\frac{\ell}{j^2}+O\!\Bigl(\frac{1}{j^3}\Bigr)\qquad(j\to\infty).
$$
Using the Stirling identity
$$
\sum_{\ell=0}^{i}(-1)^{\ell}\,\ell!\,S(i,\ell)=(-1)^i,
\tag{11}
$$
one gets
$$
a_{i,j}=\frac1j+O_i\!\Bigl(\frac{1}{j^2}\Bigr)\qquad(j\to\infty).
\tag{12}
$$
More concretely, from (5) and $|\ell/(j+\ell)-1/j|\le \ell/j^2$ for $j>\ell$,
$$
\Bigl|\,a_{i,j}-\frac1j\,\Bigr|
\le \frac{1}{j^2}\sum_{\ell=1}^{i}\ell\,\ell!\,S(i,\ell).
\tag{13}
$$
Hence for each fixed $i$ there is an $N(i)$ such that $a_{i,j}>0$ for all $j\ge N(i)$. In particular, every row has only finitely many sign changes, so only finitely many integer $j$ can be zeros in that row.
Finally, combine this with the rational‑root test. Let
$$
\tilde p_i(j):=\frac{p_i(j)}{j^{\rho_i}(j-1)^{\varepsilon_i}},
\qquad
\rho_i=\begin{cases}
2,& i\ge2\ \text{and }B_{i-1}\ne0,\\
3,& i\ge4\ \text{even},
\end{cases}
\quad
\varepsilon_i=\begin{cases}
1,& i\ge3\ \text{odd},\\
0,& \text{otherwise}.
\end{cases}
$$
Then $\tilde p_i(j)\in\mathbb{Z}[j]$ is monic with nonzero constant term. By (9), when $i$ is odd and at least $3$,
$$
\tilde p_i(0)=-\,i!\,B_{i-1}.
$$
When $i$ is even and at least $4$, dividing (8) by $j^3$ and letting $j\to0$ (using $B_{i-1}=0$) gives
$$
\tilde p_i(0)=-\,i!\,B_{i-2}-\,p'_{i-1}(1).
$$
Therefore any integer zero $j>1$ of $a_{i,j}$ must be a divisor of the corresponding integer $|\tilde p_i(0)|$. Together with the eventual positivity above, this reduces the search for integer zeros in any fixed row to a finite check. Carrying out that check in the first several rows shows that the only non‑trivial integer zero with $j>1$ that occurs is
$$
a_{4,5}=0.
$$
For completeness, here is the short derivation of (2)–(4). The base case $i=0$ is $a_{0,j}=\int_0^1 x^{j-1}\,dx$. Assume $a_{i-1,j}=(-1)^{i-1}\int_0^1 x^{j-1}Q_{i-1}(x)\,dx$. Then
$$
\begin{aligned}
a_{i,j}
&=j\int_0^1 (-1)^{i-1}\bigl(x^{j-1}-x^j\bigr)Q_{i-1}(x)\,dx\\
&=(-1)^i\int_0^1 x^{j}\,\frac{d}{dx}\bigl[(1-x)Q_{i-1}(x)\bigr]\,dx
\quad\text{(integration by parts; the boundary term vanishes)}\\
&=(-1)^i\int_0^1 x^{j-1}\,\Bigl\{x\,\frac{d}{dx}\bigl[(1-x)Q_{i-1}(x)\bigr]\Bigr\}\,dx,
\end{aligned}
$$
which is (2)–(3). Writing $Q_{i-1}(x)=\sum_{\ell}c_{i-1,\ell}x^\ell$, a direct coefficient computation using $S(i,\ell)=\ell\,S(i-1,\ell)+S(i-1,\ell-1)$ shows that the choice
$$
c_{i,\ell}=(-1)^\ell\,\ell!\,S(i,\ell)
$$
satisfies (3); this is (4). Substituting (4) into (2) yields (5) immediately. The rest follows as shown above.
$a_{9,j}$ for 
$a_{9,2^k}$ for