A test of gravitational theories including torsion with the BepiColombo radio science experiment

The Mercury Orbiter radio Science Experiment (MORE) is one of the experiments on-board the ESA/JAXA BepiColombo mission to Mercury, to be launched in October 2018. Thanks to full on-board and on-ground instrumentation performing very precise tracking from the Earth, MORE will have the chance to determine with very high accuracy the Mercury-centric orbit of the spacecraft and the heliocentric orbit of Mercury. This will allow to undertake an accurate test of relativistic theories of gravitation (relativity experiment), which consists in improving the knowledge of some post-Newtonian and related parameters, whose value is predicted by General Relativity. This paper focuses on two critical aspects of the BepiColombo relativity experiment. First of all, we address the delicate issue of determining the orbits of Mercury and the Earth-Moon barycenter at the level of accuracy required by the purposes of the experiment and we discuss a strategy to cure the rank deficiencies that appear in the problem. Secondly, we introduce and discuss the role of the solar Lense-Thirring effect in the Mercury orbit determination problem and in the relativistic parameters estimation.


I. INTRODUCTION
The Einstein's General Theory of Relativity (GR) accounts for all the experimental findings concerning gravitational interaction available at present. According to GR, the flat Minkowsky spacetime is deformed by the matter, giving rise to a non-flat and dynamic spacetime, tied by Riemannian geometry. On the other hand, it is well-known that the other fundamental interactions can be successfully described, at microscopic level, within a rigid and flat spacetime. Since GR was originally formulated as a theory involving mass distribution at macroscopic level, it is desirable to consider suitable generalizations of GR that include micro-physical processes, which could possibly induce macroscopic effects, to be constrained, in turn, by experiments (cfr., e.g., the discussion in [1]).
In the following, we will consider a particular generalization of GR which includes a non-vanishing spacetime torsion. Elementary particles, which constitute matter, are characterized by their intrinsic spin. Since spin averages out at macroscopic level, GR considers that the dynamical behavior of a macroscopic distribution of mass can be described by the energy-momentum tensor of matter alone, which is coupled to the metric g µν of a Riemann spacetime. However, at microscopic level also the spin angular momentum plays a role in the dynamics, thus it must be coupled in some way to spacetime. This fact leads to the formulation of a more general spacetime, the four-dimensional Riemann-Cartan spacetime, characterized by an additional, non-Riemannian part of the affine connection, called the contorsion tensor, which should be coupled to the spin (see, e.g., [1,2]).
Torsion as the anti-symmetric part of an asymmetric affine connection was firstly introduced one century ago bý Elie Cartan [3], who glimpsed the link between spacetime torsion and intrinsic angular momentum of matter (for an historical perspective see, e.g., [1]). If torsion is admitted, it might affect spinning particles and thus, indirectly, act on light and test particles throughout the field equations, which determine the metric. Most of the torsion theories predict a negligible amount of torsion in the solar system [4]. Indeed, beyond GR, theories of gravitation usually assume that torsion couples only to the intrinsic spin of particles and not to rotational angular momentum (see, e.g., [5]).
Recently, Mao et al. [6] considered the issue from another point of view: we can consider a given experiment in the solar system and model the problem by means of a non-standard gravitational Lagrangian, which includes a detectable torsion signal; then, whether the model is in agreement or not with the data can be experimentally tested. Later on, March et al. [7] developed a more general framework, within the scope of the parameterized post-Newtonian (PPN) approximation (see, e.g., [8]), in order to test the possible effects due to torsion around massive bodies in the solar system.
In this perspective, the radio science experiment on-board the BepiColombo mission to Mercury gives an intriguing opportunity to test possible modifications of GR in the solar system. BepiColombo is an ESA/JAXA mission launched in October 2018. The mission aims at a comprehensive exploration of the planet Mercury, thanks to two spacecraft, whose orbit insertion around the planet is expected for the end of 2025 [9]. The Mercury Orbiter Radio science Experiment (MORE) is one of the on-board experiments, devised to enable a better understanding of the geodesy and geophysics of Mercury, on one side, and of fundamental physics, on the other (see, e.g., [10][11][12]). Thanks to full on-board and on-ground instrumentation capable to perform very precise tracking from the Earth, MORE will have the chance to determine with very high accuracy the Mercury-centric orbit of the spacecraft and the heliocentric orbit of Mercury and the Earth. Taking advantage from the fact that Mercury is the best-placed planet to investigate the gravitational effects of the Sun, MORE will allow an accurate test of relativistic theories of gravitation (relativity experiment). The test consists in constraining simultaneously the value of a number of post-Newtonian (PN) parameters together with some other parameters of general interest, by means of a global non-linear least squares fit. This can be achieved, as will be described later, by means of a dedicated orbit determination software, ORBIT14, developed at the University of Pisa.
The paper is organized as follows. In Section II we briefly describe the MORE relativity experiment: we present the dynamical model and we discuss the adopted mathematical methods, together with some generalities on the orbit determination (OD) software, ORBIT14. In Section III we describe the generalized dynamical model including the possible effects due to spacetime torsion and we illustrate how it has been implemented in the software. In Section IV we present and discuss the results of the numerical simulations of the MORE relativity experiment including torsion parameters. Finally, in Section V we present a discussion and our conclusions.

