Constraints on the interior structure of Phobos from tidal deformation modeling

Phobos is one of the two moons of Mars, the only other terrestrial planet to have companions. In spite of certain common orbital characteristics between Earth’s Moon and Phobos, such as almost circular, low-inclination, and synchronous orbits, the origin of the Martian moons remains to be understood. This unsettled state of affairs reflects a scarcity of data pertinent to its interior structure and bulk composition, both of which hold the key to solving the question of the provenance of the satellites. Drawing on the importance of interior structure as a means of providing further constraints on this problem, we construct a series of models of the interior of Phobos in accordance with current observations and determine their tidal deformation. The models include a homogeneous body, ice-rock mixtures, models with positive and negative density gradients, and layered models. To compute the deformation of Phobos precisely, we solve the three-dimensional elastostatic problem to obtain the full displacement field. For this, we rely on a higher-order spectral-element method and to account for a correct representation of shape and resulting displacement field, we accurately mesh the figure of Phobos by employing the digital terrain model of Willner et al. (2014). To enable us to further distinguish between models, we also rely on currently available geophysical data (e.g., magnitude of libration in longitude, mean density, gravity field, and moments of inertia). Our results show that the homogeneous, icerock mixture, and negative density gradient models are largely degenerate and therefore difficult to separate, but that models with a density increase with depth can be differentiated from the former model families. From a purely deformational point of view, we find that for positive gradient and layered models, the largest deformation occurs on the rim of the Stickney crater, rather than immediately at the sub-Mars point on Phobos, as in the case of the other models, where tidal forces are largest. These observations will be of considerable interest and should provide important anchoring points for the future exploration of Phobos as envisaged with the Martian Moons Exploration mission.


Introduction
Phobos is the innermost satellite of Mars and is in a 1:1 spinorbit resonance as a result of which it shows, on average, the same hemisphere to Mars similarly to our Moon in its orbit about Earth. Phobos revolves around Mars in 7 h 39 min in an almost circular orbit that lies in the equatorial plane with an inclination of ∼1 • . The average distance to Mars's surface is ∼9375 km, which is well below its corotation radius (the distance at which Phobos's mean motion equals Mars's spin rate). The orbital parameters of Phobos are summarized in Table 1.
The satellite is highly irregular in shape and, as shown in Fig. 1, elongated in the direction towards Mars . The Like Earth's Moon, the origin and provenance of Phobos and Deimos remains elusive, yet provide fundamental clues on the dynamical evolution of the inner solar system. Chief among the hypotheses for the origin of the Martian moons figure asteroidal capture (Cazenave et al., 1980;Burns, 1992;Pajola et al., 2013), disruption of a common parent body (Bagheri et al., 2021), co-accretion (Safronov et al., 1986), or post-collisional re-accretion of ejecta material (Craddock, 2011;Citron et al., 2015;Rosenblatt et al., 2016;Hyodo et al., 2017;Canup and Salmon, 2018;Hesselbrock and Minton, 2017). Yet, as reviewed by e.g., Rosenblatt (2011) and Rosenblatt et al. (2019), ambiguities related to interior structure and composition make it challenging to discriminate between current hypotheses. To shed further light on this problem, additional information on the internal structure of Phobos and/or its composition is required (e.g., Usui et al., 2020).
Currently available information on the interior include satellite mass (Pätzold et al., 2014a) and shape model , which constrains the mean density of Phobos to 1.845 g/cm 3 . This relatively low value suggests a fractured and porous moon filled with either water-ice or voids or a combination thereof. Based on the estimated mean density of Phobos, a number of interior structure models, including homogeneous, fractured, or layered models of Phobos, have been investigated with the purpose of estimating bulk physical properties, including density, degree-2 gravity coefficients, libration amplitude, and moments of inertia (Borderies and Yoder, 1990;Chao and Rubincam, 1989;Pätzold et al., 2014a;Rambaux et al., 2012;Willner et al., 2014;Witasse et al., 2014;Shi et al., 2012;Tian and Zheng, 2020;Rosenblatt, 2011;Le Maistre et al., 2013Yang et al., 2019;Guo et al., 2021).
One of the principal sources of information about the interior of the satellite, which we shall rely on here, is its response to Martian tides (see also Le Maistre et al. (2019)). The tidal response of a body is determined to first order by internal properties (e.g., rigidity and density). Here, we focus on the tidal bulge height as a possible means of determining the interior structure of Phobos. The tidal bulge height is a measure of the radial displacement of a point on the surface of Phobos, and the less rigid the satellite is, the larger the resulting bulge height. If the displacement could be measured from orbit within the framework of the upcoming Martian Moons Exploration (MMX) mission (Campagnola et al., 2018;Usui et al., 2020), a relatively straightforward means of assessing interior properties could possibly be established. With the envisaged Phobos sample return as part of MMX, knowledge of its interior represents a key piece of information that is needed for improved understanding of the origin of the Martian moon(s).
In this study, we compute the gravitational potential and tidal deformation of Phobos for a series of possible models of its interior structure with a focus on locating regions of maximum displacement so as to configure a potential orbiter-satellite system for optimal return. For this purpose, we solve the full three-dimensional (3D) elastostatic problem on the deformed body to determine the displacement field and Poisson's equation for the gravitational potential. For the case of 3D models, the governing equations can only be solved numerically as attested to in previous studies (e.g., Shi et al., 2012;Le Maistre et al., 2019;Tian and Zheng, 2020). Shi et al. (2012) compared different methods including harmonic expansion, voxelization, and polyhedral approximations to compute the gravitational potential for a homogeneous density distribution. Le Maistre et al. (2019) considered a voxelized model of Phobos for a range of mass distributions and computed moments of interia, libration amplitude, and surface gravity. Tian and Zheng (2020) employed the boundary element method (BEM) for modeling the surface deformation of irregular bodies. BEM solves the governing equations on the boundaries of the body, for which the number of unknowns is reduced, but the resulting discrete systems are usually dense and parallel scalability becomes challenging. Moreover, BEM does not directly provide the volumetric fields for the gravitational potential or the displacement in the interior of Phobos.  . The big crater facing the reader is the Stickney crater.
Instead, our approach relies on the higher-order spectral-element method (Patera, 1984;Chaljub et al., 2007) and a recently developed extension to the waveform modeling package Salvus (Afanasiev et al., 2019). The spectral-element method (SEM) is particularly efficient on unstructured hexahedral meshes and suitable for iterative, matrix-free linear solvers. To properly mesh the shape of Phobos, we employ the spherical harmonics expansion of the digital terrain model of Willner et al. (2014), which results in an accurate representation of shape and resulting displacement field. Given the advantages of our numerical approach that enable us to achieve a high level of accuracy in geophysical modeling, we investigate the possibility of distinguishing between first-order, e.g., homogeneous, ice-rock mixture, and layered models using currently available observations that include the degree-2 gravitational coefficients, the libration in longitude, the principal moments of inertia, and the center-of-mass-center-of-figure offset, in addition to the aforementioned tidal displacement field.
The manuscript is arranged as follows. The theoretical background for the elasto-gravitational problem and the geophysical data pertinent to Phobos are discussed in Section 2. The numerical implementation of the spectral-element method to solve the governing equations of the elasto-gravitational problem is delineated in Section 3. In Section 4, we devise a number of interior structure models of Phobos for analysis, whose tidal, gravitational, and geophysical response is presented and discussed in Section 5.

