text:
\section{Introduction}

Planets and moons reorient with respect to the stars when internal processes or external impacts change their distribution of mass. The phenomenon was first analyzed for Earth, and is thus referred to as the ``true'' polar wander (TPW), to be distinguished from the ``apparent'' motion of the rotation pole that is perceived by an observer on a drifting continent \citep[e.g.][]{Besse2002}. Measuring the Earth's rotation pole and studying its dynamics has a long history \citep{Munk1960}, but for planets and moons with a non-negligible tidal bulge the approach to TPW is usually quite crude: only the equilibrium orientation is assessed. That is, the equilibrium inertia tensor of the body is evaluated and diagonalized, with the principal directions marking the final positions of the rotation and tidal axes \citep[for a review, see][]{Matsuyama2014}.

While the dynamics of reorientation of a tidally locked body are governed by the viscoelastic readjustment of the rotational and tidal bulges, in the equilibrium state the symmetry axes of the bulges are by definition aligned with respectively the rotational and tidal axes. Therefore, the only constituents of the inertia tensor that need to be determined in order to get the principal directions in the equilibrium state are the inertia tensor of the load, and that of the so-called remnant figure. The remnant (or fossil) figure forms when the hydrostatic shape of a planet, that is fluid across its entire depth, ``freezes'' into the growing lithosphere as the planet cools. If the centrifugal or tidal forces later change in response to TPW, the presence of an elastic lithosphere prevents the body from reaching the hydrostatic equilibrium again. Therefore, the fossil figure always reflects the primordial rotation and tides stays always aligned with the primordial directions of rotation and tides, and counterbalances the load that drives TPW \citep{Willemann1984}.

The equilibrium approach neglects the dynamics of reorientation -- it is concerned only with the final orientation of the body. Nevertheless, it can be used to estimate the TPW path when the investigated load is treated as a sequence of loads with a gradually increasing amplitude. The equilibrium orientation is then computed for each of these partial loads, yielding the evolution of reorientation, known as the fluid-limit solution \citep[e.g.][]{Keane2016}. The assumption behind this approach is that the load formation is slow when compared to the rotational and tidal bulge readjustment.

Several planets and moons are thought to have reoriented in the past, with their most striking surface features located near the poles or the equator, depending on whether the associated gravity anomaly has a negative or a positive sign \citep[e.g.][]{Keane2014,Bouley2016,Nimmo2016,Tajeddine2017}. However, the dynamic feasibility of these hypothesised, often large-angle reorientations is only rarely assessed or put in the context of the thermal and orbital history of the body. Moreover, the loading itself may depend on insolation and thus on the orientation of the body in space \citep{Ojakangas1989}. In such a case, it becomes crucial to compute the viscoelastic response to the loading as well as the TPW dynamics in a self-consistent way using a single numerical model.

Reorientation scenarios are often supported by an analysis of the surface stress patterns \citep[e.g.][]{Tajeddine2017,Keane2016}. When the direction of the centrifugal or tidal forces changes, the lithosphere is subjected to stress that is manifested by changes in the tectonic pattern. A large reorientation can thus generate a global network of extensional and compressional fractures, depending on the position of each point with respect to the old and the new centrifugal and tidal potentials \citep{Melosh1980}. In order to predict the map of tectonic stresses, the entire TPW is typically treated as a single, instantaneous event \citep[e.g.][]{Nimmo2016}, or as a progression of fracture-forming events \citep[e.g.][]{Keane2016}. Using the fluid-limit framework is somewhat paradoxical in this regard: the TPW is assumed to be slow so that the time needed for the bulge readjustment is negligible, but at the same time, reorientation is considered to be instantaneous for the purpose of surface stress evaluation (or step-wise in the case of the progression of events).

In a series of papers, \citet{Hu2017a,Hu2017b,Hu2019} developed the first method for computing the dynamics of reorientation of a tidally locked body, and tested it against the fluid-limit approach. The method was also compared with the so-called quasi-fluid approximation, which is traditionally used to evaluate the TPW on Earth \citep{Lefftz1991,Ricard1993} and in which the rotational bulge readjustment is simplified by neglecting all short-term relaxation modes. The TPW speed computed by \citet{Hu2019} lies in between the fluid-limit and the quasi-fluid solutions, and all the three methods converge when TPW is much slower than the bulge readjustment.

The algorithms developed by \citet{Hu2017a,Hu2017b,Hu2019} are relatively complex. At each time step, the problem is transformed to the ``bulge-fixed'' frame, in which the linearized Liouville equations (LLE) are solved. For a tidally locked body, the LLE are solved separately for the tidal and for the rotation vectors, and the condition of perpendicularity of the vectors is achieved by an iterative adjustment of the obtained LLE solutions \citep[see Section 2.2 and Appendix A in][]{Hu2019}.

In cases without tidal forcing, \citet{Patocka2021} obtained the same results as \citet{Hu2017a,Hu2017b} with a simple method, denoted as $\omega||\mathrm{MIA}$. The rotation axis coincides with the major (or main) inertia axis throughout the entire TPW simulation, which has been overlooked by \citet{Hu2017a,Hu2017b}, perhaps because the $\omega||\mathrm{MIA}$ assumption was originally linked to the quasi-fluid approximation. However, the two simplifications in the governing equations, namely the $\omega||\mathrm{MIA}$ approximation of the conservation of angular momentum and the quasi-fluid approximation of the viscoelastic response to loading, should be treated separately \citep{Patocka2021}. While the short-term relaxation modes must be handled with care \citep{Hu2017a,Hu2017b,Hu2019}, the $\omega||\mathrm{MIA}$ assumption is generally valid \citep{Patocka2021}.

Here, the $\omega||\mathrm{MIA}$ method is extended to encompass tidally locked bodies. We show that the solutions from \citet{Hu2019} can be reproduced when it is assumed that the tidal and rotation axes coincide with the minor (mIA) and major (MIA) axes of the total inertia tensor, respectively. This allows for a simple and robust method for computing the reorientation of synchronously rotating planets and moons, hereafter denoted as $o\omega||\mathrm{mMIA}$.

One general question to ask is whether a tidally locked body is likely to reach its equilibrium orientation along the shortest possible path, and how this path depends on the size of the host body. \citet{Hu2019} argue that the reorientation path tends to straighten as the relative size of the tidal bulge reduces, being nearly straight when the tidal/rotational bulge ratio drops to about $0.1$. The $o\omega||\mathrm{mMIA}$ assumption, on the other hand, indicates that even a tiny tidal bulge could substantially distort the TPW paths of tidally locked bodies. This is because it is easy for the load to move the minor and the intermediate inertia axes along the equatorial plane when the tidal bulge is small, making the reorientation path highly curved in effect.

This problem becomes particularly interesting when the reorientation of Pluto due to the formation of the Sputnik Planitia basin is addressed. Owing to the relatively small mass of Charon, the tidal bulge of Pluto is much smaller than its rotational bulge. Centered at 176°E 24°N, the Sputnik Planitia basin is thought to be a positive gravity anomaly despite its negative topography \citep{Keane2016}, suggesting a subsurface ocean \citep{Nimmo2016}.

Preceded by a description and validation of our method in sections \ref{sec:method} and \ref{sec:triton}, the reorientation of Pluto is investigated in sections \ref{sec:pluto} and \ref{sec:discussion}. In section \ref{sec:energy}, the $o\omega||\mathrm{mMIA}$ method is analyzed from the point of view of the energy conservation law. The main findings of the study are summarized in section 7.

\section{The $o\omega||\mathrm{mMIA}$ Approximation}\label{sec:method}

The deformation of a hydrostatically prestressed incompressible viscoelastic ice shell is calculated by integrating the conservation equations for mass and momentum and the constitutive law for a Maxwell body \citep{Tobie2008,Patocka2018,Patocka2021}. The outer boundary is treated as a quasi-free surface \citep[Eq.~7 in][]{Patocka2021}, while the bottom boundary is assumed to be in contact with an inviscid fluid that is in hydrostatic equilibrium and rotates synchronously with the shell. In section \ref{sec:pluto}, where Pluto is investigated, the interior water ocean surrounds a silicate core. In this case, we assume that the ice shell and the ocean form a coupled system that reorients simultaneously, while the core remains fixed with respect to the rotational and tidal axes and its contribution to the inertia tensor can thus be omitted \citep[cf.][]{Ojakangas1989}.

The numerical method that we use to compute the deformation is detailed and validated in \citet{Patocka2018} and \citet{Patocka2021}. The only difference here is that the centrifugal potential $\Psi$ is replaced by $\Psi + \Theta$, where $\Theta$ represents the tidal potential. While the centrifugal potential can be written as
\begin{linenomath*}
\begin{equation}\label{eqRotpot2}
    \Psi = \frac{1}{2}\left( (\bm{\omega}\cdot\bm{r})^2 - \omega^2 r^2 \right),
