This has been in my todo list for a while, but as always "it's complicated". If you can pull it off that would be a rather nice simplification.

The basic idea (just a quick write up from memory, haven't rechecked) is that the "initial" data limits of a newly created axes is x=(0, 1), y=(0, 1) (so that these are the limits of an empty axes), but this bbox is marked as "ignore_existing_data_limits"/"ignore", so that once one adds data to it for the first time the data limits of that data completely overrides the 0-1. Note that this results in a slight inconsistency with what happens if you call plot() with no data; in that case Matplotlib knows it shouldn't use the "ignored" 0-1 limits but still has no data to autoscale on and thus relies on the nonsingular expander to expand the limits to (-0.05, 0.05). mplot3d has a similar mechanism, except that of course it's different and it uses the has_data() flag to know whether to ignore_existing_data_limits (again, just from my memory, some details may be wrong).

I'm not really sure whether we can choose to change the behavior of axes() to default to (-0.05, 0.05) instead or of an empty plot() to default to (0, 1) instead (these seem rather benign API changes, but what do I know...), but I do think it would likely be simpler if the initial data limits were "empty" (thus getting rid of the ignore_existing_data_limits flag) so that even the first plot can just do a normal datalimits bbox merge and if one wants to show such an empty plot we just rely on the normal nonsingular expander.

Didn't think there would be so much overlap when I started work on #17115, but now I see these issues are fundamentally linked. I'll keep that PR draft separate as it has a different goal....but will keep an eye on this, to avoid duplicating effort.

I don't have the capacity to follow up on this. If somebody wants to take over, @anntzer's comment should be taken into account.