Find all integer pairs that produce a given Loeschian number Inspired by and drawns from Is this number Loeschian?

A positive integer \$k\$ is a Loeschian number if

\$k\$ can be expressed as \$i^2 + j^2 + i\times j\$ for \$i\$, \$j\$ integers.

For example, the first positive Loeschian numbers are: \$1\$ (\$i=1, j=0\$); \$3\$ (\$i=j=1\$); \$4\$ (\$i=2, j=0\$); \$7\$ (\$i=2, j=1\$); \$9\$ (\$i=-3, j=3\$)1; ... Note that \$i, j\$ for a given \$k\$ are not unique. For example, \$9\$ can also be generated with \$i=3, j=0\$.
Other equivalent characterizations of these numbers are:

\$k\$ can be expressed as \$i^2 + j^2 + i\times j\$ for \$i, j\$ non-negative integers. (For each pair of integers \$i, j\$ there's a pair of non-negative integers that gives the same \$k\$)

There is a set of \$k\$ contiguous hexagons that forms a tesselation on a hexagonal grid (see illustrations for \$k = 4\$ and for \$k = 7\$). (Because of this property, these numbers find application in mobile cellular communication networks.)

See more characterizations in the OEIS page of the sequence.



The first few Loeschian numbers are
0, 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 27, 28, 31, 36, 37, 39, 43, 48, 49, 52, 57, 61, 63, 64, 67, 73, 75, 76, 79, 81, 84, 91, 93, 97, 100, 103, 108, 109, 111, 112, 117, 121, 124, 127, 129, 133, 139, 144, 147, 148, 151, 156, 157, 163, 169, 171, 172, 175, 181, 183, 189, 192...

1while (\$i=-3, j=3\$) produces 9, stick to non-negative integers, so (\$i=0, j=3\$).
Loeschian numbers also appear in determining if a coincident point in a pair of rotated hexagonal lattices is closest to the origin?
The challenge
Given a non-negative integer \$k\$, output all pairs of non-negative integers \$i, j\$ such that \$i^2 + j^2 + i\times j=k\$. If none are found (i.e. \$k\$ is not Loeschian) then return nothing or some suitable flag other than \$(0, 0)\$ since that produces the first Loeschian number, \$0\$.
For reversed order pairs like \$(0, 4)\$ and \$(4, 0)\$ either include both, or one member of the pair, but it should be the same for all cases (i.e. not sometimes one and other times both).
The program or function should handle (say in less than a minute) inputs up to \$100,000\$, or up to data type limitations.
This is code golf so shortest code wins.
Test cases
 in       out
 0      (0, 0)
 1      (0, 1), (1, 0)
 3      (1, 1)
 4      (0, 2), (2, 0)
 9      (0, 3), (3, 0)
 12     (2, 2)
 16     (0, 4), (4, 0)
 27     (3, 3)
 49     (0, 7), (3, 5), (5, 3), (7, 0)
 147    (2, 11), (7, 7), (11, 2)
 169    (0, 13), (7, 8), (8, 7), (13, 0)
 196    (0, 14), (6, 10), (10, 6), (14, 0)
 361    (0, 19), (5, 16), (16, 5), (19, 0)
 507    (1, 22), (13, 13), (22, 1)
 2028   (2, 44), (26, 26), (44, 2)
 8281   (0, 91), (11, 85), (19, 80), (39, 65), (49, 56), (56, 49), (65, 39), (80, 19), (85, 11), (91, 0)
 12103  (2, 109), (21, 98), (27, 94), (34, 89), (49, 77), (61, 66), (66, 61), (77, 49), (89, 34), (94, 27), (98, 21), (109, 2)
05AB1E, 11 9 bytes
Thanks to @ovs for -2!
ÝãʒãÀ¦POQ

Oh my, I am getting in the hang of code golfing!
Prints a list of all the valid pairs (e.g [[1, 0], [0, 1]]). If there are none, the list is empty ([]). Also outputs both of any reverse integer pairs.
Try it online!
How?
You may count this as a port of the other answers, but I took a look only at the Husk answer before writing the program!
Ý             # Push a list of all numbers from 0 to the input.
 ã            # Push the cartesian power of lists. (Basically, finding all possible pairs)
  ʒ           # For each pair...
   ãÀ¦        # Find all other permutations of the pair.
      P       # Multiply each permutation.
       O      # Add the products.
        Q     # If the result is not equal to the input, yeet (throw) them from the list.
              # Automatically print the pairs not yeeted.