The lighting engine creates lightmaps for exactly those chunks that have any nearby (i.e., in a radius of 1 subchunk) non-air blocks and deletes them once the last nearby subchunk becomes empty. As a consequence, unloading a chunk, which marks all its subchunks as empty, will delete ligtmaps of still loaded neighbors if there are no other non-empty subchunks around to keep the lightmaps alive.

This can be visualized via the following steps:

• Create a redstone-ready world
• Set the render distance to 2
• Run the following commands

  /setworldspawn 1000 56 1000
  /fill 16 80 0 16 95 15 minecraft:stone
  /setblock 16 88 8 minecraft:sea_lantern
  

• Fly to chunk (-14 0)
• Fly back to chunk (0 0) and observe that the lightmap for subchunk (0 5 0) got erased

There is also a variation for skylight

• Create a redstone-ready world
• Set the render distance to 2
• Run the following commands

  /setworldspawn 1000 56 1000
  /fill 0 128 0 47 128 47 minecraft:stone
  /fill 32 80 16 32 127 31 minecraft:stone
  /fill 33 128 16 47 128 31 minecraft:air
  /fill 32 96 16 32 96 31 minecraft:air
  

• Fly to chunk (-13 0)
• Fly back to chunk (1 0) and observe that the lightmap for subchunk (1 5 1) and (1 6 1) got erased

The vanilla code contains some data retainment mechanism that puts lightmaps back to the queued lightmaps, rather then deleting them. However, this mechanism is disabled upon promoting a chunk to the light stage.
One possible solution for this issue would be to reenable this retainment mechanism once a neighbor chunk gets unloaded.

I propose an alternative to this data retainment mechanism. This will decouple the lightmap handling from the skylight-optimization distance tracking that currently controls the lightmaps. Hence this approach will naturally avoid the bug discussed here and make the vanilla data retainment mechanism obsolete.
Furthermore, this approach will naturally contain a complete fix for MC-196614 and provide a more aggressive cleanup for trivial lightmaps, compared to the current vanilla code.

• The main idea is to create lightmaps on demand, i.e., when the lighting engine wants to set a light value, and delete them once they become trivial in the appropriate sense. When there is no lightmap associated to a subchunk, then skylight is inherited from the next lightmap above or direct skylight if there is none, whereas blocklight is simply 0. In order to check when a lightmap is trivial we need to track some complexity measure that can be easily updated upon chaning light values and is 0 iff the lightmap is trivial, i.e., if it coincides with the values that would be produced with no lightmap present.
• For example, we can simply take the sum of all blocklight values and for skylight we can for each position in the subchunk take the absolute value of the difference between the light value at that position and the light value at the position above (where in case of the topmost layer the light value at the position above needs to be taken from the lightmap above or direct skylight) and then sum up all these values. Note that in case of skylight this complexity measure not only depends on the lightmap but also on its position in the world.

• One important point is now that creating and removing trivial lightmaps does not change any light values, making the lightmap handling completely transparent to the lighting propagator. This also means that newly created lightmaps will always have a complexity of 0. This does however require MC-170010 to be fixed first, as otherwise some lightmaps will not be properly initialized upon creation, causing a change of light values and hence messing around with the light propagator. This causes a bit of a cyclic dependency between the two bug reports, so they should ideally be fixed simultaneously.

• The skylight optimization is then applicable to subchunks that are not near any non-air blocks, as determined by the current distance tracking, and that don't have an associated lightmap. This second condition that having a non-trivial lightmap disables the skylight optimization basically takes care of the second part of MC-196614.
• Note that changing a light value at the bottom of a lightmap will cause subchunks below without an associated lightmap to automatically change their light value as well, since that is how missing lightmaps are handled in the light lookup. This means that before changing a light value at the bottom of a lightmap, we first have to find the next subchunk below for which skylight optimization is not applicable and create a lightmap for it, in case it does not already have one, fixing the old light value so that changes are then properly propagated through the usual skylight optimization.

• It is convenient to forcefully disable any light propagations to chunks before the pre_light stage (using terminology of MC-170012), allowing to defer the complexity initialization until the promotion to the pre_light stage. When unloading chunks, propagations should then be forcefully disabled again, so that the unloaded region is very well controlled. Ideally, this manual mechanism should not be necessary as any violation of it is always accompanied by some lighting glitch. However, in view of MC-164281 it might be a good idea to take some precautions and avoid further bugs caused by screwed complexity trackings.
• Furthermore, for the unloaded region lightmaps can be directly added to the world as nothing will be interacting with them, given that we have just forcefully disabled all such interactions. Lightmaps getting queued for already loaded chunks, i.e., on the client, might be better placed in some queue first, so they can be added at a more controlled point in time. Furthermore, when adding such a lightmap to an already loaded chunk its complexity tracking has to be reevaluated and if any value at the bottom changes, the necessary steps need to be taken for skylight, as explained above.

