Why are my grid_scores_ from a RFECV, different from the score of the same cross validated model with the optimal features from the RFECV?

Question

I'm using sklearn's RFECV to come to the optimal set of features for my classification problem. I have X with 217 numerical features, for a binary label y. I determine the optimal set of features like so:

min_features_to_select = 3  
cv = RepeatedStratifiedKFold(n_splits=2, n_repeats=50, random_state=1)
model = LogisticRegression(penalty='l1', solver="liblinear")
rfecv = RFECV(estimator=model, step=1, cv=cv,
                  scoring="roc_auc", min_features_to_select=min_features_to_select)

rfecv.fit(X, y)

To visualize the rfecv process I plot the range of numbers of features against the grid scores of the RFECV

plt.figure()
plt.xlabel("Number of features selected")
plt.ylabel("Cross validation score (nb of correct classifications)")
plt.plot(range(min_features_to_select, len(rfecv.grid_scores_) + min_features_to_select), rfecv.grid_scores_)
plt.show()

refcv_supp = rfecv.get_support()
optimal_features = list(features[refcv_supp])

Which looks like this:

The grid_scores with all feature subset of the RFECV lie around a score of 0.660. I interper this score as the performance of the model that was used in the RFECV, with each particular set of features, expressed by the roc_auc.

Now when I isolate the optimal set of features, and feed this into the exact same cross validated model, the performance does not at all look like the scores depicted in the graph...

refcv_supp = rfecv.get_support()
optimal_features = list(features[refcv_supp])
X = X.loc[:,np.array(optimal_features)]

model = LogisticRegression(penalty='l1', solver="liblinear")
cv = RepeatedStratifiedKFold(n_splits=2, n_repeats=100, random_state=1)
n_scores = cross_val_score(model, X, y, scoring='roc_auc', cv=cv, n_jobs=-1, error_score='raise')

score = np.mean(n_scores)
std = np.std(n_scores)

Which results in a mean roc_auc of 0.896 with a std of 0.046. Can anyone help me with what goes wrong here?