\end{equation}
\end{linenomath*}
where $\bm{\omega}$ is the angular velocity and $\bm{r}$ is the position vector, the tidal potential takes the following form:
\begin{linenomath*}
\begin{equation}\label{eqTidpot}
    \Theta = \frac{1}{6}o^2 r^2 - \frac{1}{2}(\bm{o}\cdot\bm{r})^2 .
\end{equation}\end{linenomath*}
For convenience, we have introduced a ``tidal vector'' $\bm{o}$ that points toward the host body and has the magnitude:
\begin{linenomath*}\begin{equation}\label{eqTidamp}
o = \frac{3\,GM_\mathrm{h}}{a^3},
\end{equation}\end{linenomath*}
where $G$ is the universal gravitational constant and $a$ is the distance from the host body of mass $M_\mathrm{h}$.

For a tidally locked planet or moon of mass $M$, the angular frequency $\omega$ is equal to the orbital frequency,
\begin{linenomath*}\begin{equation}\label{eqTidlock}
\omega = \sqrt{\frac{G(M_\mathrm{h}+M)}{a^3}} = o \sqrt{\frac{M_\mathrm{h}+M}{3M_\mathrm{h}}}.
\end{equation}\end{linenomath*}
Comparison of Eqs \eqref{eqRotpot2} and \eqref{eqTidpot} shows that the centrifugal and tidal potentials have the same degree two structure, and differ only in the sign, amplitude, and the degree zero component. The torque of the centrifugal force can be expressed as,
\begin{linenomath*}\begin{multline}
\int_{v(t)}{\bm{r}\times(-\nabla \Psi)\,\rho\,\mathrm{d}v} = \int_{v(t)}{\bm{r}\times(\bm{r}\,\omega^2 - \bm{\omega}\,(\bm{\omega}\cdot\bm{r}))\,\rho\,\mathrm{d}v} \\
= \bm{\omega}\times\int_{v(t)}{\bm{r}\,(\bm{\omega}\cdot\bm{r}) \,\rho\,\mathrm{d}v} = -\bm{\omega}\times (\bm{I}\cdot\bm{\omega}). \label{eqRotmoment}
\end{multline}\end{linenomath*}
Here $v(t)$ is the volume of the body at time $t$, $\rho$ is the density, and $\bm{I}$ is the inertia tensor,
\begin{linenomath*}\begin{equation}
    \bm{I} = \int_{v(t)}{((\bm{r}\cdot\bm{r})\, \mathbb{1} - \bm{r}\otimes\bm{r})\,\rho\,\mathrm{d}v}, \label{eqIdef}
\end{equation}\end{linenomath*}
where $\mathbb{1}$ is the identity tensor. The eigenvalues of $\bm{I}(t)$ are referred to as $C, B$, and $A$, in order of their decreasing value. In Eq.~\eqref{eqRotmoment} we employed the fact that a part of $\bm{I}\cdot\bm{\omega}$ is parallel to $\bm{\omega}$ and thus does not contribute to the cross product. Using the same procedure as above, we can express the torque of the tidal force, denoted as $\bm{M}$,
\begin{linenomath*}\begin{equation}\label{eqTidmoment}
   \bm{M} = \int_{v(t)}{\bm{r}\times(-\nabla \Theta)\,\rho\,\,\mathrm{d}v} = \bm{o}\times (\bm{I}\cdot\bm{o}).
\end{equation}\end{linenomath*}
Therefore, the Liouville equation (LE) can be expressed as follows:
\begin{linenomath*}\begin{equation}\label{eqFullLE}
    \dt{} (\bm{I}\cdot\bm{\omega}) + \bm{\omega}\times(\bm{I}\cdot\bm{\omega}) = \bm{o}\times (\bm{I}\cdot\bm{o}) .
\end{equation}\end{linenomath*}
Note that we formulate the LE in the Tisserand (or body-fixed) frame \citep{Munk1960}.

In previous studies, the time derivative in the LE was linked to the short-term viscoelastic relaxation modes of planetary mantles. The problem of TPW was solved using the so-called quasi-fluid approximation, in which the time derivative in the LE and the relaxation modes that are fast compared to TPW were omitted \citep[i.e., the quasi-fluid approximation, see][]{Lefftz1991,Ricard1993}. In the absence of tidal forcing, the LE was thus used in the form $\bm{\omega}\times(\bm{I}\cdot\bm{\omega}) = \bm{0}$. The validity of this approximation was questioned by \citep{Hu2017a,Hu2017b}, who showed that omitting the short-term relaxation modes is a potential source of error in the TPW solutions. Based on this finding, the authors assumed that the time derivative in the LE cannot be neglected, and developed a sophisticated algorithm that solves the LLE in a computational frame whose $z-$axis advances along with $\bm{\omega}$ \citep[see section 3.2 in][]{Hu2017a}. However, as shown by \citet[][Appendix A]{Patocka2021}, there is no strict connection between the time derivative in the LE and the short-term viscoelastic relaxation modes. When the shell relaxation is completely resolved, i.e.~when the employed maxwellian rheology includes all the relaxation modes, then the accuracy of the TPW solution is not harmed by employing the simple formula $\bm{\omega}\times(\bm{I}\cdot\bm{\omega}) = 0$ \citep{Patocka2021}. In other words, the source of error lies only in computing the time evolving shape of the shell inaccurately, not in dropping the time derivative term from the LE.

Motivated by the above, we write Eq.~\eqref{eqFullLE} as
\begin{linenomath*}\begin{equation}\label{eqtomap}
    \bm{\omega}\times(\bm{I}\cdot\bm{\omega}) = \bm{o}\times (\bm{I}\cdot\bm{o}).
\end{equation}\end{linenomath*}
A straightforward solution of Eq.~\eqref{eqtomap} is to place the rotation vector $\bm{\omega}$ along the major axis of inertia, MIA, and the tidal vector $\bm{o}$ along the minor axis of inertia, mIA. The orientation of the body as a function of time can thus be obtained by diagonalizing the inertia tensor $\bm{I}(t)$ at each time step. This special solution does not help in determining the magnitudes of $\bm{\omega}$ and $\bm{o}$, but note that it balances the LHS and RHS of Eq.~\eqref{eqtomap} by setting both sides to zero -- it therefore satisfies also the full LE, Eq.~\eqref{eqFullLE}, provided that $C\omega = \mathrm{const}$ (with $C$ being the time-evolving major moment of inertia and $\omega$ the spin rate). While Eq.~\eqref{eqtomap} governs the directions of the vectors $\bm{\omega}$ and $\bm{o}$ (the orientation of the body), the $C\omega = \mathrm{const}$ condition can be used to describe the changes in the spin rate. In summary, we define the $o\omega||\mathrm{mMIA}$ method as the following set of equations:
\begin{linenomath*}\begin{align}
    \bm{\omega}\times(\bm{I}\cdot\bm{\omega}) = 0, \label{eqtomap_1}\\
    \bm{o}\times (\bm{I}\cdot\bm{o}) = 0, \label{eqtomap_2}\\
    \dt{} (\bm{I}\cdot\bm{\omega}) = 0. \label{eqtomap_3}
\end{align}\end{linenomath*}

In sections \ref{sec:triton}, \ref{sec:pluto}, and \ref{sec:discussion} we perform numerical simulations of planetary reorientation. In these sections, only the directions of $\bm{\omega}$ and $\bm{o}$ are of interest. The therein reported results could be reproduced to a high degree of precision even if the changes in the spin rate were disregarded, that is, if one simply assumed $\omega = \mathrm{const}$ instead of solving Eq.~\eqref{eqtomap_3}. However, as explained in detail in section \ref{sec:energy}, Eq.~\eqref{eqtomap_3} is crucial in the energy balance of TPW, and must be included when an energetically consistent formulation is desired.

Note that Eqs \eqref{eqtomap_1} and \eqref{eqtomap_2} resemble the fluid-limit approach. The difference is that we treat $\bm{I}$ as a time dependent quantity that accounts for the existence and continuous readjustment of the rotational and tidal bulges, while in the fluid-limit approach $\bm{I}$ is just the final ($t{\rightarrow}\infty$) contribution of the load combined with that of the fossil figure. Our $o\omega||\mathrm{mMIA}$ method exploits the simplicity of the fluid-limit formulation, but it is, as shown below, equivalent to the fully dynamical method of \citet{Hu2017a}. Thus, we can capture complex, time-dependent processes, without making assumptions about the relative speed of the bulge readjustment and the TPW rate.

One way to understand the $o\omega||\mathrm{mMIA}$ method is that it filters the free oscillations of the body. This is because wobbling is the situation in which $\bm{\omega}$ and MIA have different directions, making the $\bm{\omega}\times(\bm{I}\cdot\bm{\omega})$ term non-zero, and thus triggering a non-zero rate of the angular momentum $\bm{I}\cdot\bm{\omega}$ as measured in the body-fixed frame, the increment of $\bm{I}\cdot\bm{\omega}$ being perpendicular to $\bm{\omega}$ and periodically revolving around $\bm{I}\cdot\bm{\omega}$ \citep[see e.g.~Figs~2 and 5 in][]{Spada1996b}. Avoiding the wobble, whose characteristic time scale is typically short when compared to geological time scales, is a big advantage in terms of computational efficiency, but has an obvious downside: our method cannot be applied to slowly-rotating bodies such as Venus, on which mega-wobble is thought to be the dominant reorientation mechanism \citep{Spada1996b}. On bodies other than Venus, free oscillations are typically assumed to be only a small perturbation that can be linearly combined with the secular motion of the pole \citep[e.g.][]{Martinec2014}, eventually decaying to zero \citep{Nakada2012} unless having a continuous excitation source \citep{Gross2000}.

