Received ZZZ, revised ZZZ, accepted ZZZ Published online ZZZ (Dates will be provided by the publisher.)Keywords GaN, Si, O, n-typeUsing hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster calculations, we investigate stabilization of silicon and oxygen dopants in GaN. Formation energies of Si on Ga site and O on N site are calculated at two levels of theory using conventional thermochemical and kinetic exchange and correlation density functionals (B97-2 and BB1K). We confirm the shallow donor nature of these substitutional defects. We find that the 0/1+ transition levels for both Si and O species lie well above the bottom of the conduction band, in agreement with previous supercell-based simulations. The origin of this artefact is discussed in the context of relevant experimental results.Copyright line will be provided by the publishe GaN has many industrial applications such as solid state lighting and high power microelectronics [1]. For many devices, both n-type and p-type layers are required, and oxygen and silicon are most commonly used as n-type dopants [2]. Indeed, oxygen is often assumed to incorporate in as-grown GaN unintentionally and be the source of n-type conductivity [3,4]. Although many spectroscopical data are available on intentionally O-and Si-doped GaN [5][6][7][8][9][10], indicating the presence of carrier traps and particular luminescence bands [5][6][7][8], as well as electron spin resonance signatures [9][10], to date there have been only a small number of theoretical studies on the subject [11][12][13][14][15][16][17][18], mostly performed using density functional theory within a supercell approach, which do not support the presence of such traps [13][14][15][16][17][18].Here, we report the first hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster calculations on the electronic properties of O and Si defects in GaN. The advantage of our approach is in the accurate account of polarisation effects and unambiguous definition of the vacuum reference energy. In our method, the inner cluster of 116 atoms centred on the defect is treated with density functional theory using (i) the second generation thermochemical exchange and correlation hybrid functional B97-2 [19], which is similar to those commonly used in recent plane-wave supercell calculations (21% exact exchange compared with 25% for PBE0 [20] or HSE06 [21]), (ii) the SBKJC small-core pseudopotentials on Ga [22] within the cluster and large-core refitted pseudopotentials [23][24] on the outside that provide a suitable embedding potential on the defect, and (iii) the atomic basis set of def2-TZVP quality on N, O and Si [25] and matching SBKJC basis on Ga [22]. For comparison, we use a second hybrid exchange and correlation functional employing 42% exact exchange (BB1K) [26], fitted to reproduce kinetic barriers and thermochemical data, which gives a more accurate description of electron localization than B97-2 [24,27]. This QM region is embedded in an outer cluster which is treated with an MM leve...