perft-function of chess engine is giving self-contradictory output

I am currently developing a chess engine in C++, and I am in the process of debugging my move generator. For this purpose, I wrote a simple perft() function:

int32_t Engine::perft(GameState game_state, int32_t depth)
{
    int32_t last_move_nodes = 0;
    int32_t all_nodes = 0;

    Timer timer;
    timer.start();

    int32_t output_depth = depth;

    if (depth == 0)
    {
        return 1;
    }

    std::vector<Move> legal_moves = generator.generate_legal_moves(game_state);

    for (Move move : legal_moves)
    {
        game_state.make_move(move);

        last_move_nodes = perft_no_print(game_state, depth - 1);
        all_nodes += last_move_nodes;

        std::cout << index_to_square_name(move.get_from_index()) << index_to_square_name(move.get_to_index()) << ": " << last_move_nodes << "\n";

        game_state.unmake_move(move);
    }

    std::cout << "\nDepth: " << output_depth << "\nTotal nodes: " << all_nodes << "\nTotal time: " << timer.get_milliseconds() << "ms/" << timer.get_milliseconds()/1000.0f << "s\n\n";

    return all_nodes;
}

int32_t Engine::perft_no_print(GameState game_state, int32_t depth)
{
    int32_t nodes = 0;

    if (depth == 0)
    {
        return 1;
    }

    std::vector<Move> legal_moves = generator.generate_legal_moves(game_state);

    for (Move move : legal_moves)
    {
        game_state.make_move(move);

        nodes += perft_no_print(game_state, depth - 1);

        game_state.unmake_move(move);
    }

    return nodes;
}

It's results for the initial chess position (FEN: rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - 0 1) for depths 1 and 2 match the results of stockfish's perft command, so I assume they are correct:

h2h3: 1
h2h4: 1
g2g3: 1
g2g4: 1
f2f3: 1
f2f4: 1
e2e3: 1
e2e4: 1
d2d3: 1
d2d4: 1
c2c3: 1
c2c4: 1
b2b3: 1
b2b4: 1
a2a3: 1
a2a4: 1
g1h3: 1
g1f3: 1
b1c3: 1
b1a3: 1

Depth: 1
Total nodes: 20
Total time: 1ms/0.001s

h2h3: 20
h2h4: 20
g2g3: 20
g2g4: 20
f2f3: 20
f2f4: 20
e2e3: 20
e2e4: 20
d2d3: 20
d2d4: 20
c2c3: 20
c2c4: 20
b2b3: 20
b2b4: 20
a2a3: 20
a2a4: 20
g1h3: 20
g1f3: 20
b1c3: 20
b1a3: 20

Depth: 2
Total nodes: 400
Total time: 1ms/0.001s

The results stop matching at depth 3, though:

Stockfish:

go perft 3
a2a3: 380
b2b3: 420
c2c3: 420
d2d3: 539
e2e3: 599
f2f3: 380
g2g3: 420
h2h3: 380
a2a4: 420
b2b4: 421
c2c4: 441
d2d4: 560
e2e4: 600
f2f4: 401
g2g4: 421
h2h4: 420
b1a3: 400
b1c3: 440
g1f3: 440
g1h3: 400

Nodes searched: 8902

My engine:

h2h3: 361
h2h4: 380
g2g3: 340
g2g4: 397
f2f3: 360
f2f4: 436
e2e3: 380
e2e4: 437
d2d3: 380
d2d4: 437
c2c3: 399
c2c4: 326
b2b3: 300
b2b4: 320
a2a3: 280
a2a4: 299
g1h3: 281
g1f3: 280
b1c3: 357
b1a3: 320

Depth: 3
Total nodes: 7070
Total time: 10ms/0.01s

I figured that my move generator was just buggy, and tried to track down the bugs by making a move the engine gives incorrect values for on the board and then calling perft() with depth = 2 on it to find out which moves are missing. But for all moves I tried this with, the engine suddenly starts to output the correct results I expected to get earlier!

Here is an example for the move a2a3:

When calling perft() on the initial position in stockfish, it calculates 380 subnodes for a2a3 at depth 3.
When calling perft() on the initial position in my engine, it calculates 280 subnodes for a2a3 at depth 3.
When calling perft() on the position you get after making the move a2a3 in the initial position in my engine, it calculates the correct number of total nodes at depth 2, 380:

h7h5: 19
h7h6: 19
g7g5: 19
g7g6: 19
f7f5: 19
f7f6: 19
e7e5: 19
e7e6: 19
d7d5: 19
d7d6: 19
c7c5: 19
c7c6: 19
b7b5: 19
b7b6: 19
a7a5: 19
a7a6: 19
g8h6: 19
g8f6: 19
b8c6: 19
b8a6: 19

Depth: 2
Total nodes: 380
Total time: 1ms/0.001s

If you have any idea what the problem could be here, please help me out. Thank you!

EDIT:

I discovered some interesting new facts that might help to solve the problem, but I don't know what to do with them:

For some reason, using std::sort() like this in perft():
```
std::sort(legal_moves.begin(), legal_moves.end(), [](auto first, auto second){     return first.get_from_index() % 8 > second.get_from_index() % 8; });
```
to sort the vector of legal moves causes the found number of total nodes for the initial position (for depth 3) to change from the wrong 7070 to the (also wrong) 7331.
When printing the game state after calling game_state.make_move() in perft(), it seems to have had no effect on the position bitboards (the other properties change like they are supposed to). This is very strange, because isolated, the make_move() method works just fine.

Solution

I'm unsure if you were able to pin down the issue but from the limited information available in the question, the best I can assume (and something I faced myself earlier) is that there is a problem in your unmake_move() function when it comes to captures since

Your perft fails only at level 3 - this is when the first legal capture is possible, move 1 and 2 can have no legal captures.
Your perft works fine when it's at depth 1 in the position after a2a3 rather than when it's searching at depth 3 from the start

This probably means that your unmake_move() fails at a depth greater than 1 where you need to restore some of the board's state that cannot be derived from just the move parameter you are passing in (e.g. enpassant, castling rights etc. before you made the move).