Theoretical background
The gravitational potential of Phobos can be computed as a solution to Poisson's equation where is the 3D density model of Phobos and is the gravitational constant. From the gravitational potential, we can determine the gravitational acceleration at the surface of Phobos as the magnitude of the normal derivative of . It is important to note that the spatial domain of Eq. (1) involves the full space R 3 , which requires special treatment in the case of numerical approaches (to be outlined in Section 3). For the displacement field resulting from tidal forces, we consider Phobos to be a purely elastic body that experiences periodic deformations. Consequently, the problem can be described by the momentum equation in the Lagrangian formulation (Dahlen, 1974) where is density, is displacement, is stress, is an external forcing term, and is the perturbation of the gravitational potential due to the redistribution of mass, which is determined by Hooke's law relates stress to strain where is the fourth-order elasticity tensor, ( ) = 1 2 (∇ + ∇ ) is the symmetric strain tensor and ∶ denotes contraction over adjacent indices. The constitutive relation in Eq. (4) is valid for general anisotropic media. For isotropic elastic media, which we consider in this study, the elasticity tensor reduces to the following form where and are bulk and shear moduli, respectively. Transforming Eq.
(2) to the frequency domain, we obtain 2 ( ) = ∇⋅( ∶ ( ( )))−∇( ( )⋅∇ )+∇⋅( ( ))∇ − ∇ ( )+ ( ) where is frequency. Because we only consider the purely elastic case, ( ) and are real-valued, implying no phase lag and therefore no frequency dependence. Consequently, Eq. (6) simplifies to the elastostatic problem The last term on the left-hand side ∇ is found to be at least 4 orders of magnitude smaller than the tidal force and, consequently, neglected in our simulations (cf. Fig. A.1). The force term balances the gravitational attraction of a pointlike Mars ♂ with the orbital centrifugal force . Here we do not consider the force arising from Phobos's rotation around its own axis because the impact on the resulting deformation is found to be negligibly small (∼1 ). The resultant force field can therefore be written as where ♂ is Mars's mass, ( ) is the density of a volume element d within Phobos, P is a vector connecting the centers of mass (CoMs) of Mars and Phobos that are considered to coincide with the center-offigures (CoFs) of the corresponding bodies, and is a vector connecting the volume element d within Phobos with the CoM of the satellite. Since Phobos rotates around the barycenter of the Mars-Phobos system as a solid body, the centrifugal force is the same for every unit mass inside the satellite (Chapter 4 in Murray and Dermott (2000)) and can be obtained from its orbital velocity orb for any position of Phobos along its orbit where is the distance between the CoFs of Mars and Phobos. The orbital velocity orb of Phobos can be found from the fact that Mars's gravity is always equal to the centrifugal force in the center of mass of the satellite's figure where the Mars-Phobos distance is determined from the semi-major axis ( ), orbit eccentricity ( ), and true anomaly ( ) using (e.g., Murray and Dermott, 2000) = (1 − 2 ) 1 + cos In contrast to , the gravitational attraction of Mars changes across the satellite. The tidal force tidal results from the difference between these two forces. Consequently, the magnitude of tidal increases with distance from the body's rotation axis as indicated by the length of the green arrows in Fig. 2.
Finally, while the solution to the elastostatic problem provides the full displacement field within the satellite, this static tidal deformation is not of practical interest as it is not measurable by e.g., an orbiter, because the satellite is always in a deformed state to some degree. However, the magnitude of the deformation changes over the orbit due to the changing distance to Mars.
Thus, to remove the static tidal deformation, which is caused by a static force, we compute the tidal force tidal by subtracting the force fields at periareion peri and apoareion apo from which the resultant tidal displacement field can be determined (Fig. 3

Geophysical properties of Phobos
Geophysical properties of Phobos pertinent to its interior structure that have been measured include shape, volume, mass, degree-2 gravity coefficients, and amplitude of longitudal libration (Table 2). In the following, we briefly describe the observations and how we intend to compute these.  and are the displacement fields at apoapsis and periapsis points, respectively, is the semi-major axis of the orbit, and is its eccentricity. Subtraction of one field from another allows to get rid of the static displacement and to pick out only the time-dependent displacement field. Note that the deformation of Phobos at periareion and apoareion are strongly exaggerated. The light gray Phobos illustrates the concept of longitudal libration ( ).
signifies the true anomaly and is the angle between the major orbital axis and the satellite's long axis. See main text for more details.
Volume and shape model. A body's volume is determined from its shape. Currently, there are two ways of representing the shape of a body, either through the use of continuous spherical harmonics (SPH) or through a discretized digital terrain model (DTM). The DTM requires a discretized set of surface points with known coordinates. The DTM can be obtained by different methods including various image analysis techniques such as control point analysis (Duxbury, 1991;Oberst et al., 2014;Willner et al., 2014;Gaskell, 2020), image matching, conjugate points and observations of limb profiles (Duxbury, 1974;Thomas, 1989;Simonelli et al., 1993). On the basis of the discrete satellite shape model, an analytic SPH expansion can be constructed, which describes a body's shape as an infinite set of continuous functions in spherical coordinates. In this method, the distance to the center of figure (CoF) ( , ) of a surface point at longitude and latitude is given by (Willner et al., 2010): where (sin ) are the associated Legendre polynomials of degree and order and and are expansion coefficients determined by Willner et al. (2014) up to degree and order 45. The volume based on the most recent shape model of Phobos ( Fig. 1) that derives from analysis of images from the high resolution stereo camera onboard Mars Express is 5742 ± 35 km 3 .
Mass. The mass of Phobos P can be determined from observation of the trajectory of an orbiting Mars spacecraft, as the latter feels the gravitational attraction of the satellite. The strength of the orbit perturbation can be estimated from the long-term secular drift of the satellite orbit and short-term velocity changes during close flybys of the moon, which are detected by Doppler shift (Andert et al., 2010). Radio experiment data from the MEX (MaRS Express) mission were employed for the most recent estimations of P , where P is the mass of Phobos that yields the best-fit P = 0.7072 ± 0.0013 × 10 −3 km 3 s −2 (Pätzold et al., 2014b), corresponding to a P = 1.0596 ± 0.00195 × 10 16 kg.

Physical librations and moments of inertia.
Phobos is in a 1:1 spin-orbit resonance around Mars, which implies that its mean motion = 2 ∕ , where is the orbital period, is equal to its spin rate (Lainey et al., 2007). As the orbital speed of Phobos varies along the eccentric orbit, while the rotation rate remains constant, the long axis of the satellite oscillates around the direction towards Mars's center of mass. These oscillations are called the forced longitudal librations and illustrated in Fig. 3. In addition to this, Mars raises a time-dependent gravitational torque on Phobos due to the difference between the satellite's principal moments of inertia, which causes free longitudal librations. As a result of the libration, 52% of Phobos surface is presented to Mars (Rambaux et al., 2012). The amplitude of the longitudal libration can be estimated as a difference between the true anomaly and the angle between the satellite's long axis and its major orbital axis ( Fig. 3). Because Phobos's equatorial plane almost coincides with its orbital plane, the magnitude of the oscillation in longitude is the biggest and, therefore, the easiest to detect. Latitudal librations can be ignored in the case of Phobos (Borderies and Yoder, 1990).
From a geophysical point of view, the longitudal libration is of particular interest because it is determined by the moments of inertia (MoIs) of the body. MoIs, in turn, are defined by radial moments of the density distribution. In general, the MoI is a second-order symmetric tensor with the following elements is a vector to a point within the body whose density is ( ) and are components of the symmetric inertia tensor of that same volume element. To first order, the amplitude of the main forced libration in longitude is (Borderies and Yoder, 1990): where 1 < 2 < 3 are the principal MoIs, which are the eigenvalues of the aforementioned inertia tensor . We should note that Eq. (16) really only applies to the case of a homogeneous Phobos, i.e., when the off-diagonal elements of the inertia tensor are zero. However, as the off-diagonal elements of are generally 4-5 orders of magnitude smaller than the principle elements, we can neglect the former in line with usual practice (Le Maistre et al., 2019). The forced libration amplitude of Phobos has been estimated using e.g., imagery analysis with cross-point networks (Duxbury, 1974;Duxbury and Callahan, 1989;Burmeister et al., 2018), analytical theories (Chapront-Touze, 1990;Borderies and Yoder, 1990), and numerical integration of the equations of motion (Jacobson, 2010;Rambaux et al., 2012;Yang et al., 2019). The most recent estimation based on Viking 1 Orbiter and Mars Express images (Burmeister et al., 2018) have resulted in an amplitude of 1.143 ± 0.025 • .
For comparison, Jacobson (2010) and Rambaux et al. (2012) numerically integrated Phobos's orbit, while varying the libration amplitude so as to fit the observed orbital motion. This resulted in a longitudal libration amplitude of 1.03 ± 0.22 • .
The 2-degree gravity coefficients ( 20 and 22 ), which also contain information on the interior density distribution, have been estimated by Jacobson and Lainey (2014) by numerically integrating the satellite orbit with the aim of obtaining an improved estimation of the ephemerides. A limitation of the method of Jacobson and Lainey (2014) is that 20 , 22 , and cannot be determined independently. Instead (Jacobson and Lainey, 2014) employ = 1.2 • from Willner et al. (2010) and are thus able to constrain 20 and 22 (see Table 2). From the principal MoIs, the unnormalized degree-2 gravity ( 20 and 22 ) coefficients of the body can also be determined where are the normalized principal moments of inertia and P is the radius of the sphere whose volume is identical to that of Phobos.
CoF-CoM offset. Of further interest in connection with interior structure models is the position of the CoM with respect to the CoF of the satellite. The CoF coordinates CoF are determined from the shape model, whereas CoM coordinates are determined by the density distribution within Phobos and are found by numerical integration over Phobos's volume P where ( = 1, 2, 3), , and d are Cartesian coordinates, density, and volume of an interior element, respectively. Alternatively, the displacement of the CoM can be determined from the degree-1 gravity coefficients using the following relations (Le  Maistre et  where ′ is the reference radius used for the computation of the gravity coefficients 11 , 10 , and 11 . More generally, and are the cosine and sine coefficients of the spherical harmonics expansion of the gravity field of degree and order , similarly to and in Eq. (14). For reference, the CoM of the uniform density model (to be presented in Section 4) has the predicted coordinates in the Cartesian system (−66.43 m, 63.0 m, −84.31 m), while the CoF is located at (0 m, 0 m, 0 m).

Numerical approach
We apply the spectral-element method (SEM) to solve the governing equations. SEM is a high-order finite-element method, and as such operates on the weak form of the governing equations.
(1), is defined on the full space R 3 . In order to solve it numerically with SEM, however, we need to consider a bounded domain together with a suitable boundary condition. Because the gravitational potential decays rapidly with distance away from Phobos, using a large domain surrounding Phobos and applying homogeneous Dirichlet conditions for at its boundary yield sufficiently accurate approximations. Using integration by parts, we obtain the weak form of Eq. (1) on the extended domain as where is a scalar test function satisfying = 0 on . The weak form of Eq. (7) is derived in a similar manner. Here, the computational domain is restricted to the volume of Phobos, and we impose a boundary condition at its surface, such that the stress in the direction normal to the body vanishes. By applying a test function , integrating by parts, and employing the Gauss-Green theorem to the left-hand side of Eq. (7), we obtain the weak form where the surface integral vanishes because of the free-surface boundary condition. Note that Eq. (22) involves scalar fields and for the gravitational potential and the test function, respectively, while Eq. (23) requires vector fields and . Despite this difference, however, both weak forms obey a similar structure involving a bi-linear form including first-order spatial derivatives on the left-hand-side and a volume integral over P on the right-hand-side. Hence, both equations can be treated similarly in our numerical approach, and we focus on Eq. (23) in the following.
The spectral-element method (Patera, 1984;Chaljub et al., 2007;Afanasiev et al., 2019) discretizes the weak form using an unstructured mesh that subdivides the domain into hexahedral elements. The local solution field on each element is approximated using a high-order Lagrange polynomial interpolating the Gauss-Lobatto-Legendre (GLL) collocation points. With the Galerkin projection of the weak form onto the spectral-element basis ( ) =1,…, , the bi-linear form in Eq. (23) results in a linear equation for each , as a result of which the discretized weak form can be written as a linear system where is the discretized displacement field and is the projection of the force field onto the spectral-element basis. Note that the spatial derivatives can be evaluated efficiently due to the tensorized structure of the spectral-element basis, and these integrals are computed exactly using the Gauss-Lobatto quadrature rule. The stiffness matrix in Eq. (24) depends on the elastic bulk and shear moduli. Variations in density enter on the right-hand-side in , which also holds for the discretized Poisson equation. The stiffness matrix is sparse and will not be assembled in practice. Instead, we apply an iterative preconditioned conjugate gradient solver to determine the three-dimensional displacement field. Using a diagonal Jacobi preconditioner, all computations can be carried out matrix-free, and with a very small memory footprint.
The key ingredient is to accurately mesh the shape of Phobos using the shape model from Willner et al. (2014). Starting from a cubed sphere (Ronchi et al., 1996), we utilize the spherical harmonics expansion of the digital terrain model to apply the shape deformation. Both the finite-element basis and the shape transformation of the hexahedral elements use Lagrange polynomials of order 4 defined on the GLL points. This provides a very accurate representation of the shape model and the resulting displacement field with a decent number of elements. As explained above, Poisson's equation additionally requires an extended domain . Using mesh coarsening techniques it is possible to discretize the outer domain with a relatively small number of elements, and without the need of changing the finite-element basis or quadrature rule. Examples of the resulting meshes are shown in Fig. 4.   Table 3). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) All computations were carried out with Salvus (Afanasiev et al., 2019) on a mesh that comprises 65 856 elements with 125 nodes per element representing fourth-order polynomials. The pure simulation time required to solve the elastostatic problem on this mesh and employing 12 cores is ∼9.65 min (in 1445 iterations) and therefore accomplishable on a regular laptop.

Models of Phobos's interior structure
Imaging of the surface properties of Phobos has revealed its irregular shape, cratered surface, and low albedo that suggest a certain kinship with primitive outer-belt asteroids (Thomas et al., 2011;Pieters et al., 2014;Fraeman et al., 2014;Pajola et al., 2013); yet, we lack constraints on its composition and internal structure. The low density of Phobos is generally less than that of meteoritic material (the density of carbonaceous chondritic material of type CI is 2.42 ± 0.06 g/cm 3 Macke et al., 2011), indicative of a porous body that, in addition to voids, could also contain water ice in its interior (Fanale and Salvail, 1989;Rosenblatt, 2011;Rosenblatt et al., 2019;Murchie et al., 2013). Phobos is likely to be either a moderately porous aggregate or a rubble pile. The surface of Phobos is heavily cratered, the most pronounced of which is the Stickney crater (see Fig. 1). The Stickney cratering event would most likely have shattered Phobos completely had it been a monolith or a complete rubble pile, but leave a weakly-connected Phobos intact (Asphaug et al., 1998). The extent of the porous nature of Phobos, however, remains ambiguous, although there is evidence that points to Phobos being a rubble aggregate with some type of internal strength (Pätzold et al., 2014a), similar to estimations for asteroids  Table 3). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) For the case studied here, shear modulus is 0.179 GPa and surface density is 1.5 g/cm 3 (see Table 3). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Eros, Ryugu, Bennu, and 1950 DA Sugita et al., 2019;Barnouin et al., 2019;Rozitis et al., 2014;Miller et al., 2002;Scheeres et al., 2020). Eros, which is similar in size to Phobos (volume∼2503 ± 25 km 3 ), has a bulk density of 2.67 ± 0.03 g/cm 3 , corresponding to 10%-30% porosity, and is considered to be internally fractured akin to a rubble aggregate or broken-up monolith (Zuber et al., 2000;Miller et al., 2002). The sub-kilometer-sized asteroids Ryugu and Itokawa are, on account of their low densities and distances from the Sun, considered to be highly porous bodies with up to 50% porosity in their interiors Watanabe et al., 2019). The bulk density of Ryugu is estimated to be 1.19 ± 0.02 g/cm 3 (Watanabe et al., 2019). This is below the range measured for hydrated carbonaceous asteroids of type Ch and Cgh (1.6-2.4 g/cm 3 ), but the range within which Phobos falls (Vernazza et al., 2015).    Several different interior structure models for Phobos have been considered in the literature, including a homogeneous satellite, a heavily fractured body (resulting from the Stickney impact), a locally compressed body, and a rubble pile moon consisting of voids and/or water ice, among others (Murchie et al., 1991;Britt et al., 2002;Rosenblatt, 2011;Rambaux et al., 2012;Le Maistre et al., 2013Scheeres et al., 2020). A loosely-connected body could be a consequence if Phobos was the product of a tidally disruptive or collisional event around Mars as indicated in a recent study (Bagheri et al., 2021) and/or the result of (re-)accretion from debris (e.g., Singer, 2007;Simioni et al., 2015;Craddock, 2011;Canup and Salmon, 2018;Hesselbrock and Minton, 2017;Citron et al., 2015;Rosenblatt et al., 2016). Here, we focus on four categories of models that encompass a large part of the aforementioned range of models and that satisfy the basic constraints of mass and volume. Models include: • a homogeneous model • ice-rock mixtures • models with density gradients • layered models In the following, we briefly describe the four model families. Model parameters and parameter ranges for all models are summarized in Table 3. All constructed models preserve total mass.
Homogeneous model. In this model (Fig. 5A), Phobos is assumed to be a homogeneous body with density = 1.845 g/cm 3 . We model a looselyconnected body as an effective medium using the lowest reported shear moduli ( ) for very loose sands and clays on Earth, which are of the order 10 −3 GPa (Bowles, 1997). The shear modulus for solid rocks is considered to be 10 GPa, which we employ as the maximum possible value. Bulk modulus is calculated from using Eq. (25) where is Poisson's ratio, which we assume to be 0.28 based on terrestrial experience (Ji et al., 2018).

