What you're fitting with `aov` is called a strip plot, and it's tricky to fit with `lme` because the `subject:A` and `subject:B` random effects are crossed.

Using `lmer` from the nlme package, the code would be something like

    lmer(Y ~ A*B + (1|subject) + (1|A:subject) + (1|B:subject), data=d)    

These threads from R-help may be helpful.

http://www.biostat.wustl.edu/archives/html/s-news/2005-01/msg00091.html

http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/3328