I was struggling with a similar question and found this discussion for interacting categorical smooths (r - mgcv GAM: more than one variable in `by` argument (smooth varying by more than 1 factor) - Stack Overflow) using the interaction() function, which up until now I did not know existed. Seemed to work fine for my model, I haven’t found any apparent issues so far; it did in fact generate four distinct smooths from two binary predictors with interaction.
So perhaps try something like performance ~ 1 + course*day + s(sectionIndex, by = interaction(course,day)) + (1|bib), as the specification for the population smooth. I suspect in this case that you will also want subject-specific smooths, if you have observations across time for each subject. The subject-level intercepts will not handle variation in the smooths by subject.