II. MORE RELATIVITY EXPERIMENT
It is a long-established fact that a space mission to Mercury, being the inner planet of the solar system, hence the most sensitive to the gravitational effects of the Sun, would improve significantly the determination of the PN parameters, constraining more tightly their agreement, or disagreement, with GR predictions [13]. The BepiColombo opportunity precisely suits this scope, since it is equipped with on-board instrumentation capable of very precise tracking from the Earth [14] in order to perform a comprehensive radio science experiment (MORE), consisting of a gravimetry, a rotation and a relativity experiment.
The global experiment aims at determining the following quantities of general interest: • the spherical harmonics coefficients of Mercury gravity field (see, e.g., [15] -Chap. 13) with a signal-to-noise ratio better than 10 up to, at least, degree and order 25 and the Love number k 2 [16]; • the parameters defining the model of Mercury's rotation (all the details can be found in [17]); • the "relativity" parameters, which are the PN parameters γ, β, α 1 , α 2 and the Nordtvedt parameter η, which characterize the expansion of the spacetime metric in the limit of slow motion and weak field (see, e.g., [8,18]), together with some related parameters, such as the oblateness of the Sun J 2⊙ , the solar gravitational mass µ ⊙ = GM ⊙ (where G is the gravitational constant and M ⊙ the mass of the Sun), possibly its time derivative ζ = 1/µ ⊙ dµ ⊙ /dt, and the solar angular momentum GS ⊙ which appears in the Lense-Thirring effect on the orbit of Mercury (see [19] for a detailed discussion).
To achieve the challenging scientific goals of MORE, it is mandatory to perform a very precise determination of the orbit of the spacecraft around Mercury and of the orbit of Mercury and the Earth (replaced, in fact, with the Earth-Moon barycenter (EMB) for the purposes of our analysis). This is enabled, in turn, by state-of-the-art onboard and on-ground instrumentation [14,20]. The on-board transponder will be able to collect the radio tracking observables (range, range-rate) up to a goal accuracy (in Ka-band) of about σ r = 15 cm at 300 s for one-way range and σ rr = 1.5 × 10 −4 cm/s at 1000 s for one-way range-rate.

II.1. The heliocentric dynamics of Mercury and the EMB
To perform the OD of the spacecraft, Mercury and the EMB at the required level of accuracy, their dynamics need to be modeled very carefully. In practice, although we are dealing with a unique global experiment and a unique set of measurements, we can conceptually separate the gravimetry-rotation experiments from the relativity experiment (see, e.g., [21]). As a matter of fact, gravimetry and rotation mainly involve short period phenomena, being preferentially related to the orbital motion of the spacecraft around Mercury (which has a periodicity of 2.3 hours), whereas relativistic phenomena take place over longer time scales, of the order of months or years, thus they can be analyzed mainly by studying the heliocentric motion of Mercury and the Earth, taking place over 88 and 365 days, respectively. Furthermore, the chance to perform the relativity experiment independently from the others, for the purpose of simulations, is even more legitimate if we consider the goal accuracies of the observations: properly comparing σ r and σ rr over the same integration time according to Gaussian statistics, we find that σ r /σ rr ∼ 10 5 s. Thus, we can infer that range observations are more accurate when observing phenomena with periodicity longer than 10 5 s, like relativistic phenomena, while the opposite holds for gravimetry and rotation, which are, in turn, performed mainly by means of range-rate observations.
All the details concerning the Mercury-centric dynamical model of the spacecraft can be found in [21,22]. In the following, we will focus on the heliocentric dynamics of Mercury and the EMB, which represent, in fact, the core of the MORE relativity experiment. Adopting a Lagrangian formulation, the equation of motion for Mercury and EMB are given by the Eulero-Lagrange equations: where r i , v i are position and velocity, respectively, of the i-th body (i = 1 for Mercury, i = 2 for EMB) in the Solar System Barycenter (SSB) reference frame and L is the Lagrangian of the problem. The Lagrangian can be decomposed, in turn, as where L N ew is the Lagrangian of the Newtonian N -body problem, L GR is the correction due to GR in the PN limit and L P P N includes the contribution due to PN and related parameters. In particular, the term L P P N describes how each parameter individually affects the dynamics and it assumes the form: The explicit expression of each term can be found in [21]. Let us recall that γ and β are expected to be equal to unit in GR, while η, α 1 and α 2 are all zero in GR. From Equation (1), the total acceleration acting on the i-th body can be derived, taking into account that the main term is the N -body Newtonian acceleration, a N ew i , while the other terms are small perturbations, and it takes the following approximated expression: whereB accounts for the acceleration of the SSB (see [21] for details).