• Removal of trivial lightmaps should only be done once every update cycle or upon saving or something similar, in order to avoid rapidly removing and reallocating lightmaps

One concern that might come up is regarding the interaction between skylight optimization and unloaded chunks. More concretely, note that the last accessible chunk is in the full state but its neighbors are only guaranteed to be in pre_light state. Hence there can be light updates into chunks that are only in pre_light state and their neighbors are not guaranteed to be loaded at all. So the issue one might see here is that there are light updates into chunks that do not have complete information about their neighbors, so the skylight optimization might be applied in situations where it shouldn't. We will in the following sketch a proof that the extra condition that the skylight optimization is disabled for subchunks with an associated lightmap is in fact sufficient to ensure correct results. Note that some similar discussion could have been already placed under MC-170012 or MC-196614, but I think it fits here nicely as well. Of course one could simply avoid this whole discussion by just increasing the chunk loading distance by 1, so that light updates only take place in chunks that are in full state; however that would be rather inefficient given that this is actually not an issue.

For the following discussion we will assume MC-170012 to be fixed and use its terminology. We will assume from the chunk loading mechanism that whenever a block change in a chunk c1 could potentially cause a light update into chunk c2 (or some boundary touching it) then c2 should be loaded in pre_light stage. Note that any violation of this will result in a lighting glitch simply because light updates would get stuck at chunk boundaries, so this would be a completely unrelated issue on its own. More technically, we can state this condition as follows: If there occurs a block update in some chunk c1 and we are given two neighbors c and c2 of it (where we also count diagonal neighbors and also say that a chunk is a neighbor of itself) such that they are also neighbors to one another and such that c has its initial lighting done, i.e., was generated into the light stage but is not necessarily loaded currently, then c2 should be loaded in pre_light stage.
Note that this assumption might not be true on the client. However there are other mechanisms taking care of such effects for the client light syncing, so we don't want to worry about this here.

We want to show that the end result of our lighting model is correct, for which it suffices to show local correctness, i.e., each block has a light value that is precisely the maximum over all contributions of its neighbors, where the contributions are specified in the ideal lighting model where no skylight optimization takes place and all chunks are loaded. It is always true that subchunks for which no skylight optimization is applied are locally correct and also unloaded chunks retain their local correctness since by the assumption above no block directly adjacent to an unloaded chunk changes its light value. Hence it is enough to show local correctness for subchunks for which skylight optimization is applied.

So, for sake of concreteness assume that skylight optimization is applied to subchunk (0 0 0) for which we want to show local correctness. We may assume inductively that all subchunks at y-level >= 1 are locally-correct.
Consider some horizontal chunk border of it, e.g., the border to chunk (-1 0), consisting of the blocks (0 0 0)..(0 15 15). If no light value for a block adjacent to this border changed, i.e., for no block in the region (-1 -1 -1)..(1 16 16) then it retains its local correctness and there is nothing to check for this border. Otherwise, such a light change must come from a block change in a chunk from the region (-1 -1)..(0 1) and the skylight source propagated through this block change must lie in the same chunk region and furthermore this skylight source and the block change must lie in chunks that are neighbors of one another. Hence we can conclude by the assumption about chunk loading that some quadrant containing the border is loaded in pre_light stage, i.e., the quadrant (-1 -1)..(0 0) or (-1 0)..(0 1). By the specification for the skylight optimization we then know that the subchunks of this quadrant that lie in the region (-1 -1 -1)..(1 1 1) contain only air blocks. Similarly, we conclude that a corner-column of subchunk (0 0 0) either retains local correctness or the unique quadrant containing the column in its interior is loaded in pre_light stage and the intersection with (-1 -1 -1)..(1 1 1) consists purely of air blocks.
Consider now the region consisting of all such quadrants intersected with (-1 -1 -1)..(1 1 1). This region then has the following properties which allow to deduce quite a lot of information about the lighting model:

• It is star-shaped with respect to chunk (0 0)
• It contains only air blocks
• Its boundary can be decomposed into 3 parts:
  • The top faces
  • The faces that are one subchunk away from (0 0 0) (excluding the top faces)
  • The faces touching (0 0 0)
• The faces and corner-columns of (0 0 0) that also belong to the boundary of this region retain local correctness, so there is nothing to check for those