\section{Validation: Test Case}\label{sec:triton}

First, we compare the $o\omega||\mathrm{mMIA}$ solution over the algorithm of \citet{Hu2019}. We perform a series of simulations for a model Triton, described in Table 1 and Fig.~1 of \citet{Hu2019}. The tidally locked moon is loaded with a point mass of respectively 1.5, 3.0, and 6.0 $\cdot 10^{17}$ kg, placed at 15° colatitude and -15° longitude. The ice shell includes a 10 km thick elastic lithosphere that formed after the body had reached the hydrostatic shape. The time evolution of colatitude and longitude of $\bm{\omega}(t)$ and $\bm{o}(t)$ computed with the $o\omega||\mathrm{mMIA}$ method (Fig.~\ref{fig:Triton}a,b, dotted lines) matches perfectly the original solution published by \citet{Hu2019} (solid lines). Despite its simplicity, Eq.~\eqref{eqtomap} is thus equivalent to the algorithms developed by \citet{Hu2019}.

In Fig.~\ref{fig:Triton}c, the solutions shown in panels a and b are plotted on the surface of the globe. This graphic representation allows the evolution of $\bm{\omega}(t)$ and $\bm{o}(t)$ to be shown with a single line each. On the other hand, it does not provide a clear view of how the reorientation progresses in time.

The solutions in panels a-c are plotted in the body-fixed (Tisserand) frame, i.e., in the frame in which the governing equations are formulated and solved. A more convenient way to visualize reorientations of tidally locked bodies was proposed by \citet{Hu2017a}, who used the so-called bulge-fixed frame. In this frame, the rotation axis always intersects the surface of the body at colatitude 0°N, while the sub-host point (tidal axis) is always located at longitude 0° and colatitude 90° (Fig.~\ref{fig:Triton}d). In this representation, the lines show how the colatitude and longitude of the load change in time (i.e., the lines show the time-evolving geographic position of the load).

\begin{figure*}[t]
    \centering
    \includegraphics[width=1.0\textwidth]{figs/Fig1ab.png}
    \includegraphics[width=1.0\textwidth]{figs/Fig1cd.png}
    \caption{The test case of Triton. In panels a) and b), the colatitude and longitude (dotted lines) of respectively $\bm{\omega}$ and $\bm{o}$ as obtained in the body-fixed frame are plotted as a function of time and compared over the solutions from the top panel of Fig.~1 in \citet{Hu2019} (solid lines). In panel c), the simulation results are depicted in polar coordinates. In panel d), the same results are plotted in the bulge-fixed frame, in which the rotational axis intersects the surface of the body at colatitude 0°N at all times (black star), while the sub-host point is at all times located at longitude 0° and colatitude 90° (black cross). The red square indicates the initial position of the loads, while the coloured circles mark the equilibrium positions. Temporal evolution is illustrated by the small black dots that evenly sample the trajectory of each load in time intervals of 250 ky.}
    \label{fig:Triton}
\end{figure*}

In the three simulations presented in Fig.~\ref{fig:Triton}, Triton is assumed to orbit a planet much heavier than itself ($M_\mathrm{h}/M\rightarrow \infty$ in Eq.~\eqref{eqTidlock}). Following the exercise in Fig.~2 of \citet{Hu2019}, we perform additional simulations with a point load of mass $6\cdot 10^{17}$ kg, this time varying the $M_\mathrm{h} / M$ ratio, to investigate a system in which the mass of the host body is comparable with or even smaller than the mass of the studied body ($M_\mathrm{h} / M = 1, 0.1,$ and $0.01$, see Fig.~\ref{fig:MhMinfluence}).

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{figs/Fig2.png}
    \caption{The influence of $M_\mathrm{h}/M$ on the path of a positive load, depicted in the bulge-fixed frame (cf.~Fig.~\ref{fig:Triton}d). Temporal evolution is illustrated by the small black dots that are evenly sampled in time at intervals of 250 ky. The load path straightens when $M_\mathrm{h}/M$ or the load size decrease.}
    \label{fig:MhMinfluence}
\end{figure}

If $M_\mathrm{h}/M\rightarrow \infty$, the position of the rotation axis changes by almost 90° in less than 2 Myr for the heaviest load, while the sub-host point takes longer to readjust (green lines in Fig.~\ref{fig:Triton}a and b). As a result, the load path is highly curved. This is not obvious in the body-fixed frame (Fig.~\ref{fig:Triton}c), but becomes apparent when the solution is plotted in the bulge-fixed frame (green line in Fig.~\ref{fig:Triton}d). As the $M_\mathrm{h} / M$ ratio is decreased, the time evolution of $\bm{\omega}$ and $\bm{o}$ becomes more balanced, and the equilibrium is reached in a more straightforward manner. Similarly, the load path straightens when $M_\mathrm{h} / M$ is constant and the load amplitude is decreased, but this is mainly because the final (equilibrium) position of the load changes significantly (dotted lines in Fig.~\ref{fig:MhMinfluence}).

In the next section, we apply our method to Pluto to simulate the ancient reorientation due to the Sputnik Planitia basin. The ``host'' body in this case is Charon, corresponding to $M_\mathrm{h} / M \approx 0.1$. For this mass ratio, the motion of a positive load occurs along a nearly straight line in the bulge-fixed frame (Fig.~\ref{fig:MhMinfluence}, brown curve). Since the sign of the load on Pluto is still a subject of debate, we also investigate the case where the load is negative and then compare the results of our simulations with observations.

\section{Reorientation of Pluto}\label{sec:pluto}

The topography of Pluto is dominated by a deep basin known as Sputnik Planitia. Its elliptic shape and multiring structure indicate an impact origin \citep{McKinnon2017}, but the basin could also have formed by the weight of accumulated nitrogen ice \citep{Hamilton2016}.

The geographic location of Sputnik Planitia (24°N, 176°E) suggests that it is a positive gravity anomaly, balanced ca.~25° away from the anti-Charon point by the remnant figure of Pluto \citep{Keane2016}. In order to explain how such a deep basin could increase the local gravity, \citet{Nimmo2016} speculate that there must be a dense interior ocean below the ice shell. One problem with this scenario is that impact simulations indicate that the post-impact ocean uplift is not sufficiently large to compensate the negative surface topography, and it would soon disappear due to the low viscosity of ice near the ice-ocean boundary \citep{Johnson2016}. Both \citet{Keane2016} and \citet{Nimmo2016} focus on the proximity of the basin center to the tidal axis, but they take little account of the fact that the center is very close to the 180° meridian, i.e., to the plane which contains both the tidal and rotational axes.

In the study of \citet{Hamilton2016}, on the other hand, Sputnik Planitia is shown to have reached 180° longitude already during the tidal despinning of Pluto caused by its large satellite Charon, and the latitude of the nitrogen deposit that presumably formed the basin (24°N) is shown to agree well with the band of low insolation of the highly tilted Pluto. However, the reorientation of the dwarf planet is addressed in a qualitative rather than a quantitative way, and it is not clearly explained how the basin could stay away from the equator without being balanced by a fossil rotational bulge. Since the authors argue that Sputnik Planitia was formed early in Pluto's history, it is unlikely that such a fossil bulge could have existed.

Finally, it was proposed by \citet{Kihoulou2022} that, regardless of its origin, Sputnik Planitia could have reoriented Pluto's primordial, thin ice shell such that its center was at the anti-Charon point. As Pluto cooled and the crustal thickness increased, the sign of the basin gravity anomaly changed, resulting in a further reorientation stage, during which the load began to drift northward. Assuming a 50 km thick elastic lithosphere (i.e., comparable to that assumed by \citet{Keane2016} and \citet{Nimmo2016}), \citet{Kihoulou2022} showed that this motion stops when the load reaches a latitude of about 25°, corresponding to the present position of Sputnik Planitia. It should be mentioned, however, that dynamics of the second reorientation stage were simplified in that the initial position of the load was fixed at the longitude of 180° and the stability of the solution with respect to small perturbations of the initial position was not tested.

The TPW is controlled by the viscoelastic relaxation of the ice shell occurring in response to changes in the rotational and tidal potentials. We assume that the shell is 150 km thick \citep{Johnson2016} and its viscosity varies with radius as follows \citep{Goldsby2001}:
\begin{linenomath*}\begin{equation}\label{eqicevisc}
    \eta_\mathrm{ice} = \frac{T(r)\, d_\mathrm{grain}^2}{ 3 A_\mathrm{act} } \exp\left(\frac{E_\mathrm{act}}{R_\mathrm{gas} T(r)}\right) ,