Ice-rock mixtures.
In this type of model (Fig. 6A), Phobos is assumed to be a conglomerate of ice and rocky material. If Phobos originates from the outer asteroid belt, which is beyond the solar system's ''snow line'', ice is expected to make up part of the body. For example, significant amounts of ice have been observed on several asteroids from the main belt, such as 1 Ceres and 433 Eros (Snodgrass et al., 2017). Although the presence of water on Phobos is yet to be fully confirmed, a weak OH − absorption peak has been observed in spectra from Mars Reconnaissance Orbiter (Fraeman et al., 2014).
The density of ice is taken to be 1.0 g/cm 3 which lies within the range indicated in Zuber et al. (2007), while the density of the rocky material ( rocky ) is computed so as to satisfy Phobos's mass. We further assume that the icy material is in the form of spheres with a radius ranging between 200 m and 1000 m that are randomly distributed in the interior and do not intersect. Since the ice-rock fraction is unknown, we considered four sub-classes of the model with differing mass proportions of ice and rock: Ice-rock I: ice fraction = 15%, rocky = 2.0 g/cm 3 .
The densities of rocky material of the first two models correspond to densities of CI-and CM-type carbonaceous chondritic material (Britt et al., 2002), respectively. CI and CM have reflectance spectra that are similar to certain units on Phobos (labeled the ''blue'' units) (Rivkin et al., 2002), respectively. The densities of the third and fourth models are bracket those of mafic and ultramafic rocks. For the purpose of considering a larger range of models, we constructed an additional 10 ice-rock mixtures for each subclass by randomly varying the distribution of ice and rock. These models will serve for improving statistics when discussing geophysical properties.
Models with density gradients. This model category represents a variant of the layered models, in that density increases continuously with depth ( Fig. 7A) or, as an alternative hereof, decreases in density from the surface to the center of the body (Fig. A.2A). The latter model has been suggested as a possibility for the sub-kilometer-sized asteroid Bennu based on an analysis of the difference of the measured zonal gravity coefficients 1 through 4 with those predicted from a constant-density model (Scheeres et al., 2020) that pointed to surface and interior density heterogeneities. The physical mechanism responsible for creating a density inversion with depth could be related to a past period of rapid rotation while Bennu sojourned in the main asteroid belt that was either YORP-induced or resulted from its accretion following the disruption  of a progenitor body (Scheeres et al., 2016(Scheeres et al., , 2020Hergenrother et al., 2019). Since surface and center density are unknown, we considered two sub-classes of density gradient models with density varying linearly between the surface ( s ) and the center ( c ): Positive gradient model: s = 1.5 g/cm 3 and c = 2.882 g/cm 3 . Negative gradient model: s = 2.0 g/cm 3 and c = 1.382 g/cm 3 .
Again, to consider a larger model range, we constructed 2 additional models with higher and lower gradients by changing the surface density (±0.1 g/cm 3 ).
Layered models. While Phobos, unlike the much larger asteroids Vesta and Ceres, has probably not undergone differentiation, we nevertheless consider a layered model of Phobos to enlarge the model range. The layered models represent intermediate cases between homogeneous and positive gradient models in terms of density. To simulate this case, we divided the body into a number of separate layers with constant density and shear modulus (Fig. A.2E). To achieve a semi-continuous change in properties from the surface to the center, we found that six layers provide a satisfactory balance between a simple ''core-mantle'' model and a continuous density increase. We considered two types of density variations: Layer I: from s = 1.6 g/cm 3 to c = 2.4 g/cm 3 . Layer II: from s = 1.6 g/cm 3 to c = 3.0 g/cm 3 .
All layer boundaries in this model type are the downsized DTM surfaces centered in the CoF forming a Matryoshka-like, or ''object-withinsimilar-object'', structure of the interior. The layers are characterized by variable thickness, since the ratio between the equatorial and polar thicknesses of any layer in Fig. A.2E is the same as the ratio between the largest and smallest dimensions of the satellite's figure. As in the case of the ice-rock mixtures, we constructed an additional set of perturbed cases of both Layer I and II models by varying the density of each layer by up to 0.1 g/cm 3 . In all, we considered 100 models for each of the layered family models. Fig. A.1. Comparison of the tidal force field (top) and perturbation in the tidal gravity field due to redistribution of mass (bottom) for the case of a homogeneous model of Phobos with shear modulus = 0.01 GPa. Since the contribution from the distortion in the gravity field is ∼4 orders of magnitude smaller than the unperturbed tidal force field, we neglect the former effect in our simulations.

Results
In the following, we discuss the model classes indicated by boldfaced values in Table 3. Because of the differences in material parameters between models, particularly rigidity, which is less well constrained, we focus on the form of the surface displacement distribution rather than amplitude.
Surface gravity. Computed surface gravity across the body for each of the main Phobos models is shown in Figs. 5-7 panels B (and Fig. A.2B,F) and is seen to be relatively consistent across models, with the exception of the ice-rock mixture which shows ''spotty'' features. Yet, the geographical locations of maximum and minimum gravitational acceleration ( ) are similar across models. The largest occur at the North and South poles, since the latter are closest to the center of the satellite, whereas the smallest is located on the Eastern rim of Stickney crater (furthest from the center). From a geophysical perspective, surface gravity could potentially be measured by a gravimeter, which would provide information on localized density anomalies related to subsurface porosity as discussed in more detail by Le Maistre et al. (2019).
Tidal deformation. The tidally-induced full 3D surface and interior displacements for each of the Phobos models are shown in Figs. 5-7 panels C and D (and Fig. A.2C-D,G-H) in the form of the amplitude of the displacement projected normal to the surface (panels C), which for irregular bodies generally differs from the radial component of displacement, and displacement magnitude (panels D). The importance of the former projection is apparent for the displacement in the region of the Stickney crater, where the relatively large deformation occurs almost in the plane of the surface, resulting in almost lateral rather than radial mass movement. This is illustrated in Fig. A.3.
From the results, the following general features can be extracted: the hemispheres opposite to and facing Mars experience, as expected, the largest deformation. We also notice that the location of maximum deformation shifts from the center of the Mars-facing hemisphere (homogeneous, ice-rock, and negative gradient) toward the rim of the Stickney crater closest to the sub-Mars point on Phobos (indicated by the colored dot) for the positive gradient and layered models. This can be understood as follows. As illustrated in Figs. 5D and 6D (and Fig. A.2D), the deformation amplitude in the interior of the homogeneous, ice-rock, and negative gradient body is predominantly determined by the tidal force field, which is strongest toward the surface (cf. Fig. 2). Consequently, the outer layers in these three models experience the largest deformation. In contrast, in the positive gradient and layered models (Figs. 7D and A.2H), the core is endowed with strength relative to the other three models and resists being deformed. Accordingly, the total displacement immediately below the sub-Mars point suffers less deformation. In comparison, there is more outer-layer material underneath the Stickney crater that undergoes deformation and which turns out to be the largest in the positive gradient and layered models. The deformation pattern along the equator to the East of the sub-Mars point clearly contains the signature of the Stickney crater in all models G). The smallest surface deformations are, in line with expectations, concentrated in an annulus around the center of the moon that lies perpendicular to its long axis.
In the 3D tidal deformation calculations performed so far, we considered a range of models each with fixed and (cf. Table 3). To understand the deformation behavior more generally for a wider range of ( has a smaller impact on the deformation and is ignored in the following), we varied for all nine models continuously within the ranges specified in Table 3 and recomputed the tidal deformation for each of the ''new'' models. The results of this experiment are summarized in Fig. 8, which shows the maximum deformation experienced by the body as a function of effective (hereinafter eff ). The latter replaces since we compare homogeneous models (single ) with models consisting of different 's (as in the ice-rock mixtures and gradient and layered models) and represents an appropriately-averaged for the more complex models. 1 In addition to these models, we also consider the responses of a spherically symmetric model of Phobos (with a radius as defined in Eq. (19)) and the 1D-equivalent of that, the radially symmetric model of Le Maistre et al. (2013), which relies on the analytical solution for the surface deformation where is tidal potential, is surface gravitational acceleration, is Newton's gravitational constant, and , P , and P are density, radius, and mass of a spherically symmetric Phobos, respectively.
From Fig. 8 we can make a number of observations. Firstly, the maximum tidal deformation for all models increases in an exponential fashion with decreasing rigidity. Secondly, the spherically symmetric 3D solution (bold orange line) is equivalent to the response of the 1D radial model (dashed green line). Thirdly, the response of homogeneous Phobos (blue line), but including proper shape, is slightly larger than the two aforementioned cases. This is reasonable since the maximum deformation occurs along the long axis of Phobos, which is pointing in the direction of Mars. The other 3D models (stars and inverted triangles) are located around the response line of homogeneous Phobos and vary slightly with eff because of the averaging scheme applied. In summary, this plot provides a first-order means of determining the effective shear modulus from observations of displacement.
1 The most commonly used schemes for averaging material properties are the Voigt (arithmetic mean), the Reuss (harmonic mean), and the geometric mean (Berryman, 2013). With the exception of the homogeneous model, we rely on the geometric mean geom = ∏ with ∑ = 1, where is the number of nodes in the mesh, rigidity of node , and the relative weight of the node.   Table 3). (E) Layered density model of Phobos. (F-H) as (B-D). For the model shown in (E), shear modulus is 0.3558 GPa and core density is 3.0 g/cm 3 (see Table 3). The green dots indicates the location of the sub-Mars point. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Geophysical predictions: comparison with observations. Following Le Maistre et al. (2019), here we compare the observed geophysical properties (libration, gravity, and MoI) with those computed based on all of the nine Phobos models. Computed and observed libration amplitude and degree-2 gravity coefficients are shown in Fig. 9. The computed libration magnitude (Fig. 9A), which ultimately depends on the principal MoIs (discussed below) and therefore density distribution, is similar for most models (around 1.14 • ), with maybe the exception of some of the ice-rock mixtures (I and IV), and in overall agreement with the observed values of Willner et al. (2010), Oberst et al. (2014), Burmeister et al. (2018) and Lainey et al. (2021). This implies that from libration magnitude only, it is currently not possible to distinguish between the various interior structure models; yet, as shown by Le Maistre et al. (2019), if local density anomalies are present underneath Stickney crater, for example, these can result in changes in libration magnitude that are sufficiently distinct from the observations that such models can be singled out.
The 2-degree gravity coefficients reflect the radial distribution of mass in the interior of the satellite. Predicted and observed 20 and 22 values are shown in Fig. 9B. Most models are located outside of the error bounds of the current observations of Jacobson and Lainey (2014) and Yang et al. (2019) (shown in the inset). While this is also the case for the gradient and layered models, the latter nevertheless result in larger 20 values relative to Jacobson and Lainey (2014), reflecting the increasing density towards the center of the satellite. The distribution of the gravity coefficients for the gradient and layered models shows a linear trend that intersects the homogeneous model, while the other models are distributed around the line. The gradient models are found to be located on either side of the layered models with the positive gradient models fitting the 22 of Jacobson and Lainey (2014), but only slightly out of range with regard to 20 , whereas the negative gradient models are unable to explain the 22 observation of Jacobson and Lainey (2014). An explanation for the discrepancy between predictions and observations could possibly reside in the estimation of the observed gravity coefficients. The latter were obtained by Jacobson and Lainey (2014) based on observations of the satellite orbit. As mentioned earlier Jacobson and Lainey (2014)  We also compute the principal moments of inertia for all models, which are shown in Fig. 10A-C. There are currently no observations with which to compare the computed MoIs. Nevertheless, relative to the principal MoIs of the homogeneous body, we find, not unsurprisingly, little variation for the ice-rock mixtures, since their bulk density distributions are almost homogeneous. For the layered models, however, a difference in principal MoIs of up to 20%-30% is observed; the more the mass is concentrated toward the center, the lower the MoI and the larger the difference relative to the homogeneous model. This indicates that ''random'' ice-rock mixtures are generally indistinguishable from a homogeneous model, but that a layered model is potentially distinguishable. The uncertainties associated with the ice-rock mixtures and layered models (error bars in Fig. 10A-C) reflect the range in MoI for the randomly constructed models of each of these families. The range in ice-rock mixtures induces relatively small variations in MoIs, which is to be expected given that these mixtures ''behave'' as homogeneous bodies. The gradient models, in contrast, are distinct to the homogeneous case, particularly the positive gradient model, which shows a decrease of up to ∼8%, whereas the negative gradient model is characterized by increased MoIs relative to the homogeneous model.
Another means of deriving information on density variations in the interior of the body is through the displacement of the CoM relative to the CoF. The CoM can in principle be determined from orbit through the degree-1 gravity coefficients (Eq. (21)). However, as illustrated in Fig. 10D-F, all presently considered models result, within model uncertainties, in effectively the same CoF-CoM offsets. This results from the fact that all these models consist of concentric layers centered in the CoF with constant density within a given layer of a given model. There is thus little to be deduced from the latter, although more complicated models with local density variations underneath the Stickney crater appear to provide a stronger signal (Le Maistre et al., 2019). However, given the few constraints currently available, there is little point in constructing more complicated models that encompass e.g., lateral density variations.

