JavaScript (ES7),  129 123  122 bytes

Saved 1 byte thanks to a suggestion from @Coder on my answer to a previous challenge

f=(x,y,n=0,m=[...s=1/0+''].map(_=>[...s]))=>m.map((r,Y)=>r.map((V,X)=>V<=n|((x-X)*(y-Y))**2^4||f(X,Y,n+1,m)),m[y][x]=n)&&m

Try it online!

Commented

f = (                // f is a recursive function taking:
  x, y,              //   the current position (x, y)
  n = 0,             //   a move counter n
  m = [...           //   the output matrix m[], initialized to an array
    s = 1 / 0 + ''   //   of 8 arrays made of the string "Infinity"
  ].map(_ => [...s]) //   i.e. ["I", "n", "f", "i", "n", "i", "t", "y"]
) =>                 //
m.map((r, Y) =>      // for each row r[] at position Y in m[]:
  r.map((V, X) =>    //   for each entry at position X in m[]:
    V <= n |         //     do nothing if V is less than or equal to n
    (                //     or if the squared product of
      (x - X) *      //       x - X and
      (y - Y)        //       y - Y
    ) ** 2 ^ 4 ||    //     is not equal to 4
    f(               //     otherwise, do a recursive call:
      X, Y,          //       pass the new position (X, Y)
      n + 1,         //       increment the move counter
      m              //       pass m[] unchanged
    )                //     end of recursive call
  ),                 //   end of inner map()
  m[y][x] = n        //   set the current square to n
) && m               // end of outer map(); return m[]

Charcoal, 59 58 bytes

⊞υ⟦⁰ＮＮ⟧Ｆυ«≔⊟ιη≔⊟ιθＪηθＩΣιＦ⁸Ｆ⁸«Ｊλκ¿›⁼⁵⁺Ｘ⁻λη²Ｘ⁻κθ²℅ＫＫ⊞υ⟦⊕Σικλ

Try it online! Link is to verbose version of code. Explanation:

⊞υ⟦⁰ＮＮ⟧

Start a breadth-first search by placing a 0 at the input coordinates.

Ｆυ«

Loop over the squares as they are discovered.

≔⊟ιη≔⊟ιθＪηθＩΣι

Set this square to the current number of steps.

Ｆ⁸Ｆ⁸«Ｊλκ

Loop over the whole board.

¿›⁼⁵⁺Ｘ⁻λη²Ｘ⁻κθ²℅ＫＫ

If this is an empty square that is a Euclidean distance of √5 from the current square, then...

⊞υ⟦⊕Σικλ

... save this square as being one more step away.

gbz80 machine code on Better Game Boy, 128 127 124 123 118 117 bytes

11 03 04 21 80 FF 06 08 0E 08 7C 22 0D 20 FC 3E
0D 22 3E 0A 22 05 20 F0 70 26 2F CD 7A 01 52 18
08 64 64 01 00 80 FF 00 00 76 24 3E 37 BC 28 43
3E 07 BA 38 3E BB 38 3B 3E 80 2E 0A 83 2D 20 FC
82 4F F2 BC 38 2D 7C E2 D5 14 1C 1C CD 7A 01 14
1D CD 7A 01 1D 1D CD 7A 01 15 1D CD 7A 01 15 15
CD 7A 01 15 1C CD 7A 01 1C 1C CD 7A 01 14 1C CD
7A 01 D1 25 C9

Eight more bytes can be saved from severe compromises to printed data, read the assembly source to learn more.

As a detail of the Game Boy platform, it must be placed at 0x150, and be jumped to from the entrypoint at 0x100. I opted to not count that because of it being standard GB dev practice (the entrypoint gives you only four bytes to work with).

I know that this answer sucks in terms of length (all those calls!), but I was overcome by the urge to do it for a problem after seeing someone else do x86 real-mode machine code elsewhere, and for this problem to get it to 0x80 bytes. I have no idea why I did this, I'm an ASM noob, but here we are.

BGB is specified because of the desire to print the answer. BGB has a debug print function that is only usable by its debugger (and that of no$gmb and a few others).

The initial position is encoded in the second and third byte and zero-indexed: second is y, third is x.

As for checking my answer, I suppose the only way is to have you download my ROM, or for you to assemble my assembly using rgbds. I have put both on GitHub. My assembly is incredibly undercommented.