II.2. Mathematical methods
The parameters of interest for the MORE relativity experiment are determined simultaneously by means of a global non-linear least squares (LS) fit. Following, e.g., [15] -Chap. 5, the non-linear LS fit aims at determining a set of parameters, u, which minimizes the target function: where m is the number of observations, W is the matrix containing the observation weights and ξ(u) = O −C(u) is the vector of the residuals, that is the difference between the observations O (i.e., the tracking data) and the predictions C(u), resulting from the light-time computation as a function of all the parameters u (see [23] for all the details). The procedure to compute the set u ⋆ of the parameters which minimizes Q is based on a modified Newton's method called differential correction method. First, we define the design matrix B and the normal matrix C as The stationary points of the target function are the solution of the normal equation: where ∆u ⋆ = u ⋆ − u. The method consists in applying iteratively the correction until convergence. From a probabilistic point of view, the inverse of the normal matrix, Γ = C −1 , can be interpreted as the covariance matrix of the vector u ⋆ (see, e.g., [15] -Chap. 3), carrying information on the attainable accuracy of the estimated parameters.

II.3. The ORBIT14 software
Since 2007, the Celestial Mechanics Group of the University of Pisa has developed (under an Italian Space Agency agreement) an orbit determination software, ORBIT14, dedicated to the BepiColombo and Juno radio science experiments (see, e.g., [24,25]). All the code is written in Fortran90.
In the case of BepiColombo, the software includes two separated stages: • the data simulator: awaiting for real data, it generates the simulated observables and the nominal value for the orbital elements of the Mercury-centric orbit of the spacecraft and the heliocentric orbits of Mercury and the EMB; • the differential corrector: it is the core of the code, solving for the parameters of interest by means of a global non-linear LS fit, within a constrained multi-arc strategy ( [26]).
A comprehensive view of the structure of the software can be found in, e.g., [21]. The software was successfully used for the analysis of the real Doppler data of the NASA Juno mission [27].

III. DYNAMICAL MODEL WITH TORSION
Following [7], we will consider a class of gravitational theories allowing for non-vanishing torsion based on a Riemann-Cartan spacetime. This spacetime is a four-dimensional manifold endowed with a Lorentzian metric g µν and an affine connection Γ λ µν such that ∇ λ g µν = 0 (mathematical details can be found, e.g., in [1]). The connection is uniquely determined by the metric and by the torsion tensor: In this case, the connection departs from the Levi-Civita connection (which is symmetric in the first two indices and holds in the Riemann spacetime of GR) by an additional anti-symmetric term, called the contorsion tensor, K µν λ , which cancels out in case of vanishing torsion tensor.