\end{equation}\end{linenomath*}
where $T(r)$ is the conductive temperature profile, $R$ is the universal gas constant, $d$ is the grain size and $A$ and $E$ are diffusion creep parameters (see Table 1). Temperature of the outer surface is 47 K and the ice/water interface is at 265 K (accounting for the pressure dependence of the melting temperature). A viscosity cutoff of $10^{24}$ Pa s is imposed for numerical reasons. Due to the low surface temperature of Pluto, most of the shell is governed by the cut-off viscosity. While the characteristic time scale of TPW strongly depends on the viscosity of ice and the ice shell thickness, it can be demonstrated that the load path is much less affected by the choice of these parameters. Since the relaxation time is not of primary importance for the discussion of Pluto's reorientation, we prescribe $\eta_\mathrm{max}$ rather arbitrarily, assess the plausibility of the above hypotheses by analyzing the load paths in the bulge-fixed frame, and discuss the robustness of our results only later in section \ref{sec:discussion}.

We investigate the TPW induced by a disc load, whose thickness is increasing linearly from zero to the listed value over the time $t_\mathrm{load}$ (Table 1). As in the case of the point-load studied in section \ref{sec:triton}, the spherical disc is treated as a fixed contribution to the inertia tensor, i.e., no compensation of the load over time is considered. Note than when an impact basin is formed, it is a priori unclear whether the sign of the gravity anomaly it creates is positive or negative \citep[see, e.g., the Extended Data Figure 1 in][]{Keane2016}. In the absence of data, choosing a representative load amplitude is problematic. We investigate spherical discs of thicknesses ranging from -300 to 300 m, allowing us to capture the possible scenarios.

We begin the study of Pluto's reorientation with a demonstration of how the $M_\mathrm{h}/M$ ratio affects the paths of positive and negative loads, initially placed near the north pole and the sub-host point, respectively (Fig.~\ref{fig:Pluto_thereback}). Such an exercise is instructive, since most of the previous studies of icy bodies have focused on the case where the host body is much heavier than the orbiting planet or moon ($M_\mathrm{h}/M\rightarrow\infty$).

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.45\textwidth]{figs/Fig3a_ref.png}
    \includegraphics[width=0.45\textwidth]{figs/Fig3b_ref.png}
    \caption{Load paths on a satellite (Table \ref{tabPluto}), viewed in the bulge-fixed frame. The size of the host body, $M_\mathrm{h}$, is varied, and for $M_\mathrm{h}=0.1212$ also the load amplitude is varied. Dotted lines represent models without a fossil bulge. Solid lines show simulations in which a 50 km thick elastic lithosphere has formed only after the hydrostatic shape was reached -- the equilibrium position of the load is shown by solid circles in these cases. The black cross is the sub-host point and the black star is the north pole. a) Positive loads, located initially at 1°E, 89°N. b) Negative loads, initially placed at 1°E, 1°N.}
    \label{fig:Pluto_thereback}
\end{figure*}

In this exercise, the loads are placed far away from their equilibrium positions, with the goal to address complete overturns of the body. In each simulation, Pluto is first spinning at angular velocity $\omega_0$ (Table 1) until it reaches an initial equilibrium state. In the first set of simulations (solid lines in Fig.~\ref{fig:Pluto_thereback}), the initial state is assumed to be hydrostatic and the elastic lithosphere to be created only after the equilibrium was established, leading to the formation of a "fossil" bulge \citep[e.g.][]{Matsuyama2014}. In the second set of simulations (dotted lines), the elastic lithosphere is included from the beginning of the simulation so that no fossil bulge is created. Fig.~\ref{fig:Pluto_thereback} shows the load paths predicted for positive (a) and negative (b) loading and different values of $M_\mathrm{h}/M$. The duration of each simulation is 2 Gyr.

\begin{table}
\centering
\begin{tabular}{c c c}
\hline
   \multicolumn{3}{c}{Internal structure of Pluto} \\
   Outer radius & 1188 & km \\
   Ice/water radius & 1038 & km \\
   Water/core radius & 858 & km \\
   Ice density & 950 & kg/m$^3$ \\
   Water density & 1000 & kg/m$^3$ \\
   Core density & 3360 & kg/m$^3$ \\
   Ice shear modulus & 3.49 & GPa \\
   Surface temperature, $T_\mathrm{surf}$ & 47 & $K$ \\
   Ice/water temperature, $T_\mathrm{bot}$ & 265 & $K$ \\
   Activation energy, $E$ & 59 & kJ \\
   Exponential prefactor, $A$ & $9\cdot 10^{-8}$ & m$^2$K/(Pa s) \\
   Grain size, $d$ & 10 & mm \\
   Cut-off viscosity, $\eta_\mathrm{max}$ & $10^{24}$ & Pa s \\
   Lithosphere thickness, $D_\mathrm{EL}$ & 50 & km \\
   Initial spin rate, $\omega_0$ & $1.14074 {\cdot} 10^{-5}$ & rad/s \\
\hline
   \multicolumn{3}{c}{Parameters of the disc load} \\
   Lateral extent & 25 & ° \\
   Thickness, $h$ & -300 to 300 & m \\
   Density & 1000 & kg/m$^3$ \\
   Growth time, $t_\mathrm{load}$ & 10 & Myr  \\
\hline
\end{tabular}
\caption{Model of Pluto, the parameters of ice shell are taken from \citet{Johnson2016} and \citet{Kihoulou2022b}.}
\label{tabPluto}
\end{table}

A comparison of panels a and b in Fig.~\ref{fig:Pluto_thereback} shows that the load paths for negative loads (Fig.~\ref{fig:Pluto_thereback}b) significantly differ from those for positive loads (Fig.~\ref{fig:Pluto_thereback}a). Even though the initial perturbation in longitude is only 1°, negative loads never move directly towards the north pole, regardless of the value of $M_\mathrm{h}/M$ and of whether the body has a remnant figure or not.

In case of models without a fossil bulge (dotted lines), the load paths straighten as $M_\mathrm{h}/M$ decreases when the load is positive, while the opposite is true for a negative load.

When a fossil bulge is considered, the load paths obtained for the positive loads are nearly straight. The negative loads first move along the equator and turn to the north only after approaching the center of the trailing hemisphere. This peculiar behaviour is analyzed below. The dependence of the solution on the load size is illustrated for the Charon/Pluto mass ratio ($M_\mathrm{h}/M=0.1212$), showing that the extent of reorientation decreases non-linearly with the load magnitude.

To understand the observed behaviour, one must focus on the way in which the inertia tensor of the load combines with that of the initial figure. While the inertia contribution of the positive load combines with the initial figure such as to overturn the major and the intermediate axes of the total inertia tensor $\bm{I}$, the negative load acts such as to overturn the intermediate and the minor axes. This is expressed by the evolution of the $C-B$ and $B-A$ differences during the growth time of the load (10 Myr, Table 1). When the load is negative, it is $B-A$ that decreases during the first few Myr, while for the positive load the $C-B$ difference narrows. It is thus the value of $B-A$ that stabilizes the longitude of the negative load. When the mass of the host body $M_\mathrm{h}$ is small, the initial figure approaches that of a tidally undeformed body, and the initial values of $A$ and $B$ are close to each other (cf.~the purple and orange vertical line segments in Fig.~\ref{fig:MIA}a). As a result, the negative load is more likely change its longitude as $M_\mathrm{h}/M$ decreases (cf.~the solid brown, purple, orange, and blue circles in Fig.~\ref{fig:Pluto_thereback}b).

When the load is positive, TPW occurs preferentially along 90° longitude, because reorientation in the plane perpendicular to $\bm{o}$ does not require any readjustment of the tidal bulge. As $M_\mathrm{h}/M$ decreases, this stabilization becomes less important, allowing for the load to reach its equilibrium position in a more straightforward manner.

The curves in Fig.~\ref{fig:MIA}a are quite complicated, because the evolution of the moments of inertia is affected by several processes: i) the linear rise of the load (full onset is marked by the vertical black dashed line), ii) the progress of TPW itself (Fig.~\ref{fig:MIA}b), and iii) readjustment of the tidal and rotational bulges in response to the TPW.

A more graphic explanation of the behaviour obtained in Fig.~\ref{fig:Pluto_thereback} is shown in Appendix \ref{sec:Janal}, where the principal directions of the equilibrium inertia tensor are plotted with respect to the initial figure of the body. In Fig.~\ref{fig:appendix}, we select a few load amplitudes and analyze how the individual contributions combine to form the total inertia tensor $\bm{I}$.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{figs/MIA_PN.png}
    \caption{a) Moments of the deviatoric part of the inertia tensor $\bm{I}$ during the TPW simulations plotted by the dotted orange and purple lines in Fig.~\ref{fig:Pluto_thereback} (i.e.~without a fossil bulge). Solid lines are for the positive load of 300 m thickness, dashed lines are for the negative load of the same magnitude. The values are normalized by the diagonal element of the spherical part of the inertia tensor, i.e.~by one third of the trace of $\bm{I}$. Full onset of the load, $t_\mathrm{load}$, is marked by the vertical dashed line. b) The temporal progress of the reorientation. We show the angular distance of the load from its final position for the four cases depicted in panel a.}
    \label{fig:MIA}