Edit 1: Shortened corrective increments and decrements into a push and pop, saving one byte.

Edit 2: Moving various things around allowed saving 3 bytes, as well as much better pruning.

Edit 3: Remove decrement that was never necessary, saving 1 byte. That's embarrassing.

Edit 4: Removing a useless load and moving things around in the init loop to allow use of postincrement indexing saved 5 bytes, and the ability to remove 9 8 bytes with severe compromises was identified.

Edit 5: thejonymyster reminded me that I may assume valid input, moving one byte from compromises to saved.

Perl 5, 155 bytes

sub{$_=$"x64;($x,$y,$i,@_)=@_,"$x$y"=~/^[0-7]{2}$/&&s/(?<=^.{@{[$x+$y*8]}}) /0+$i/e&&push@_,map{/./;$x+$&-3,$y+$'-3,$i+1}21,41,12,52,25,45,14,54while@_;$_}

Try it online!

Somewhat ungolfed:

$f=

sub {
  $_ = $"x64;                      #init board string to 64 spaces
  ($x, $y, $i, @_) = @_            #move x, y and i from head of @_ param array
  ,
  "$x$y"=~/^[0-7]{2}$/             #if x,y within board
  &&                               #and
  s/(?<=^.{@{[$x+$y*8]}}) /0+$i/e  #space found+replaced by i+1 in current x,y
  &&                               #and
  push @_,                         #push possible jump to cells to todo-queue
    map {/./; $x+$&-3, $y+$'-3, $i+1 }
    21, 41, 12, 52, 25, 45, 14, 54 #knight pos cells digits x-3,y-3
        while @_;                  #...all that while there's more to do
  $_                               #return filled up 64 char board string
}

;print &$f( 1 , 0 ) =~ s/.{8}/$&\n/gr,"\n";

Output:

30323234
23212343
12143234
23232343
32323434
43434345
34343454
45454545

Python (without numpy), 142 135 bytes

lambda x,y,k=[9]*64,r=range(64):[k:=[(x*8+y!=i)*k[i]and-~min(k[j]for j in r if(i%8-j%8)**2^(i//8-j//8)**2==5)for i in r]for _ in r][-1]

Attempt This Online!

Takes in x and y coordinate as separate inputs, outputs a flattened square list. The code for finding if two positions are a knight's move from each other is taken from Arnauld's answer to "Where can the knight be in N moves?".

-7 bytes from @CursorCoercer by moving the code onto a single line & into a lambda

C (clang), 153 bytes

12 bytes thanks to ceilingcat!

e;r(*b,x,y,l,i){for(;i--&&x>0&y>0&x<9&y<9&l<9;b[e=x+y*8]>l?b[e]=l:0)e="XVTPJFDB"[i]-65,r(b,x+e%5-2,y+e/5-2,l+1,8);}f(*b,x,y){r(memset(b,1,324),x,y,0,8);}

Try it online!

Python Numpy, 126 bytes

lambda*z:H[z]
from numpy import*
x=mat(tri(8)).I
x-=x*x
x+=x.T
H=8-inner(*([kron(x,x)-eye(64)<0,1]**c_[:8]).T).A
H.shape=4*[8]

Attempt This Online!

How

The entire solution is precomputed, the function itself does no more than a simple lookup.

The finalised lookup table is 8x8x8x8, but while building it is reduced to 64x64.

We start by setting up the single step matrix and then take powers to cover sequences of moves.

1 1D considerations

The 1D case is easy but useful to build on. We can go 1 or 2 units up or down. If we represent that in a 2D source-target matrix x this means that the two diagonals just above the main diagonal must be set and the two below.

2 Go 2D

The two dimensions must be combined in an AND like way because the restriction 1 or 2 up or down holds along both axes. The Kronecker product does that, but there is the additional constraint that if it is 1 or 2 along one axis it must be the other along the other axis. This can be enforced by giving the 1. and 2. off diagonals different signs. After Kroneckering the valid bits will all be negative.

3. Construct 1D

There are many ways to construct x. The shortest I found was: Get the upper or lower triangle. Take the inverse. This is almost what we want only shifted one unit to the centre. This can be fixed by taking the difference with the square. Finally, we must add (or OR) the transpose.

4. Putting the bits together

After taking the Kronecker product we need to add the identity matrix so taking nth power gives the indicator of what's reachable in up to n steps. Conveniently, as the matrix at this point is boolean valued multiple paths to the same target are only counted as one. To get the right distance indices we must count in how many powers a target does not occur. Therefore, the type of the matrices must be changed to int after taking powers but before adding them together. We do this by taking the inner product with a vector of ones which forces coercion to int and then sums.

Previous Python Numpy, 131 bytes

lambda y,x:65-((D:=Z==y*8+x)+sum(D:=D|M@D for _ in Z)).reshape(8,8)
import numpy
Z=numpy.c_[:64];X=Z//8+Z%8*1j
M=abs((X-X.T)**2)==5

Attempt This Online!

Conceptually straight-forward approach. Form the full 64-by-64 indicator matrix of knight's moves (by checking for Euclidean distance \$=\sqrt 5\$). Repeatedly multiply by it the indicator vector of the starting position.

Python 3 with numpy, 138 136 133 bytes

lambda*a:(S:=sum(D:=-abs(mgrid[:8,:8].T-a).T,0))-(amin([*D//2,S//3,2//S],0)+S&-2)+equal(*D)*(S&-~amin(S)//6==-4)*2
from numpy import*

Attempt This Online!

Example, with arguments 1,2 (a=[1,2]):

mgrid[:8,:8]
[[[0 0 0 0 0 0 0 0]   [[0 1 2 3 4 5 6 7]
  [1 1 1 1 1 1 1 1]    [0 1 2 3 4 5 6 7]
  [2 2 2 2 2 2 2 2]    [0 1 2 3 4 5 6 7]
  [3 3 3 3 3 3 3 3]    [0 1 2 3 4 5 6 7]
  [4 4 4 4 4 4 4 4]    [0 1 2 3 4 5 6 7]
  [5 5 5 5 5 5 5 5]    [0 1 2 3 4 5 6 7]
  [6 6 6 6 6 6 6 6]    [0 1 2 3 4 5 6 7]
  [7 7 7 7 7 7 7 7]]   [0 1 2 3 4 5 6 7]]]

This produces two matrices that hold the co-ordinates.

D:=-abs(mgrid[:8,:8].T-a).T
[[[-1 -1 -1 -1 -1 -1 -1 -1]  [[-2 -1  0 -1 -2 -3 -4 -5]
  [ 0  0  0  0  0  0  0  0]   [-2 -1  0 -1 -2 -3 -4 -5]
  [-1 -1 -1 -1 -1 -1 -1 -1]   [-2 -1  0 -1 -2 -3 -4 -5]
  [-2 -2 -2 -2 -2 -2 -2 -2]   [-2 -1  0 -1 -2 -3 -4 -5]
  [-3 -3 -3 -3 -3 -3 -3 -3]   [-2 -1  0 -1 -2 -3 -4 -5]
  [-4 -4 -4 -4 -4 -4 -4 -4]   [-2 -1  0 -1 -2 -3 -4 -5]
  [-5 -5 -5 -5 -5 -5 -5 -5]   [-2 -1  0 -1 -2 -3 -4 -5]
  [-6 -6 -6 -6 -6 -6 -6 -6]]  [-2 -1  0 -1 -2 -3 -4 -5]]]

Transpose, subtract a, and transpose again: this subtracts the arguments from the two matrices, respectively. Also take the absolute value and negate.

This produces the negative vertical and horizontal distances from the given square.

S:=sum(D,0)
[[ -3  -2  -1  -2  -3  -4  -5  -6]
 [ -2  -1   0  -1  -2  -3  -4  -5]
 [ -3  -2  -1  -2  -3  -4  -5  -6]
 [ -4  -3  -2  -3  -4  -5  -6  -7]
 [ -5  -4  -3  -4  -5  -6  -7  -8]
 [ -6  -5  -4  -5  -6  -7  -8  -9]
 [ -7  -6  -5  -6  -7  -8  -9 -10]
 [ -8  -7  -6  -7  -8  -9 -10 -11]]

Add the two matrices for the negative taxicab distance from the given square.

The reason for the negation is that we really want rounding-up division, but only have easy access to rounding-down division with //.

\$\lfloor\frac{-n}{d}\rfloor = -\lceil\frac{n}{d}\rceil\$

O:=amin([*D//2,S//3,2//S],0)
[[-1 -1 -2 -1 -1 -2 -2 -3]
 [-1 -2  0 -2 -1 -2 -2 -3]
 [-1 -1 -2 -1 -1 -2 -2 -3]
 [-2 -1 -1 -1 -2 -2 -2 -3]
 [-2 -2 -2 -2 -2 -2 -3 -3]
 [-2 -2 -2 -2 -2 -3 -3 -3]
 [-3 -3 -3 -3 -3 -3 -3 -4]
 [-3 -3 -3 -3 -3 -3 -4 -4]]

The *D//2,S//3 part produces concentric octagons corresponding to the furthest points the knight can reach in a number of moves. 2//S is 0 with a RuntimeWarning for S=0, −2 for S=1, −1 for S≥2; it makes a difference only in the case S=1, changing the values of the squares adjacent to the given square from −1 to −2.

S-(O+S&-2)
[[1 2 3 2 1 2 3 4]
 [2 3 0 3 2 3 2 3]
 [1 2 3 2 1 2 3 4]
 [2 1 2 1 2 3 2 3]
 [3 2 3 2 3 2 3 4]
 [2 3 2 3 2 3 4 3]
 [3 4 3 4 3 4 3 4]
 [4 3 4 3 4 3 4 5]]

A correction for parity: the number of moves should have the same parity as the the taxicab distance S. Add S, round down to the nearest multiple of 2 (by zeroing the low bit), and subtract S; also negate and simplify.

Finally, some more adjustments:

+equal(*D)*(S&-~amin(S)//6==-4)*2

• The squares two steps away along a diagonal should have value 4 instead of 2.

• If the given square is a corner square, the diagonally adjacent square should also have value 4 instead of 2.

• equal(*D) checks for squares sharing a diagonal with the given square.

• amin(S), the negated highest taxicab distance, is −14 if the given square is a corner square, and otherwise ranges from −8 (for a centre square) to −13.

• -~amin(S) is one higher: −13 for corner squares, −7 to −12 otherwise.

• Floor-dividing by 6 produces −3 for corner squares and −2 otherwise.

• On the diagonals, S is even. ANDing with −2 = ...1110₂ leaves it unchanged, whereas ANDing with −3 = ...1101₂ zeroes the twos bit, mapping both −2 and −4 to −4.

The result:

[[1 2 3 2 1 2 3 4]
 [2 3 0 3 2 3 2 3]
 [1 2 3 2 1 2 3 4]
 [4 1 2 1 4 3 2 3]
 [3 2 3 2 3 2 3 4]
 [2 3 2 3 2 3 4 3]
 [3 4 3 4 3 4 3 4]
 [4 3 4 3 4 3 4 5]]

Knight, 272 264 193 190 188 bytes

;=nF;=xP;=yP;=sS*'9'64+*8y xT0;W>9=n+1n;=r~1W>8=r+1r;=c~1W>8=c+1c;=a~3W>4=a+1a;=b~3W>4=b+1b&?5+*a a*b b&&&&>8=x+a r<~1x&>8=y+b c<~1y>E Gs+x*8yT=qGs+r*8cT=sSs+x*8yT+1q;=i~8W>64=i+8iO Gs i 8

Input is 0-based horizontal and vertical coordinates, each on separate lines.

A fitting language for this challenge... or maybe not...

Try it Online!

Here's the Knight code translated to Python: Try It Online!

C (gcc), 191 179 176 bytes

-15 bytes thanks to @ceilingcat

f(a,x,y)int*a;{a[x=x*8+y]=!memset(a,1,256);r(a,x);}r(a,x,b,y)int*a;{for(b=0;++b<9;x/8==y>>3&a[x]+1<a[y-=8*~(~b&1)*~-(b&2)]&x>>6==y>>6&&r(a,y,a[y]=a[x]+1))y=x-~(b&1)*~-(b/2&2);}

Try it online!

K (ngn/k), 48 46 bytes

8 8#{x&1+&/'x@&'4=*/4#i-\:'i:!8 8}/9*~(!64)=8/

Try it online!

Jelly,  24 23  22 bytes

W8p¤ạṢ€ċؽʋƇ$ƬŒṬĖP€o/’

A monadic Link that accepts a pair of integers and yields a list of lists of integers.

Try it online! (The footer pretty prints the grid)

How?

Repeatedly finds coordinates reachable from the current list of coordinates starting with just the given coordinate. Then creates grids from each of these lists with their grid number at the identified coordinates and reduces these with logical OR to keep the first non-zero value at each coordinate. This has the initial position as \$1\$ rather than \$0\$, so all numbers are then reduced by one.

W8p¤ạṢ€ċؽʋƇ$ƬŒṬĖP€o/’ - Link: integers [x, y]
W                      - wrap -> [[x, y]]
             Ƭ         - collect up while distinct applying:
            $          -   last two links as a monad:
   ¤                   -     nilad followed by link(s) as a nilad:
 8                     -       eight
  p                    -       ([1..8]) Cartesian power ([1..8]) -> all coordinates
           Ƈ           -     keep those for which:
          ʋ            -       last four links as a dyad:
    ạ                  -         absolute difference (vectorises)
     Ṣ€                -         sort each
        ؽ             -         [1,2]
       ċ               -         count occurrences
              ŒṬ       - convert to grids of zeros with ones at the coordinates 
                Ė      - enumerate -> [[1, 1st grid], [2, 2nd grid], ...]
                 P€    - product of each -> grids with 1s replaced by grid number
                    /  - reduce by:
                   o   -   logical OR (vectorises)
                     ’ - decrement

Alternative 22, same method just identifying the existence of [1,2] and [2,1] by filtering for those which, when sorted, are invariant under the operation of getting the range of their length ([1,2]):

W8p¤ạṢJƑ$Ƈ¥Ƈ$ƬŒṬĖP€o/’
            

Retina, 121 bytes

.+
*#8*$(8*_ ¶
^(#)*((?<-1>.)*)_
$2$#1
s)+`(?<=(\d)(....|.{11}|.{13})?.{7})_|_(?=.{7}(....|.{11}|.{13})?(\d))
$.(_$1*_$4*

Try it online! Takes input as a pair of 0-indexed digits. Explanation:

s)`

Run the whole script in single-line mode where . also matches newlines.

.+
*#8*$(8*_ ¶

Treat the input as base 10 and insert that many #s followed by an empty chessboard of 8 rows of 8 _s terminated by spaces and newlines to bring each row up to 10 characters.

^(#)*((?<-1>.)*)_
$2$#1

Replace the nth character of the chessboard with a 0, removing n. ($#1 is just a weird way of saying 0; I can't use 0 literally because it will tokenise with the previous 2.)

+`

Repeat until no more changes can be made.

(?<=(\d)(....|.{11}|.{13})?.{7})_|_(?=.{7}(....|.{11}|.{13})?(\d))

Search for all empty squares that are a knight's move away from a square that already contains a digit (which will always be the largest digit so far).

$.(_$1*_$4*

Replace them with the next digit.

Wolfram Language (Mathematica), 55 bytes

GraphDistance[8~KnightTourGraph~8,#+8#2+1]~Partition~8&

Try it online!

Using the built-in KnightTourGraph function.

If the strict output format is required:

Wolfram Language (Mathematica), 69 bytes

-6 bytes and fixed a bug thanks to @att and @Steffan

StringRiffle[GraphDistance[8~KnightTourGraph~8,#+8#2+1]~Partition~8]&

Try it online!

MATL, 34 32 31 29 bytes

8t&Oli(9:"tt3Zv&+3=Z+g<@Q*+]q

Input is a cell array with the vertical and horizontal coordinates, starting at top left, 1-based.

Try it online!

How it works

8t&O      % Push 8×8 array of zeros
li(       % Write 1 at the position taken as input
9:"       % For each k in [1, 2, ..., 9]
  t       %   Duplicate. This is the current pattern, to be extended. A 0 in an
          %   entry means position not yet reached. A non-zero value indicates
          %   the required number of moves plus 1
  t       %   Duplicate again
  3Zv     %   Symmetric range: [1 2 3 2 1]
  &+      %   5x5 matrix of pair-wise additions
  3=      %   True for entries that equal 3. This gives a 5×5 matrix containing
          %   the 8 possible moves for a knight
  Z+      %   2-dimensional convolution of the two above matrices, maintaining
          %   the size of the former (8×8) in the output
  g       %   Convert the result to logical: non-zeros become 1
  <       %   Less than (element-wise)? This indicates *new* available positions
  @Q*     %   Multiply element-wise by k+1
  +       %   Add, element-wise. This updates the pattern with the new values
]         % End
q         % Subtract 1 element-wise. Implicit display