A search is presented for additional neutral Higgs bosons in the τ τ final state in proton-proton collisions at the LHC. The search is performed in the context of the minimal supersymmetric extension of the standard model (MSSM), using the data collected with the CMS detector in 2016 at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 35.9 fb−1. To enhance the sensitivity to neutral MSSM Higgs bosons, the search includes production of the Higgs boson in association with b quarks. No significant deviation above the expected background is observed. Model-independent limits at 95% confidence level (CL) are set on the product of the branching fraction for the decay into τ leptons and the cross section for the production via gluon fusion or in association with b quarks. These limits range from 18 pb at 90 GeV to 3.5 fb at 3.2 TeV for gluon fusion and from 15 pb (at 90 GeV) to 2.5 fb (at 3.2 TeV) for production in association with b quarks, assuming a narrow width resonance. In the mhhod + scenario these limits translate into a 95% CL exclusion of tan β > 6 for neutral Higgs boson masses below 250 GeV, where tan β is the ratio of the vacuum expectation values of the neutral components of the two Higgs doublets. The 95% CL exclusion contour reaches 1.6 TeV for tan β = 60.[Figure not available: see fulltext.].