Extract Features at Multiple Image-Scales
Not sure I have understood well your second question but I am gonna try to see if I can help you.
Question 1
Yes, the authors meant just that. I see were the confusion might come. So:

The authors say: &quot;[...] yields a field of 512-dimensional descriptor vectors&quot;
You say: &quot;[...] i get back 512 individual image-patches of size 14x14&quot;
You get: &quot;i get a feature-tensor of shape (1, 14, 14, 512)&quot;

So when the authors say a field of K-dimensional vector, with the word field they implicitly say that is a 2D vector (so $[H,W]$) of K dimensions, so $[H,W,K]$ which is exactly what your are getting (only with the batch dimension first, that Keras ignores in its input, that's why I prefer pytorch, it is more explicit).
So being explicit, from an input image $I:[B,H,W,C]=[1, 224, 224, 3]$ you get a field (or patch, or feature) vector, a 2D vector, of $K=512$ dimensions (or features, or descriptors) so $F:[B,H_f, W_f, K]=[1,14, 14, 512]$
What you are doing in your snippet is taking 1D vectors of $K=512$ dimensions, so you are extracting K-dimensional vectors from a K-dimensional feature map (or patch or field).
The flatten() goes in the second question.
Question 2
There are a lot of ways of doing feature-extraction at multiple scales. Say for example you extract features at 3 levels so you get:

$P_1:[H_1, W_1, C_1]=[56, 56, 100]$
$P_2:[H_2, W_2, C_2]=[26, 26, 200]$
$P_3:[H_3, W_3, C_3]=[13, 13, 300]$

This vectors are not concatenable. You indeed could flatten() them over the $[H,W]$ dimensions and concat over the channel dimension. But you would loose the 2D spatial relationship of images. What is normally done (from FPN, Retinanet EfficientNet, EfficientDet...) is the following:

You extract the previous feature maps at different scales
You process them over a FCN head of your network (so last dimension is bounding box coordinates or class probabilities)
You concatenate over last dimension

After extracting $P_1, P_2, P_3$ and passing by a FCN classification head on your network you would get (given you have 10 classes):

$P_1':[H_1, W_1, C_1]=[56, 56, 10]$
$P_2':[H_2, W_2, C_2]=[26, 26, 10]$
$P_3':[H_3, W_3, C_3]=[13, 13, 10]$

Now you have a classification vector per each &quot;pixel&quot; on each feature dimension. So now it makes more sense to flatten them to get 1D vectors of 10 class probs:

$P_{concat}=[N, 10]$ where $N$ is the flattened number of dimensions from $H,W$ of each scale feature map

If you want the implementation details go to RetinaNet in papers with code. I would link one very nice implementation but it is in pytorch, and I see you use Keras.