  NGC7538 IRS1 is considered the best high-mass accretion disk candidate around
an O-type young star in the northern hemisphere. We investigated the 3D
kinematics and dynamics of circumstellar gas with very high linear resolution,
from tens to 1500 AU, with the ultimate goal of building a comprehensive
dynamical model for this YSO. We employed four different observing epochs of
EVN data at 6.7 GHz, spanning almost eight years, which enabled us to measure,
besides line-of-sight (l.o.s.) velocities and positions, also l.o.s.
accelerations and proper motions of methanol masers. In addition, we imaged
with the JVLA-B array highly-excited ammonia inversion lines, from (6,6) to
(13,13), which enabled us to probe the hottest molecular gas very close to the
exciting source(s). We found five 6.7 GHz maser clusters which are distributed
over a region extended N-S across ~1500 AU and are associated with three peaks
of the radio continuum. We proposed that these maser clusters identify three
individual high-mass YSOs, named IRS1a, IRS1b, and IRS1c. We modeled the maser
clusters in IRS1a and IRS1b in terms of edge-on disks in centrifugal
equilibrium. In the first case, masers may trace a quasi-Keplerian thin disk,
orbiting around a high-mass YSO, IRS1a, of up to 25 solar masses. This YSO
dominates the bolometric luminosity of the region. The second disk is both
massive (<16 Msun within ~500 AU) and thick, and the mass of the central YSO,
IRS1b, is constrained to be at most a few solar masses. In summary, we present
compelling evidence that NGC7538 IRS1 is not forming just one single high-mass
YSO, but consists of a multiple system of high-mass YSOs, which are surrounded
by accretion disks, and are probably driving individual outflows. This new
model naturally explains all the different orientations and disk/outflow
structures proposed for the region in previous models.
