Python, 84 bytes

def f(n):
 z,*y=n*[0,-1],
 while z:yield map([*range(2*n)].pop,y+z);p,*z,_=z;y+=-p,0

Attempt This Online!

Previously:

Python, 89 bytes

lambda n:[map([*range(2*n)].pop,([0,0,1,0]*n)[:2*j]+[j%-2,j%2-1]*(n-j))for j in range(n)]

Attempt This Online!

Returns flattened generators of nx2 snakes. Test code beautifies the output.

How?

Keeps a snake in memory and pops its bits into the right places as it traverses the grid. This works well because the pop indices are either 0,0,1,0,0,0,1,0,0,0,1,... (wiggly bit) or ...,0,-1,0,-1,0,... (long final loop).

Python 3, 87 bytes

lambda n:[[[k+[r:=x//2,n*2+~r][k+x&1],x^r%2][r<k]for x in range(n*2)]for k in range(n)]

Try it online!

Outputs a flattened list representing a two-column snake with n rows. The idea is to make a snake with k rows of "coils" and n-k rows for the "tail", for k each from 0 to n-1.

I tried a few ways to replace the inner list selection k+[r:=x//2,n*2+~r][k+x&1] that handles the snake's "tail" being in the first or second column. The best I got were 1 byte longer, but I feel like I'm missing something.

lambda n:[[[n+k+((r:=x//2)-n^(k+x)%-2),x^r%2][r<k]for x in range(n*2)]for k in range(n)]
lambda n:[[[k+(r:=x//2)^-(k+x&1),x^r%2][r<k]%(2*n)for x in range(n*2)]for k in range(n)]

K (ngn/k), 30 bytes

{?<'i{(y#x),|y_x}\x,&2!i:!2*x}

Try it online!

Generates flat snake matrices - first row concatenated with second row.

Snakes ara permutations of 2*N items. Observe their inverse permutations, for instance for N=5 we have:

(0 1 2 3 4 9 8 7 6 5
 0 5 6 7 8 9 4 3 2 1
 0 5 6 1 2 3 4 9 8 7
 0 5 6 1 2 7 8 9 4 3
 0 5 6 1 2 7 8 3 4 9)

To get from the first row to the second, we reverse everything starting from index 1 to the end. From the second to the third row - reverse everything from index 3. Then, index 5, etc, for all odd integers below 2*N.

It's easy to reverse a permutation in k - that's monadic <, also known as "grade".

This function:

{(y#x),|y_x}

accepts a permutation x and an index y, and returns x with reversed items after index y. It takes us from one snake to the next.

Generating the initial snake (0 1 2 3 4 9 8 7 6 5) is tricky. Its first half goes up and the second goes down. Luckily, we can reuse the function above, passing permutation i:0 1..(2N-1) and index N.

x,&2!i:!2*x builds a list of start indices for flipping, and i{..}\ goes through them using i (from the previous paragraph) as the initial value. The last iteration will create a duplicate result - we fix that by using distinct (monadic ?) at the end.

J, 46 bytes

i.((_2]`|.\,)@{.,0|:[|.(-@#]`|.\,)@}.)"{i.@,&2

Try it online!

Uses ovs's observation that each snake is made up of a two parts: a section of "switchbacks" where we go back and forth, followed by a single long "all the way to the end of the aisle and back" loop. Each section may be 0 length, in which case the other section is the entire thing.

Rather than using recursion, though, we just explicitly list all possible splits of the input. Eg, for n=3 all possible splits are:

┌───┬───┐
│   │0 1│
│   │2 3│
│   │4 5│
├───┼───┤
│0 1│2 3│
│   │4 5│
├───┼───┤
│0 1│4 5│
│2 3│   │
├───┼───┤
│0 1│   │
│2 3│   │
│4 5│   │
└───┴───┘

Then we apply the "switchback" and "loop" transformations to each piece:

• The "switchback" xform reverses odd rows.
• The "loop" xform flattens, breaks into two equal pieces, and reverses the 2nd piece, and transposes. It also reverses the columns as needed so the numbers match up with the first piece.

K (ngn/k), 33 bytes

Returns the snakes as a list of 2×N arrays.

2_{((0 1,'|2+)'x),,1(|y+)\!y}/!1+

Try it online!

All snakes have the same structure: On the left is a "wiggly" part, where they move up and down at every position, then the tail is "straight" (going to the very right, then turning around at the end).

The snakes can be constructed recursively: For a given size N, there is the fully "straight" snake:

1(|y+)\!y         / (0..N-1; 2*N-1..N)

And all snakes of size N-1 with an extra wiggle at the head:

(0 1,'|2+)'x
(        )'x      / for each snake of size N-1:
       2+         /   increment every position by 2,
      |           /   reverse the snake as the extra wiggle swaps top and bottom,
 0 1,'            /   prepend new wiggly head.

Jelly, 11 bytes

(brute-force - meant to do this earlier but got distracted golfing my Python port of the below 20 byte method!)

The idea to head the results, saving 2 bytes, is taken from DLosc's Pip -p answer.

2pŒ!ạƝ§$ÐṂḣ

A monadic Link that accepts \$n\$ and yields a list of the \$n\$ snakes, each as a list of 1-indexed coordinates in visiting order.

Try it online! (Gets slow, quicky!)

How?

2pŒ!ạƝ§$ÐṂḣ - Link: positive integer, N
2p          - two Cartesian product {[1..N]} -> AllCoordinates
  Œ!        - all permutations (starts with all those starting with `[1,1]` )
        ÐṂ  - keep those {P} minimal under:
       $    -   last two links as a monad - f(P):
     Ɲ      -     for neighbouring pairs of coordinates in {P}:
    ạ       -       absolute difference (vectorises) -> [dr, dc]
      §     -     sums -> ManhattanDistances
          ḣ - head to index {N}

Jelly, 20 bytes

(non-brute-force)

ḤRW;ḶḤo‘Ʋ’ḣ;Ṛ{Qʋ\ḊỤ€

A monadic Link that accepts \$n\$ and yields a list of the \$n\$ snakes, each as a flattened list of 1-indexed snake locations.

Try it online! Or see the test-suite (converts to 0-indexing to align with the question text).

How?

ḤRW;ḶḤo‘Ʋ’ḣ;Ṛ{Qʋ\ḊỤ€ - Link: positive integer, N
Ḥ                    - double {N} -> 2N
 R                   - range {that} -> [1..2N]
  W                  - wrap {that} -> [[1..2N]]
        Ʋ            - last four links as a monad - f(N):
    Ḷ                -   lowered range {N} -> [0..N-1]
     Ḥ               -   double {that} -> [0,2,4,...,2N-2]
       ‘             -   increment {N} -> N+1
      o              -   {doubles} logical OR {N+1} -> [N+1,2,4,...,2N-2]
   ;                 - {wrapped} concatenate {that} -> [[1..2N],N+1,2,4,...,2N-2]]
         ’           - decrement {that} -> [[0..2N-1],N,1,3,...,2N-3]]
                \    - cumulative reduce {that} with:
               ʋ     -   last four links as a dyad - f(Current, V):
          ḣ          -     head {Current} to index {V}
            Ṛ{       -     reverse {Current}
           ;         -     {head} concatenate {reversed}
              Q      -     deduplicate {that}
                         -> i.e. f(a-z, 3) -> abczyx...d
                 Ḋ   - dequeue {that} (to remove the [0..2N-1])
                   € - for each {of those}:
                  Ụ  -   grade up

05AB1E, 19 bytes

·Lœʒн}ʒ.ā{€θI‰üαOP

Inspired by DLosc's Pip answer, so make sure to upvote that answer as well.

Outputs a list of the 1-based flattened \$2\times n\$ snakes.
Brute-force, and times out for \$n\geq5\$.

Try it online or verify the first few smaller test cases.

Explanation:

·               # Double the (implicit) input-integer
 L              # Pop and push a list in the range [1,2*n]
  œ             # Pop and get all permutations of this list
   ʒн}          # Only keep the permutations that start with a 1:
   ʒ            #  Filter it by:
    н           #   Pop and keep the first item
                #   (only 1 is truthy in 05AB1E)
     }          #  Close the filter
   ʒ            # Filter further by:
    .ā          #  Enumerate: pair each value with its 0-based index
      {         #  Sort the pairs based on the first values
       €        #  Map over each pair:
        θ       #   Leave just the last item of the pair (the index)
         I‰     #  Divmod those indices by the input-integer to get 2D-indices
           ü    #  For each overlapping pair of 2D-indices:
            α   #   Get the absolute difference of both the x and y
             O  #  Sum each inner pair
              P #  Take the product of all those sums
                #  (again, only 1 is truthy in 05AB1E)
                # (after which the resulting list is output implicitly)

Charcoal, 57 36 bytes

ＮθＥθ⭆¹Ｅ²Ｅθ⎇‹πι⁺⊗π﹪⁺νπ²⁺ι⎇﹪⁺ιν²⁻⊖⊗θππ

Try it online! Link is to verbose version of code. Pretty-prints each pair of lists on its own line. Explanation: Directly calculates which segment of snake appears in each square of each grid.

Ｎθ                                      Input `N` as a number
   θ                                    Input `N`
  Ｅ                                     Map over implicit range
      Ｅ²Ｅθ                              Map over the grid
            π                           Current column
           ‹                            Is less than
             ι                          Current snake
          ⎇                             If true then
                   ν                    Current row
                  ⁺                     Plus
                    π                   Current column
                 ﹪                      Modulo
                     ²                  Literal integer `2`
              ⁺                         Plus
                π                       Current column
               ⊗                        Doubled
                       ι                Otherwise current snake
                      ⁺                 Plus
                           ι            Current snake
                          ⁺             Plus
                            ν           Current row
                         ﹪              Modulo
                             ²          Literal integer `2`
                        ⎇               If nonzero then
                                 θ      Input `N`
                                ⊗       Doubled
                               ⊖        Decremented
                              ⁻         Minus
                                  π     Current column
                                   π    Else current column
    ⭆¹                                   Pretty-print

Previous 57-byte solution:

ＮθＦθ«≔Ｅ²⟦⟧η≔ηζＦι«Ｆζ⊞λＬΣζ≔⮌ζζ»Ｆ⊗⁻θι⊞§ζ⁰ＬΣζＦ⁻θι⊞§ζ¹⊟§ζ⁰⟦⭆¹η

Attempt This Online! Link is to verbose version of code. Pretty-prints each pair of lists on its own line. Explanation:

Ｎθ

Input N.

Ｆθ«

Loop over each snake.

≔Ｅ²⟦⟧η

Start with two empty list of snake indices.

≔ηζ

Start the snake in the first list.

Ｆι«

For the ith snake, make it wiggle down and up i times.

Ｆζ⊞λＬΣζ

Add the next two indices, one to each list.

≔⮌ζζ

Reverse the order of the lists, so that the snake wiggles on each loop.

»Ｆ⊗⁻θι⊞§ζ⁰ＬΣζ

Add the rest of the snake to the first list for now.

Ｆ⁻θι⊞§ζ¹⊟§ζ⁰

Remove the excess snake from the first list and add it to the second list, making it double back on itself.

⟦⭆¹η

Output the lists in their original order.

Python 3.8, 90 bytes

def f(n):x=r=[*range(2*n)];return[map((x:=x[:i]+x[:i-1:-1]).index,r)for i in[n]+r[1:-1:2]]

A function that accepts n and returns a list of snake generators, each of which can produce a flattened, 0-indexed, snake.

Try it online!

How?

Port of my Jelly answer.

Start with [0..2N-1] and repeatedly take the first i elements followed by the rest in reverse for i in [n,1,3,5,...,2n-3] to find the next ranking of snake locations.

To get the snake locations from a ranking we map across [0..2N-1] using the index function of the ranking.

Pip -P, 28 bytes

1=M{aMP$+:BAD_}FIPMFL2CGaH:a

Brute-force approach; the time complexity is pretty BAD. Any input higher than 3 times out. Attempt This Online!

Explanation

1=M{aMP$+:BAD_}FIPMFL2CGaH:a
                     2CGa     ; 2-by-a coordinate grid
                   FL         ; Flatten into a list of 2*a coordinate pairs
                 PM           ; Get all permutations
               FI             ; Filter by this function:
   {a         }               ;  To a given permutation
     MP                       ;  map this function to all consecutive pairs:
             _                ;   First pair
          B                   ;   Second pair
           AD                 ;   Absolute difference (element-wise)
       $+:                    ;   Sum
                              ;   This will be 1 if the coordinate pairs are
                              ;   adjacent, greater than 1 otherwise
  M                           ;  Max of all results
1=                            ;  is 1 (i.e. all consecutive pairs are adjacent)
                              ; PM includes permutations that don't start at
                              ; [0;0], but the ones that do are at the beginning,
                              ; so to select them, we simply...
                         H:a  ; Take the first a elements

JavaScript (V8), 94 bytes

Prints one grid per line in flattened transposed format, with the snakes starting at 1.

f=(w,i=0,n=0,...r)=>i<0|i/w/2|r[i]||++n<w*2|[-2,0,2].map(d=>f(w,i+d^!d,r[i]=n,...r))||print(r)

Try it online!

Commented

f = (         // f is a recursive function taking:
  w,          //   w = width
  i = 0,      //   i = position (as y + x * 2)
  n = 0,      //   n = current cell (1-indexed)
  ...r        //   r[] = flattened transposed grid
) =>          //
i < 0 |       // abort if i is negative,
i / w / 2 |   // or i = w*2,
r[i] ||       // or r[i] is not zero'ish
++n < w * 2 | // increment n and test if it's less than w*2
[-2, 0, 2]    // list of directions
.map(d =>     // for each direction d:
  f(          //   recursive call:
    w,        //     pass w unchanged
    i + d     //     column update: add d to i
    ^ !d,     //     change the row if d=0
    r[i] = n, //     pass n and update r[i]
    ...r      //     pass a copy of r[]
  )           //   end of recursive call
) ||          // end of map()
print(r)      // if n < w*2 was false, print the grid

JavaScript (Node.js), 92 bytes

n=>[...Array(n)].map((_,i,A)=>A.map((_,j)=>[(t=[j<i?++k:n*2+i-j,++k])[j=j+k&1],t[j^1]],k=0))

Try it online!

n=>[...Array(n)].map((_,i,A)=>A.map((_,j)=>[(t=[j<i?k++:--y,k++])[j=j+k&1],t[j^1]],k=0,y=n*2))

Python, 100 bytes

f=lambda n,k=0:(m:=range(k,n))and[[m,range(n,n+n-k)[::-1]]]+[[[k,*j],[k+1,*i]]for i,j in f(n+1,k+2)]

Attempt This Online!

Returns the snakes as a list of 2×N matrices.

Python, 108 bytes

f=lambda n,k=0,t=1,a=[]:n*[1]and[a+[[k+i,k+2*n-i-1][::t]for i in range(n)]]+f(n-1,2+k,t^-2,a+[[k,k+1][::t]])

Attempt This Online!

Returns the snakes as a list of Nx2 matrices.

Python3, 177 bytes

def f(n):
 q=[(0,0,[(0,0)])]
 for x,y,p in q:
  if len(p)==2*n:yield p;continue
  q+=[(*T,p+[T])for X,Y in[(0,1),(0,-1),(1,0),(-1,0)]if n>y+Y>=0<=x+X<2and(T:=(x+X,y+Y))not in p]

Try it online!