import arcpy

#Set these as the name of the layer in ArcMap's TOC
counties = "County_layer"
hospitals = "Hospital_layer"
#Set this as the name of the field to find max value for
field = "Doctors"

object_id = []
hospital_fl = arcpy.MakeFeatureLayer_management(hospitals, "hospital_fl")

with arcpy.da.SearchCursor(counties, ("SHAPE@")) as rows:
    for row in rows:
        polygon_geom = row[0]
        arcpy.SelectLayerByLocation_management(hospital_fl, "COMPLETELY_WITHIN", polygon_geom)
        doctors = {row[0]: row[1] for row in arcpy.da.SearchCursor(hospital_fl, ("OID@", field))}
        max_doctors = max(doctors.values())
        for k, v in doctors.items():
            if v == max_doctors:
                object_id.append(k)

arcpy.SelectLayerByAttribute_management(hospital_fl, "NEW_SELECTION", "OBJECTID IN {0}".format(tuple(object_id)))
arcpy.SpatialJoin_analysis(counties, hospital_fl, "in_memory\Counties_output", "JOIN_ONE_TO_MANY", "KEEP_ALL", "#", "CONTAINS")
arcpy.Delete_management(hospital_fl)