Attitude Determination (AD) constitutes an important navigation component for vehicles that require orientation information, such as spacecraft or ships. Global Navigation Satellite Systems (GNSS) enable resolving the orientation of a vehicle in a precise and absolute manner, by employing a setup of multiple GNSS antennas rigidly mounted onboard the tracked vehicle. To achieve high-precision attitude estimation based on GNSS, the use of carrier phase observations becomes indispensable, with the consequent added complexity of resolving the integer ambiguities. The use of inertial aiding has been extensively exploited for AD, since it enables tracking fast rotation variations and bridging short periods of GNSS outage. In this work, the fusion of inertial and GNSS information is exploited within the recursive Bayesian estimation framework, applying an Error State Kalman Filter (ESKF). Unlike common Kalman Filters, ESKF tracks the error or variations in the state estimate, posing meaningful advantages for AD. On the one hand, ESKF represents attitude using a minimal state representation, in form of rotation vector, avoiding attitude constraints and singularity risks on the covariance matrix estimates. On the other hand, second-order products on the derivation of the Jacobian matrices can be neglected, since the errorstate operates always close to zero. This work details the procedure of recursively estimating the attitude based on the fusion of GNSS and inertial sensing. The GNSS attitude model is parametrized in terms of quaternion rotation, and the overall three-steps AD procedure (float estimation, ambiguity resolution and solution fixing) is presented. The method performance is assessed on a Monte Carlo simulation, where different noise levels, number of satellites and baseline lengths are tested. The results show that the inertial aiding, along with a constrained attitude model for the float estimation, significantly improve the performance of attitude determination compared to classical unaided baseline tracking.Traditionally, the GNSS compass model relates the double difference observations to the baseline vectors relating the antennas, expressed in the global frame. The use of this model in combination with the standard LAMBDA method has been often applied for AD [6][7][8][9]. However, classical LAMBDA is not designed for AD, due to the nonlinear constraints related to the prior knowledge of the distance and relative orientation between the antennas [10]. Aiming to improve LAMBDA for GNSS carrier phase-based AD, Teunissen and Giorgi introduced the Constrained-LAMBDA (C-LAMBDA) and the Multivariate Constrained-LAMBDA (MC-LAMBDA) methods. First, C-LAMBDA modifies the original LAMBDA method introducing the information on the baseline length as a constraint, exploited during the process of AR and solution fixing [11][12][13]. Although C-LAMBDA is shown to provide significant improvement against classical LAMBDA, additional information on the relative orientation between antennas is not included. In MC-LAMBDA, the...