I will check up almost all numbers up to 100 digits (except some numbers) that there is no other repdigit which is the area of a Pythagorean Triangle except $666666.$ I will use the CAS Maxima. Its documentation can be found here. I will also explain why I do not expect that there are more numbers with this property.
Let us summarize the facts already stated in your post and in comments:
the cathetus of all Pythagorean Triangle can be generated by
$$f(m^2-n^2),\, 2fmn$$
where $m,n,f \in \mathbb{Z^+},$
the corresponding hypothenuse is $f(m^2+n^2)$, the corresponding area $A$ is $$A=f^2(m^2-n^2)mn$$
We are looking for numbers $m,n,f,t \in \mathbb{Z^+}, n<m, d \in \{1,2,3,4,5,6,7,8,9\}$ such that
$$\frac{10^t-1}{9}d=f(m-n)(m+n)mn$$
At least one of $n$, $m$ of $m+n$ is even, so $2|A$ and therefore $d \notin \{1,3,5,7,9\}.$
Either $n$ or one of $(m-n),m,(m+n)$ is divisible by $3$, so $3|A$. Therefore $d=6$ or $d \in \{2,4,8\}$ and $3|t$. $t$ is the number of places of the repdigit.
If we check the area $A$ we can see that $A \bmod 10$ is in $\{0,4,6\}$ That can easily be done by Maxima by iterating over all possible values for $m,n,f \mod {10}$
(%i) unique(
create_list(mod(m*n*(m^2-n^2),M),
f,0,M-1,
m,0,M-1,
n,0,M-1)),M=10;
(%o) [0,4,6]
But repdigit numbers with digit $0$ we can ignore, so we have $d=6$ or $d=4$ and $3|t$
There was another property that I discovered by accident:
(%i) unique(
create_list(mod(m*n*(m^2-n^2),M),
f,0,M-1,
m,0,M-1,
n,0,M-1)),M=17;
(%o) [0,1,2,4,6,7,8,9,10,11,13,15,16]
So $A \bmod 17$ is never in $\{3,5,12,14\}$
This has the consequence that for repdigit numbers with digit $4$ the numbers where the number of digits cannot be $9,13,14,15 \mod 16$ and for repdigit numbers with digit $6$ the number of digits cannot be $3,7,10 \bmod 16$. But we will not use this in the following investigations.
We can also estimate the size of the factors of $A$:
$$
1 \le n \le cm \\
m <=m+n \le m(1+c) \\
m(1-c) \le m-n <=m \\
\implies \\
(1-c)m^3\le A \le c(1+c)m^4
$$
$$
cm \le n \le m \\
(1+c)m \le m+n \le 2m \\
1 \le m-n \le (1-c)m \\\
\implies \\
c(1+c)m^3 \le A \le 2(1-c)m^4
$$
We have $$\min\{c(1+c),(1-c)\} m^3 \le A \le \max\{c(1+c),2(1-c)\} m^4$$
for arbitrary $c$. We choose $c=\sqrt 2-1$ on the left side and $c=\frac 3 2 $ on the right side to get
$$(2-\sqrt 2) m^3 \le A \le \frac 3 2 m^4$$
We also know that not all of the four factors $n,m-n,m,m+n$ of $A$ can be smaller than $\sqrt[4]A$ because otherwise their product would be less than $A$. The smallest of these four numbers is either $n$ or $m-n$, the largest is $m+n$, the second largest is $m$. So
$$\min\{n,m-n\}\le \sqrt[4]A$$
Now how can we check if a given repdigit number $N$ satisfies the equation $N=nm(m-n)(m+n)$? We can check for each factor $b$ of $N$ if there is an appropriate $m$. This means, we solve for each factor $b$ of $N$ the equation
$$x^3b-xb^3=N$$
If this equation has an integer solution then we have found a Pythagorean Triangle generated by $(m, n)=(x, b). $
We have ignored that a factor $f^2$ can still be part of the area formula. So we have to repeat this procedure for all equations
$$x^3n-xn^3=\frac N {f^2}$$
where $f^2$ is a square factor of $N$. But if we do not find an integer solutions for these equations than $N$ ist not the area of a Pythagorean Triangle.
Actually we do not have to test all $n$ that divide $N$. From the inequalities above we know that $n<frac 1 {2-\sqrt(2)}\sqrt[3]A$. But in my implementation I use the fact that $\min\{n,m-n\}\le \sqrt[4]A$ So we have to solve the equation $x^3b-xb^3=N$ for the divisor $b$ of $N$ or we have $b=m-n$, then we substitute $n$ by $m-b$ in $m^3n-mn^3=N$ and get the equation
$$2bx^3-3b^2x^2+b^3x=N$$
Here a integer solution gives the tuple $(m, n)=(x, x-b)$
Here is a sketch of the algorithm in a pseudo programming language. The function squarepart(n) retrieves the largest square number dividing an integer
R is a repdigit number
for f in divisors of squarepart(R)^(1/2):
for d in divisors of R/f^2 less than (R/f^2 )^(1/4):
find root of x^3*d-x*d^3=R/f^2
if the root is a positive integer then
m:=x
n:=d
(f,m,n) is a solution
find root of 2*x^3*d-3*d^2*x^2+d^3*x=R/f^2
if root x is a positive integer then
m:=x
n:=m-d
(f,m,n) is a solution
Maxima has the function divisors(n)? that returns the list of divisors for a number. And it has the function $realroots(f,bound)$ that retrieves all real roots of a polynomial. If bound is less than $1$ and the $root$ is an integer than Maxima returns this integer.
One can check with the function integerp(number) if a number is an integer
My procedure in Maxima looks like this
get_pythbase4(R):=block([N, m, found:[], programmode:true,
result],
for f in divisors(apply("*",map(lambda([w],w[1]^floor(w[2]/2)),ifactors(R)))) do (
N:R/f^2,
for d in sublist(listify(divisors(N)),lambda([t],t<=isqrt(isqrt(N)+1)+1)) do
for equation in [x^3*d-x*d^3=N,2*x^3*d-3*d^2*x^2+d^3*x=N] do (
for solution in realroots(equation,1/2) do (
m:rhs(solution),
if integerp(m) then (
if (equation=(x^3*d-x*d^3=N)) then
result:['m=m,'n=d, 'f=f]
else
result:['m=m,'m-'n=d, 'f=f],
found:cons(result,found)
)
))
),
return(found)
);
The first for loop loops over all $f$ where $f^2$ is a divisor of the repdigit number $R$. The second one loops all divisors of $\frac R{f^2}$. The third loops constructs from d and N both equations we mentioned above. The fourth loop loops over the at most three roots of equations. The if-clause checks if the real solution which always exists, is an integer solution. If so, the result will be displayed.
we get the following result:
(%i) get_pythbase4(666666);
(%o) [[m = 37,n = 26,f = 1],[m = 37,m-n = 11,f = 1]]
for all other numbers I tried, e.g. ((10^60-1)/9*6) I got
(%i) get_pythbase4((10^60-1)/9*6);
(%o) []
This last calculation took about 50 seconds on my notebook.
Finally I executed the following function
run_loop3():= block([all_solutions:[],B,res],
for t:3 thru 100 do (
if (not member(t, [ 71,82,83,84,89,90,93,94,95,96,97,98 ])) then (
B:(10^t-1)/9,
print('t=t, 'M=6),
res:get_pythbase4(B*6),
if (not emptyp(res)) then (
print([B*6,res])
),
if (mod(t,3)=0) then (
print('t=t, 'M=4),
res:get_pythbase4(B*4),
if not emptyp(res) then (
print([B*4,res])
)
)
)
)
);
This loops over all repdigit numbers of length less or equal $100$ except that of length $\{71,82,83,84,89,90,93,94,95,96,97,98\}$ and calls the function get_pythbase4 for that number.
This gave me the following output
(%i) run_loop3();
t = 3 M = 6
t = 3 M = 4
t = 4 M = 6
t = 5 M = 6
t = 6 M = 6
[666666, [[m = 37, n = 26, f = 1], [m = 37, m - n = 11, f = 1]]]
t = 6 M = 4
t = 7 M = 6
t = 8 M = 6
t = 9 M = 6
...
t = 91 M = 6
t = 92 M = 6
t = 99 M = 6
t = 99 M = 4
t = 100 M = 6
So only for $666666$ is displayed a solution.
Note 1: The run_loop function that I actually executed still contained checks for M=2 and M=8, as I only realized later that this was not necessary. Therefore my output is a little bit "beautified".
Note 2:
The reason why I excluded some numbers for the loop was that Maxima was not able to factorize them in reasonable time or, for 84,90, 96, that the number of different prime factors was larger than 20 and therefore the numbers of divisors was larger than $10^6$ which would take a long time to process.
Note 3: The procedure factorizes the same number B more for each f and for each M. That is not very performant, but this makes the source code simpler and I have no interest in optimizing this procedure.
This value converges to $0$ as $N$ goes to infinity.
Note 4: The prime factorization of the missing number can be calculated at Dario Alpern's Web site. My algorithms can be modified that it uses the prime factorization supplied by this page to find the missing results.
I want to (over-)estimate the number of Pythagorean Triangles with an area less or equal $N$. So how many tuples $(f,m,n)$ are there such that the area they define is smaller $N$. We have
$$\sum_{f=1}^{\sqrt N} \sum_{m=1}^{\sqrt[3]{\frac N {f^2}}}\sum_{n=1}^{m}1 \\
=\sum_{f=1}^{\sqrt N} \sum_{m=1}^{\sqrt[3]{\frac N {f^2}}} m \\
\le \sum_{f=1}^{\sqrt N} \frac12 (\sqrt[3]{\frac N {f^2}} )^2 \\
= \frac12 N^{\frac 23} \sum_{f=1}^{\sqrt N} f^{-\frac 43}\\
\le\frac12 N^{\frac 23} (1+\int_1^{\sqrt N}f^{-\frac 43}df)\\
\frac12 N^{\frac 23} ( 1-(\sqrt N)^{-\frac13}+1)\le N^{\frac 23}$$
So not more than $N^{\frac 23}$ number of the $N$ numbers les or equal $N$ are the area of a Pythagorean Triangle. $N$ has about $\log_{10}N$ digits. Now assume we draw randomly $\log_{10}N$ numbers from these $N$ numbers. Then we expect $$\left(\frac {N^{\frac 23}}N \cdot \frac{\log_{10}N}N \right)N= \frac{\log_{10}N}{N^{\frac 23}}$$ of the drawn numbers are equal to the area of a Pythagorean Triple. This value converges to $0$ if $N$ goes to infinity. So if the set of repdigit numbers less or equal $N$ behaves like a set of the same size of random numbers then I would not expect to find more repdigit numbers that are equal to the area of a Pythagorean Triangle. But of course it is possible that there is a relation between repdigit numbers and the Pythagorean Triangles and they do not behave like a set of random number.
Bonus: Another digit pattern popular in recreational mathematics are palindromes. To find palindrome Pythagorean Triangles is much easier. Here is a small Python program that checks if Pythagorean Triangle with $1\le n<m\le 1000, f=1$ if its area is a palindrome number
for m in range(1,1000):
for n in range(1,i):
A=m*n*(m**2-n**2)
s=str(A)
if (s==s[::-1]):
print(m,n,A)
Its output is
2 1 6
21 8 63336
37 26 666666
58 26 4053504
62 49 4383834
78 1 474474
103 51 42066024
103 101 4244424
126 22 42666624
158 11 43177134
236 208 610262016
407 101 6390000936
$666666$ is in this list because it is a palindrome, too. An we see there are much more palindromes than repunit numbers that meet our requirement. The number of palindromes less than $N$ is $\sqrt N$, because we can choose the digits of the left halve of a digit an then get the right part my mirroring the left part. The expected number that the elements of a random set of size $\sqrt N$ are Pythagorean Triangle area is therefore $N^{\frac 16}$. So I expect there are Palindrome Numbers of size $N$ that are equal to the area of Pythagorean Triangles even for large $N$. But nevertheless they may be hard to find. But actually I only calculated an upper bound of the number of Pythagorean Triangles and therefore of the expected numbers.
