EDIT : What follows was written as an answer to the initial redaction of the question, where B was supposed to be nilpotent (i. e. B^k-B=\mathbf{0} for some prime k, with no specification of k). Furthermore, it uses an wrong definition of nilpotency (first-class thinko...). Therefore, it answers a different problem, with very weak relation to the real question. Kept for what it's worth...
Not an answer, but a contribution. Furthermore, I won't try to solve the (hard) problem of solving in \mathbb{Z}, but, at least as a first step, to solve it in \mathbb{Q}.
The problem has two parameters :
- The matrices' dimension n
- the index (degree) of B's nilpotency.
Each of the matrix equations translates to n^2 equations involving 2n^2 unknown equations ; similarly, the nilpotency of B translates to n^2 equations of maximal dregree k involving n^2 parameters.
However, these equations are far from being independent. For example, the nilpotency of B translates to
n=2\, k=2 : 4 equations of 4 unknowns ; the resulting ideal has dimension 2, which means that this system has two degrees of freedom.
n=3\\,k=2 ; 9 equations of 9 unknowns, the resulting system still has 4 degrees of freedom.
The "obvious" brute-force way (a polynomial system in a_{i,j}\,b_{i,j}) doesn't work. However, one may try to solve for B given A, i. e. build a system of polynomials in b_{i,j} whose coefficients are taken in (the fraction field of) the ring of polynomials in a_{i,j}. For example, after running :
NN=2
k=2
R1=PolynomialRing(QQ, ["a%d%d"%(u, v) for v in range(NN) for u in range(NN)])
# R1.inject_variables()
A=matrix(NN,R1.gens())
R2=PolynomialRing(FractionField(R1), ["b%d%d"%(u, v) for v in range(NN) for u in range(NN)])
# R2.inject_variables()
B=matrix(NN,R2.gens())
S1=(A^2*B+A*B*A+B*A^2).list()
S2=(B^2*A+B*A*B+A*B^2).list()
# Nilpotency
S3=(B^k-B).list()
J1=R2.ideal(S1+S2+S3)
one gets :
sage: R2.ideal(S1+S2+S3).variety()
[{b11: 0, b01: 0, b10: 0, b00: 0}]
sage: J1.dimension()
0
sage: J1.variety()
[{b11: 0, b01: 0, b10: 0, b00: 0}]
This result is easily repeated for k\in\,(2,3,5,7) ; for k=11, the system failed to compute J1.dimension()
after about 10 minutes (I threw the towel...).
Similarly, in the case n=3,\,k=2, the system has been searching the same J1.dimension()
for about 3 hours.
FWIW, the (gratis-but-not-free) Wolfram engine's Reduce
exhibits a similar behaviour., However, Solve
seems to be able to get the (unique) solution (null B) for n=3,\,k=3 in a reasonable time ; getting the (same) answer for n=4,\,k=2 needs about a minute. Going beyond n=4 or k=7 leads to unreasonable (=exceeding my patience) times.
This brute-force symbolic approach using only the polynomials machinery (or whatever the Wolfram engine uses) seems bound to fail. An approach using more results of linear algebra is needed.
HTH,
You don't specify the index (degree) of B (i. e. the smallest (prime) k for which Bk=B).
@Emmanuel Charpentier, I have updated. Thanks
If B2=0, the second equation reduces to simply BAB=0.
@Max Alekseyev, yeah, you are correct. But I have not been able to solve the two equations simultaneously. Clearly B=0 is a solution. But Does there exist a non zero matrix B such that B2=0 and satisfies the two matrix equation simultaneously.
The first equation translates into (A+B)3=A3 and so it may be work to look for B being the difference between two cubic roots from the same matrix.