New powerful thermal modelling for high-precision gravity missions with application to Pioneer 10/11

The evaluation of about 25 years of Doppler data has shown an anomalous constant deceleration of the deep space probes Pioneer 10 and 11. This observation became known as the Pioneer anomaly (PA) and has been confirmed independently by several groups. Many disturbing effects that could cause a constant deceleration of the craft have been excluded as possible source of the PA. However, a potential asymmetric heat dissipation of the spacecraft surface leading to a resulting acceleration still remains to be analysed in detail. We developed a method to calculate this force with very high precision by means of finite element (FE) modelling and ray tracing algorithms. The elaborated method is divided into two separate parts. The first part consists of the modelling of the spacecraft geometry in FE and the generation of a steady state temperature surface map of the craft. In the second part, this thermal map is used to compute the force with a ray-tracing algorithm, which gives the total momentum generated by the radiation emitted from the spacecraft surface. The modelling steps and the force computation are presented for a simplified geometry of the Pioneer 10/11 spacecraft including radioisotope thermoelectric generators (RTG), equipment/experiment section and the high gain antenna. Analysis results how that the magnitude of the forces to be expected are non-negligible with respect to the PA and that more detailed investigations are necessary. The method worked out here for the first time is not restricted to the modelling of the Pioneer spacecraft but can be used for many future fundamental physics (in particular gravitational physics) and geodesy missions like LISA, LISA Pathfinder or MICROSCOPE for which an exact disturbance modelling is crucial.


Introduction
The Doppler data available for Pioneer 10/11 have shown that both spacecraft are subject to an anomalous acceleration of a Pioneer = 8.74 × 10 −10 m s −2 , which has become known as the Pioneer anomaly (PA) [1]. This anomaly was detected in the early 1980s as a small, constant frequency drift in the Pioneer Doppler signal and published first in 1998 by Anderson et al [3,4]. It has been independently confirmed by other groups [5,6]. Until now no satisfactory explanation for the effect could be found. Several systematic effects have been ruled out as the cause of the anomaly [2,7,8]. However, the influence of the heat distribution of the craft still remains to be analysed in detail. If the heat dissipation pattern is asymmetric, a resulting recoil acts on the spacecraft. According to a first estimate done by the JPL, only 60 W of directed power could be sufficient to explain the anomaly [2]. As the Pioneers had a total power of 2580 W [2] provided by four identical radioisotope thermoelectric generators (RTG), a small 3 asymmetry in the radiation pattern can cause a disturbance acceleration in the magnitude of the PA. In addition, all electrical payloads aboard generate waste heat, which is radiated by the main compartment and the louver system thus adding to an anisotropic radiation pattern. Approximation [9] and estimation methods [10] support the theory of thermal effects as a major contribution to the overall residual acceleration. This demands a precise thermal analysis of the Pioneer craft.

The Pioneer anomaly
The anomalous acceleration of the Pioneer spacecraft was observed the first time after their last flyby which turned the spacecraft into an escape trajectory. This acceleration is nearly constant, that is, it is constant to within 3%. No convincing explanation related to additional masses in the Solar system [11] has been found. Despite of the remarkable order of magnitude equality a Pioneer ∼ cH , where c is the velocity of light and H the Hubble constant, the influence of the cosmic expansion has been excluded as origin of the anomalous acceleration [7,12,13]. A drag force due to dust in the Solar system also has been ruled out [14].
One of the main uncertainties in the modelling of the spacecraft is a possible anisotropic heat radiation. Even this potential cause might possibly be ruled out from the constancy of the observed acceleration: the main part of power that is radiated away from the spacecraft originates from the RTGs which obtain their energy from the decay of plutonium with a halflife time of 87.74 years. Therefore, after less than three years the decrease of the available energy should lead to a corresponding decrease of the anisotropic radiation and should have been noticed by a decrease of the resulting acceleration by more than 3%. After 10 years the power, and thus also any anisotropic thermal radiation, decreases by more than 10%. Since the PA is constant to within approximately 3% for more than 10 years it has been concluded that an anisotropy of the thermal radiation cannot be the cause of the PA. One may, however, speculate that a degradation of the surfaces of the spacecraft (an effect that has been observed for the Wilkinson Microwave Anisotropy Probe (WMAP) [15] to be much larger than expected) may lead to a time-varying emissivity or reflectivity and, thus, to an additional change in the force acting on the spacecraft. This change might be such that the decrease of the radiation force may be compensated by an increase in the emissivity. These loopholes in the conventional explanation of the non-thermal cause of the PA can be accounted for only if one has a complete thermal model of the spacecraft at hand. In this case, one can use the housekeeping data 2 to obtain information of the behaviour of the properties of surfaces and/or electrical components. This is an additional reason that proves that a complete understanding of the PA requires a complete thermal model of the Pioneer spacecraft.

