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.

Your first attempt is equivalent to `aov(Y ~ A*B + Error(subject), data=d)`, which doesn't include all the random effects; your second attempt is the right idea, but the syntax for crossed random effects using lme is very tricky. 

Using `lme` from the `nlme` package, the code would be

    lme(Y ~ A*B, random=list(subject=pdBlocked(list(~1, pdIdent(~A-1), pdIdent(~B-1)))), data=d)

Using `lmer` from the `lme4` 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 (and to give credit, that's where I got the `nlme` code from).

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

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

http://www.mail-archive.com/[email protected]/msg10843.html

This last link refers to p.165 of Pinheiro/Bates; that may be helpful too.


EDIT: Also note that in the data set you have, some of variance components are negative, which is not allowed using random effects with lme, so the results differ.  A data set with all positive variance components can be created using  a seed of 8.  The results then agree. See [this answer](http://stats.stackexchange.com/questions/14088/why-do-lme-and-aov-return-different-results-for-repeated-measures-anova-in-r/14185#14185) for details.  

Also note that `lme` from `nlme` does not compute the denominator degrees of freedom correctly, so the F-statistics agree but not the p-values, and `lmer` from `lme4` doesn't try too because it's very tricky in the presence of unbalanced crossed random effects, and may not even be a sensible thing to do.  But that's more than I want to get into here.