\end{figure}

The results in Fig.~\ref{fig:Pluto_thereback}b cast a doubt on the hypothesis of \citet{Kihoulou2022}, who have suggested that Sputnik Planitia first moved to the anti-Charon point, where it resided until the ice shell grew thicker and a fossil bulge was formed, and only then moved towards the north due to the disappearance of the ocean uplift for the thicker ice shell (and thus change of sign of the gravity anomaly). In fact, already a slight perturbation in longitude would trigger a rapid motion of the load in the direction perpendicular to the plane spanned by $\bm{\omega}$ and $\bm{o}$, making it unlikely for the basin to reach its present geographic location. Note that some perturbation is always required to start a reorientation process from the sub-host (or the anti-host) point, because these points represent the equilibrium states (stable or unstable, depending on the sign of the load). As shown in Fig.~\ref{fig:Pluto_thereback}b, a negative load can move directly from the anti-Charon point to the north, as suggested by \citet{Kihoulou2022}, only if the latitudinal position of the load is disturbed but its longitudinal position remains fixed at exactly 180° (cf.~also the blue line in Fig.~\ref{fig:Pluto_Ham}b).

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.45\textwidth]{figs/ham_pos.png}
    \includegraphics[width=0.45\textwidth]{figs/ham_neg.png}
    \caption{Same as Fig.~\ref{fig:Pluto_thereback}, only here it is the initial position and amplitude of the disc load that is varied. The ratio $M_\mathrm{h}/M$ is 0.1212 for all cases. Black pentagon marks the mirror image of Sputnik Planitia's center, i.e., its position in the first quadrant of polar coordinates. a) Disc thickness is 300 m. Initial longitude of the positive anomaly is taken as 0°, 20°, 40°, 60°, 80°, and 90°E. b) Negative anomaly of a varying magnitude, initially placed at 24°N, 4°E, with the exception of the blue line, for which the disc is placed at 24°N, 0°E.}
    \label{fig:Pluto_Ham}
\end{figure*}

In Fig.~\ref{fig:Pluto_Ham}, we explore the concept of \citet{Hamilton2016}. Note that we again employ the planes of symmetry of the problem and place Sputnik Planitia at 24°N, 4°E instead of 24°N, 176°E (black pentagon), because we choose to work in the first quadrant of polar coordinates. In Fig.~\ref{fig:Pluto_Ham}b, we investigate what would have happened if the load had first reached its present location and then changed from positive to negative. Inspection of the figure shows that the answer depends on the presence or absence of a fossil bulge and on the size of the load. In the absence of fossil figure (dotted lines), Pluto would undergo a large reorientation, moving Sputnik Planitia far away from its present position within a few hundred Myr. Regardless of its size, the negative load always tends to move to the east, making it unlikely that Sputnik Planitia could remain as close to the 180° meridian as it is at present. A similar trend is found for models with a fossil bulge (solid lines), but in this case, the final position of the basin depends on the disc thickness. The basin is displaced by less than 3° if the load amplitude is 30 m (solid purple line), while the distance between the initial and final positions of the basin is more than 25° for loads with magnitudes exceeding 100 m.

In Fig.~\ref{fig:Pluto_Ham}a, we impose a positive load at the latitude of 30°N (Pluto's coldest region) and calculate its path for different initial longitudes. In the absence of fossil figure (dotted lines), the load always ends up at the anti-Charon point, i.e., 24° away from the present center of Sputnik Planitia. The present position of the basin can be better predicted by models that include the effect of a fossil bulge. For a 300 m thick disc the equilibrium position has a lower latitude than 24°N regardless of the initial longitude, but for a thinner disc the present position of Sputnik Planitia can be reached (cf.~also Fig.~2a in \citet{Keane2016}). Our results suggest that either i) the elastic lithosphere had already been formed when the nitrogen ice began to accumulate, or ii) Pluto has been experiencing a continuous reorientation, while the geographic location of the nitrogen ice does not change due to an ongoing redeposition.

\section{The rate of reorientation}\label{sec:discussion}

The TPW rate depends on the load amplitude, the internal structure and the shape of the fossil figure. As to the internal structure, two parameters are of particular importance: the thickness of the ice shell and the cut-off viscosity $\eta_\mathrm{max}$, representing the effective viscosity of low-temperature ice. To a first approximation, these two parameters determine the characteristic time of viscoelastic readjustment of the tidal and equatorial bulges and thus control the rate of reorientation.

In case of Pluto, neither of these parameters is known with certainty. The ice thickness is estimated to be between 100 and 180 km \citep[][]{Johnson2016, Denton2020, Kihoulou2022b}, while $\eta_\mathrm{max}$ depends on a number of factors such as the overall stress level, composition and deformation history of ice, and is not well constrained by observations.

The simulations presented in section \ref{sec:pluto} (Figs \ref{fig:Pluto_thereback}-\ref{fig:Pluto_Ham}) are computed for $\eta_\mathrm{max}=10^{24}$ Pa s and an ice thickness of 150 km. In this case, the time required to reach the equilibrium orientation varies from a few Myr to 2 Gyr.  When $\eta_\mathrm{max}$ is decreased by two orders of magnitude (Fig.~\ref{fig:robustness}a,b) the rate of reorientation speeds up, with the final orientation being reached within a few tens of Myr. However, the qualitative behaviour does not change. When the ice shell thickness is varied the load paths are even less affected (Fig.~\ref{fig:robustness}c,d).

Note that the changes in load's latitude and longitude can be highly non-linear in time. For instance, the curved segment of the load path depicted in blue in Fig.~\ref{fig:Pluto_thereback}a takes considerably more time than the preceding drop in latitude. When the shell thickness is increased to 300 km, this last stage of reorientation is not completed within the simulation time of 2 Gyr. Due to non-linearities, reaching the equilibrium orientation may require time that is comparable to the age of the solar system in some cases. With a dynamic method, such peculiarities can easily be revealed, but within the fluid-limit framework they would remain unnoticed.

In all the presented simulations, loading is imposed only gradually in time, with the load magnitude rising linearly over the time $t_\mathrm{load}=10$ Myr (Table 1). If $t_\mathrm{load}$ was significantly reduced, the values of $B$ and $A$ would intersect at some point in time, and thus the two principal directions corresponding to the intermediate and minor moments of inertia would instantaneously revert (``flip''). Such an event would turn the rotational dynamics unstable, triggering a 90° wobble that would disturb the body considerably (similarly, the principal directions of inertia could flip due to a crater formed by a sudden impact). The dynamics of wobble on a tidally locked body are outside the scope of the present paper.

\section{Energy Balance}\label{sec:energy}

In this section, we derive an energy conservation law for a tidally locked planet or moon. The derived formula allows a detailed examination of the different types of energy (thermal, elastic, kinetic, rotational, tidal, and gravitational) and we present it in a form that is suitable for testing the accuracy of any TPW solver.

Conservation of energy in a non-inertial, rotating frame can be written as:
\begin{linenomath*}\begin{multline}
\dt{}\int_{v(t)}\left(\rho\epsilon+\frac{_1}{^2}\rho\bm{v}\cdot\bm{v}\right)\, \mathrm{d}v \\
=\int_{v(t)}\left(-\nabla U-\nabla\Theta-\nabla\Psi-2\bm{\omega}\times\bm{v}- \dt{\bm{\omega}}\times\bm{r}\right)\cdot\bm{v}\,\rho \,\mathrm{d}v,
\label{eq_bal}
\end{multline}\end{linenomath*}
where $v(t)$ is the volume of the body at time $t$. The terms in the integral on the left-hand side (LHS) correspond to the specific internal energy \citep[Eqs (20) and (21) in][]{Patocka2018} and the kinetic energy associated with the deformation, respectively, while the integral on the right-hand side (RHS) represents the power of the body forces acting in a rotating system.

It is easy to see that the power delivered by the Coriolis force is zero, as the force is always perpendicular to the velocity vector $\bm{v}$. Since, by definition, the relative angular momentum, $\bm{h}$, in the Tisserand frame vanishes, the total power of the Euler force is also zero:
\begin{linenomath*}\begin{multline}
\int_{v(t)}\left(-\rho\dt{\bm{\omega}}\times\bm{r}\right)\cdot\bm{v}\,\mathrm{d}v = -\dt{\bm{\omega}}\cdot\int_{v(t)}\bm{r}\times(\rho\bm{v})\,\mathrm{d}v \\
= -\dt{\bm{\omega}}\cdot\bm{h}\ = \bm{0}. \label{eqEulerPo}
\end{multline}\end{linenomath*}
The power of the centrifugal force is \citep[see Eq.~(A7) in][]{Patocka2021}:
\begin{linenomath*}\begin{multline}
   \int_{v(t)}{-\nabla\Psi \cdot\bm{v} \rho\,\mathrm{d}v} = \int_{v(t)}{\left[ (\bm{r}{\cdot}\bm{v})\,\omega^2 -(\bm{\omega}{\cdot}\bm{v})\cdot(\bm{\omega}{\cdot}\bm{r})\right] \rho\,\mathrm{d}v} \\
   = \frac{1}{2}\bm{\omega}\cdot \dt{\bm{I}}\cdot\bm{\omega}. \label{eqCentrPo}
\end{multline}\end{linenomath*}
Using Eq.~\eqref{eqCentrPo} and taking into account that the power $-\nabla U \cdot \bm{v}$ that is associated with the deformation of the body in its own gravitational potential, $U$, is provided the gravitational energy $\int \rho U/2 \,\mathrm{d}v$ \citep[Eq.~(A21) in][]{Patocka2018}, Eq.~\eqref{eq_bal} can be rewritten as follows:
\begin{linenomath*}\begin{multline}
\dt{}\int_{v(t)}\left(\epsilon+\frac{_1}{^2}\bm{v}\cdot\bm{v} + \frac{1}{2} U\right)\,\rho\, \mathrm{d}v \\
=\frac{1}{2}\bm{\omega}\cdot \dt{\bm{I}}\cdot\bm{\omega} - \int_{v(t)} \nabla \Theta\cdot\bm{v} \rho\,\mathrm{d}v . \label{eq_balrot}
\end{multline}\end{linenomath*}
Adding the rate of the rotational energy, $E_\mathrm{rot}=(\bm{\omega}{\cdot} \bm{I}{\cdot}\bm{\omega})/2$, to both sides, we get:
\begin{linenomath*}\begin{align}
\dt{}\int_{v(t)}\left(\epsilon+\frac{1}{2}\bm{v}\cdot\bm{v} +\frac{1}{2} U - \Psi\right)\,\rho\,\mathrm{d}v \nonumber\\ = \frac{1}{2}\bm{\omega}\cdot \dt{\bm{I}}\cdot\bm{\omega} + \dt{} \left( \frac{1}{2}\bm{\omega}\cdot \bm{I}\cdot\bm{\omega} \right) - \int_{v(t)} \nabla \Theta\cdot\bm{v} \rho\,\mathrm{d}v. \label{eq_balrot2}
\end{align}\end{linenomath*}
Note that the rotational energy is obtained by integrating $-\Psi$, while in the case of the gravitational energy the integrand is $+U/2$. Combining the first two terms on the RHS and using the symmetry of $\bm{I}$, we get:
\begin{linenomath*}\begin{multline}
\dt{}\left(E_\mathrm{int} + E_\mathrm{kin} + E_\mathrm{grav} + E_\mathrm{rot}\right) \\ = \dt{( \bm{I} \cdot \bm{\omega} )} \cdot\bm{\omega} - \int_{v(t)} \nabla \Theta\cdot\bm{v} \rho\,\mathrm{d}v. \label{eq_balrot3a}
\end{multline}\end{linenomath*}
The first term on the RHS can be expressed in terms of the torque $\bm{M}$, because the projection of the LE onto $\bm{\omega}$ reads:
\begin{linenomath*}\begin{equation}\label{eq_LEdotw}
\bm{M} \cdot \bm{\omega} = \left( \dt{( \bm{I} \cdot \bm{\omega} )} + \bm{\omega}\times(\bm{I}\cdot\bm{\omega})\right) \cdot\bm{\omega} = \dt{( \bm{I} \cdot \bm{\omega} )} \cdot\bm{\omega}.
\end{equation}\end{linenomath*}
Eq.~\eqref{eq_balrot2} then takes the form
\begin{linenomath*}\begin{multline}
\dt{}\left(E_\mathrm{int} + E_\mathrm{kin} + E_\mathrm{grav} + E_\mathrm{rot}\right) \\ = \bm{M} \cdot\bm{\omega} - \int_{v(t)} \nabla \Theta\cdot\bm{v} \rho\,\mathrm{d}v. \label{eq_balrot3b}
\end{multline}\end{linenomath*}
Since Eq.~\eqref{eq_balrot3b} is valid if and only if the LE is satisfied, it can be used to test the accuracy of the TPW solution: For a tide-free body (i.e., $\bm{M}=0$ and $\Theta=0$), the sum of $E_\mathrm{int} + E_\mathrm{kin} + E_\mathrm{grav} + E_\mathrm{rot}$ must be conserved. As shown by \citet{Patocka2021}, this is satisfied only if the changes in $\Psi$ are treated in a self-consistent manner and to a high degree of precision in the course of deformation.

If the body is subject to tidal forcing, the RHS in Eq.~\eqref{eq_balrot3b} is non-zero, implying that the sum of the energies on the LHS is no more conserved. Analogously to Eq.~\eqref{eqCentrPo}, the second term on the RHS can be expressed as
\begin{linenomath*}\begin{equation}\label{eqTidPo}
    \int_{v(t)} \nabla \Theta\cdot\bm{v}\, \rho \,\mathrm{d}v  = \frac{1}{2}\bm{o}\cdot \dt{\bm{I}}\cdot\bm{o} + \dt{}\int_{v(t)} \frac{1}{3} r^2 o^2 \rho \, \,\mathrm{d}v.
\end{equation}\end{linenomath*}
For convenience, we define the ``tidal'' energy
\begin{linenomath*}\begin{equation}\label{eqTidEn}
    E_\mathrm{tid} = \int_{v(t)} \Theta \rho\,\,\mathrm{d}v = \frac{1}{2} \bm{o}\cdot\bm{I}\cdot\bm{o} - \int_{v(t)} \frac{1}{3} r^2 o^2 \,\rho\,\,\mathrm{d}v.
\end{equation}\end{linenomath*}
The tidal energy is the potential energy possessed by the body due to the gravity of the host body. Therefore, $\Theta$ is taken with a positive sign in Eq.~\eqref{eqTidEn}, but the usual factor of 1/2 is omitted. Recalling that $\bm{M}=\bm{o}\times(\bm{I}\cdot\bm{o})$ and using the identity $\bm{a}\cdot(\bm{b}\times\bm{c}) = \bm{c}\cdot(\bm{a}\times\bm{b})$, the energy balance, Eq.~\eqref{eq_balrot3b}, can be rearranged as follows:
\begin{linenomath*}\begin{multline}
\dt{}\left( E_\mathrm{int} + E_\mathrm{kin} + E_\mathrm{grav} + E_\mathrm{rot} + E_\mathrm{tid} \right) \\ =(\bm{I}\cdot\bm{o})\cdot \left(\bm{\omega}\times\bm{o} + \dt{\bm{o}} \right) . \label{eq_balrot3}
\end{multline}\end{linenomath*}
The RHS of Eq.~\eqref{eq_balrot3} describes the energy that is exchanged with the host body throughout the reorientation process.

Note that in deriving Eq.~\eqref{eq_balrot3}, no specific assumptions have been made about the form of the LE. Under the $o\omega||\mathrm{mMIA}$ approximation, Eq.~\eqref{eq_balrot3} reduces to:
\begin{linenomath*}\begin{equation}\label{eq_balwMIA}
\dt{}\left( E_\mathrm{int} + E_\mathrm{kin} + E_\mathrm{grav} + E_\mathrm{rot} + E_\mathrm{tid} \right)= Ao\dt{o},
\end{equation}\end{linenomath*}
where $A$ is the minimum moment of inertia and the RHS can be evaluated from the solution of Eqs~\eqref{eqtomap_1}--\eqref{eqtomap_3} by using the relationship between $\omega$ and $o$, Eq.~\eqref{eqTidlock}.

In reality, reorientation may lead to small deviations from the $o\omega||\mathrm{mMIA}$ assumption, but these deviations are likely to be small because they are effectively dampened by dissipative processes.

The use of the energy balance, Eq.~\eqref{eq_balwMIA}, is illustrated in Fig.~\ref{fig:energy}.
Inspection of the figure shows that most types of energy included in Eq.~\eqref{eq_balwMIA} are of the same order of magnitude and the LHS (dark blue line) is equal to the RHS (thin red line), confirming that the numerical solution is correctly implemented. The relative changes of the spin rate in this simulation are of the order $10^{-5}$, which corresponds to a change in the Pluto-Charon distance of $\sim 100$ m. Note that the rotational and tidal energies are minimized during the reorientation, while the gravitational energy of the disc load in the gravitational potential of the studied body is increased, because the disc is lifted against the surface gravity by the readjustment of the rotational and tidal bulges.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{figs/enbal_P300_hfP.png}
    \caption{Contributions to the total energy and their evolution for model with $M/M_\mathrm{h}=0.1212, h=300$ m and no fossil bulge. The terms in Eq.~\eqref{eq_balwMIA} are integrated over time and shown in different colors. The initial time of the plot $t=10$ Myr corresponds to the full onset of loading. $E_\mathrm{kin}$ is negligible and thus not shown for figure clarity.}
    \label{fig:energy}
\end{figure}

\section{Summary}\label{sec:Sum}

We have developed a new method to investigate the dynamics of reorientation of a tidally locked body. Our method gives the same results as that of \citet{Hu2017a}, but is much simpler, requiring only the evaluation of the eigenvalues and eigenvectors of the inertia tensor. The theoretical framework of the method is completed by an analysis of the energy balance, which can be used to test the numerical solution. A new release of the code LIOUSHELL that was used to perform the simulations is freely available on GitHub \citep{Patocka_LIOUSHELL20}.

Unlike most studies that only investigate the final (equilibrium) orientation of the body \citep[recently, e.g.,][]{Schenk2020,Matsuyama2021,Johnson2021}, the $o\omega||\mathrm{mMIA}$ method presented here can be used to predict the wander of the north pole and the sub-host point, i.e., the evolution of the rotation and tidal vectors in the body-fixed frame, or the path of the surface load in the geographic frame. Since the rate of reorientation depends on the internal structure and material parameters of the body, the modeling of reorientation dynamics can be used to explore the properties of the body and to test different scenarios of its evolution.

We investigate the path of a positive load placed near one of the rotation poles and the path of a negative load placed near the tidal axis, varying the mass ratio $M_\mathrm{h}/M$ where $M$ is the mass of the body under consideration and $M_\mathrm{h}$ is the mass of the body whose gravitational pull generates the tidal force. When $M_\mathrm{h}/M$ decreases, the load path straightens for the positive load and becomes more curved for the negative load. The negative load tends to rapidly change its longitude when $M_\mathrm{h}/M$ is small, because then the equatorial moments of inertia are close in value.

The asymmetric response of a tidally locked body to positive and negative loads has important consequences for the reorientation of Pluto ($M/M_\mathrm{h}=0.12$) and the present-day location of the Sputnik Planitia basin. In particular, we show that it is unlikely that the load associated with Sputnik Planitia changed sign from positive to negative after reaching the anti-Charon point, as proposed by \citet{Kihoulou2022}. A more plausible scenario is that the basin represents a positive load and its center stays away from the equator due to a fossil bulge that existed already at the time of Sputnik Planitia's formation \citep{Keane2016}. If it is the uplift of a dense interior ocean that compensates the negative topography of the basin \citep{Nimmo2016}, then the ice shell of Pluto must be less than 200 km thick and insulated by a layer of clathrates at its bottom \citep{Kamata2019,Kihoulou2022b}. Our results do not support the hypothesis that Sputnik Planitia was formed early in Pluto's history \citep{Hamilton2016}, because in the absence of fossil figure the basin would move directly to the anti-Charon point.

\section*{Acknowledgements}

V.P.~and M.K.~acknowledge support by the Czech Science Foundation through project nr.~22-20388S.

\printcredits

redpajama_set_name: RedPajamaArXiv


Authored by Abdullah Mohammadi, Ruta Nimkar and Emily Savage, this publication in the Migration Research Series analyses the perceptions that Afghan smugglers have of themselves and of their relationships with their communities in Afghanistan. The paper is based on interviews with 23 smugglers in three sites in Afghanistan and considers community dynamics and low-level smugglers rather than high-level organizers of smuggling networks. It highlights that smuggling networks have a long-standing and respected place in Afghan culture. The paper also provides an analysis of factors that affect perceptions of and trust in smugglers among Afghan society. It concludes with some implications to support policy responses and programming concerning migrant smuggling and migration in Afghanistan.


text:
Trending City Guides
Cannabis + Drugs
Everything to know about legal cannabis consumption in Barcelona
Photo: content_creator/Shutterstock
Lena Papadopoulos
Since the early 2000s, the Catalan region of Spain, particularly Barcelona, has been growing as a hotspot for cannabis tourism. You may have heard about the city's somewhat open approach to marijuana regulation, but note that cannabis consumption there isn't the same as in Canada or in US states where recreational weed is legal. It has similarities to the marijuana scene in Amsterdam, but there are some major key differences.
If you plan on exploring Barcelona's cannabis culture on your next visit, this is everything you need to know.
Is cannabis legal in Spain?
The short answer is no — cannabis is technically illegal. Trafficking, selling, or carrying cannabis in Spain is a crime. Smoking or displaying it in public will cost you a large fine, any cannabis you have will be seized, and you may be subject to jail time. The personal use of marijuana, however, has been decriminalized.
It's acceptable to have or smoke marijuana in moderate amounts in private spaces. This includes private residences and private cannabis clubs. In Barcelona, you're allowed to grow up to two plants in your own private residence, or you can forfeit your right to grow marijuana and opt to join a cannabis club instead. As a visitor to the city, it's not likely you're going to grow your own cannabis, so if you'd like to enjoy a blunt or two while in Barcelona, your best (and safest) bet is to join a club.
Everything you need to know about Barcelona's cannabis clubs
The only decriminalized way to acquire and/or consume cannabis outside of a private residence in Barcelona is by way of a cannabis club. These clubs should not be confused with medical marijuana dispensaries. Because cannabis is officially illegal, Spain does not offer medical marijuana, doctors cannot prescribe cannabis, and there are no dispensaries in the city.
Cannabis clubs in Barcelona are private entities committed to the protection, security, and anonymity of their members. The personal information of cannabis club members is not shared with any other entity, including the government. Authorities have no access to member information of any kind, and they can't enter a cannabis club without specific authorization to do so. For this reason, clubs make a concerted effort to remain discreet and avoid causing any disruption to nearby residents.
Because it's a private, secure entity, one cannot just casually walk into a cannabis club. You must have both an ID and a membership in order to enter. You can only acquire membership to a club through sponsorship or by referral through a current club member. You may be thinking it won't be easy to obtain membership as a tourist, but it's really not complex or difficult to do so.
How to join a cannabis club in Barcelona
Joining a cannabis club in Barcelona can be as simple as sending an email. Some online organizations, such as Marijuana Games, will offer to act as your sponsor for membership provided you meet their basic requirements. All you have to do is send some preliminary information to the email address listed, and they will set up a registration appointment with a club on your behalf.
You will receive an email reply from your "sponsor" letting you know when your appointment is and at which club (you don't get to choose which club you'll join through this method).
At your appointment, you'll fill out some paperwork provided by the club. Again, this information is protected and secure. You'll also be required to present valid, government-issued identification (passport, driver's license, state ID, etc.). It's important to note that you cannot register without a residential address in Spain. This doesn't mean you have to be a resident, but the address of a hostel or hotel will not suffice. This can be the address of an Airbnb or a local friend's residence. You won't be asked to authenticate this address, and nothing will be sent to the location you provide.
In order to join the club, you must also pay an annual membership fee. This fee is meant to offset the club's expenses and does not cover the cost of any cannabis. The membership fee tends to run between 25 and 50 euros per year, but it can sometimes be lower or higher depending on the club. You'll need to pay this fee in cash.
How to get cannabis from a Barcelona cannabis club
Photo: VTT Studio/Shutterstock
For legal reasons, you do not "buy" cannabis from a club. By joining a club, you become part of a group that grows and disperses cannabis among club members. In this sense, you don't pay for the marijuana itself. Instead, whatever you pay to acquire cannabis is used to cover the cost of growing the marijuana and maintaining the club's facilities.
Additionally, clubs place limits on how much cannabis you can acquire. Because it's only meant for personal use, you're typically restricted to five grams per day and can't exceed a certain amount per month. As this isn't regulated by any law, the exact number depends on each club.
Because you can face consequences for possessing marijuana in public, the club will advise you not to take it outside the property. Most members do take cannabis from the club, however. If you choose to take cannabis from the club, it's recommended to keep it in your undergarments, as that's the one place police can't search (unless they have a warrant to do so). If you're caught with marijuana, you should say as little as possible about where you got it, and do not mention the name of your club.
While consuming cannabis in Barcelona may sound risky, there's really nothing to worry about as long as you join a club (which isn't difficult) and follow the basic guidelines.
More like this: This guy smoked cannabis at every single coffeeshop in Amsterdam so you don't have to
What did you think of this article?
More on Barcelona
8 rooftop bars in Barcelona to take in the city's nightlife
Tara Tadlock
The 8 best places to meet travelers in Barcelona
Dayana Aleksandrova
14 awesome things to do on a Friday night in Barcelona when you're stone cold sober
12 things moving to Barcelona taught me about sex
This is what the Sagrada Familia will look like when completed in 2026
Morgane Croissant
15 side effects I had as a traveler in Barcelona
8 things Barcelona locals want you to know
9 best outdoor spots to hike or chill around Barcelona
20+ spots for the best nightlife in Barcelona
Top 10 things you can do for free in Barcelona
Everything you always wanted to know about tapas but were afraid to ask
Ana Bulnes
Barcelona's 4 coolest neighborhoods to stay, eat, and party
More on Cannabis + Drugs
Everything to know before going to a state with recreational marijuana
Nickolaus Hines
Oakland is the second US city to legalize magic mushrooms
Eben Diskin
The biggest 4/20 parties in the US
This guy smoked cannabis at every single coffeeshop in Amsterdam so you don't have to
Can you fly with marijuana in states that it's legal?
Cannabis dispensary etiquette: dos and don'ts to know before you go
Amanda Gabriele
Denver decriminalized magic mushrooms. Here's why this is huge.
Tim Wenger
Denver becomes the first city in the US to decriminalize magic mushrooms
Everything you need to know about doing magic mushrooms
Jessica Devnani
Australian cafe reportedly sold a pot brownie to a family and the kids were not alright
Kosher cannabis is real, but it's not as simple as getting it blessed by a rabbi
Matthew Meltzer
Everything to know about driving with cannabis in your car in legal recreational states
Be a local. Wherever you go.
Fearless travel in your inbox.
Join our tribe of millions.
Submit your work to Matador!
Matador's growing Creators Community is the place to connect with fellow travel journalists.
Matador Network © 2006-2020
stoked by
We use cookies to understand how you use our site and to improve your experience. This includes personalizing content and advertising. For more information on how we use cookies consult our revised Privacy Policy and Terms of Service.
Disable Cookies I accept

redpajama_set_name: RedPajamaCommonCrawl


For beginner cooks with hectic schedules, Do-It-Yourself pizza is a quick and tasty weeknight dinner solution. In the time that it takes to have pizza delivered, you can actually make your own!
and whatever toppings you like.
The nice thing about DIY pizza is that you can have fun personalizing it. You can incorporate your favorite flavors and textures, or whatever happens to be in the fridge. There's no fixed recipe to follow! Let's get started.
Ready-made crust: The bread section of the grocery store typically offers several brands of ready-made crusts that you simply top with whatever you like and bake according to the instructions. Some of these crusts even come with a package of tomato sauce! They're usually sold in rounds that make enough for several people.
If you like thin-crust pizza, you can't get thinner than this! It's a good idea to pre-heat the tortilla in a toaster oven or regular oven at 350F just for a few minutes to make it a bit more solid, then turn the oven up to 450F. Put the wrap(s) on a foil or parchment paper-lined baking sheet.
Tortillas are the base for Humbecue Pizzas — dressed with hummus, barbecue sauce, spinach, ham and feta cheese. Photo courtesy of Jenn at Eating Bender.
Then add whatever sauce, cheese and toppings you like. Go easy on toppings so that the whole thing doesn't get soggy. Put the tortilla in the oven for 10 minutes and check to see when it's done – the edge of the tortilla should become brown and crispy. Here's a recipe for Humbecue Tortilla Pizzas, which are made with hummus, barbecue sauce, spinach, ham and feta cheese.
Pita bread is another great option for those who like thin, crispy pizzas. There's no need to toast the pita before putting on the toppings; follow the same baking directions as with the tortilla pizza.
Here, plain old pita gets dressed up with Italian sausage, peppers and carrots.
This Pita Pizza (shown above) topped with Italian sausage, peppers and carrots looks simple and scrumptious.
For a more substantial crust, use focaccia bread, which is often sold in big rectangular pieces. When it comes to baking the pizza, preheat the oven to 450F and cook for 10 to 15 minutes.
Here at startcooking.com we have already posted a video on English Muffin Pizzas, which are topped with tomato sauce, ground beef and cheese.
These mini pizzas are irresistible, and if you don't have any ground beef on hand, you can make them with tuna.
If you would prefer to substitute crumpets for English muffins, go right ahead, old chap!
A lot of bagel shops now have bagel pizzas.
Making your own is easy as long as you have the necessary ingredients on hand! Buy your bagels ahead of time and store them (sliced!) in the freezer. Bagels freeze beautifully as long as they are wrapped tightly in plastic wrap and then put in a zip lock bag.
Pizza traditionally has a layer of tomato sauce between the crust and the cheese. You can buy ready-made pizza sauce in a can or jar, or use tomato sauce (the kind you would use for pasta). If you don't have either of these on hand, you could try tomato paste or salsa.
Some pizza lovers prefer their pizzas "white" – with no tomato sauce. White pizzas may be topped with pesto, ricotta or Parmesan cheese and pine nuts. Or with alfredo sauce (which you can buy ready-made) and other toppings.
In this recipe for Spinach, Mushroom and Goat Cheese Pizza , from Cooking By the Seat of my Pants, the sauce is a combination of sour cream and Greek vinaigrette.
Spinach, mushroom and goat cheese pizza — not something you see on a takeout menu!
Photo courtesy of Jerry D. Russell at Cooking by the Seat of my Pants.
If you love basil, try Pioneer Woman's version of pizza using pesto rather than tomato sauce, fresh mozzarella, tomatoes, basil and parmesan.
Mild: Most pizza places use mozzarella cheese, which turns nice and gooey when it melts, and has a very mild flavor. Other mild cheeses that can substitute for mozzarella are colby, cheddar, gouda, edam and Monterey jack. You can grate the cheese yourself or buy bags of pre-grated cheese. Ricotta, bocconcini and fresh mozzarella are other mild options that offer a gourmet twist because of their rich, milky textures.
Strong: Try mixing mozzarella with stronger-tasting cheeses, or going for something really distinctive. You'll find that you can use a smaller amount of cheese if it's strong-tasting. In our video on English Muffin Pizzas, startcooking.com suggests blue cheese with chopped fresh sage and walnuts. Some other strong cheeses to try are: parmesan, provolone, goat cheese, aged cheddar and gorgonzola.
This is where you can really make a pizza your own. Here's a list of topping ideas to get you started, as well as some great-tasting combinations.
Tip: Any meat toppings should be cooked before being added to the pizza.
Good luck with your pizza creations!
What a great post! I just posted about pizzas using the new Pillsbury thin crust, but I have tried Boboli for breakfast pizza and english muffins for a lunch pizza. Great inspiration using tortillas as a very thin crust!!
Hi Kathy….is there a substitute for parchment paper? We always have all sizes of tortillas at our house but alas NO paper? HELP!!
English Muffin Pizzas were an occasional Sunday- night-supper-family-treat when I was growing up. (My mom would always cook a roast with all the trimmings for the big Sunday lunch.) My Dad was usually in charge of making the Pizza and we all loved it! To this day I can't make them without thinking of my Dad!
Cynthia, it all depends on what you are cooking. For a tortilla, use tin foil instead.
Mmm…pizza!! This post has some really great tips, thank you for sharing! And an even bigger thanks for including one of my recipes I just sent you an email!
Jenn, your recipe is a big hit – thanks for sharing!
I have made Pizza using the Boboli crust before and yes it's super simple. Now I would like to advance and try to make my own crust. I gave it a shot today, but the written instructions online are confusing if you have nevermade dough before. Would you mind to post a video on how to make pizza crust?
Making one's own pizza crust is a great idea BUT we do have to remember this site is called "startcooking"!
Working with yeast can be tricky for a beginner, but I will add it to the ever growing list of requests!
Pizza crust is easy to make for beginners if you have a bread machine. You still have to stretch the dough, but you'd have to do that if you bought pizza dough from a pizza parlor. And with a bread machine, all the kneading is done in the machine.
Thanks for the great tip! Unfortunately most beginner cooks do not own a bread machine. In fact, I have discovered that a blender is usually the most common piece of kitchen equipment.
Thanks 4 the great idea- pizza is my favorite junk food!! After reading this, i've made my own pizza at least 3times already using ragu pizza sauce & pita bread with various toppings – i luv em! & so do my friends! Who needs 2order out? Very easy to make – Wow!
Even got my husband, who has NEVER had pizza in his life, taste it & now eating it too!! (Never succeeded in getting him to try pizza until now for 7yrs! – I guess actually seeing it being made did the trick).
This sight is great – simple, clear, easy-to-follow instructions & videos! Thank you!
One more thing: I brush olive oil (just because this is what I buy) on both sided of the pita bread b4 adding the toppings. Not sure if it makes a difference with how the crust turns out, as i've not tried it without it.
That's great news! Thanks for the olive oil tip – I'll have to give it a try.


text:
Shabbat Dinner with
Professor Jeffrey Hoffman
Join us for Shabbat Dinner on Friday, November 30 at 7:30 with
special guest, Dr. Jeffrey Hoffman.
Dr. Jeffrey Hoffman made five flights as a space shuttle astronaut, including the first mission to repair the Hubble Space Telescope in 1993, when the orbiting telescope's flawed optical system was corrected. Trained as an astrophysicist, he also flew on the 1990 Spacelab shuttle mission that featured the Astro-1 ultraviolet astronomical observatory in the shuttle's payload bay. Over the course of his five missions he logged more than 1,211 hours and 21.5 million miles in space. He was also NASA's first Jewish male astronaut, and the second Jewish man in space after Soviet cosmonaut Boris Volynov.
He joined the MIT Aeronautics and Astronautics faculty in 2001 as a senior lecturer, and since 2002 has been a Professor of the Practice in that department. His research specialties include human space flight operations, space flight technology, human-machine interactions, extravehicular activity, and conducting laboratory research in space. His teaching interests include space systems design and space policy. Dr. Hoffman instructed a course in systems engineering on the space shuttle that is available for free in video format from academic earth.
Shabbat Dinner is open to MIT affiliates by RSVPing below.