Thermal modelling
Our approach for the calculation of the disturbance force resulting from thermal dissipation is based on the finite element (FE) method and methods of ray tracing. For the calculation a thermal FE simulation of the satellite to be analysed is needed. With this analysis at hand the resulting dissipation force can be computed. This paper describes the principle of the method and shows the force computation in detail for a test case model of the Pioneer spacecraft 4 with simplified geometries. The results show that a part of the PA can possibly be explained with the anisotropic heat radiation. Up to now only the RTG fuel and a general heat load on the equipment/experiment section have been considered as heat sources in the model. Future analysis will have to take into account all heat sources and the complete detailed Pioneer geometry in order to calculate a more exact result. As will be shown, the method is not restricted to the Pioneer spacecraft but can be used for many other fundamental physics or geodesy missions like LISA, LISA Pathfinder or MICROSCOPE. For this a thermal FE model of the spacecraft has to be created for each new mission as input to the force computing algorithm. Together with a further input giving the actual environmental conditions, the exact force acting on a satellite in a certain environment can be calculated with very high precision.
In our model, we calculate the resulting force directly from the discretised satellite surface elements without the need to use a control volume or radiation point sources as in other approaches. By this method the actual geometry of the craft is included in the radiation calculations which considerably increases the obtainable accuracies with respect to other modelling approaches (e.g. thermal modelling of LISA [16], where only a limited part of the total radiation pattern is used in the computations). The first step in the calculation of the disturbance force consists of the generation of a thermal map of the spacecraft surface. This surface temperature distribution can be acquired from a thermal FE analysis. The geometry of the spacecraft has to be modeled in discrete FE where the actual shape is idealized to develop a sound balance between needed modelling detail and available computation resources. This modelling step has to be carried out for each new spacecraft to be analysed. The creation of the FE model will be shown in detail for a test model of the Pioneer 10/11 spacecraft in section 5. After application of the thermal constraints and loads to the thermal model, an FE solver can be used to obtain the steady state temperature node solutions. Using the FE solution and the housekeeping data, characteristic parameters of the model can be calibrated such that the resulting thermal map fits the available temperature data. Furthermore parameter and case studies can be made by changing material or geometrical parameters of the model and acquiring the solution again. Thus, the influence of different parameters on the resulting dissipation force can be evaluated which enables a detailed analysis of, e.g. surface degradation effects or failure scenarios. The temperatures computed for the model surface nodes are the input parameters for the force computation algorithm. Based on the node coordinates, the temperatures and the material parameters, the resulting force is obtained using ray tracing methods to compute the influence of absorption, reflection and shadowing. The algorithm is realized as a set of C-routines with graphic interfaces in MATLAB. Section 6 discusses the steps needed to compute the heat dissipation force for the Pioneer spacecraft. Other approaches for thermal modelling and force computation based on ray tracing are discussed in [17]- [19].

Theoretical background
A radiating surface generates a recoil force which is proportional to the radiated power. The fundamentals needed for the computation of this force will be presented in this section. The energy flux emitted by a grey radiator with the emissivity ε A for a specified wavelength λ is defined as the spectral radiance L λ [20]. where h is Planck's constant, c is the speed of light, k is the Boltzmann constant and T is the surface temperature. Integration over all wavelengths of the spectrum gives the radiance L: If the radiating surface is an ideal radiator, the radiation pattern is hemispheric and the distribution of the intensity I over the hemisphere can be expressed by Lambert's cosine law: where I n = L A is the intensity of radiation in normal direction of the emitting surface. The hemisphere above the radiating surface is determined by the angles φ and β with 0 φ 2π and 0 β π/2. Each specific solid angle element d can be identified by specific values for φ and β as displayed in figure 1. The integration of the radiance L over the complete hemisphere results in the area specific energy flux E dA which is emitted by the area dA, where d = sin β dβ dφ. The energy flux emitted into a specific solid angle element can be expressed as: The total energy flux emitted into the hemispherical surface integrated over the area of the radiating surface leads to the total power output into the hemisphere: 6 which is Boltzmann's law for radiating grey bodies. The total power output P tot can be formulated in terms of the intensity as P tot = I n π.
(7) Now we have to consider that a specific energy flux E dA, of E dA, = L cos β sin β dβ dφ dA (8) generates the recoil force for each solid angle element [21]. Due to the symmetry of the hemisphere only components normal to the emitting plate will contribute to a resulting recoil force. Thus, the solid angle force components evolve to Integration over the whole hemisphere leads to

