Assuming that each subject is only in one group, you have a nested design. Conceptually, it makes more sense to treat group as a fixed effect. As you have only three groups, it wouldn't make much sense statistically speaking anyways because you'd be asking the software to estimate a variance for
group assuming a normal distribution from only three observations.
Your sample size is quite low which might be a problem. Nevertheless, I suggest starting with a simple random intercept model:
lmer(outcome~day + group + (1|subject), data = dat)
This model fits a global intercept which is simply the intercept for the reference group, deviations from that intercept for the remaining groups, a single slope for the effect of
day and a random intercept for
subject. Hence, this model assumes that each group has the same trajectory over time, the same slope but different intercepts. To allow each group their own slope, you could fit the following model with an interaction between
lmer(outcome~day*group + (1|subject), data = dat)
Lastly, you could also allow for random slopes:
lmer(outcome~day*group + (day|subject), data = dat)
All these models assume a linear relationship between the
day which might be unrealistic. If you suspect nonlinear relationships, you could easily accomodate those by using polynomials or (my recommendation) restricted cubic splines (aka natural splines) with, say, 3 knots. Finally, have a look at this post which goes into more details about such longitudinal models.