One solution is to decouple each chunk from its neighbours.
Let's suppose you have nxn chunks, with a 1-tile wide border strip between them, forming a grid lattice.
We'll generate the corners of this grid first, then use the neighbouring corner information to generate the border edges, then use the completed border frame to fill the chunk.
So, starting with the corners. You know you have a corner at every pair (i * (n+1), j * (n+1))
where i and j are integers.
Run this coordinate through a hash function, and use the resulting scrambled bits to select a tile to place here. You can also fold your world seed into the hash to get different corner patterns for different seeds.
This way, the tile selected at each corner depends only on its position and the world seed, not its neighbours or the order in which they were generated. You can use other generation strategies here too, say if you want to create large scale structure to your world, you could sample a noise function at each corner to decide whether it's sea or lowland or mountain, rather than selecting them in an uncorrelated random fashion.
So, we now have the 4 corners of the chunk we want to generate picked out. Next we generate the edges.
First, seed a random number generator with the coordinates of the middle tile in this edge (and your world seed). We'll use that for all the random choices ahead, so they depend only on the edge's position in the grid, not the history of what you generated before.
Take the 1 x (n+2) rectangle of tiles spanning from one corner to its next neighbour, and run a modified wave function collapse on that strip. This modified version will match only neighbours to the left and right, and ignore neighbours above and below (or vice versa when generating the vertical edges).
You might need to constrain the tiles you can select at/near the corners just to make sure you can't create a situation that's impossible to fill in later, like....
A B
C ?
where A is the corner of a chunk, A & B are valid neighbours, A & C are valid neighbours, but no tile exists that can match with a B above and a C to the left. If your tileset is reasonably permissive then this situation shouldn't require major constraints to avoid.
After repeating this process for each edge surrounding our chunk, we have neighbours on all sides, and we're ready to fill it in.
Again we'll spin up a seeded random number generator with the coordinates of the middle of this chunk, so we don't leak previous generator state into our subsequent choices. Then we'll run the big wave function collapse pass on this inner square, using the border of already-collapsed neighbours we have on all four sides.
Since the interior depends only on the edges, and the edges depend only on the corners, and the corners depend only on position, we'll generate the same corners/edges/chunk interior in a given location, no matter what order we generated them in.
The method I presented above will work with many tilesets, but not all, and it has a risk of showing a bias. The pattern of tiles you get from the 1-dimensional strip WFC we use on the edges between chunks might not have the same statistical properties as a randomly chosen 1-tile-wide strip carved out of 2D WFC with the same tile set. Interactions with tiles chosen on adjacent rows can create/propagate additional constraints that change the distribution we observe.
Here's an example of a tileset that would cause trouble for this algorithm (assuming no tile rotation/flipping is allowed). The constraints of the corner tiles in the top and bottom row mean that yellow "north" roads and red "south" roads can only appear in alternation, and the number of blue gaps we expect to see between them is a function of the transition probabilities of the orange horizontal road. But when looking at a single strip at a time, WFC doesn't know that. It knows an empty blue tile can go next to a yellow road or a red road, and a yellow road or a red road can go next to a blue tile, so it might place blue-yellow-blue-yellow with no red in-between, or even if adjusted to avoid this (such as by distinguishing "left of road blue" from "right of road blue"), may still be prone to creating improbable spacings between the reds and yellows.
Even in less drastic scenarios than this, there might be multi-tile formations that span across rows and skew the observed probability of certain left-right pairs in a single row, compared to WFC's baseline probabilities.
To solve this issue, Boris the Brave suggests a layered block-based generation strategy. It proceeds like so:
Initialize your world to an empty tile or default pattern of your choosing. (Could even be a large-scale Perlin noise autotile, etc.)
Split the world into blocks, with 1-tile gaps between them. Run WFC on the interior of those blocks, using the tiles in the gaps as the initial constraints.
Because of the separation, each of these blocks can be computed independently. Seed your randomizer with the coordinates of the block to ensure each one is deterministically generated, regardless of order.
Now shift your block grid by half a block space, so that the new blocks cross what had been the vertical edge between blocks in the previous pass. Repeat. (So we erase most of the tiles we just generated in step 2, keeping only the tiles in the gaps)
Shift the block grid again, to span the horizontal edge between blocks in the previous pass. Repeat.
Shift the block grid once more, to span the vertical edge from the previous pass. Repeat.
So we've run the generator 4x times, re-doing a lot of work. But now every tile that's left was generated from the interior of a 2D WFC at some point in the chain, and so any slice we take of the resulting map should have all the statistical properties we expect, without bias between edges / chunk centers.
There's a great animation here demonstrating how this ensures you only ever need to evaluate a finite number of blocks to get a deterministic output block.
Based on this, I do wonder if you could get "good enough" anti-biasing in fewer passes, by getting rid of the 1-tile gaps between blocks/chunks, and populating them in a checkerboard pattern. Say all "black" squares get populated first, with no edge constraints at all. Then all "white" squares get populated based on their four adjacent black squares, ensuring connectivity. You could then follow up with a third pass at an offset to the original checkerboard in case it shows visible biases along its grid. This would mean something like 8-9 block evaluations for a single output block, rather than 13 for the scheme above. Something to experiment with!