Ray tracing algorithm
The ray tracing algorithm is used to compute the forces generated by a given surface heat distribution. This heat distribution is obtained from a thermal FE analysis which will be discussed in section 5 in more detail. For a complex geometry such as the Pioneer RTGs the surface elements not only radiate into space but can also radiate into other elements of the model thus causing absorption and reflection effects. In order to compute the amount of radiation that is absorbed or reflected it is necessary to check if other elements are visible or shadowed with respect to the radiating element. In addition the reflected radiation can cause further reflections on other surfaces. This section will explain how these effects are modeled in the ray tracer and how the resulting force is computed taking into account the radiation losses due to absorption and reflection effects. The input files for the ray tracer, which are exported from the thermal FE analysis, include the following data: • node coordinates in Cartesian frame, • element node list (allocation of nodes to elements), • material parameters (emissivity, reflectivity and absorptivity), • nodal surface temperatures (steady state solution). Based on the input data a mathematical model of the craft is generated and processed. The radiation fluxes resulting from emission, absorption and reflection can be converted to equivalent force contributions where the total recoil force is Here n is the total number of surface elements, F e (i) is the force contribution of unblocked emission from element i, F abs (i) is the loss resulting from absorption of radiation at element i and F ref (i) is the force gain at element i from reflection. The following sections will discuss the modelling of each of the three force contributions. Allocation of solid angle elements and force (see [32]).

Emission
The first step is the calculation of the dissipation force generated by emission. From section 2, we know that the resulting force for a radiating plate element is normal to the element surface. The normal vector e n (i) in the global frame is computed from the element node positions and stored in D for each surface element i. With (11) we can compute the resulting force vector for each surface element i with Boltzmann's law for grey body radiation and the effective factor 2/3. The emission coefficient ε A enables the modelling of different optical properties for each radiating surface in the model. Thus, the resulting force vector can be expressed as

Absorption
For the modelling of the influence of absorption and reflection, methods of ray tracing have to be applied. The algorithm we developed computes the equivalent force contributions for each model surface separately. In each computation step an element is considered as actively emitting while all other elements in the model are considered as passive. Here, the passive elements are used to model absorption and reflection resulting from the emission of the active element. After the flux computations have been finished, the next element is set to active and the computation is repeated until all surface elements in the model have been processed. The definition of the solid angle element is used to generate a pattern of rays for each model surface, which specifies the fraction of radiation emitted into each specific direction. The hemisphere into which the radiation is emitted is represented by a two-dimensional array of elements which are characterized by the angles β and φ. Figure 2 (left) shows the resulting allocation of solid angle elements S φ,β in a tesseral division over the hemispherical surface. The vector from the radiating element centre e c (i) to the centre of a solid angle element e c,Omega (φ, β) is defined as the ray direction R(i, φ, β): For the computation of the force all elements that are visible to the active element have to be determined. The visibility of any receiving element is characterized by two requirements. An element j visible to the active element i has to be inside the hemisphere of the sending element. This first requirement can be formulated as where j is the index of the receiving element. The second requirement for visibility is that the receiving element surface has to be orientated towards the surface of the active element. As we consider only two-dimensional elements for the force computation no radiation into the back side of an element is allowed. This second criterion can be formulated as A receiving element is visible to the sending element if both criteria in equations (15) and (16) are met. After checking the visibility criteria for each receiving element, the results are stored into the visibility matrix U for further processing. For the modelling of shadowing the receiving elements that are visible (flagged by a 1 in the visibility matrix U) are sorted by distance from the active element and stored again in the sorted order. While processing the hit detection the rays that are closest to the sending elements are checked for hits first and shut down if a hit has been detected. Thus, the rays will never reach elements which are shadowed by other elements. Although an element may be visible from the radiating element not all rays will go into the direction of the receiving element. Therefore, only rays pointing roughly in the direction of the receiving element will be considered for the intersection calculation. This criterion can be formulated as Radiation can only be absorbed or reflected by the receiving element if the ray hits the receiving element surface. Thus, the intersection point of the receiving element plane and each ray is computed using where N 1 , N 2 and N 3 are the node coordinates of the receiving element. The solution of this system for r , s, t gives the coordinates of the intersection point P with The four node vectors of the receiving element are now used to generate two adjacent triangle surfaces on the element surface. The triangles 4 ) are now considered as vertex vectors for the use of barycentric coordinates. The node vectors and the intersection point vector have to be projected into the receiving element plane. For this, we need the angle between the global coordinate x y-plane (with normal vector e N ,z ) and the element plane as well as the rotation axis E with = arccos The transformation into the receiving element plane is realized using the quaternion rotation matrix A [22] which is with q 1 = e 1 sin 2 , q 2 = e 2 sin 2 , q 3 = e 3 sin 2 , q 4 = cos 2 (22) and The vectors in the new coordinate system are given by The element nodes and the intersection point are in the same plane and share the same z-coordinate in the new coordinate system. Therefore, the dimension of the vertexes can effectively be reduced to two. The coordinates of the point P(x, y) in barycentric frame [23] are defined as where the barycentric coordinates c 1 , c 2 and c 3 are the three unknowns of this system. All the other values are given by the transformed triangle vertex points and the transformed intersection point coordinate. With the solutions c 1 , c 2 and c 3 of the system (25) it can be checked whether the intersection point lies within the triangle or not. The criterion for the intersection point P being inside the triangle is If Pis in either one of the triangles the surface element has been hit and fluxes between the active and receiving element have to be computed. The ray and all following rays that hit the element are deactivated, thus modelling the shadowing effect of the target surface. For the hit element the so-called view factors are computed, which denote the fraction of radiation emitted by the active element 1 which reaches a receiving element 2. Let 1 and 2 be the angles from the connection vector of the two element middle nodes e c (i, j) to each element normal vector, respectively (see figure 3). The view factor from surface A 1 to A 2 is defined as [24] (27) and the radiation flux from element 1 to element 2 can be approximated to with P tot as the total power emitted from element 1 and r = | e c ( j) − e c (i)|. Now the absorption contributions can be calculated for each surface element using (28):

