Why `max_features=n_features` does not make the Random Forest independent of number of trees?
Interesting puzzle indeed.

First things first. The DecisionTreeClassifier has some stochastic behavior. For instance, the splitter code iterates through the features at random:

        f_j = rand_int(n_drawn_constants, f_i - n_found_constants,
                       random_state)


Your data is small and comes from the same distribution. What this means is that you'll have a lot of identical purity scores depending on how iteration is done. If you (a) increase your data, or (b) make it more separable, you'll see the problem should ameliorate.

To clarify: if the algorithm computes the score for feature A and then computes the score for feature B and it gets score N. Or if it computes first the score for feature B and then for feature A and it gets the same score N, you can see how each decision tree will be different, and have different scores during test, even if the train test is the same (100% if max_depth=None of course). (You can confirm this.)

During my exploration of your question, I have produced the following code with my own implementation of a random forest. Since it took me some time, I figured I might as well paste it here. :) Seriously, it can be useful. You can try to disable random_state from my implementation to see what I mean.

from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
import numpy as np


class MyRandomForestClassifier:
    def __init__(self, n_estimators):
        self.n_estimators = n_estimators

    def fit(self, X, y):
        self.trees = [DecisionTreeClassifier(random_state=1).fit(X, y)
                      for _ in range(self.n_estimators)]
        return self

    def predict(self, X):
        yp = [tree.predict(X) for tree in self.trees]
        return ((np.sum(yp, 0) / len(self.trees)) &gt; 0.5).astype(int)

    def score(self, X, y):
        return accuracy_score(y, self.predict(X))


for alpha in (1, 0.1, 0.01):
    np.random.seed(1)
    print('# alpha: %s' % str(alpha))
    N = 1000
    X = np.random.random((N, 10))
    y = np.r_[np.zeros(N//2, int), np.ones(N//2, int)]
    X[y == 1] = X[y == 1]*alpha
    Xtr, Xts, ytr, yts = train_test_split(X, y)

    print('## sklearn forest')
    for n_estimators in (1, 10, 100, 200, 500):
        m = RandomForestClassifier(
            n_estimators, max_features=None, bootstrap=False)
        m.fit(Xtr, ytr)
        print('%3d: %.4f' % (n_estimators, m.score(Xts, yts)))

    print('## my forest')
    for n_estimators in (1, 10, 100, 200, 500):
        m = MyRandomForestClassifier(n_estimators)
        m.fit(Xtr, ytr)
        print('%3d: %.4f' % (n_estimators, m.score(Xts, yts)))
    print()


Summary: Each DecisionTreeClassifier is stochastic, data such as yours, which is small and comes from the same distribution, are bound to produce slightly different trees, even if the random forest itself is deterministic. You can fix this by passing the same seed to each DecisionTreeClassifier which you can do using random_state=something. RandomForestClassifier also has a random_state parameter which it passes along each DecisionTreeClassifier. (This is slightly incorrect, see the edit.)

EDIT2: While this removes the stochasticity component of the training, the decision trees would still be different. The thing is that sklearn ensembles generate a new random seed for each child based on the random state they are given. They do not pass along the same random_state.

You can see this is the case by checking the _set_random_states method from the ensemble base module, in particular this line, which propagates the random_state across the ensembles' children.