This answer is using pieces from my related answer to the related question MSE question 4779869,
where many details to the computations below are shown. Some simple blocks of sage code collect
experimentally all $b$-polynomials involved. The results are tabulated - as much as the MSE site permits to do.
This answer gives proofs. Some details are missing, in order to have a complete proof, but the path / the story should be clear.
We start with the equation
$$v^2 =(a^3-b)(b^3-a)\ ,$$
considered over the function field $\Bbb Q(a)$.
Instead of furthermore using a specialization $a=n$. I will work with $a=A$ to make this specialization
"fixed" in notation. Also, i will use $u$ instead of $b$. So we deal with the curve:
$$
C=C(A)=C(A)_{u,v}\ :\qquad v^2 =(u^3-A)(A^3-u)\ .
$$
As explained in loc. cit. (MSE), we move birationally
$(u,v)\leftrightarrow (x,y)$, and in it in fact $u\leftrightarrow x$ only, via
$x=c/(A^3 - u) - A^6$, $u=A^3-c/(x+A^6)$, $c:=A^9-A$,
to the elliptic curve
$$
\bbox[yellow]{\qquad
E=E(A)=E(A)_{x,y}\ :\qquad\qquad y^2 = x^3 - 3A^4x - A^2(A^8+1)\ .
\qquad
}
$$
Dividing by $A^6$, the obvious substitutions lead to the form:
$$
Y^2 = X^3 - 3X +\left(A^4+\bar A^4\right)\ .
$$
(A bar is used for the inverse.)
Here we have the reciprocal symmetry explicitly shown, and the particular solution
$X=A^2+1+\bar A^2$, resp. $x=A^4 + A^2 +1$ can be checked. It corresponds to $u=a$.
For $x$ consider the lift $G=G(A)$ on $E(A)$ given by
$$
G(A) = (A^4 + A^2 + 1, \ A^6 + A^4 + A^2 +1)\ .
$$
Then the many rational functions $b_k(A)$ from the question (with $A=n$) occur in the following manner.
Start with $G(A)$. Build $kG(A)=(x_k(A),y_k(A))\in E(\ \Bbb Q(A)\ )$, extract $x_k(A)$ and pass from
$E$ to $C$ to obtain $u_k(A)=b_k(A)$ and the corresponding $v_k(A)$.
In the sequel we need the addition formula for points of the above curve.
Given $P_1,P_2$, the sum is $P_3$ - all points being of the shape $P_k=(x_k,y_k)\ne \infty$ - we have:
$$
\tag{$*$}
$$
$$
\begin{aligned}
x_3 &= -x_1-x_2 +\lambda^2\ ,\\
y_3 &= -y_1 -(x_3-x_1)\lambda\,\qquad\text{ where $\lambda$ is the slope }\lambda=\frac{y_2-y_1}{x_2-x_1}
\end{aligned}
$$
in the case of $x_1\ne x_2$, and else we pass to the limit, $x_3=(3x_1^2-3A^4)/y_1$.
It may be useful to record some relations for the passage from one world to the other one, in between notations as in MSE:
$$
\begin{aligned}
u &=A^3 -\frac 1{x_1}=A^3-\frac c{x_2}=A^3-\frac c{x+A^6}=\frac {A^3x+A}{x+A^6}
=\begin{bmatrix} A^3 & A \\ 1 & A^6\end{bmatrix}\cdot x\ ,\\
x &= x_2-A^6 =cx_1-A^6=\frac c{A^3-u}-A^6=\frac{A^6u-A}{-u+A^3}
=\begin{bmatrix} A^6 & -A \\ -1 & A^3\end{bmatrix}\cdot u\ ,\\
v &=\frac{y_1}{x_1^2}=\frac {y_1}c\cdot\frac{c^2}{x_2^2}=\frac{cy}{(x+A^6)^2}\ ,\\
y &= y_2=cy_1=cvx_1^2=\frac {cv}{(A^3-u)^2}\ .
\end{aligned}
$$
Questions:
(i) Consider $b_k(A)/A$. Is this expression of the shape $B_k(A, 1)/B_k(1,A)$
for some homogeneous polynomial $B_k$? If yes, we may and do assume below that it is normed so that $B_k(A,1)$ is monic, i.e. $B_k(1,0)=1$.
(ii) If yes, what is $B_k(0, 1)$, the free coefficient of $B_k(A,1)$?
(iii) If yes, what is $B_k(1, 1)$, the sum of the coefficients of $B_k$?
Answers to the questions are following.
Let us assume (i). Then we have
$$
\begin{aligned}
x_k(A)
&=
\frac c{A^3-u_k(A)}-A^6
=
\frac {A^9-A}{A^3-u_k(A)}-A^6
=
\frac {A^8-1}{A^2-\frac{B_k(A,1)}{B_k(1,A)}}-A^6
\\
&\overset{A\to 0}\longrightarrow
\frac {0-1}{0-\frac{B_k(A,1)}{B_k(1,A)}}-0=\frac 1{B_k(0,1)}
\end{aligned}
$$
To show now the formula claimed on experimental data in the OP
$$
\bbox[yellow]{
\qquad
B_k(0,1)={\color{red}{k^2}}
\qquad
}
$$
it is enough to show:
Lemma: Consider $x_k(A)$, $y_k(A)$ as rational functions of $A$.
Then zero is not a pole and not a zero for them, and we have moreover:
$$
x_k(0)=\frac 1{\color{red}{k^2}}\ ,\qquad y_k(0)=\frac 1{k^3}\ .
$$
The above can be easily checked for the first few values of $k$,
let us have explicit samples:
$$\tag{$\dagger$}
$$
$$
\small
\begin{aligned}
(x_1,y_1) &= G=(A^4 + A^2 + 1,\ A^6 + A^4 + A^2 + 1)\ ,\\
(x_2,y_2) &= 2G=\left(
\frac 1{\color{red}4}(A^4 + 10A^2 + 1),\ \frac 18(A^6 - 17A^4 - 17A^2 + 1)\right)
\ ,\\
(x_3,y_3) &= 3G
\\
&= \left(
\frac{A^{12} + 101 A^{10} + 291 A^8 + 238 A^6 + 291 A^4 + 101 A^2 + 1}
{{\color{red}9} (A^2 - 1)^4 }\ ,\right.
\\
&\qquad\left.
\
\frac{
(A^{12} - 214 A^{10} - 2481 A^8 - 2804 A^6 - 2481 A^4 - 214 A^2 + 1)
(A^4 + 1)
(A^2 + 1)}
{27 (A^2 - 1)^6}
\right)\ ,
\end{aligned}
$$
and so on.
Proof of the lemma: For $k=1,2$ things are checked. Assume by induction the formula holds for some $k\ge 2$.
Then the addition law (*) on the curve computed for $(k+1)G=G+kG$, then evaluated in $A=0$ gives:
$$
\small
\begin{aligned}
\lambda_k(0) &:= \frac{y_1-y_k}{x_1-x_k}(0)=\frac{1-\frac 1{k^3}}{1-\frac 1{k^2}}=1 + \frac 1{k(k+1)}\ .
\\
x_{k+1}(0) &= -x_1(0)-x_k(0)+\lambda_k^2(0)=-1-\frac 1{k^2}+ \left(1 + \frac 1{k(k+1)}\right)^2=\frac 1{\color{red}{(k+1)^2}}\ ,\\
y_{k+1}(0) &= -y_1(0)-(x_{k+1}(0)-x_1(0))\lambda_k(0)
=-1-\left(\frac 1{(k+1)^2-1}\right)\cdot\frac{k^2+k+1}{k(k+1)}
\\
&=-1+\frac{(k+2)(k^2+k+1))}{(k+1)^3}=\frac 1{(k+1)^3}
\ .
\end{aligned}
$$
$\square$
We have so far modulo the mirror property (i) the formula (ii), which was the target as i started to type.
(Soon (iii) also was in there.)
So let us show (ii). The sample values for $x_k,y_k$ above motivate the following
Lemma: The rational functions $x_k(A)$, $y_k(A)$ use symmetric polynomials in their numerator
and denominator (for suitable representations). Their total $A$-degree (at $\infty$ is $4$, respectively $6$.
Proof: Use the symmetry in the $(X,Y)$ world, then come back to the $(x,y)$ world.
$\square$
The main step is:
Lemma: For all $k\ge 1$ we have $x_k=p_k/q_k$ with reciprocal polynomials
$p_k,q_k=d_k^2$. Let $P_k,Q_k$ be the homogenized versions of $p_k,q_k$.
Then there exist further homogeneous polynomials $B_k,V_k$ with
$$
\bbox[lightyellow]{\qquad
\begin{aligned}
u_k(A) &= A\cdot \frac {B_k(A,1)}{B_k(1,A)}=\frac{A^2P_k(A,1) + Q_k(A,1)}{A^6Q(A,1) +P(A,1)}\ ,\\
v_k(A) &= A\cdot \frac {V_k(A,1)}{B^2_k(1,A)}\ ,
\end{aligned}\qquad}
$$
Explicitly:
$$
\bbox[lightyellow]{
\qquad
b_k(A)=B_k(A,1)$ \text{ is a factor of }a^2p_k(A)+ q_k(A)\ ,\qquad p_k,q_k\text{ reciprocal.}
\qquad
}
$$
The above factor is obtained after the simplification of the fraction $(A^2p_k+q_k)/(A^6q_k+p_k)$.
A known simplification occurs for $k=3j$, when both $(A^2p_k+q_k)$ and
$A^6q_k +p_k$ are divisible by the
$$\color{blue}{
\text{Known factor: }(A^6+A^4+A^2+1)
}
\ .
$$
Moreover, $V_k$ is antisymmetric, $V_k(s,t)=-V_k(t,s)$.
As an illustration let us plot the involved polynomials for first small $k$ values.
(For bigger values use the linked code.)
For $k=1$ we have $u_1=b_1=A$, and $v_1=A^3-A$, $B_1=1$, $V_1(A,1)= A^2-1$.
And indeed, $V_1(s,t)=s^2-t^2=-(t^2-s^2)=-V_1(t,s)$.
For $k=2$ computations deliver
$x_2=\frac 14(A^4+10A^2+1)$,
$u_2=b_2=A\cdot B_2(n,1)/B_2(1,n)$ with
$$
B_2(n,1)= n^6 + 10n^4 + n^2 + 4\ ,
$$
and $v_2=A\cdot V_2(n,1)/B_2(1,n)^2$
with
$$V_2(n,1)=2n^{14} - 34n^{12} - 34n^{10} + 2n^{8} - 2n^{6} + 34n^{4} + 34n^{2} - 2\ .$$
Further experimental data are on parallel answers, and/or in my linked MSE-post.
Proof: We show only what we need, the formula for $u_k$. (Olny $B$-computations, no $V$-computation.)
Recall that $x_k$ has a representation as $P(A,1)/Q(A,1)$ with degree $4$,
so $\deg P=4+d$ for some $d$, and $\deg Q =d$.
One can see inductively, that only even powers occur.
Then
$$
u_k =
\begin{bmatrix} A^3 & A \\ 1 & A^6\end{bmatrix}\cdot \frac PQ
=A\cdot \frac{A^2P(A,1)+Q(A,1)}{A^6Q(A,1) + P(A,1)}
\ .
$$
To have a better image, let us put the coefficients of $A^2P, Q; A^6Q, P$ on parallel lines.
The coefficients of $P$ are $p_0,p_1=0,p_2,\dots,p_2,p_1=0,p_0$. They are corresponding to degree $4+d$.
The coefficients of $Q$ are similarly $q_0,q_1=0,q_2,\dots,q_2,q_1=0,q_0$.
Then the numerator has coefficients:
a0 a1 a2 a3 a4 a5 a6 ......... a7 a6 a5 a4 a3 a2 a1 a0 0 0 --- from A^2 P
0 0 0 0 0 0 q0 ......... q8 q7 q6 q5 q4 q3 q2 q1 q0 --- from Q
And the denominator has the coefficients:
q0 q1 q2 q3 q4 q5 q6 q7 q8 .... q0 0 0 0 0 0 0 --- from A^6 Q
0 0 a0 a1 a2 a3 a4 a4 a6 .... a6 a5 a4 a3 a2 a1 a0 --- from P
It is clear that the claimed reciprocity holds,
the numerators makes a polynomial $B(A,1)=A^2P(A,1)+Q(A,1)$, and
the denominator is exactly $B(1,A)=A^6Q(A,1) +P(A,1)$.
This shows (i), and thus the formula for the free coefficient
$\color{red}{k^2}$ in the numerator of $b_k$ is proven.
$\square$
It remains (iii).
I will invest some typing effort to show the path.
For this, note that we need in the last notation
$B(1,1)=P(1,1)+Q(1,1)$.
So we want to specialize $A=1$. The corresponding "limit curve"
$y^2=x^3 -3x-2$
is no longer elliptic,
so the group law degenerates and can be easily examined.
I will say some words in this EDITED answer on this.
It may be also of interest to obtain an explicit recursion for the polynomials $B$.
Computations need specific data.
We define homogeneous polynomials $P_k,D_k,R_k$ with dehomogenized versions $p_k,d_k,r_k$,
that appear in $kG=(x_k(A),y_k(A))$, $k\ge 1$,
$$
x_k(A)=\frac {p_k(A)}{d_k(A)^2}\ ,\qquad
y_k(A)=\frac {r_k(A)}{d_k(A)^3}\ ,
$$
where the base factor $d_k(A)$ of the denominators is relatively prime to the monic numerators
$p_k,q_k$
of $x_k,y_k$.
The recursion starts with the values:
$$
\begin{aligned}
p_1 &= A^4 + A^2 +1\ , & r_1 &= A^6 + A^4 + A^2 + 1\ , & d_1 = 1\ ,\\
p_2 &= A^4 + 10A^2 +1\ , & r_2 &= (A^2 + 1)(A^2 - 4A - 1)(A^2 + 4A - 1)\ , & d_2 = 2\ ,\\
\end{aligned}
$$
and the recursion step is determined by applying the formula $(*)$, with an in between step of computing
the slopes $\lambda_k$. We will decide later which is the recursion step $s$. (It may be one, two, three, six,
depending on the joy of typing arguments and the specific purpose.)
Formulas are easily written:
$$
\begin{aligned}
(x_{k+s}, y_{k+s}) &= (k+s)G=kG\oplus sG=(x_k,y_k)\oplus(x_s,y_s)\ ,
\\
\lambda_k=\lambda_{k,s}
&:=
\frac {y_k-y_s}{x_k-x_s}
=
\frac 1{d_k}\cdot\frac{r_k -y_s d_k^3}{p_k - x_s d_k^2}\ ,
\\
\frac{p_{k+s}}{d_{k+s}^2}
&=x_{k+s}
\\
&=-x_s-x_k+\lambda_k^2
\\
&=-x_s-\frac{p_k}{d_k^2}+
\frac 1{d_k^2}\cdot\frac{(r_k -y_sd_k^3)^2}{(p_k-x_s d_k^2)^2}
\ ,
\\
\frac{r_{k+s}}{d_{k+s}^3}
&=y_{k+s}
\\
&=-y_s-(x_{k+s}-x_s)\lambda_k\ .
\end{aligned}
$$
From here extract the recursion formula for $d_{k+s}$, $p_{k+s}$, $r_{k+s}$.
Finally consider the homogenized versions $D_k,P_k,R_k$, and set:
$$
B_k = \underbrace{(A^2P_k + D_k^2)}_{B_k^*}/\text{factor}\ .
$$
Here the factor is the gcd of $B_k^*:=(A^2P_k+Q_k)$ and $(A^2Q_k+P_k)$, as in the Main Lemma.
We have a known factor.
No further multiplicative cancellation occurs in the experiments with small numbers $k$,
but if they occur we do not implement them here, sorry.
This is the only missing point, and let us address it as a "gap".
The table for the degrees is:
$$
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|}
\hline
k & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline
\deg d_k & 0 & 0 & 4 & 6 & 12 & 16 & 24 & 30 & 40 & 48 \\\hline
\deg p_k = 2\deg d_k + 4 & 4 & 4 & 12 & 16 & 28 & 36 & 52 & 64 & 84 & 100 \\\hline
\deg r_k = 3\deg d_k + 6 & 6 & 6 & 18 & 24 & 42 & 54 & 78 & 96 & 126 & 150 \\\hline
\deg B_k & 1 & 6 & 8 & 18 & 24 & 38 & 48 & 66 & 80 & 102 \\\hline
\end{array}
$$
A one-step recursion would not take care of the fact that the degree drops
by chance for $k=1$, and this propagates now with step two. So we may want to use the step two
while giving the one or the other proof.
Let us break into
two cases for the above reason, odd and even values of $k$.
$$
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|}
\hline
k & 1 & 3 & 5 & 7 & 9 & 11 & 13 & 15 & 17 & 19 \\\hline
\deg d_k & 0 & 4 & 12 & 24 & 40 & 60 & 84 & 112 & 144 & 180 \\\hline
\deg p_k = 2\deg d_k + 4 & 4 & 12 & 28 & 52 & 84 & 124 & 172 & 228 & 292 & 364 \\\hline
\deg r_k = 3\deg d_k + 6 & 6 & 18 & 42 & 78 & 126 & 186 & 258 & 342 & 438 & 546 \\\hline
\deg B_k = 2\deg d_k + \delta_{k1} & 1 & 8 & 24 & 48 & 80 & 120 & 168 & 224 & 288 & 360 \\\hline
\end{array}
$$
The extracted formulas are:
$$
\tag{$k$ odd $\ge1$}
$$
$$
\begin{aligned}
\deg d_k &= \frac 12(k^2-1)\ ,\\
\deg p_k &= (k^2+3)\ ,\\
\deg r_k &= \frac 32(k^2+3)\ ,\\
\deg B_k &= k^2-1=2\deg d_k\ ,
\end{aligned}
$$
which can be tested inductively starting from $k=3$ with step two modulo the "gap".
$$
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|}
\hline
k & 2 & 4 & 6 & 8 & 10 & 12 & 14 & 16 & 18 \\\hline
\deg d_k & 0 & 6 & 16 & 30 & 48 & 70 & 96 & 126 & 160 \\\hline
\deg p_k = 2\deg d_k + 4 & 4 & 16 & 36 & 64 & 100 & 144 & 196 & 256 & 324 \\\hline
\deg r_k = 3\deg d_k + 6 & 6 & 24 & 54 & 96 & 150 & 216 & 294 & 384 & 486 \\\hline
\deg B_k = 2\deg d_k + 6 & 6 & 18 & 38 & 66 & 102 & 146 & 198 & 258 & 326 \\\hline
\end{array}
$$
And the pattern is:
$$
\tag{$k$ even $\ge2$}
$$
$$
\begin{aligned}
\deg d_k &= \frac 12(k^2-4)\ ,\\
\deg p_k &= k^2\ ,\\
\deg r_k &= \frac 32 k^2\ ,\\
\deg B_k &= k^2+2\ ,
\end{aligned}
$$
which can also be tested inductively.
All the above modulo the "gap".
This is a line break in the proof of (iii), recall what we have so far.
Some chain of polynomials $d,p,r;B$ have been constructed, indexed by $k$,
so that $p/d^2$ and $r/d^3$ reproduce the multiples of $G$ as elements in $E(A)(\ \Bbb Q(A)\ )$,
they are defined recursively using only ring operations. The situation
so far allows to evaluate all the polynomials in the point one.
Then $k\cdot G(1)$ makes no sense any longer, since $E(1)$ is not an elliptic curve,
but the operations involving the polynomials $d,p,r;B$ do make sense, when evaluated in one.
Now i have to say something, and the best claim i can make in this framework is to
use experimental data, based on it make short sentences.
We simply plot the values of the $d,p,r,b$ polynomials, when computed in one.
I am providing the sage code, since it is simple:
def EA(A, field=QQ):
return EllipticCurve(field, [-3*A^4, -A^10 - A^2])
R0.<n> = PolynomialRing(ZZ)
F = FractionField(R0)
E = EA(n, field=F)
G = E.point([n^4 + n^2 + 1, n^6 + n^4 + n^2 + 1])
for k in [1..20]:
xk, yk = (k*G).xy() # xk, yk are in F, rational functions of n
dk = R0( sqrt(xk.denominator()) )
pk = R0( xk.numerator() )
rk = R0( yk.numerator() )
Bk = (n^3 - (n^9 - n)/(xk + n^6)).numerator()
dk1, pk1, rk1, Bk1 = dk(1), pk(1), rk(1), Bk(1)
Bk1_formula = pk1 + dk1^2
info = ' & '.join([latex(factor(expr)) if expr else str(expr)
for expr in (dk1, pk1, rk1, Bk1, Bk1_formula)])
print(f"{k} & {info}\\\\\\hline")
And we obtain:
$$
\begin{array}{|r|r|r|r|l|}
\hline
k & d_k(1) & p_k(1) & r_k(1) & b_k(1)=B_k(1,1) & p_k(1) +d_k^2(1)\\\hline
1 & 1 & 3 & 2^{2} & 1 & 2^{2}\\\hline
2 & 2 & 2^{2} \cdot 3 & -1 \cdot 2^{5} & 2^{4} & 2^{4}\\\hline
3 & 0 & 2^{10} & 2^{15} & \color{blue}{2^{8}} & 2^{10}\\\hline
4 & -1 \cdot 2^{7} & 2^{14} \cdot 3 & -1 \cdot 2^{23} & 2^{16} & 2^{16}\\\hline
5 & -1 \cdot 2^{12} & 2^{24} \cdot 3 & 2^{38} & 2^{24} & 2^{26}\\\hline
6 & 0 & 2^{36} & 2^{54} & 2^{36} & 2^{36}\\\hline
7 & 2^{24} & 2^{48} \cdot 3 & 2^{74} & 2^{48} & 2^{50}\\\hline
8 & 2^{31} & 2^{62} \cdot 3 & -1 \cdot 2^{95} & 2^{64} & 2^{64}\\\hline
9 & 0 & 2^{82} & -1 \cdot 2^{123} & \color{blue}{2^{80}} & 2^{82}\\\hline
10 & -1 \cdot 2^{49} & 2^{98} \cdot 3 & -1 \cdot 2^{149} & 2^{100} & 2^{100}\\\hline
11 & -1 \cdot 2^{60} & 2^{120} \cdot 3 & 2^{182} & 2^{120} & 2^{122}\\\hline
12 & 0 & 2^{144} & 2^{216} & 2^{144} & 2^{144}\\\hline
13 & 2^{84} & 2^{168} \cdot 3 & 2^{254} & 2^{168} & 2^{170}\\\hline
14 & 2^{97} & 2^{194} \cdot 3 & -1 \cdot 2^{293} & 2^{196} & 2^{196}\\\hline
15 & 0 & 2^{226} & -1 \cdot 2^{339} & \color{blue}{2^{224}} & 2^{226}\\\hline
16 & -1 \cdot 2^{127} & 2^{254} \cdot 3 & -1 \cdot 2^{383} & 2^{256} & 2^{256}\\\hline
17 & -1 \cdot 2^{144} & 2^{288} \cdot 3 & 2^{434} & 2^{288} & 2^{290}\\\hline
18 & 0 & 2^{324} & 2^{486} & 2^{324} & 2^{324}\\\hline
19 & 2^{180} & 2^{360} \cdot 3 & 2^{542} & 2^{360} & 2^{362}\\\hline
20 & 2^{199} & 2^{398} \cdot 3 & -1 \cdot 2^{599} & 2^{400} & 2^{400}\\\hline
\end{array}
$$
All the above, for it is easier for me to make sentences now.
The last two columns almost agree for $k\ge 2$, this is just a check of the claimed formula for $B_k$
in terms of $P_k$ and $D_k$. The blue differences come from the "known factor" (from the Main Lemma),
it is $\color{blue}{(A^6+A^4+A^2+1)_{\text{in }A=1} = 1+1+1+1=4}$.
For $k$ a multiple of $3$ the denominator vanishes in one.
(So for some arguments the induction / recursion has a convenient step three. Or maybe six.)
For the proof of this fact, note that the expression $(\dagger)$ for $3G$ comes with denominators
which are powers of
$d_3(A)=3(A^2-1)^2$. Simple code gives for $d_6(A)$ the expression
$$
\small
2 \cdot 3 \cdot (A - 1)^{2} \cdot (A + 1)^{2} \cdot (A^{12} - 214A^{10} - 2481A^{8} - 2804A^{6} - 2481A^{4} - 214A^{2} + 1)\ .
$$
By induction, $d_{3k}(A)$ has a zero of order two in $A=1$.
To have a better argument for the vanishing, note that the specialization $A=0$
leads to multiples of $G$ with non-vanishing $y$-component. So let switch the view in the curve
$E(A)$ from the afine $(x,y)$ perspective with $z=1$ to the $(X,Z)$ perspective with $Y=1$,
so that in the projective space we have $[x:y:1]=[X:1:Z]$.
Then we have in parallel:
$$
\begin{aligned}
&E(A)_{x,y}\ :\ & y^2 &= x^3 -3A^2x - A^2(A^8+1) \ , \\
&E(A)_{X,Z}\ :\ & Z &= X^3 -3A^2XZ^2 - a^2(A^8+1)Z^3 \ , \\[2mm]
&E(1)_{x,y}\ :\ & y^2 &= x^3 -3x - 2 = (x-2)(x+1)^2 \ , \\
&E(1)_{X,Z}\ :\ & Z &= X^3 -3XZ^2 -2Z^3 = (X-2Z)(X+Z)^2 \ , \\[2mm]
&& \infty(A)_{x,y} &= [0:1:0]= \infty \ ,\\
&& \infty(A)_{X,Z} &= [0:1:0]=(0,0)_{X,Z} \ ,\\[2mm]
&& G(1)_{x,y} &= (3,4) = [3:4:1]\ ,\\
&& G(1)_{X,Z} &= [3:4:1]=\left(\frac 34,\frac 14\right)_{X,Z}\ ,\\
&& (2G)(A=1)_{x,y} &= (3,-4) = [3:-4:1]\ ,\\
&& (2G)(A=1)_{X,Z} &= [3:-4:1]=\left(-\frac 34,-\frac 14\right)_{X,Z}\ ,\\
&& (3G)(A=1)_{x,y} &= \infty = [0:1:0]\ ,\\
&& (3G)(A=1)_{X,Z} &= [0:1:0]=(0,0)_{X,Z}\ ,\\
\end{aligned}
$$
and from here, we can continue the computation of $(kG)$ computed in $A=1$ using the same
geometric definition of multiples on the degenerated elliptic curve $E(1)$, as we do on $E(A)$ over $\Bbb Q(A)$.
For instance, in order to compute $(4G)_{X,Z}$ in $A=1$ we intersect the line through $G(A=1)=(3/4,1/4)_{X,Z}$
and $(3G)(A=1)=(0,0)_{X,Z}$ -- which is $X=3Z$ -- with the cubic $E(1)$.
The substitution of $X=3Z$ into the equation of $E(1)$ gives in $Z$ the equation $Z=(27-9-2)Z^3=16Z^3$, so $Z$ is among
$0,\pm \frac14$. The third point on the line and the curve has $Z=-1/4$, corresponding to $(-3/4,-1/4)_{X,Z}$,
and it remains to draw through this last point and the infinity point, intersect once more to get $G$ again.
We have the repetition after three steps, so the $X$-component always vanishes for $3kG(A=1)$, so the
corresponding $x$-component vanishes.
This shows that the point $(x,y)=(p/d^2,r/d^3)=[pd:r:d^3]$ is always among $[3:4:1]$, $[3:-4:1]$, $[0:1:0]$.
The period for the pattern is three.
An inductive argument shows that - modulo the "gap" - we have:
$$
\begin{aligned}
k &= 6j + 0 & d_k(1) &= 0\ ,\\
& & p_k(1) &= 2^{k^2}\ ,\\
& & r_k(1) &= 2^{3k^2/2} = p_k(1)^{3/2}\ ,\\
& & b_k(1) &= p_k(1)+d_k(1)^2=2^{k^2}\ ,
\\[2mm]
k &= 6j + 1 & d_k(1) &= 2^{(k^2-1)/2}\ ,\\
& & p_k(1) &= 3\cdot 2^{k^2-1}=3d_k(1)^2\ ,\\
& & r_k(1) &= 2^{(3k^2+1)/2}\ ,\\
& & b_k(1) &= p_k(1)+d_k(1)^2=2^{k^2+1}\ ,
\\[2mm]
k &= 6j + 2 & d_k(1) &= 2^{(k^2-2)/2}\ ,\\
& & p_k(1) &= 3\cdot 2^{k^2-2}=3d_k(1)^2\ ,\\
& & r_k(1) &= -2^{(3k^2-2)/2}\ ,\\
& & b_k(1) &= p_k(1)+d_k(1)^2=2^{k^2}\ ,
\\[2mm]
k &= 6j + 3 & d_k(1) &= 0\ ,\\
& & p_k(1) &= 2^{k^2+1}\ ,\\
& & r_k(1) &= -2^{3(k^2+1)/2} = -p_k(1)^{3/2}\ ,\\
& & b_k(1) &= 2^{k^2-1}=(p_k(1)+d_k(1)^2)\color{blue}{/4} = p_k(1)\color{blue}{/4}\ ,
\\[2mm]
k &= 6j + 4 & d_k(1) &= -2^{(k^2-2)/2}\ ,\\
& & p_k(1) &= 3\cdot 2^{k^2-2}=3d_k(1)^2\ ,\\
& & r_k(1) &= -2^{(3k^2-2)/2}\ ,\\
& & b_k(1) &= p_k(1)+d_k(1)^2=2^{k^2}\ ,
\\[2mm]
k &= 6j + 5 & d_k(1) &= -2^{(k^2-1)/2}\ ,\\
& & p_k(1) &= 3\cdot 2^{k^2-1}=3d_k(1)^2\ ,\\
& & r_k(1) &= 2^{(3k^2+1)/2}\ ,\\
& & b_k(1) &= p_k(1)+d_k(1)^2=2^{k^2}\ .
\end{aligned}
$$
The above shows how the arithmetic information is
covered in $d_k$, in case it is not vanishing,
then $p_k,r_k$ are constrained to be $p_k=3d_k^2$, and $(-1)^{k-1} 4d_k^3$.
Some simple sage code may be of interest, it computes the $B$-polynomials (and implicitly all data in between. Please adapt to see other wanted variables with corresponding prints.)
Similar pieces of code that compute "more" and for different purposes
are on my MSE answer(s).
def EA(A, field=QQ):
return EllipticCurve(field, [-3*A^4, -A^10 - A^2])
R0.<n> = PolynomialRing(ZZ)
F = FractionField(R0)
E = EA(n, field=F)
G = E.point([n^4 + n^2 + 1, n^6 + n^4 + n^2 + 1])
def B(k):
a = n
x = (k*G)[0]
u = a^3 - (a^9 - a)/(x + a^6)
b = u
# v = sqrt( (a^3 - b)*(b^3 - a) )
B = R0( (b/n).numerator() ).homogenize()
if b == n * B(n, 1) / B(1, n):
return B(n, 1)
# we may want to compute some B(k) for some relatively big k
print(f"B(13) = {latex(B(13))}")
And the above delivers:
$$
\tiny
B(13) = n^{168} + 643636n^{166} + 5416986934n^{164} - 6595169725244n^{162} + 20823931764302085n^{160} + 2875645846603778256n^{158} + 1667712145664558474856n^{156} + 159719934229962104174544n^{154} - 16418235274964244829099494n^{152} + 954887218043767314805499448n^{150} + 119731846710387314592686553492n^{148} + 3545291054165619315617657176152n^{146} + 121504163937978834746972642259210n^{144} + 3483488952052995289250579056814736n^{142} + 83193822695689099812768800490620328n^{140} + 963669940858723246664288414820239376n^{138} + 5691469353214542915269951938210105161n^{136} + 34636022418791988733051310409079244676n^{134} + 772541505470326563462542035514293522878n^{132} + 12825827388403891535299751323799386857108n^{130} + 138375990668894146984438432219206293217189n^{128} + 1133788329625647602603370172507874631269952n^{126} + 7592111666637440834997403231308573289580064n^{124} + 43764898682113299066393550990233916516372032n^{122} + 220209983529714884603307386851596893018428056n^{120} + 964750413127091863839044880196438110773385120n^{118} + 3687661350705896058636957906339919636823931888n^{116} + 12401212124649612962711080846137647902633076512n^{114} + 37170447403151095547819375350251357999929628184n^{112} + 100558209997856791867475940285188752412545861696n^{110} + 248161480118927171797911195972622236332506749088n^{108} + 563180918961626841685377680475438225636927001152n^{106} + 1182200905381161473183913287890740359064114659202n^{104} + 2303517567013892455193606912279403544735275537192n^{102} + 4176619501453304059258890215236080772082409113292n^{100} + 7054672921795819076752476136496891377159205078088n^{98} + 11110486539818301367775609188302841995113100887274n^{96} + 16320842530963850967065747412813453323987341073760n^{94} + 22369933935840320524337471821250807854472109461424n^{92} + 28613368251106838111952897091956658826075046771552n^{90} + 34161418782068102116676283742638185275584785620668n^{88} + 38067666977807856406875432807516827053194535910736n^{86} + 39591896044137661357489144863004740708845299012344n^{84} + 38416428184074426024353311084072898577233512604688n^{82} + 34758529974089937437794383641742859184577458327196n^{80} + 29295902751838084379673739856980382438415938704608n^{78} + 22971124728627050111639034597162544229637487183536n^{76} + 16724561176091393920903041897975104152397111166432n^{74} + 11277526309480866434237969013007318587784647334746n^{72} + 7019642047412867390482205533910867351330296121064n^{70} + 4015198416796069582405560580736721363005646871116n^{68} + 2098714389423006258602615266138135065936380182920n^{66} + 994739957392593453251893263108577583908969756626n^{64} + 423353217820373867605601641881099107708577721152n^{62} + 159535773297002729507424474343153022492340924064n^{60} + 52182270960305643285335423308984192212446749504n^{58} + 14348595244558616092247378416556523542802658648n^{56} + 3127138039859359323247812486544851441415556512n^{54} + 467690447726443685407311684343560305168336880n^{52} + 21997995433324589748997186648535586518194464n^{50} - 8114465233693039008333355772108118158268072n^{48} - 684223605949013017419151498549204480745664n^{46} + 1077917101896619148096733705656236255590432n^{44} + 570615057602475835936996367920813091645760n^{42} + 162025245456359277677940854988272484389277n^{40} + 30607069413717696641244864966871424783268n^{38} + 4174844335887760703273324515919511090174n^{36} + 436477806016997976391827054641725899252n^{34} + 35548978009085705306173155297537417201n^{32} + 1927492454500235987135304419497529424n^{30} + 61328816631688764307329243196871976n^{28} - 402516608559576880843779649143984n^{26} - 212640181426567207064206638104262n^{24} - 7095123653702930401996019826888n^{22} + 73964684163965477297514101652n^{20} + 3136114410853974725886870360n^{18} + 97009175247955676705825514n^{16} + 711325305233993038133520n^{14} + 5666753838034609840872n^{12} - 33680411373627173232n^{10} - 20879330408465859n^{8} + 107551370047636n^{6} + 84462003190n^{4} - 7758620n^{2} + \color{red}{169}\ .
$$
It is of course simple to compute the value in $1$ or in $i$ or in a primitive $8$.th root $u$.
Why these values? Since the discriminant of $E(A)$ is (up to factor)
$$
4\Big( -3A^4\Big)^3 + 27\Big(-A^2(A^8+1)\Big)^2
=27A^4(A^8-1)^2\ ,
$$
so specializing to a root of unity of order (dividing) $8$ leads to a degeneration.
Here is the list of the first few specializations:
u = QQbar( cos(2*pi/8) + i*sin(2*pi/8) )
R.<j> = QuadraticField(-1, latex_name='i')
for k in [2..9]:
Bk = B(k)
Bk1, Bki, Bku = Bk(1), ZZ(Bk(i)), R(Bk(u))
info = ' & '.join([latex(val.factor()) if val else str(val) for val in (Bk1, Bki, Bku)])
print(f"{k} & {info} \\\\\\hline")
In a table:
$$
\tiny
\begin{array}{|r|r|r|l|}
\hline
k & b_k(1) & b_k(i) & b_k(u)\\\hline
2 & 2^{4} & 2^{2} \cdot 3 & \left(i\right) \cdot (i + 1)^{2} \cdot 3 \\\hline
3 & 2^{8} & 2^{4} \cdot 3^{2} & (i + 1)^{4} \cdot 3^{2} \cdot (4 i + 5) \\\hline
4 & 2^{16} & -1 \cdot 2^{8} \cdot 3^{4} & \left(i\right) \cdot (i + 1)^{8} \cdot 3^{4} \\\hline
5 & 2^{24} & 2^{12} \cdot 3^{6} & \left(-i\right) \cdot (i - 4) \cdot (i + 1)^{12} \cdot 3^{6} \cdot (-8 i + 13) \\\hline
6 & 2^{36} & 2^{18} \cdot 3^{9} & \left(i\right) \cdot (i + 1)^{18} \cdot 3^{9} \\\hline
7 & 2^{48} & 2^{24} \cdot 3^{12} & (i + 1)^{24} \cdot 3^{12} \cdot (440 i + 441) \\\hline
8 & 2^{64} & -1 \cdot 2^{32} \cdot 3^{16} & \left(i\right) \cdot (i + 1)^{32} \cdot 3^{16} \\\hline
9 & 2^{80} & 2^{40} \cdot 3^{20} & (i + 1)^{40} \cdot 3^{20} \cdot (4360 i + 4361) \\\hline
\end{array}
$$
