As a complement to other answers, I present here a very general method to see if a polynomial is irreducible of not. It needs a little more computation but it is not specific to the polynomial we consider.
There is an algorithm for absolute (i.e. over an algebraically closed field) factorisation of bivariate polynomials. It has been discovered by Picard (1906) and has been recently used by Ruppert (1999) and Gao (2003) to design efficient algorithms.
Let $f$ in $k[x,y]$ of bidegree $(m,n)$. Let assume that $\gcd(f, \partial_xf) = 1$, if it is not the case, we have an obvious factor.
Let $V$ be the vector space
$$ \left\{ (A,B)\in k[x,y]^2 \ \middle| \ \deg A \leq (m-1,n),\ \deg B \leq (m,n-1) \right \},$$
and $E$ the subspace of all $(A,B)$ in $V$ such that
$$ \frac{\partial}{\partial y}\left(\frac A f\right) = \frac{\partial}{\partial x}\left(\frac B f\right).$$
There is an obvious element in $E$: the couple $\left(\partial_x f,\partial_y f\right)$. More interestingly, if $g$ is a polynomial which divides $f$, then the couple $\left(\frac fg \partial_x g, \frac fg \partial_y g\right)$ is in $E$ as well. If $f = gh$, you can check that
$$ \left(\partial_x f,\partial_y f\right) = \left(\tfrac fg \partial_x g, \tfrac fg \partial_y g\right) + \left(\tfrac fh \partial_x h, \tfrac fh \partial_y h\right).$$
Now, the following fact should not surprise you : the dimension of E over $k$ is exactly the number of irreducible factors of $f$ over the algebraic closure of $k$.
So just with linear algebra over the rational numbers, I can prove that $xy+x^6+y^6$, as well as the two other polynomials you gave, are irreducible just by computing the rank of a matrix.
For the sake of completeness, here is a little piece of Maple code to perform the computation.
f := x*y+x^6+y^6:
# A and B are generic elements of V.
A := add(add(a[i,j]*x^i*y^j,i=0..degree(f,x)-1),j=0..degree(f,y)):
B := add(add(b[i,j]*x^i*y^j,i=0..degree(f,x)),j=0..degree(f,y)-1):
# This is the system of equation defining E.
collect(numer(diff(A/f,y)-diff(B/f,x)),[x,y],distributed):
eqs := [coeffs(%,[x,y])]:
# We compute a basis a generic element of E.
sol:=solve(eqs, indets(eqs)):
# The number of free variables in the generic
# element gives the dimension of E.
nops(indets(subs(sol, indets(eqs))));
This code give you the number of irreducible factors of $f$ !
- Picard, É. and Simart, G. (1906). Théorie des fonctions algébriques de deux variables indépendantes II
- Ruppert, W. M. (1999). Reducibility of polynomials $f(x,y)$ modulo $p$, J. Number Theory, 77(1), 62--70
- Gao, S. (2003). Factoring multivariate polynomials via partial differential equations, Math. Comp., 72(242), 801--822 (electronic)