III.1. Spacetime with torsion in a PPN framework
In order to develop a model to test torsion within a PPN framework consistent with the dynamical model described in Section II.1, the metric and the torsion tensor need to be parameterized in a region of space at a distance r from the Sun such that the quantity ǫ µ⊙ = µ ⊙ /r ≪ 1, i.e., at large distances compared with the Schwarzschild radius. Since the PN framework we are adopting includes terms up to second order in the small quantity ǫ µ⊙ (see, e.g., [6,7]), the metric and the torsion tensor must be expanded up to the same degree of accuracy.
In general, the torsion tensor has 24 independent components, each being a function of time and position, but, adopting symmetry arguments and a perturbative approach up to O 2 (ǫ µ⊙ ), March et al. [7] showed that ultimately there are only 2 independent components of the torsion tensor, which can be parameterized by means of 3 independent parameters, the torsion parameters t 1 , t 2 , t 3 .
In a Riemann-Cartan spacetime two different classes of curves can be considered, autoparallels and geodesics curves, which reduce both to the geodesics of the Riemann spacetime when torsion vanishes [1]. In particular, autoparallels are curves along which the velocity vector is transported parallel to itself by the connection Γ λ µν , while along geodesics the velocity vector is transported parallel to itself by the Levi-Civita connection. In GR the two types of trajectories coincide while, in general, they may differ in the presence of torsion. Since geodesics curves in this generalized framework turn out to be the same as in the standard PPN framework (see [7] for details), new predictions related to torsion may arise only when considering autoparallel trajectories, which explicitly depend on torsion parameters.
After some mathematical manipulations and simplifications and accounting only for terms up to the second order in the small quantity ǫ µ⊙ , March et al. [7] achieved the final equation of motion for a test body moving along an autoparallel trajectory, written in rectangular coordinates: where α ∈ {1, 2, 3}, v 2 = 3 α=1 (ẋ α ) 2 and m is the mass of the central massive body (G = c = 1 has been assumed) [28]. We also recall that, differently from our approach, in [7] all the PN parameters other than γ and β are assumed to be zero, hence they are not included in the development and parameterization of the metric and the connection.

III.2. Implementation of torsion in ORBIT14
To properly account for the possible dynamical effects due to spacetime torsion, we need to add the torsion contribution to the global acceleration acting on the i-th body, described by Equation (2). In principle, we could directly use Equation (4), but, for the purposes of the MORE relativity experiment, we need to take into account the possible dynamical effects due to all the PN and related parameters we are interested in, since they all need to be determined simultaneously in the LS fit. Nevertheless, assuming that the only contribution to the Lagrangian L P P N is due to γ and β, Equations (2) and (4) must coincide apart from the terms due to torsion. This, in fact, is the case; thus, after some manipulation, we can isolate the contribution due to torsion on the motion of the i-th body, which reads [29]: where r i,⊙ ≡ r i − r ⊙ is the position of the body i with respect to the SSB and r i,⊙ = |r i,⊙ | is its module (the same notation holding for the velocity vector). Two issues stand out by looking at Equation (5). The first one concerns the fact that the contribution to the acceleration due to each of the three torsion parameters depends on the gravitational mass of the Sun, µ ⊙ , which is one of the parameters to be determined as well; in the case of t 1 , also a coupling with γ occurs. Thus, we need to properly take into account this dependency in the computation of the design and normal matrices, which must contain the partial derivatives of the new acceleration term with respect to µ ⊙ and γ. As a consequence, we can expect that a considerable correlation between µ ⊙ , γ and the torsion parameters could show up by solving simultaneously for all the parameters.
The second issue concerns specifically β and t 3 . Looking at Equation (4), it turns out that the dynamical effect due to these two parameters shows the same proportionality to the term r i,⊙ /r 4 i,⊙ . Thus, we can define an overall acceleration term due to the combined effect of both parameters, given by: As a consequence, in principle we cannot solve simultaneously for the two parameters within the LS fit, since the signatures due to each of the two effects on the observations cannot be distinguished in any way. From a computational point of view, this means that the design matrix would have two linearly dependent columns, so that the normal matrix would be degenerate, and therefore not invertible. This circumstance is usually known as rank deficiency and a general description of the issue in the case of an OD problem can be found in [15] -Chap. 6. Similar issues have been already dealt with in the past in the framework of the MORE relativity experiment (a detailed discussion on the approximate rank deficiency between β and the oblateness of the Sun, J 2⊙ , can be found in [11], while the general issue in the case of the MORE relativity experiment has been widely discussed by the authors in [19,25]). Rank deficiencies in an OD problem can be cured in two ways: solving for fewer parameters (this technique is known as descoping) or using more observations. The second way can be pursued by adding independent observations (e.g., using data from other experiments [30]) or, similarly, by adding a number of constraints equal to the order of the rank deficiency, which represent the available apriori knowledge of the problem (see, e.g., [19]). The first approach (descoping) has been applied, in fact, in the framework of the relativity experiment until now: assuming a vanishing torsion, as in GR, is equivalent to assuming that t 3 = 0; thus, looking at Equation (6), the actual effect on the dynamics is all due to the parameter β. This will be the strategy applied in the scenarios of the reference simulation and of simulation (a) described in Section IV.2. On the other hand, if we consider that a non-vanishing torsion may exist, and, in particular, that t 3 could be not null, the descoping approach cannot be adopted anymore. In this case we cannot separate the effects of β and t 3 and we can only estimate the linear combination (β − t 3 ). In we want to determine the two parameters, we need to add some further information on, at least, one of the two parameters. This can be done by adding an apriori constraint on the parameter β, given by the present knowledge on its value. If the apriori constraint is sufficiently tight, then the normal matrix can be inverted again and the problem can be solved. This approach will be adopted in the case of simulation (b) in Section IV.2.