We need to show local correctness for those blocks in (0 0 0) that are strictly inside this region, i.e., excluding those borders for which we already know it. We now consider two lighting models: The real model applying the skylight optimization according to our specification, where we replace propagations at y-levels >= 16 with ideal ones, as we already know local correctness in that area by induction. And the semi-ideal model that additionally treats subchunk (0 0 0) as non-optimizable, i.e., uses ideal propagations. We will show that both models give the same solution with boundary conditions given by the boundary of the region. This suffices since the solution to the semi-ideal model is locally correct on subchunk (0 0 0) by construction. In order to show this it is in turn sufficient to show that both solutions agree on subchunk (0 0 0) since we can then add it to the boundary conditions and both models are identical outside this subchunk, hence have the same solution.
For this we show that for any propagation path from the boundary of the constructed region into subchunk (0 0 0) in the real model there exists a corresponding propagation path from the boundary to the same block in the semi-ideal model that has at least the same contribution, and conversely for every propagation path in the semi-ideal model there is a corresponding path in the real model that has at least the same contribution.
Note that the distinction between changed and unchanged borders was necessary as the unchanged borders are now treated as sources in this model, hence proving local correctness in this model does not provide useful information about them (which is no problem as we already know that these are fine).
We consider the following cases according to the decomposition of the boundary as above:

• Both in the real and semi-ideal model there are no propagation paths from the faces that are one subchunks away from (0 0 0), excluding the top faces, so this case is easy
• Both in the real and semi-ideal model the only propagation paths starting from the top faces start with direct skylight access, i.e., with a skylight value of 15. An optimal propagation path in the real model first goes straight down until y=16, then propagates in the xz-plane and then propagates for the remaining y-level using the skylight optimization, hence not reducing the light value. An optimal path in the semi-ideal model only does part of its xz-movement at y=16 and does the rest at the final y-level, since subchunk (0 0 0) is treated as non-optimizable. This gives a 1:1 correspondence between optimal real and optimal semi-ideal paths giving the same light contributions, as desired.
• Finally consider propagation paths starting from a face touching (0 0 0). Both for real and semi-ideal propagation paths we can find an alternative one giving at least the same contribution that first travels along the boundary of the region, then travels on the xz-plane inside a single non-optimizable subchunk (including (0 0 0) in case of the semi-ideal model) and then possibly travels downwards through an optimizable subchunk using skylight optimization. To prove this consider such a path appended by a single-step propagation into a non-optimizable subchunk. Using that the region consists purely of air blocks, the propagations can then be reordered in such a way that the subchunk is crossed along the boundary of the region, so that this part can be absorbed into the first part of the path and the resulting path is again of the desired form.
  Note that the boundary of the region is "locally correct" with respect to the real model, i.e., each light level on the boundary is at least the maximum of all real contributions from neighbors lying on the boundary, or in other words, the propagated value of a real path along the boundary is at most the given boundary value. This is true simply because the boundary values are chosen from the global solution of the real model. The same is also true for the semi-ideal model, since it only differs from the real model on subchunk (0 0 0) and for blocks of this subchunk that belong to the boundary of the region we know that they are locally correct in the ideal model and hence also in the semi-ideal model. Hence we can leave out this first part of the propagation paths travelling along the boundary as the new starting point will then still lie on the boundary and has at least the same contribution by the argument just given.
  So we can restrict our attention to those optimal paths first travelling on the xz-plane of a single non-optimizable subchunk and then possibly travelling downwards through an optimizable subchunk using skylight optimization.
  In the real model an optimal path now starts at y=16 (which is treated as a non-optimizable subchunk in both models) then travels on the xz-plane and then travels down into subchunk (0 0 0) using skylight optimization, so that the value doesn't change. An optimal semi-ideal path starts at the correct y-value and then travel on the xz-plane. The important observation is now that since subchunk (0 0 0) is optimizable, it does not have an associated lightmap and hence its light values are constant along columns (including y=16). Hence we can transform an optimal real path into an optimal semi-ideal path by starting at the correct y-value, which then has the same starting value and gives the same contribution. And also conversely an optimal semi-ideal path can be transformed into a real path by starting at y=16 with the same starting value and same contribution. This produces again a 1:1 correspondence between optimal real and semi-ideal propagation paths giving rise to the same contributions, as desired.

This finishes this rather lengthy proof.

Note that this argument requires a slightly inefficient specification for the skylight optimization, namely that non-air blocks also make the neighbor subchunks above non-optimizable. This is the one implemented by Vanilla. (This was used when transforming paths to first travel along the boundary and then only inside a single chunk. Without this assumptions the subchunks at y=-1 might not be empty and there would be more complicated optimal propagation paths).
In the absence of unloaded chunks this is not needed and it would be sufficient that a non-air block only makes neighbor subchunks that have less or equal y-value non-optimizable. In the presence of unloaded chunks the argument given above requires the slightly inefficient version and in fact one can (easily) construct examples showing that the more efficient specification would produce wrong results and is hence not sufficient.
So if one would like to implement this optimization, extra care needs to be taken in the presence of unloaded chunks and the optimization needs to be disabled in those cases or information about unloaded neighbors needs to be provided through other means.

Best,
PhiPro