Discussion and conclusion
Comparison of geophysical predictions with current observations have shown that it is generally difficult to discriminate between homogeneous models and ice-rock mixtures given current uncertainties. Degree-2 gravity coefficients and MoIs, on the other hand, appear to provide a means of differentiating gradient and layered models from the other two families. The main reason for this is that the models considered here do not impart strong enough changes to the principal MoIs that these are measurably different. Le Maistre et al. (2019), for example, considered the case where the structure immediately beneath the Stickney crater was more fractured than material elsewhere in the body. Such models involve changes to the principal MoIs and therefore impart perturbations in the libration amplitude that in principle are detectable. Here, we abstained from considering such models given the current dearth in observations that would allow us to move beyond first-order differences.
To infer interior structure, in particular a density increase with depth (i.e., positive gradient and layered models), would be through measurements of the surface deformation, preferably along the equator as shown in this study. For these two configurations, the largest deformation occurs toward the rim of the Stickney crater ( Fig. 7C and Fig. A.2G). These signals are relatively robust and in combination with the determination of the degree-2 gravity coefficients could potentially provide a detection of a density increase with depth. According to our modeling results, failure to detect this Stickney crater-related signal could indicate either of the following possibilities: (1) a more rigid satellite; (2) a porous or rubble aggregate moon; (3) a density decrease with depth; (4) any combination thereof; or (5) presence of lateral heterogeneities. If rigid enough, detection of surface deformation would be suppressed to the extent that it may be undetectable, even at the sub-Mars point (in case of a rubble aggregate). Yet, as argued earlier, an entirely solid satellite is very difficult to reconcile with the dynamical evidence (cf. Bagheri et al., 2021), leaving a rubble aggregate as the most likely alternative with maximum deformation at the sub-Mars point. The latter is also detectable from a polar orbit (cf. Fig. 6C).
Based on our modeling, we have seen that tidal surface deformation measurements on Phobos is an important source of information to constrain the effective shear modulus of the body. If we consider rigidities similar to those of Martian regolith simulants (0.01-0.1 GPa) (Delage et al., 2017), then the tidal deformation amplitude can achieve tens of centimeters or even up to a 1 m if is an order of magnitude smaller (cf. Fig. 8). In Fig. A.4 several locations (yellow dashed lines), where the maximum tidal surface deformation is expected to occur on the surface of Phobos, are shown as possible areas of special focus to be considered in future missions. The results presented here are, as noted, based on a purely elastic tidal response. Planetary bodies, however, consist of anelastic materials that dissipate energy, as a result of which, the tidal response, i.e., the tidal bulge raised on Phobos by Mars, will not be in line with the acting force but will lag by a phase that is proportional to the amount of tidal dissipation inside Phobos (Efroimsky, 2012). For rubble aggregate or complete rubble pile bodies, the tidal response is likely to be influenced by its rheological behavior. While the rheological behavior of rubble pile-like bodies is not yet well understood, studies have shown that considerable tidal dissipation could occur within such objects (Nimmo and Matsuyama, 2019;Goldreich et al., 2009;Brasser, 2020). Consequently, our tidal deformation model needs to be augmented with an appropriate viscoelastic model as implemented in e.g., Renaud and Henning (2018), Khan et al. (2018), Bagheri et al. (2019), Lau and Faul (2019). Generally, the tidal deformation of such bodies increases as these become less viscous, i.e., as viscosity decreases, to the extent that the resultant bulge height on Phobos could reach values in excess of 1 m for the configurations considered here (Dmitrovskii, 2020). The bulge heights shown here therefore represent lower bounds. The implementation of anelastic effects in our physical modeling scheme is currently under consideration.