IV. NUMERICAL ANALYSIS
In this Section, we describe the results of the numerical simulations of the MORE relativity experiment, including the dynamical effects due to torsion. In particular, we aim at checking whether the torsion parameters can be estimated by means of the experiment and, if so, at which level of accuracy.

IV.1. Simulation scenario and assumptions
First of all, we define a reference scenario. In this case, the LS fit aims at determining the relativity parameters listed in Section II, without accounting for possible effects due to torsion in the solution (see, e.g., [19]). We recall that the relativity parameters are: the PN parameters β, γ, α 1 , α 2 , the Nordtvedt parameter η, the oblateness of the Sun J 2⊙ , the gravitational mass of the Sun µ ⊙ , its time derivative ζ and the angular momentum of the Sun GS ⊙ .
The starting epoch of the simulation is set to 14 March 2026 and the nominal duration is set to 1 year. In Section IV.4 we will consider also possible benefits due to an extension of the mission duration up to an additional year. We assume that the radio tracking observables are affected only by random effects with a standard deviation of σ r = 15 cm at 300 s and σ rr = 1.5 × 10 −4 cm/s at 1000 s, respectively, for Ka-band observations. The software is capable of including also a possible systematic component to the range error model and to calibrate for it, but, since we are interested in a covariance analysis, we do not account for this detrimental effect here (this issue has been partially discussed in [31]).
For each relativity parameter we are interested in, we set an apriori constraint given by the present knowledge of the parameter, as listed in Table I. We point out that the adopted apriori on γ is a conservative estimate derived by a full set of simulations, carried out by the authors in [32], of the Superior Conjunction Experiment (SCE), planned during the cruise phase of BepiColombo [33]. Concerning the torsion parameters, they have not been experimentally estimated so far. In [7], referring to the present knowledge of β and γ, the authors derived the following constraints on the values of t 2 and t 3 : Moreover, we make use of an important assumption: we link the PN parameters by the Nordtvedt equation [34] which means that we are considering only metric theories of gravitation. The addition of the Nordtvedt equation in our model is motivated by the fact that β and J 2⊙ are expected to show an almost 100% correlation, since their dynamical effect on the orbit of Mercury is comparable (see, e.g., [11] for an extensive discussion of this symmetry in the case of the MORE relativity experiment). Thus, the addition of the constraint removes the degeneracy between β and J 2⊙ , but this result is, in turn, obtained at the cost of forcing an almost 100% correlation between β and η [35].

IV.2. Simulation results
In the following, we denote as "reference simulation" the case in which the solution is computed in the standard reference scenario described in Section IV.1, which consists in estimating the state vectors of Mercury and the EMB at the central epoch of the mission (6 + 6 parameters) and the relativity parameters, for a total of 21 parameters. Then, we define as simulation (a) the case of the reference scenario with the addition of the torsion parameters t 1 and t 2 in the solve-for list, for a total of 23 estimated parameters, and we estimate (β − t 3 ) in place of β. Finally, we label as simulation (b) the case of a simulation in which all the three torsion parameters, t 1 , t 2 , t 3 , are simultaneously determined by means of the LS fit together with the relativity parameters and the state vectors, for a total of 24 parameters. In addition, in each scenario the state vector of the spacecraft (position and velocity, 6 parameters) is estimated at the central epoch of each of the 365 observed arcs, for a total of 2190 parameters. These parameters are handled as local parameters (i.e., they belong only to a single arc: see, e.g., [26] for details) and they result uncorrelated with the relativity parameters, which are global parameters. For the sake of clarity, the global parameters estimated in each of the three scenarios are summarized in Table II. The results in terms of formal accuracies (1-σ) determined in the three different scenarios are compared in Table  III for what concerns the state vectors (position and velocity in cartesian coordinates) of Mercury (labeled with the index M ) and the EMB (labeled with the index E), while in Table IV the corresponding results for the relativity parameters are shown.
Concerning the determination of the state vectors at the central epoch, shown in Table III, we can immediately observe that, as expected, the results are the same in all scenarios, suggesting that the there is no correlation between the torsion parameters and any of the components of the orbital conditions of Mercury and the EMB, i.e., solving for the torsion parameters does not weaken the determination of the state vectors. The two orbits can be determined TABLE IV. Formal accuracies expected for the relativity parameters in the case of: reference simulation, simulation (a) and simulation (b), respectively. In the last column, we introduce a reference number, N , which labels each parameter, that will be adopted in Section IV.3. Note that σ(µ⊙) is expressed in cm 3 /s 2 , σ(ζ) in y −1 and σ(GS⊙) in cm 3 /s 2 .