Reflection
The surface element receiving the radiation partially absorbs and partially reflects the incoming ray. The fraction of radiation emitted by element i reflected at element j (with subsequent reflection k) is characterised by the reflection coefficient γ : Considering specular reflection the direction of the reflected flux P i, j,k,ref can be computed for each incoming (absorbed flux) by where e N is the normal direction of the reflecting element and r inc is the incoming flux direction. The reflection force contribution can then be computed for each surface element with where n ref is the total number of considered reflections. With this the recoil force computation method is completely characterized.

Method verification and performance
In order to test the performance and to verify the introduced thermal recoil force computing method a simple test case is processed analytically and compared with the numerical results. For this test case the geometry shown in figure 4 has been selected where two parallel surfaces with the same surface a · b are facing each other with the distance l. Surface 1 emits radiation, which is absorbed by surface 2 thus leading to a deviation F abs from the total recoil force caused by unblocked emission. This resulting force is computed analytically and numerically with varying element sizes A e . The analytical result can be acquired by introducing the geometry where The absorbed force can be computed using (28) and (9) assuming a surface temperature of 300 K, an emissivity coefficient of 0.9 and surface dimensions of a/l = b/l = 1 as a test case configuration. This leads to an analytical value of F abs,analytical = 2.745 × 10 −7 N . With the same configuration and parameters the test case is processed with the ray tracing method using a fixed pattern of 100 × 100 emitted rays per FE model surface. For this the resulting absorbed flux and the respective force loss is computed for different element sizes thus effectively changing the ratio χ = √ (A e )/l. The numerical results for the resulting force losses for different χ are plotted in figure 5. As can be seen in figure 5 (left) the numerical approximation differs from the theoretical value for low ratios of 1/χ and converges to the exact result for values of 1/χ > 4. This implies that precise force computations are possible for FE models with a fine mesh structure. By computing the maximum χ for a given model the quality of the numerical approximation can be quantified. A higher number of elements raises the number of ray tracing equations to be solved and thus the computation time considerably as can be seen in figure 5 (right). Therefore one has to find a compromise between the desired accuracy and the computational power/time available. In practice, one can find an optimal value for χ by processing a number of thermal recoil force computations while reducing the element size. If the difference between two consecutive simulations is sufficiently small no further reduction of χ is needed.
Besides the size of the FE surfaces the accuracy of the solution highly depends on the ray tracing quality. As explained in section 3, each surface element in the model emits a set of rays in order to check possible further surfaces that absorb or reflect emitted radiation. If the total number of emitted rays is too small, some elements which can be seen from the emitting element will not be detected as hit (because no ray intersects) and therefore the force result will be erroneous. For a better understanding of this effect, the test case geometry is processed with a fixed element size and a varying number of emitted rays. The results are plotted in figure 6 (left). For a small number of rays per FE surface the resulting computed force is far below the theoretical value of 2.745 × 10 −7 N. This results from the fact that only one hit can be detected per traced ray. If the total number of rays is too small, elements that are visible to the currently active element can be misjudged as non-hit because no ray intersects with them due to a too coarse angular ray distribution. Thus, the influence of these elements is not included in the total force loss sum. Looking at the characteristics of figure 6 one infers that the computed force loss (which corresponds to the lost flux) approaches the theoretical value while the number of rays per element is increased. For a higher number of rays more surface elements visible to the currently emitting element are detected, thus reducing the overall error. The computation time for different number of rays clearly shows that the processing time grows linearly with the number of rays. This is not surprising as the intersection computation which has to be processed for each ray in the model is the most demanding task within the ray tracing algorithm. In order to validate the quality of the solution and to ensure that enough rays have been included in the simulation, the computation has to be performed with a fixed element size and a varying number of rays. The optimal number of rays has been found when the solution converges.
Using the two surface test case, reflection behaviour can also be analysed. For this the reflectivity of the absorbing surface is set to γ = 0.9 and the orientation angle is varied from 0 • (facing) to 90 • (perpendicular). Figure 7 (left) shows the characteristics of the resulting recoil and the absorption/reflection corrections. As expected, the starting configuration results in a recoil force aligned with the z-direction, which is reduced (from the free emission recoil) by a factor of roughly 2 due to absorption and reflection. As the orientation angle grows, more radiation is reflected to the sides and F abs/ref,x grows until a maximum is reached at an angle of 45 • . For growing angles the radiation components reflected back into the z-direction, as well as the fraction of absorbed energy decrease. Thus, at an angle of 90 • the effective recoil force is aligned with the surface normal and the resulting z-component reaches a maximum value, which resembles free emission. This follows intuitively as the receiving surface is no longer blocking  any of the emitted radiation. In figure 7 (right) the characteristics of the force components (aligned with the surface normal, orientation angle of 0 • ) lost to absorption and reflection are plotted. Without considering reflections the result matches the theoretical absorption value (dashed line). Taking into account, reflection an oscillating behaviour of the result becomes visible. This oscillation results from the change of direction for each subsequent ray reflection and the decreasing available reflection energy (reduced by absorption and emission of reflected radiation into free space).