Parameter
Reference with the MORE experiment at the level of some meters, for what concerns the position, and of some m/s, concerning the velocity, with a general trend of a better estimate in the x − y plane, while the z-direction seems weaker. We point out that no constraint has been imposed on the two orbits (contrarily to what the authors have done in [19]): the issue of possible symmetries affecting the geometry of our experiment is very critical and a dedicated paper will be devoted by the authors to this specific topic in the near future.
Moving the attention to Table IV, first of all some observations on the results achievable in the case of the reference scenario are in order. Comparing the estimated accuracy with the present knowledge, shown in Table I, we can expect a significant improvement in the knowledge of γ, α 1 , α 2 and J 2⊙ thanks to the MORE experiment. A minor improvement could be achieved in the estimate of the gravitational mass of the Sun, µ ⊙ , while no significant improvement is expected, at the moment, for the knowledge of ζ (see [21] for a discussion about this parameter) and GS ⊙ (see [19] for a further discussion).
Concerning the torsion parameters, which are the focus of this work, we can observe that in the case of simulation (a) we achieve the relevant result that the torsion parameters t 1 and t 2 can be estimated by means of the MORE relativity experiment at the level of 10 −5 , which could place a significant constraint on the validity or not of torsion theories. In the case of simulation (b), i.e. estimating also the parameter t 3 , an order of magnitude is lost in the determination of the torsion parameters. This is mainly a consequence of the strong correlation between β and t 3 : since the three torsion parameters are expected to be correlated one with the other (see Equation (5)), the worsening affects the estimate of all the three parameters. Anyway, they can be estimated at the level of some parts in 10 −4 and this is still a significant result.
The addition of the torsion parameters in the solve-for list worsens the determination of most of the relativity parameters, slightly more in the case of simulation (b) with respect to simulation (a). The worsening is more pronounced in the case of the parameter µ ⊙ , whose accuracy deteriorates by a factor 3 in simulation (a) and by a factor 4 in simulation (b). This result is not surprising since the dynamical effect due to torsion depends on µ ⊙ , as pointed out in Section III.2. In particular, in the case of simulation (b), the accuracy of the parameter µ ⊙ is comparable to the present knowledge and this fact suggests that, in the absence of an apriori on µ ⊙ , it would be difficult to solve for this parameter simultaneously with the torsion parameters. In the case of the parameters β and η, the worsening is limited by less than a factor 2. This fact may seem contradictory in the case of simulation (b). Indeed, in Section III.2 we have seen that the dynamical effect due to β and t 3 is the same and, as a consequence, we could expect that the addition of t 3 in the solve-for list should cause a significant worsening in the determination of β. This is, in fact, the case: if we do not constrain the value of β in some way and we simultaneously estimate t 3 , we would find that the normal matrix is not invertible. The problem is solved by adding the present knowledge on β as an apriori constraint on the parameter: we are not able to improve the present knowledge of β (by looking at Table IV -simulation (b), the obtained formal accuracy is equal to the present knowledge), but at least we are able to separate its effect from that of t 3 .
A special discussion deserves the issue of the determination itself of β and η. We recall that the two parameters are highly correlated since they are linked through the Nordtvedt relation. In the past, we observed that the determination of β and η strongly depends on the level at which the experiment is capable to determine the orbits of Mercury and the EMB. In particular, referring to the results shown in [19], we observed that the determination of the orbits of Mercury and the EMB with an accuracy below the meter and m/s for the position and velocity, respectively, leads to a great benefit in the estimate of β and η (cfr. Tables 2 and 3 in [19]). The problem is that an accuracy at the level of 10 cm or less in the determination of the two orbits seems unrealistic, due to a number of detrimental effects, in particular due to the uncertainty in the knowledge of the masses and state vectors of the bodies in the solar system, which influences in turn the position and velocity of Mercury and the EMB (mainly in the case of Jupiter, see [41] for an extensive discussion). Indeed, we point out that, for the purpose of the experiment, the mass and position of the bodies included in the N -body Newtonian acceleration (see Equation (2)) are taken from the JPL ephemerides, thus they are assumed as perfectly known. As a consequence, we cannot expect to significantly improve our knowledge on β and η with the relativity experiment set up in the way described here. Different choices can be made, considering a wider set of independent observations due to different space missions to be fitted simultaneously, as proposed in [30].