The Pioneer test case
As a first application for the presented thermal recoil force computing method a test case based on the Pioneer 10/11 geometry is processed. This model enables a check of the order of magnitude of thermal perturbations for the Pioneer spacecraft.
For the presented case the input FE model is generated with a set of connected ANSYS Parametric Design Language (APDL)-macros [26]. This macro consists of five independent major modelling steps as displayed in figure 8: • Creation of the geometry and premesh with primitives and cutting operations, • definition and assignment of the material parameters, • meshing, The test case model includes the three main geometrical elements of the Pioneer spacecraft, namely high gain antenna, equipment/experiment sections as well as the RTGs.
A global startup macro coordinates the interfaces between the submacros and can be used to apply different sets of parameters to the simulation (e.g. different material parameters, different heat loads), which enables parametric studies. In a global solution macro the system of equations defined by the geometrical model and the resulting FE mesh are processed by the ANSYS FE solver. A set of simulation parameters including time-stepping, solver type and solution convergence criteria can be specified.

The RTG model
The geometry of the Pioneer 10/11 RTG is described in detail in [27]. It consists of the fuel assembly with the plutonium fuel in the centre and different shielding layers for radiation containment. The RTG has six radiation fins, which are sloped at the ends. For better radiation properties the whole surface is painted white. In order to generate the RTG FE model the geometry has to be discretized and simplified. For the test case, the outer geometry of the RTGs are modelled exact while the interior composition is modelled as two solid cylinders thus simplifying the composition of different capsule shielding layers. The interior cylindrical volume represents the fuel assembly where the heat load is applied. The surrounding outer cylinder models the heat shield, which smoothes the heat flow between the heat source and the radiating outer surfaces.
The considered material parameters are listed in table 1. The resulting FE mesh is shown in figure 9 on the right side. Radiation processing contains a precise view factor calculation between all elements in the model for the implementation of radiation exchange between the model surfaces and into space. Here, a global radiation matrix including the view factors between every finite surface combination in the model is computed in preprocessing and included as a single superelement containing the radiation constraints for the FE solution. As boundary conditions an outer temperature of 3 K was chosen to model the environmental conditions in deep space. This value is defined as a thermal constraint for the so-called  Figure 9. RTG premesh(left) and mesh(right) (see [32]).
spacenode, which acts as a heat sink for all radiating surfaces. The outer RTG surfaces (fins, body and closures) are defined as heat radiators. The thermal load on the fuel cylinders is modelled as a uniform heat generation on the FE assigned to the fuel cylinder volumes following the approach of Scheffer [9,28], where the total power provided by the four Pioneer RTGs at a time d (given in years) can be computed with where and P start = 2580 W, P el,ref = 68 W, p = 2.6 W.
Here, P el (d) is the electrical power usable for the electrical subsystems and P RTG (d) the total waste power which has to be dissipated. This implies that for different simulated mission times the heat load on the fuel cylinder FE is adjusted according to the decay state of the radio isotopic fuel. In addition to the decay of the fuel, the conversion rate between available energy and provided electrical energy drops over time by ageing of the generator hardware. In order to model this effect, available RTG fin root temperature measurements [29] are used as thermal boundaries for the nodes at the model fin roots.