IV.3. Analysis of the correlations
To carry out a complete description of the results, the correlations showing up between the estimated parameters need to be analyzed. We recall that the correlation of the parameter A with the parameter B, ρ(A, B), where A and B are, respectively, the i-th and the j-th parameters in the solve-for list, is proportional to the out-of-diagonal Γ ij term of the covariance matrix Γ (see Section II.2) and it is normalized to 1 by means of the formal accuracy of the two parameters, σ(A) and σ(B).
A global view of the behavior of correlations between the parameters is given in Figures 1, 2 and 3 for the cases of the reference simulation, simulation (a) and simulation (b), respectively. The estimated parameters are labeled  Tables III  and IV. with a number N = 1, ..., 24 as in Tables III and IV: indices from 1 to 12 refer to the state vectors of Mercury and the EMB (position and velocity components), from 13 to 21 to the relativity parameters, in the order indicated in Table IV, and, finally, from 22 to 24 to the three torsion parameters. The color-bar on the right of each figure spans from the value 0 (dark blue) in the case that no correlation is found to the value 1 (yellow) in the case that an exact correlation shows up.
In the case of the reference simulation (Figure 1) a two-block structure can be clearly recognized, which can be found also in the cases of simulation (a) and (b). This structure reveals the fact that the subset of parameters 1-12, i.e., the components of the state vectors, and the subset of the parameters 13-21, i.e., the relativity parameters, are almost not correlated one with the others, while one parameter can be highly correlated only with the other parameters of the same subset. In particular, we can observe that the correlation between β and η is ρ(β, η) = 0.99 in all the three cases, as expected. Moreover, GS ⊙ shows a non-null correlation with both J 2⊙ and µ ⊙ , equal to ρ(J 2⊙ , GS ⊙ ) = 0.58 and ρ(µ ⊙ , GS ⊙ ) = 0.55. The issue of the correlation showing up between J 2⊙ and GS ⊙ is described in [19]. Finally, J 2⊙ shows a residual correlation with β, and consequently with η, equal to ρ(β, J 2⊙ ) = 0.59 and ρ(η, J 2⊙ ) = 0.59 and also µ ⊙ turns out to be correlated with these two parameters, at the level of ρ(β, µ ⊙ ) = 0.65 and ρ(η, µ ⊙ ) = 0.65.
In the case of simulation (a), shown in Figure 2, the behavior of correlations is similar to the reference case. Now that the parameters t 1 and t 2 are included in the solve-for list, it shows up a non-null correlation between t 1 and µ ⊙ , equal to ρ(µ ⊙ , t 1 ) = 0.85, as it was expected by looking at Equation (5). In the case of t 2 , instead, the correlation with µ ⊙ is only ρ(µ ⊙ , t 2 ) = 0.46. Conversely, almost no correlation shows up between the torsion parameters and 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Relativity parameters   24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6 Tables III  and IV. γ: this can be ascribed to the presence of a tight apriori on γ, as expected from the SCE during the cruise phase. Moreover, t 1 is also partially correlated with β and η, at the level of ρ(β, t 1 ) = 0.63 and ρ(η, t 1 ) = 0.63.
Finally, in the case of simulation (b), shown in Figure 3, the full correlation between the three torsion parameters shows up. If we look at Equation (5), we can see that, as long as we assume t 3 = 0 (simulation (a)), the effects due to t 1 and t 2 are, at least partially, independent, thus a low correlation between t 1 and t 2 is expected; indeed, we find that their correlation is below 0.1. When also t 3 comes into play, since its effect is proportional to one of the terms due to t 1 (see Equation (5)), we find that ρ (b) (t 1 , t 3 ) = 0.99 and, as a consequence, ρ (b) (t 3 , t 2 ) = 0.99 and ρ (b) (t 1 , t 2 ) = 0.98. The correlation between the torsion parameters and the gravitational mass of the Sun, µ ⊙ , persists also in this case, at the level of ρ(µ ⊙ , t 1 ) = 0.53, ρ(µ ⊙ , t 2 ) = 0.66 and ρ(µ ⊙ , t 3 ) = 0.62, respectively.

IV.4. Possible benefits from an extended mission
The nominal duration of the BepiColombo orbital phase is fixed to one year, with the scientific operations in orbit starting on March 2026. Anyway, the possibility of an extension up to one further year is envisaged. Thus, to check the possible benefits due to a 2-year duration of the experiment, we re-ran the simulations in the three scenarios described in Section IV.2. The results on the achievable 1-σ accuracies for the relativity experiment are shown in Table V, where for each case we show also the improvement factor, R 1-to-2 , defined as the ratio of the 1-year accuracy over the 2-year accuracy. We can observe that the main benefits of an extended mission concern the determination of β, η and α 1 , whose accuracies can be improved almost by a factor 5. This fact can be easily understood considering that the effects of these parameters on the planetary dynamics are more visible over long time scales, hence the benefit of a longer mission. On the other hand, γ and µ ⊙ would not benefit in any way from a duration longer than one year. Regarding the torsion parameters, their knowledge could be improved by a factor 2 by considering a 2-year long orbital phase.

V. DISCUSSION AND CONCLUSIONS
In this paper we have shown how to account for a possible non-vanishing spacetime torsion in the framework of the MORE relativity experiment. In particular, we have described the implementation of the updated dynamical model within the ORBIT14 software, developed by the Celestial Mechanics Group at the University of Pisa. The implementation has been done, in turn, by parameterizing the torsion effects by means of three parameters, t 1 , t 2 , t 3 , to be estimated in a global LS fit in the same way as the PN and related parameters, which characterize the "classical" MORE relativity experiment of BepiColombo.
We have shown that the torsion parameters can be estimated with MORE, at least, at the level of some parts in 10 −4 , which is a remarkable result in order to constrain torsion theories of gravitation. Our analysis showed that the main limitations in estimating the torsion are due to the correlation of the torsion parameters with the gravitational mass of the Sun, µ ⊙ , and of t 3 , in particular, with the PN parameter β. This issue can be overcome by adding an TABLE V. Formal accuracies expected for the relativity parameters considering a 2-year long orbital experiment, in the case of: reference simulation, simulation (a) and simulation (b), respectively. For each case we show the 1-σ accuracy (columns 2, 4, 6) and the improvement factor, R1-to-2, defined as the ratio of the 1-year over the 2-year accuracy (columns 3, 5, 7). Note that σ(µ⊙) is expressed in cm 3 /s 2 , σ(ζ) in y −1 and σ(GS⊙) in cm 3 /s 2 .

Parameter
Reference R1-to-2 Simulation (a) R1-to-2 apriori constraint on µ ⊙ and β. The consequence is that, when all the three torsion parameters are estimated in the global LS fit (i.e., simulation scenario (b)), β and µ ⊙ cannot be determined at a better level with respect to the apriori value. Thus, we can conclude that the torsion parameters can be estimated at a significant level of accuracy by the MORE relativity experiment, but this can be done, in turn, by relaxing the requirements on the estimate of β and µ ⊙ . On the other hand, any independent improvement in the knowledge of β and µ ⊙ could be adopted as an apriori on these parameters within our fit and it could allow furher improvement on the determination of the torsion parameters. To quantify this point in the case of β, we ran two illustrative simulations starting from the scenario of simulation (b) but we tightened the apriori constraint on β up to 1 × 10 −5 (case I) and 5 × 10 −6 (case II). The results are shown in Table VI, where the case of simulation (b) is compared with case I and case II for what concerns the achievable accuracies on β, t 1 , t 2 , t 3 . Comparing the results of simulation (b) with respect to case I and II, we can observe that a possible independent improvement on the knowledge of β up to 5 × 10 −6 would not allow any significant improvement in the determination of the torsion parameters. This result suggests that an accuracy of some parts in 10 −4 is the ultimate level of accuracy by which the MORE relativity experiment itself is capable to constrain the torsion parameters.

V.1. Acknowledgments
The results of the research presented in this paper have been performed within the scope of the contract "addendum n. 2017-40-H.1-2020 all'accordo attuativo n. 2017-40-H" with the Italian Space Agency (ASI).