The equipment/experiment section model
The equipment/experiment section model geometry was taken from available technical drawings given in [30]. The model includes the outer shapes of the main compartments, the launch adapter and the louver system arranged radially around the launch adapter. The louver opening angle can be adjusted thus enabling the modelling of the passive heat control of the main compartment. The louver blades are highly reflective while the second surface mirrors located beneath them possess a high emissivity. All other outer surfaces are modelled with the optical properties of the thermal insulation. The considered material parameters for the different subcomponents are given in table 1. For this first study the interior of the craft is modelled as a solid homogeneous heat conducting medium where the total heat load is distributed homogeneously within the equipment/experiment section centre volumes. Future models will also include the exact geometry of the interior compartment, thus modelling the heat produced by each electrical appliance separately. For the simulation of radiation emitted by the craft three superelements containing the radiation matrices of louver blade surfaces, second surface mirrors and equipment/experiment section outer surfaces, respectively, are defined.

The high gain antenna model
The geometry of the high gain antenna model is based on technical drawings given in [30] and is shown in figure 10. The shape is approximated by a set of coordinates on the antenna cross-section. A spline curve that is based on the cross-section coordinates is rotated around the global z-axis to generate the antenna volume. The antenna surface is painted white on the front and bare aluminium at the rear; all relevant material parameters are given in table 1. For this first assessment of the thermal recoil force magnitude the antenna itself is considered to be thermally neutral and is only included as a reflecting surface for the computation of the thermal recoil force resulting from emissions of the equipment/experiment section and the RTGs. After the definition of the FE models for RTG, equipment/experiment section and antenna the ANSYS FE solver is used to process the system of differential equations by iteration until a steady state is reached.

Computation of the thermal recoil force
The scope of this analysis is a first assessment of the magnitude of recoil forces resulting from anisotropic heat radiation for the Pioneer 10/11 satellite, the dependence of these forces on different input parameters (thermal loads) and the evaluation of these results with respect to the PA. Due to many simplifications of the model geometry, heat loads and boundary conditions the force results will not deliver an exact solution for the thermal recoil forces but can be used to estimate the magnitude of the effects and to derive the characteristics of the dependence of the forces on the applied heat loads. For a first assessment, one can identify four main contributions to anisotropic heat radiation which may lead to a resulting recoil directed against flight direction. Figure 11 shows the radiation sources and the qualitative flux paths, respectively.

RTG radiation reflected at high gain antenna
The RTGs, being the main source of waste heat radiation, are situated at a distance of roughly 3 m from the main compartment centre. Thus, the resulting fluxes absorbed and reflected by the high gain antenna dish surfaces will be small and only a fraction of the total emitted RTG power contributes to a resulting force. Nevertheless, owing to the large available power (2580 W at BOL) a small fraction of that power may already result in accelerations in the magnitude of the PA.

RTG radiation reflected at the equipment/experiment section panels
The high gain antenna and the equipment/experiment section side panels are also within the field of view of the RTGs. Here the distances of the sources and the absorbing/reflecting surfaces are quite small, leading to considerable radiation fluxes. The main part of these radiation fluxes will be reflected in radial direction and do not contribute to a force aligned with the direction of flight. Nevertheless, resulting from the mounting position of the RTGs (not aligned with panel centres) and secondary reflection at the high gain antenna a small fraction of the total power radiated into the panels may still contribute to a resulting recoil.

Louver and equipment/experiment section radiation emitted into space
The louver system dissipates the main part of the electrical waste power into space. The resulting recoil force acts against flight direction and is directly proportional to the louver blade opening angle. For a completely opened louver the power is radiated directly into flight direction, converting the main part of the radiation into a force against flight direction. For lower surface temperatures the louver angle decreases, reflecting more radiation back to the craft, which leads to a decreasing recoil. The radiation emitted by the free surfaces of the equipment/experiment section rear wall adds to the louver emission. Although the temperatures of the outer surfaces are quite low in comparison to the RTG surfaces, the fact that the emitted radiation is not blocked by other surfaces and the big surface size may lead to a considerable contribution to the recoil force.

Equipment/experiment section radiation reflected at high gain antenna
Although the major part of equipment/experiment section is radiated by the louver system, a fraction of the total power is emitted by the surface pointing to the high gain antenna and then reflected several times until the radiation can be emitted into space. This effect may cause a minor contribution to the total recoil force. For each of the identified effects a test scenario is formulated and processed with the introduced force computation method. The thermal equilibrium states are simulated in thermal FE analyses for different boundary conditions (variation of RTG and equipment/experiment section heat load, corresponding fin root temperature constraint). Figure 12 shows the complete input model and the actual Pioneer 10/11 geometry for comparison. As can be seen, the dominant geometrical features of the Pioneer satellite have been included into the FE model. All computed forces are delivered in the global Cartesian frame with the z-axis pointing against the direction of flight.

Analysis results
The ANSYS sparse matrix FE solver is used to solve the resulting system of equations. After the solution has been acquired, a superimposed mesh of rectangular shell elements is generated on

Test case I: RTG radiation reflected at high gain antenna
The contribution to the total recoil force caused by RTG radiation reflected at the high gain antenna has been computed for the mission years 1974 to 2002 with resolution of 2 years. The corresponding boundary conditions are acquired from (35) and fin root temperature measurements from Pioneer 10. The nodal loads (temperatures) of the high gain antenna are constrained in this respect so that only the RTG surface temperatures will be included in the subsequent force computation. Following the argument in the previous section, the equipment/experiment section surfaces are not exported, thus only allowing reflection of emitted radiation at the antenna surface. A resulting exemplary surface temperature distribution derived from an FE analysis for thermal conditions in 1976 is shown in figure 13. The highest temperatures are obtained in the fuel assembly where the heat load has been applied. The main RTG body has a moderate temperature with only minor gradients at the body surface. The fins have the lowest temperatures with a high gradient to the outer fin tip. A difference in temperature for both RTGs is clearly visible. This may result from an insulating influence of the main compartment on the RTG closer to the craft and is introduced by the inclusion of fin root temperature data in the FE simulation. For each processed mission year the surface temperature distribution has been computed in an FE analysis and then processed in the ray tracer to compute the corresponding recoil force. Here an angular resolution of 150 × 150 rays per surface element and a maximum number of 3 reflections per ray have been processed. The resulting forces are converted to effective disturbance accelerations assuming a reference mass of 246.5 kg (Pioneer 10 dry mass of 233 kg and half of the propellant mass of 13.5 kg [30]). Figure 14 (left) shows  RTG fin root temperature data [29] (left) and resulting disturbance acceleration acting against flight direction, test case I (right).
smoothed Pioneer 10 RTG temperature measurements for one of the RTG assembly consisting of two RTGs. As can be seen, the temperatures of each RTGs fin root differ by about 5 K during the mission. The temperature data used as boundary for the FE solution is indicated with an arrow for each specific computation at the lower end of the graph. Figure 14 right shows the computed Pioneer 10 RTG recoil acceleration against the direction of flight for different mission times. Here, the recoil mainly results from RTG radiation absorbed and reflected back from the high gain antenna surface into the flight direction. The computed acceleration decreases almost linearly with time from 7.21 × 10 −10 m s −2 to 5.59 × 10 −10 m s −2 by about 22% within the analysed time frame. A decrease of this magnitude, if present, should be visible in the now prepared 30-year interval of Doppler data as was also stated by [31]. Looking at the obtained values, the computed RTG acceleration is clearly significant with respect to the PA (slightly below 2/3 of PA at mission end). Due to the test case geometry (only RTGs and antenna as emitting /reflecting surfaces) this value might exceed the RTG contribution obtained with a full model because parts of the main compartment may shield antenna surfaces from RTG radiation. Due to the r 2 -dependency of the magnitude of radiation fluxes between RTGs and antenna this shielding influence is not dominant because those antenna fractions which are shielded by the main compartment are also those with the biggest distance from the RTGs, thus contributing less to the resulting acceleration. A detailed evaluation of this effect will be performed in future investigations.

Test case II: RTG radiation reflected at side panels
The computed acceleration results for test case II show that the major part of RTG radiation emitted to the equipment/experiment section side panels is reflected radially into space and has only a minor influence on the acceleration aligned with flight direction. The computed acceleration components in the z-direction for maximum RTG power are well below 10 −12 m s −2 . The small influence of this effect is quite understandable, looking at the optical properties and the geometric configuration of the craft. Looking from the position of one of the side panels, the angles to the visible antenna surfaces are quite small, resulting in small view factors for the antenna-panel surface pairs. Therefore, only a small amount of the RTG radiation reflected at the side panels will reach the antenna surface. Due to the small influence on the total acceleration, even for maximum RTG power, this effect will be neglected in the further analysis.

Test case III: equipment/experiment section rear and louver system radiation emitted into space
The louver emissivity exceeds the effective emissivity of the MLI surface considerably. Thus the main fraction of equipment section power is emitted by the rear. Figure 15 (left) shows the input heat load from [9], which has been distributed evenly among the equipment/experiment section volumes. The resulting recoil accelerations acting against the direction of flight are plotted in figure 15 (right). Again a Pioneer 10 mass of 246.5 assuming a half full fuel tank has been used as reference. As can be seen, the acceleration decreases linearly with decreasing input heat load from a maximum of 3.65 × 10 −10 to 1.75 × 10 −10 m s −2 . With this the computed equipment/experiment section recoil contribution is about a factor of 2 smaller than the RTG contribution, but is still considerably high. The current model lacks the detailed interior composition of the craft as well as important parts of the thermal control subsystem (shunt radiator). In addition, in a more realistic approach a part of the total power has to be radiated by the antenna and counteracts the compartment contribution. The obtained results have therefore to be viewed as order of magnitude estimates, which may decrease as more 22 model detail is added. A more precise result can and will be obtained with a complete model including interior and exterior details as well as a temperature sensor fit comparable with the one performed for test case I.

Test case IV: equipment section front radiation reflected at high gain antenna
For the assessment of the acceleration resulting from the equipment/experiment section front radiating into the high gain antenna, thermal FE simulations ranging from equipment/ experiment section heat loads of 59 to 132 W have been performed. The results show that this test case only contributes little to a resulting recoil. Due to the low emissivity of the MLI, the main fraction of waste heat will be emitted by the louvers, even when all louvers are closed.
Looking at the geometry one infers that radiation emitted by the equipment section front only reaches open space at low angles. From the total emitted radiation a considerable amount is absorbed from either the antenna or equipment section surface (multiple reflections). From the rest of the radiation a fraction is emitted in flight direction, another fraction is emitted against flight direction. Due to the size of the antenna the fraction contributing to a force directed against flight direction is slightly larger but the overall effect is weak and can be neglected.

Summary of thermal recoil force calculation
The results obtained for the introduced test cases show that the acceleration caused by anisotropic heat radiation has to be included in PA investigations. The magnitudes of the computed accelerations, although obtained for simplified model geometries, are in agreement with results and estimates obtained by other groups [9,10] and underline the importance of more detailed analyses of the thermal effects. The characteristics and dependence of the accelerations on the applied heat loads have been shown in a qualitative analysis. It has been emphasized that the models have to be improved with more geometric detail before an exact solution for the thermal recoil acceleration can be provided. Nevertheless the qualitative analysis of the time evolution of the accelerations is a considerable improvement in the understanding of the importance of accurate thermal modelling for the Pioneer 10/11 spacecraft.

Conclusion
With the method presented in this paper a powerful tool has been created, which allows the heat dissipation force calculation for satellite surfaces in unprecedented detail. The method has been fully applied to a test case model which includes the main geometrical shape of the Pioneer 10/11 spacecraft using simplified interiors in this first approach for the evaluation of the magnitude of thermal recoil forces. Due to simplifications in the structural, material and load models the results cannot be interpreted as a final solution for the exact value of the thermal recoil force for Pioneer 10/11 but they can give a first estimate to the importance of accurate thermal modeling with respect to ongoing PA investigations. The acquired results show that the heat dissipation acceleration may explain a significant part of the PA, which makes further investigations in terms of a thorough thermal analysis with detailed modelling of the Pioneer craft necessary. The presented analysis delivers qualitative characteristics of the dependence of the thermal recoil force directed against flight direction for different heat load states. With this it has been possible to resolve the change of the resulting recoil force during the mission. All computed values presented in this study are of course only valid for the processed test case geometries, sets of selected material parameters and corresponding assumed loads. In order to evaluate more accurate values for the acceleration originating from the asymmetric heat dissipation, the recoil force has to be computed for the whole detailed geometry of the craft, considering all heat sources onboard (Radio isotopic Heating Units (RHU), louvers, shunt radiator, payloads) in the analysis. This also implies that a more detailed FE model, including the interior composition of equipment/experiment section, the exact antenna shape and external/internal payloads for the Pioneer 10/11 spacecraft, will have to be developed. Although it is a strong guess that the interior composition of the RTG has only a weak effect on the resulting thermal force, this has not yet been proven in an analysis. Therefore a more detailed RTG model, taking into account the interior composition with different insulation sheets, capsule supports and radiative heat transfer within the RTG assembly, will be needed. The models will be generated in coordination and cooperation with the members of the Pioneer collaboration in order to define baseline geometries to be able to compare the methods and assumptions of the different groups working on thermal modelling. Furthermore the new detailed force computation will benefit from the many contributions regarding thermal modelling and thermal force computation which have been done so far by our colleagues in the Pioneer collaboration. The presented ray tracing algorithm delivers accurate results for the thermal dissipation force and shows good computation performance for complex spacecraft geometries. With a thermal FE model of the spacecraft and the proper importing methods, the ray tracer can be used for many other satellite missions, which will be a necessity, especially for the fundamental physics missions like LISA, LISA pathfinder and MICROSCOPE where the exact modelling of the disturbances is crucial for the mission goal. With minor changes it will also be possible to use the ray tracer for the modelling of other surface forces such as solar wind or solar radiation.