# A freely precessing magnetar following an X-ray outburst

Gregory Desvignes, Patrick Weltevrede, Yong Gao, David Ian Jones, Michael Kramer, Manisha Caleb, Ramesh Karuppusamy, Lina Levin, Kuo Liu, Andrew G. Lyne, et al.

► To cite this version:

Gregory Desvignes, Patrick Weltevrede, Yong Gao, David Ian Jones, Michael Kramer, et al. A freely precessing magnetar following an X-ray outburst. Nature Astronomy, 2024, 8, pp.617-627. 10.1038/s41550-024-02226-7. insu-04611879

HAL Id: insu-04611879

https://insu.hal.science/insu-04611879

Submitted on 14 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A freely precessing magnetar following an X-ray outburst

Gregory Desvignes, Patrick Weltevrede, Yong Gao, David Ian Jones, Michael Kramer, Manisha Caleb, Ramesh Karuppusamy, Lina Levin, Kuo Liu, Andrew G. Lyne, et al.

► To cite this version:

Gregory Desvignes, Patrick Weltevrede, Yong Gao, David Ian Jones, Michael Kramer, et al. A freely precessing magnetar following an X-ray outburst. Nature Astronomy, 2024, 8, pp.617-627. 10.1038/s41550-024-02226-7. insu-04611879

HAL Id: insu-04611879

https://insu.hal.science/insu-04611879v1

Submitted on 14 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A freely precessing magnetar following an X-ray outburst

Received: 16 November 2022

Accepted: 15 February 2024

Published online: 8 April 2024

<figure><img src="image_3.png" /></figure>

Check for updates

Gregory Desvignes Ⓓ<fnref n="1" />,²✉, Patrick Weltevrede Ⓓ<fnref n="3" />, Yong Gao Ⓓ⁴,⁵,

David Ian Jones⁶, Michael Kramer Ⓓ<fnref n="1" />,<fnref n="3" />, Manisha Caleb Ⓓ⁷,⁸,

Ramesh Karuppusamy Ⓓ<fnref n="1" />, Lina Levin Ⓓ<fnref n="3" />, Kuo Liu Ⓓ<fnref n="1" />, Andrew G. Lyne Ⓓ<fnref n="3" />,

Lijing Shao Ⓓ<fnref n="1" />,⁵,⁹, Ben Stappers<fnref n="3" /> & Jérôme Pétri Ⓓ¹⁰

Magnetars—highly magnetized neutron stars—are thought to be the most likely progenitors for fast radio bursts (FRBs). Freely precessing magnetars are further invoked to explain the repeating FRBs. We report here on new high-cadence radio observations of the magnetar XTE J1810−197 recorded shortly after an X-ray outburst. We interpret the polarization variations of the magnetar radio emission as evidence for the magnetar undergoing free precession following the outburst while its magnetosphere slowly untwists. The observations of precession being damped on a timescale of months argue against the scenario of freely precessing magnetars as the origin of repeating FRBs. Using free-precession models based on relaxing ellipticity with a decay of the wobble angle, we find the magnetar ellipticity to be in good agreement with theoretical predictions from nuclear physics. Our precise measurement of the magnetar’s geometry can also further help in refining the modelling of X-ray light curves and constrain the star’s compactness.

Magnetars are rare, typically slowly spinning neutron stars whose emission is thought to be powered by the decay of their large (≥10¹² G) magnetic fields<fnref n="1" /> in opposition to normal rotation-powered pulsars. They occasionally undergo bright X-ray outburst phases, supposedly originating from a sudden quake in the star’s crust, resulting in a twisted magnetosphere that fuels the outburst². Out of the 24 confirmed magnetars currently known, only 6 have been shown to emit in the radio band<fnref n="3" />.

The discovery of a bright fast radio burst (FRB) from the galactic magnetar SGR 1935+2154 recently gave strong credence to radio magnetars being the progenitors of at least some extra galactic FRBs⁴ (see, for example, ref. 5 for a review on FRBs). Models of freely precessing

magnetars were then put forward to explain the periodicities in the activity observed in the train of pulses from some repeating FRBs⁶,⁷.

The first recorded X-ray outburst of the magnetar XTE J1810−197 (position 18h09m05.17s − 19°43′51.8″, J2000) happened around late 2002 (ref. 8) and bright radio pulsations were detected 3 years later⁹ with a periodicity of $P_s = 5.54$ s. The radio emission slowly decayed with time until it ceased to be detected in late 2008 (ref. 10). Radio pulsations from XTE J1810−197 were again detected on 8 December 2018 (Modified Julian Date (MJD) 58,460.6)¹¹, following a second X-ray outburst known to have occurred between 20 and 26 November 2018 (ref. 12).

We present a set of 62 polarimetric observations of XTE J1810−197 recorded with the Lovell and Effelsberg radio telescopes at

<footnote n="1">¹Max-Planck-Institut für Radioastronomie, Bonn, Germany. ²Laboratoire d'Études Spatiales et d'Instrumentation en Astrophysique, Observatoire de Paris, Université Paris-Sciences-et-Lettres, Centre National de la Recherche Scientifique, Sorbonne Université, Université de Paris, Meudon, France.</footnote>

<footnote n="3">³Jodrell Bank Centre for Astrophysics, The University of Manchester, Manchester, UK. ⁴Department of Astronomy, School of Physics, Peking University, Beijing, China. ⁵Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing, China. ⁶Mathematical Sciences and STAG Research Centre, University of Southampton, Southampton, UK. ⁷Sydney Institute for Astronomy, School of Physics, The University of Sydney, Sydney, New South Wales, Australia. ⁸ASTRO3D: ARC Centre of Excellence for All-sky Astrophysics in 3D, Canberra, Australian Capital Territory, Australia. ⁹National Astronomical Observatories, Chinese Academy of Sciences, Beijing, China. ¹⁰Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, Strasbourg, France. ✉e-mail: gdesvignes@mpifr-bonn.mpg.de</footnote>
<figure><img src="image_2.png" /></figure>

# A freely precessing magnetar following an X-ray outburst

Received: 16 November 2022

Accepted: 15 February 2024

Published online: 8 April 2024

<figure><img src="image_3.png" /></figure>

Check for updates

Gregory Desvignes Ⓓ<fnref n="1" />,²✉, Patrick Weltevrede Ⓓ<fnref n="3" />, Yong Gao Ⓓ⁴,⁵,

David Ian Jones⁶, Michael Kramer Ⓓ<fnref n="1" />,<fnref n="3" />, Manisha Caleb Ⓓ⁷,⁸,

Ramesh Karuppusamy Ⓓ<fnref n="1" />, Lina Levin Ⓓ<fnref n="3" />, Kuo Liu Ⓓ<fnref n="1" />, Andrew G. Lyne Ⓓ<fnref n="3" />,

Lijing Shao Ⓓ<fnref n="1" />,⁵,⁹, Ben Stappers<fnref n="3" /> & Jérôme Pétri Ⓓ¹⁰

Magnetars—highly magnetized neutron stars—are thought to be the most likely progenitors for fast radio bursts (FRBs). Freely precessing magnetars are further invoked to explain the repeating FRBs. We report here on new high-cadence radio observations of the magnetar XTE J1810–197 recorded shortly after an X-ray outburst. We interpret the polarization variations of the magnetar radio emission as evidence for the magnetar undergoing free precession following the outburst while its magnetosphere slowly untwists. The observations of precession being damped on a timescale of months argue against the scenario of freely precessing magnetars as the origin of repeating FRBs. Using free-precession models based on relaxing ellipticity with a decay of the wobble angle, we find the magnetar ellipticity to be in good agreement with theoretical predictions from nuclear physics. Our precise measurement of the magnetar’s geometry can also further help in refining the modelling of X-ray light curves and constrain the star’s compactness.

Magnetars are rare, typically slowly spinning neutron stars whose emission is thought to be powered by the decay of their large (≥10¹² G) magnetic fields<fnref n="1" /> in opposition to normal rotation-powered pulsars. They occasionally undergo bright X-ray outburst phases, supposedly originating from a sudden quake in the star’s crust, resulting in a twisted magnetosphere that fuels the outburst². Out of the 24 confirmed magnetars currently known, only 6 have been shown to emit in the radio band<fnref n="3" />.

The discovery of a bright fast radio burst (FRB) from the galactic magnetar SGR 1935+2154 recently gave strong credence to radio magnetars being the progenitors of at least some extra galactic FRBs⁴ (see, for example, ref. 5 for a review on FRBs). Models of freely precessing

magnetars were then put forward to explain the periodicities in the activity observed in the train of pulses from some repeating FRBs⁶,⁷.

The first recorded X-ray outburst of the magnetar XTE J1810–197 (position 18h09m05.17s – 19°43′51.8″, J2000) happened around late 2002 (ref. 8) and bright radio pulsations were detected 3 years later⁹ with a periodicity of $P_s = 5.54$ s. The radio emission slowly decayed with time until it ceased to be detected in late 2008 (ref. 10). Radio pulsations from XTE J1810–197 were again detected on 8 December 2018 (Modified Julian Date (MJD) 58,460.6)¹¹, following a second X-ray outburst known to have occurred between 20 and 26 November 2018 (ref. 12).

We present a set of 62 polarimetric observations of XTE J1810–197 recorded with the Lovell and Effelsberg radio telescopes at

<footnote n="1">¹Max-Planck-Institut für Radioastronomie, Bonn, Germany. ²Laboratoire d'Études Spatiales et d'Instrumentation en Astrophysique, Observatoire de Paris, Université Paris-Sciences-et-Lettres, Centre National de la Recherche Scientifique, Sorbonne Université, Université de Paris, Meudon, France.</footnote>

<footnote n="3">³Jodrell Bank Centre for Astrophysics, The University of Manchester, Manchester, UK. ⁴Department of Astronomy, School of Physics, Peking University, Beijing, China. ⁵Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing, China. ⁶Mathematical Sciences and STAG Research Centre, University of Southampton, Southampton, UK. ⁷Sydney Institute for Astronomy, School of Physics, The University of Sydney, Sydney, New South Wales, Australia. ⁸ASTRO3D: ARC Centre of Excellence for All-sky Astrophysics in 3D, Canberra, Australian Capital Territory, Australia. ⁹National Astronomical Observatories, Chinese Academy of Sciences, Beijing, China. ¹⁰Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, Strasbourg, France. ✉e-mail: gdesvignes@mpifr-bonn.mpg.de</footnote>
1.5 and 6 GHz, respectively, between 8 December 2018 and 18 June 2020 (MJD 59,018). These observations provide calibrated polarimetric pulse profiles with almost daily cadence during the first month of monitoring. An overview of the observations and data processing can be found in the Methods section.

The polarimetric pulse profiles exhibit a high degree of linear polarization *L* at both frequencies, consistent with previous observations of this magnetar<sup>13–15</sup>. Thanks to our observing cadence, we noticed rapid and systematic changes of the position angle (PA) of *L* with time (Supplementary Figs. 1–7), including two reversals in the sign of the gradient of the PA, first between MJDs 58,464–58,466, also previously reported in ref. 15, and a second one between MJDs 58,551–58,589.

The standard model to interpret the polarization of a pulsar's profile is the rotating vector model (RVM)¹⁶. Assuming a dipolar field geometry, it describes the observed PA sweep of the polarized radio emission as it crosses our line of sight. Under this theory, the exact shape of the PA sweep is determined by the viewing geometry, that is the angle between the pulsar spin axis and the pulsar-observer line of sight $\zeta = \alpha + \beta$ where $\alpha$ is the magnetic inclination angle and $\beta$ is the angle between the magnetic axis and our line of sight at its closest approach¹⁷, hereafter referred to as the impact parameter. Although deviations from the RVM can be observed (due to, for example, propagation effects throughout the pulsar magnetosphere), the validity of the geometrical nature of the RVM has been demonstrated, at least for rotation-powered pulsars¹⁸.

An extension to the RVM for the suggested twisted magnetic fields of magnetars<sup>2,19,20</sup> predicts that the untwisting of the magnetic field can produce a vertical shift in PA, and to a lesser extent a horizontal shift in rotational phase, of the PA sweep while retaining the same PA gradient under the pulse. The off-centred dipole model for polarization<sup>21</sup> does predict some small changes in the gradient of the PA but none of the theories can account for the reversal in the sign of the PA gradient without a change in the viewing geometry. Furthermore, at the radio emission height of tens of thousands of km above the magnetar's surface, inferred from our polarization data (Methods), multipolar magnetic fields are unlikely to be dominant over the dipolar magnetic field. The sudden apparent change in the gradient of the PA of the magnetar Swift J1818.0–1607 on MJD 59,062 (ref. 22), together with a depolarization of the linear polarization, can be explained by an apparent orthogonal polarization mode<sup>23</sup>. Precession is therefore the only known physical process that could cause the observed systematic variations in the PA and its PA gradient sign reversal. In this work, we consider both forced precession and free precession.

In the case of forced precession, the direction of the neutron star (NS) spin vector changes significantly, due to, for example, its motion in a curved spacetime caused by a massive companion star. In this case, $\alpha$ remains fixed with time and $\beta$ is modulated by the relativistic spin-precession period as predicted by Einstein's theory of General Relativity<sup>18,24</sup>.

In contrast, in the case of free precession (see, for example, ref. 25), precession originates from the aspherical deformation of the NS expressed as its ellipticity $\epsilon$. The angle between this deformation axis and the angular momentum of the NS is defined as the wobble angle $\theta$. The direction of the NS spin vector is assumed to be approximately fixed in space, that is $\zeta$ remains constant with time. The magnetic axis precesses around the symmetry axis misaligned with an angle $\xi$, causing periodic modulations in $\alpha$ and $\beta$. This precession occurs at the frequency $\omega_p = \epsilon\omega \cos\theta$ where $\omega = \frac{2\pi}{P_s}$ is the NS angular frequency.

Figure 1 illustrates the free-precession geometry.

Despite initial claims of detections of free precession in the pulsars B1642–03 and B1828–11 through the timing of their pulse arrival times<sup>26,27</sup>, these pulsars have since then been shown to exhibit magnetospheric mode changes causing the observed profile variation<sup>28,29</sup>, complicating the interpretation of their timing behaviour. Free precession

<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Geometry of a freely precessing biaxial neutron star in the inertial frame.** The NS spin vector **ω** of period $P_s$ rotates around the angular momentum vector **L** with a very small angle $\theta_1$. The deformation axis **n**<sub>d</sub> is misaligned with respect to **L** by $\theta$. The NS magnetic axis **B** rotates around **n**<sub>d</sub> at an angle $\chi$ with the free-precession period $P_f \approx \frac{P_s}{\epsilon\cos\theta}$ where $\epsilon$ is the ellipticity of the NS. **L**, **ω** and **n**<sub>d</sub> always lie in the same plane. As the NS precesses, the magnetic inclination angle $\alpha$ and $\beta$, the closest approach between the NS magnetic axis and our line of sight (LOS), varies with time while $\zeta = \alpha + \beta$ is assumed to be constant.</figcaption></figure>

has also been invoked to explain the pulse phase modulation with period <1 day observed in the hard X-ray emission of the magnetars 4U 0142+61 (ref. 30), 1E 1547.0–5408 (ref. 31) and SGR 1900+14 (ref. 32).

Although XTE J1810–197 is not known to orbit around a massive companion, we nonetheless attempted to compare the cases of forced and free precession. We performed a simultaneous fit of the RVM to all PA data from our 62 epochs, with a common $\alpha$ parameter between epochs in the case of forced precession. For the free-precession scenario, we replaced the common $\alpha$ parameter with a common $\zeta$ parameter. The remaining RVM parameters ($\beta$ and two offset parameters, per epoch) are set free. Finally we applied the nested sampling tool POLYCHORD<sup>33</sup> to explore the 189-dimension parameter space and compute the Bayesian log-evidence $\log Z$ for each case. We take the 68% confidence levels on the one-dimensional marginalized posterior as our 1$\sigma$ uncertainties (more details on the analysis can be found in Methods).

The log-Bayes factor, the difference in 'log-evidence' between the two models, is -75 in favour of the free-precession model and therefore unarguably supports, as expected, free precession against forced precession as the origin of the observed PA changes. When including the effects of a decaying, eastward-twisted magnetic field (Methods), the log-Bayes factor increases to ~244 in favour of the free-precession model. This analysis gives $\zeta = 167.26^\circ \pm 0.22^\circ$ and we show in Fig. 2 the temporal variations of $\alpha$ and $\beta$ according to this model, with $160^\circ < \alpha < 175^\circ$ and $|\beta| < 8^\circ$. From Fig. 2 we can see that the rapid rate of change of $\alpha$ and $\beta$ (-0.4$^\circ$ day$^{-1}$ at its observed maximum near MJD 58,470) is decreasing with time arguing for damped precession. We can constrain the time relaxation of the twisted magnetic field $\tau_t > 730$ days (95% confidence level) with a large initial twist parameter $n_0 \approx 0.04$ at the time of the outburst (Extended Data Fig. 1). This evidence for the untwisting of the magnetosphere is consistent with the
1.5 and 6 GHz, respectively, between 8 December 2018 and 18 June 2020 (MJD 59,018). These observations provide calibrated polarimetric pulse profiles with almost daily cadence during the first month of monitoring. An overview of the observations and data processing can be found in the Methods section.

The polarimetric pulse profiles exhibit a high degree of linear polarization *L* at both frequencies, consistent with previous observations of this magnetar<sup>13–15</sup>. Thanks to our observing cadence, we noticed rapid and systematic changes of the position angle (PA) of *L* with time (Supplementary Figs. 1–7), including two reversals in the sign of the gradient of the PA, first between MJDs 58,464–58,466, also previously reported in ref. 15, and a second one between MJDs 58,551–58,589.

The standard model to interpret the polarization of a pulsar's profile is the rotating vector model (RVM)¹⁶. Assuming a dipolar field geometry, it describes the observed PA sweep of the polarized radio emission as it crosses our line of sight. Under this theory, the exact shape of the PA sweep is determined by the viewing geometry, that is the angle between the pulsar spin axis and the pulsar-observer line of sight $\zeta = \alpha + \beta$ where $\alpha$ is the magnetic inclination angle and $\beta$ is the angle between the magnetic axis and our line of sight at its closest approach¹⁷, hereafter referred to as the impact parameter. Although deviations from the RVM can be observed (due to, for example, propagation effects throughout the pulsar magnetosphere), the validity of the geometrical nature of the RVM has been demonstrated, at least for rotation-powered pulsars¹⁸.

An extension to the RVM for the suggested twisted magnetic fields of magnetars<sup>2,19,20</sup> predicts that the untwisting of the magnetic field can produce a vertical shift in PA, and to a lesser extent a horizontal shift in rotational phase, of the PA sweep while retaining the same PA gradient under the pulse. The off-centred dipole model for polarization<sup>21</sup> does predict some small changes in the gradient of the PA but none of the theories can account for the reversal in the sign of the PA gradient without a change in the viewing geometry. Furthermore, at the radio emission height of tens of thousands of km above the magnetar's surface, inferred from our polarization data (Methods), multipolar magnetic fields are unlikely to be dominant over the dipolar magnetic field. The sudden apparent change in the gradient of the PA of the magnetar Swift J1818.0–1607 on MJD 59,062 (ref. 22), together with a depolarization of the linear polarization, can be explained by an apparent orthogonal polarization mode<sup>23</sup>. Precession is therefore the only known physical process that could cause the observed systematic variations in the PA and its PA gradient sign reversal. In this work, we consider both forced precession and free precession.

In the case of forced precession, the direction of the neutron star (NS) spin vector changes significantly, due to, for example, its motion in a curved spacetime caused by a massive companion star. In this case, $\alpha$ remains fixed with time and $\beta$ is modulated by the relativistic spin-precession period as predicted by Einstein's theory of General Relativity<sup>18,24</sup>.

In contrast, in the case of free precession (see, for example, ref. 25), precession originates from the aspherical deformation of the NS expressed as its ellipticity $\epsilon$. The angle between this deformation axis and the angular momentum of the NS is defined as the wobble angle $\theta$. The direction of the NS spin vector is assumed to be approximately fixed in space, that is $\zeta$ remains constant with time. The magnetic axis precesses around the symmetry axis misaligned with an angle $\xi$, causing periodic modulations in $\alpha$ and $\beta$. This precession occurs at the frequency $\omega_p = \epsilon\omega \cos\theta$ where $\omega = \frac{2\pi}{P_s}$ is the NS angular frequency.

Figure 1 illustrates the free-precession geometry.

Despite initial claims of detections of free precession in the pulsars B1642–03 and B1828–11 through the timing of their pulse arrival times<sup>26,27</sup>, these pulsars have since then been shown to exhibit magnetospheric mode changes causing the observed profile variation<sup>28,29</sup>, complicating the interpretation of their timing behaviour. Free precession

has also been invoked to explain the pulse phase modulation with period <1 day observed in the hard X-ray emission of the magnetars 4U 0142+61 (ref. 30), 1E 1547.0–5408 (ref. 31) and SGR 1900+14 (ref. 32).

Although XTE J1810–197 is not known to orbit around a massive companion, we nonetheless attempted to compare the cases of forced and free precession. We performed a simultaneous fit of the RVM to all PA data from our 62 epochs, with a common $\alpha$ parameter between epochs in the case of forced precession. For the free-precession scenario, we replaced the common $\alpha$ parameter with a common $\zeta$ parameter. The remaining RVM parameters ($\beta$ and two offset parameters, per epoch) are set free. Finally we applied the nested sampling tool POLYCHORD<sup>33</sup> to explore the 189-dimension parameter space and compute the Bayesian log-evidence $\log Z$ for each case. We take the 68% confidence levels on the one-dimensional marginalized posterior as our $1\sigma$ uncertainties (more details on the analysis can be found in Methods).

The log-Bayes factor, the difference in 'log-evidence' between the two models, is -75 in favour of the free-precession model and therefore unarguably supports, as expected, free precession against forced precession as the origin of the observed PA changes. When including the effects of a decaying, eastward-twisted magnetic field (Methods), the log-Bayes factor increases to ~244 in favour of the free-precession model. This analysis gives $\zeta = 167.26^\circ \pm 0.22^\circ$ and we show in Fig. 2 the temporal variations of $\alpha$ and $\beta$ according to this model, with $160^\circ < \alpha < 175^\circ$ and $|\beta| < 8^\circ$. From Fig. 2 we can see that the rapid rate of change of $\alpha$ and $\beta$ (-0.4$^\circ$ day$^{-1}$ at its observed maximum near MJD 58,470) is decreasing with time arguing for damped precession. We can constrain the time relaxation of the twisted magnetic field $\tau_t > 730$ days (95% confidence level) with a large initial twist parameter $n_0 \approx 0.04$ at the time of the outburst (Extended Data Fig. 1). This evidence for the untwisting of the magnetosphere is consistent with the

<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Geometry of a freely precessing biaxial neutron star in the inertial frame.** The NS spin vector **ω** of period $P_s$ rotates around the angular momentum vector **L** with a very small angle $\theta_1$. The deformation axis **n**<sub>d</sub> is misaligned with respect to **L** by $\theta$. The NS magnetic axis **B** rotates around **n**<sub>d</sub> at an angle $\chi$ with the free-precession period $P_f \approx \frac{P_s}{\epsilon \cos \theta}$ where $\epsilon$ is the ellipticity of the NS. **L**, **ω** and **n**<sub>d</sub> always lie in the same plane. As the NS precesses, the magnetic inclination angle $\alpha$ and $\beta$, the closest approach between the NS magnetic axis and our line of sight (LOS), varies with time while $\zeta = \alpha + \beta$ is assumed to be constant.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**Fig. 2 | Temporal evolution of the geometry of XTE J1810–197.** The top and bottom panels show the magnetic inclination angle $\alpha$ and the impact parameter $\beta$, respectively. Assuming the model independent free-precession scenario (constant $\zeta$ across all epochs), the orange and blue points represent the measurements from the Jodrell 1.5 GHz and Effelsberg 6 GHz data, respectively. We take the 95% confidence levels on the one-dimensional marginalized posterior as the uncertainties for the $\alpha$ and $\beta$ measurements. The dot-dash line at $\beta = 0^\circ$ represents the crossing of the line of sight over the magnetic pole where the gradient of the PA is predicted to change sign. The period of time when the second X-ray outburst occurred is marked yellow. The dashed and dotted lines delimit the 95% confidence levels on $\alpha$ and $\beta$ for the model of crust–core coupling with relaxing ellipticity and the phenomenological model of wobble angle decay with relaxing ellipticity model, respectively.

observation of the decrease in pulse width³⁴ and the variable spin-down rate of the star and its delay after the outburst³⁵.

In light of these new results, we reanalysed four epochs (Extended Data Figs. 2 and 3 and Methods) recorded with the Effelsberg radio telescope during the summer of 2006 (ref. 13), more than 3 years after the first recorded X-ray outburst, and applied the RVM to each of these epochs to revisit the viewing geometry of the magnetar in 2006. Taking the averaged $\zeta$ value, we find $\zeta = 155^\circ \pm 10^\circ$, indicating as expected a stable viewing angle between 2006 and 2018. These data could however not constrain the twisting of the magnetosphere. Assuming that free precession also occurred following the 2002 outburst of the magnetar (as we suggest for the 2018 outburst), it would probably have been damped by the time of the first radio detection in 2006, suggesting that $\alpha$ and $\beta$ should have remained constant. Indeed, no variations in the shape of the PA were detected between 2006 and 2008 (ref. 10), confirming the stable geometry.

We also considered different models of biaxial free precession of a NS to interpret the observed PA variations, including relaxation of $\epsilon$ (to either a zero or non-zero value) and relaxation of $\epsilon$ combined with either a phenomenological or a frictional crust–core coupling model for the decay of $\theta$. All models included the effects of a decaying and twisted magnetic field. As the geometry of XTE J1810–197 remained stable during 2006–2008, we assume here that free precession was set in motion at the time $T_0$ of the 2018 X-ray outburst. We again used POLYCHORD to perform the model comparison (see Methods for more details on these models and our analysis). Our results (Extended Data Table 1) show with a Bayes factor of 45 that the preferred model (Table 1) is based on relaxing ellipticity combined with a phenomenological description for the decay of $\theta$. The second preferred model is the crust–core coupling model with relaxation of $\epsilon$. The predicted ellipticity $\epsilon$ (Extended Data Fig. 4) drops from a value of ~$(2.4 \pm 0.1) \times 10^{-6}$ and ~$(-1.7 \pm 0.03) \times 10^{-6}$ at $t = T_0$ to a constant value of $(1.25 \pm 0.07) \times 10^{-7}$ and $(1.0 \pm 0.02) \times 10^{-7}$ about 5 months later for the preferred phenomenological and crust–core coupling models, respectively. The waveforms

of the predicted geometry for these two models are shown in Fig. 2. Both predict a large positive $\beta$ at the time of the X-ray outburst and a relatively stable negative $\beta$ from 2020 onward. The discrepancy in $\beta$ between 2006 and 2008 and the value predicted by our model at the epoch $T_0$ could be attributed to either chaotic events at the time of the outburst that cannot be reproduced by our model or by the decay of ellipticity deviating from a simple exponential form.

The upper limit for the ellipticity of a NS derived from crustal elasticity³⁶⁻³⁸ is on the order of $\epsilon_{ela} \approx 10^{-6} \left(\frac{\theta_{br}}{0.1}\right)$. The parameter $\sigma_{br}$ is the braking strain of the elastic crust. Molecular dynamics simulations for high pressure Coulomb crystals suggest that $\sigma_{br} \approx 0.1$ (ref. 39). Strong internal magnetic fields B in the NS also create deformation with an ellipticity on the order of $\epsilon_{mag} \approx 10^{-6} \left(\frac{B}{10^{15} \text{ G}}\right)^2$, which is a crude estimation but consistent with more rigorous calculations (for example, refs. 40–43). The maximum value of $\epsilon = 1.7 \times 10^{-6}$ to ~$2.4 \times 10^{-6}$ at $t = T_0$ depending on the model considered is consistent with the maximum ellipticity $\epsilon_{ela} \approx 10^{-6}$ that a NS could sustain. Assuming the surface magnetic field strength $B_{surf} \approx 2.7 \times 10^{14}$ G of XTE J1810–197 (ref. 44) is representative of the internal magnetic field, then $\epsilon_{mag} \approx 7 \times 10^{-8}$ is also broadly consistent with the constant $\epsilon \approx 1 \times 10^{-7}$ derived from the two preferred models. A possible interpretation is that the constant contribution to $\epsilon$ is sourced from $\epsilon_{mag}$ while the crustal strain causing the decaying contribution to $\epsilon$ relaxed in the first 3 months after the outburst.

XTE J1810–197 is not the only magnetar observed in radio prior and after an X-ray outburst. PSR J1622–4950 was first discovered by the High Time Resolution Universe South survey in 2009 (ref. 45) and monitored in radio until its emission ceased to be detected in 2014. Radio emission from this source was again detected on 5 April 2017, with an X-ray outburst happening two weeks before the radio redetection. However, the first polarimetric pulse profiles used for that study were recorded about 5 months after the outburst. A comparison of the polarimetric pulse profiles pre and post outburst did not reveal any large changes in the magnetar’s viewing geometry⁴⁶. We can speculate that either the outburst of PSR J1622–4950 did not excite the free precession or the precession damped on a timescale similar to or shorter than that of XTE J1810–197 (for example, due to a smaller initial wobble angle).

The precession of XTE J1810–197 can potentially lead to changes in the observed polarimetric pulse profiles as our line of sight cuts through different parts of the beam. Additionally, free precession is expected to cause variations in the spin-down rate of XTE J1810–197, as the spin-down torque depends upon the angle $\alpha$, varying approximately as $\sin^2\alpha$ (ref. 47). No large variations are observed in the shape of the pulse profile in the first week of observations, except for some intensity fluctuations, despite the large change in $\beta$, including a change of sign. However, while we only observe large amplitude variations in the left- and right-handed circular polarization mixed with propagation effects as the PA gradient first changes sign, we observe a reversal of the sign of the circular polarization at the time of the second crossing of the magnetic pole (Extended Data Fig. 5), as expected⁴⁸ and observed in the case of the relativistic binary pulsar J1906+0746 (ref. 18). With our measured range of change of $\alpha$ (Fig. 2), we estimate the spin-down rate (Fig. 2 and equation 9 of ref. 47) to vary by about 700% to 30% for a plasma conductivity in the pulsar magnetosphere ranging from zero to infinity, respectively. However, we find no correlation between the spin-down rate presented in Fig. 3 of ref. 35 and $\alpha$. This indicates that other physical processes such as the untwisting of magnetic field lines², high-energy emission, the magnetar’s wind, magnetospheric flares or seismic activities⁴⁹ must contribute more to variations in the spin-down rate of this magnetar than changes in the magnetic inclination angle.

The geometry of XTE J1810–197 has been independently constrained by the modelling of its X-ray pulse profile⁵⁰⁻⁵². Assuming the X-ray spot axis is aligned with the magnetic axis (as valid for a dipolar
<figure><img src="image_5.png" /></figure>

**Fig. 2 | Temporal evolution of the geometry of XTE J1810–197.** The top and bottom panels show the magnetic inclination angle $\alpha$ and the impact parameter $\beta$, respectively. Assuming the model independent free-precession scenario (constant $\zeta$ across all epochs), the orange and blue points represent the measurements from the Jodrell 1.5 GHz and Effelsberg 6 GHz data, respectively. We take the 95% confidence levels on the one-dimensional marginalized posterior as the uncertainties for the $\alpha$ and $\beta$ measurements. The dot-dash line at $\beta = 0^\circ$ represents the crossing of the line of sight over the magnetic pole where the gradient of the PA is predicted to change sign. The period of time when the second X-ray outburst occurred is marked yellow. The dashed and dotted lines delimit the 95% confidence levels on $\alpha$ and $\beta$ for the model of crust–core coupling with relaxing ellipticity and the phenomenological model of wobble angle decay with relaxing ellipticity model, respectively.

observation of the decrease in pulse width³⁴ and the variable spin-down rate of the star and its delay after the outburst³⁵.

In light of these new results, we reanalysed four epochs (Extended Data Figs. 2 and 3 and Methods) recorded with the Effelsberg radio telescope during the summer of 2006 (ref. 13), more than 3 years after the first recorded X-ray outburst, and applied the RVM to each of these epochs to revisit the viewing geometry of the magnetar in 2006. Taking the averaged $\zeta$ value, we find $\zeta = 155^\circ \pm 10^\circ$, indicating as expected a stable viewing angle between 2006 and 2018. These data could however not constrain the twisting of the magnetosphere. Assuming that free precession also occurred following the 2002 outburst of the magnetar (as we suggest for the 2018 outburst), it would probably have been damped by the time of the first radio detection in 2006, suggesting that $\alpha$ and $\beta$ should have remained constant. Indeed, no variations in the shape of the PA were detected between 2006 and 2008 (ref. 10), confirming the stable geometry.

We also considered different models of biaxial free precession of a NS to interpret the observed PA variations, including relaxation of $\epsilon$ (to either a zero or non-zero value) and relaxation of $\epsilon$ combined with either a phenomenological or a frictional crust–core coupling model for the decay of $\theta$. All models included the effects of a decaying and twisted magnetic field. As the geometry of XTE J1810–197 remained stable during 2006–2008, we assume here that free precession was set in motion at the time $T_0$ of the 2018 X-ray outburst. We again used POLYCHORD to perform the model comparison (see Methods for more details on these models and our analysis). Our results (Extended Data Table 1) show with a Bayes factor of 45 that the preferred model (Table 1) is based on relaxing ellipticity combined with a phenomenological description for the decay of $\theta$. The second preferred model is the crust–core coupling model with relaxation of $\epsilon$. The predicted ellipticity $\epsilon$ (Extended Data Fig. 4) drops from a value of $\sim (2.4 \pm 0.1) \times 10^{-6}$ and $-(1.7 \pm 0.03) \times 10^{-6}$ at $t = T_0$ to a constant value of $(1.25 \pm 0.07) \times 10^{-7}$ and $(1.0 \pm 0.02) \times 10^{-7}$ about 5 months later for the preferred phenomenological and crust–core coupling models, respectively. The waveforms

of the predicted geometry for these two models are shown in Fig. 2. Both predict a large positive $\beta$ at the time of the X-ray outburst and a relatively stable negative $\beta$ from 2020 onward. The discrepancy in $\beta$ between 2006 and 2008 and the value predicted by our model at the epoch $T_0$ could be attributed to either chaotic events at the time of the outburst that cannot be reproduced by our model or by the decay of ellipticity deviating from a simple exponential form.

The upper limit for the ellipticity of a NS derived from crustal elasticity³⁶⁻³⁸ is on the order of $\epsilon_{ela} \approx 10^{-6} \left(\frac{\theta_{br}}{0.1}\right)$. The parameter $\sigma_{br}$ is the braking strain of the elastic crust. Molecular dynamics simulations for high pressure Coulomb crystals suggest that $\sigma_{br} \approx 0.1$ (ref. 39). Strong internal magnetic fields $B$ in the NS also create deformation with an ellipticity on the order of $\epsilon_{mag} \approx 10^{-6} \left(\frac{B}{10^{15} G}\right)^2$, which is a crude estimation but consistent with more rigorous calculations (for example, refs. 40–43). The maximum value of $\epsilon = 1.7 \times 10^{-6}$ to $-2.4 \times 10^{-6}$ at $t = T_0$ depending on the model considered is consistent with the maximum ellipticity $\epsilon_{ela} \approx 10^{-6}$ that a NS could sustain. Assuming the surface magnetic field strength $B_{surf} \approx 2.7 \times 10^{14}$ G of XTE J1810–197 (ref. 44) is representative of the internal magnetic field, then $\epsilon_{mag} \approx 7 \times 10^{-8}$ is also broadly consistent with the constant $\epsilon \approx 1 \times 10^{-7}$ derived from the two preferred models. A possible interpretation is that the constant contribution to $\epsilon$ is sourced from $\epsilon_{mag}$ while the crustal strain causing the decaying contribution to $\epsilon$ relaxed in the first 3 months after the outburst.

XTE J1810–197 is not the only magnetar observed in radio prior and after an X-ray outburst. PSR J1622–4950 was first discovered by the High Time Resolution Universe South survey in 2009 (ref. 45) and monitored in radio until its emission ceased to be detected in 2014. Radio emission from this source was again detected on 5 April 2017, with an X-ray outburst happening two weeks before the radio redetection. However, the first polarimetric pulse profiles used for that study were recorded about 5 months after the outburst. A comparison of the polarimetric pulse profiles pre and post outburst did not reveal any large changes in the magnetar’s viewing geometry⁴⁶. We can speculate that either the outburst of PSR J1622–4950 did not excite the free precession or the precession damped on a timescale similar to or shorter than that of XTE J1810–197 (for example, due to a smaller initial wobble angle).

The precession of XTE J1810–197 can potentially lead to changes in the observed polarimetric pulse profiles as our line of sight cuts through different parts of the beam. Additionally, free precession is expected to cause variations in the spin-down rate of XTE J1810–197, as the spin-down torque depends upon the angle $\alpha$, varying approximately as $\sin^2\alpha$ (ref. 47). No large variations are observed in the shape of the pulse profile in the first week of observations, except for some intensity fluctuations, despite the large change in $\beta$, including a change of sign. However, while we only observe large amplitude variations in the left- and right-handed circular polarization mixed with propagation effects as the PA gradient first changes sign, we observe a reversal of the sign of the circular polarization at the time of the second crossing of the magnetic pole (Extended Data Fig. 5), as expected⁴⁸ and observed in the case of the relativistic binary pulsar J1906+0746 (ref. 18). With our measured range of change of $\alpha$ (Fig. 2), we estimate the spin-down rate (Fig. 2 and equation 9 of ref. 47) to vary by about 700% to 30% for a plasma conductivity in the pulsar magnetosphere ranging from zero to infinity, respectively. However, we find no correlation between the spin-down rate presented in Fig. 3 of ref. 35 and $\alpha$. This indicates that other physical processes such as the untwisting of magnetic field lines², high-energy emission, the magnetar’s wind, magnetospheric flares or seismic activities⁴⁹ must contribute more to variations in the spin-down rate of this magnetar than changes in the magnetic inclination angle.

The geometry of XTE J1810–197 has been independently constrained by the modelling of its X-ray pulse profile⁵⁰⁻⁵². Assuming the X-ray spot axis is aligned with the magnetic axis (as valid for a dipolar
<table><caption>**Table 1 | Parameter values for the two preferred models of free precession**</caption><thead><tr><th>Parameter</th><th>Value of model A</th><th>Value of model B</th></tr></thead><tbody><tr><td>Viewing angle, ζ (deg)</td><td>169.07±0.22</td><td>168.50<sup>+0.26</sup><sub>-0.21</sub></td></tr><tr><td>Initial wobble angle, θ<sub>0</sub> (deg)</td><td>19.52±0.38</td><td>30.52<sup>+0.66</sup><sub>-0.74</sub></td></tr><tr><td>Angle between the magnetic and symmetry axis, χ (deg)</td><td>173.34±0.13</td><td>171.37<sup>+0.20</sup><sub>-0.16</sub></td></tr><tr><td>Initial phase of the precession, Φ<sub>0</sub> (deg)</td><td>45<sup>+5</sup><sub>-4</sub></td><td>108±3</td></tr><tr><td>Constant ellipticity of the NS, ε<sub>0</sub></td><td>(1.24±0.03)×10<sup>-7</sup></td><td>(9.17±0.14)×10<sup>-8</sup></td></tr><tr><td>Initial ellipticity of the NS, ε<sub>1</sub></td><td>(2.37±0.05)×10<sup>-6</sup></td><td>(1.58±0.03)×10<sup>-6</sup></td></tr><tr><td>Ellipticity relaxation timescale, τ<sub>ε</sub> (days)</td><td>19.55±0.35</td><td>36.43±0.46</td></tr><tr><td>Wobble angle decay timescale, τ<sub>δ</sub> (days)</td><td>74.30±0.26</td><td>—</td></tr><tr><td>Frictional coupling timescale, τ<sub>c</sub> (s)</td><td>—</td><td>2.49±0.01</td></tr><tr><td>Ratio between the moment of inertia of the crust and the core, κ</td><td>—</td><td>&lt;0.01</td></tr><tr><td>Start time of the precession, T<sub>0</sub> (MJD)</td><td>58,444.5±0.4</td><td>58,444.9±0.5</td></tr><tr><td>Initial twist parameter, n<sub>0</sub></td><td>0.040<sup>+0.008</sup><sub>-0.040</sub></td><td>0.062<sup>+0.015</sup><sub>-0.062</sub></td></tr><tr><td>Twisted magnetic field relaxation timescale, τ<sub>t</sub> (days)</td><td>&gt;1,500</td><td>&gt;1,500</td></tr></tbody></table>

Model A denotes the phenomenological model of decaying wobble angle with relaxing ellipticity and model B denotes the model of frictional crust-core coupling with relaxing ellipticity. Both models include the effects of an decaying, eastward-twisted magnetic field.

magnetic field), we can compare the radio and X-ray derived geometry. We find our measurement of α is consistent with the results from the different analyses of X-ray data (if we take into consideration we use the historic convention of the PA increasing clockwise on the sky⁵³ or if we refer to the opposite magnetic pole, that is α = 180° - α). The X-ray analyses, however, provide different and inconsistent estimates for ζ originating from the different assumptions behind their complex modelling. Our results rule out the modelling of the 2003 outburst by ref. 50 and strongly favour the analysis of ref. 51 based on, for example, thermal quiescent emission, a dipolar model for the magnetar surface temperature (instead of the uniform temperature assumed in ref. 50) and inclusion of light deflection from General Relativity. Combined with the Very Long Baseline Array estimate of the magnetar's distance of 2.5⁰±⁰.⁴ kpc (ref. 54), our accurate measurement of ζ could eventually constrain the compactness (mass M over radius R) of the magnetar. Based on the modelling by ref. 51 assuming M = 1.4 M<sub>☉</sub> our constraint on ζ indicates R ≳ 15 km.

Neutron stars with asymmetric deformation along an axis different from the rotation axis are expected to emit gravitational waves (GW) at both once and twice their rotation frequency⁵⁵. Numerous attempts have been made at detecting the GW emitted by normal and fast-spinning pulsars using ground-based interferometers such as the Laser Interferometer Gravitational-wave Observatory and Virgo. For slowly spinning NSs like magnetars and XTE J1810–197, the frequencies of the GW correspond to 0.18 and 0.36 Hz, falling in the frequency range (10⁻⁴ Hz to 1 Hz) of the Laser Interferometer Space Antenna gravitational wave observatory. However, the characteristic strain amplitude of the GW at the Earth is a function of P<sub>s</sub>⁻² (see, for example, ref. 56) and several orders of magnitude below the expected sensitivity of the Laser Interferometer Space Antenna gravitational wave observatory.

Freely precessing magnetars were recently proposed⁶,⁷ to explain the periodicity observed in the activity of some repeating FRBs⁵⁷,⁵⁸. We have shown here that free precession of XTE J1810–197, probably excited at the time of its 2018 X-ray outburst, has been damped within a timescale of a few months, with a corresponding decay in both the wobble angle and the ellipticity. These results argue against free precession

of magnetars as the mechanism behind repeating FRBs with activity periods of the order of months, at least for magnetars with B<sub>surf</sub> similar to XTE J1810–197.

Using our flux-calibrated single-pulse observations of XTE J1810–197 from Effelsberg, we investigated if the brightest single pulses (SPs) detected with peak flux density >200 Jy (still several orders of magnitude below the typical FRB luminosity⁵⁵) were occurring at some specific impact parameter β. However, we find no link between β and the brightest SPs (Extended Data Fig. 6).

In summary, high-cadence radio and X-ray observations of magnetars, especially shortly after the detection of an outburst, are key to understanding the physics of free precession and could provide important information for testing the internal and magnetic structure of NSs.

## Methods

### Observations and data reduction

Following the detection of the radio revival of XTE J1810–197 on 8 December 2018, the Lovell Telescope at the Jodrell Bank Observatory (JBO) was used for a regular monitoring of the magnetar, at almost daily cadence during the first weeks of observations¹¹. The JBO data were acquired with the ROACH pulsar backend⁵⁹ tuned to a central frequency of 1,532 MHz (L-band). A 384 MHz bandwidth split over 768 channels was calibrated using a matrix template matching technique⁶⁰ and observations of the pulsars B0540+23, B0611+22 and B1737–30. A total of 36 observations, folded with the ephemeris from ref. 11, were recorded between 8 December 2018 and 23 January 2019, before the telescope went down for maintenance. One additional observation with the Lovell Telescope was made on 9 March 2019, during an extended period of maintenance. Most JBO observations are between 30 to 60 min long.

At the Effelsberg observatory, we observed XTE J1810–197 with the S60 receiver of the 100 m telescope tuned to a frequency of 4.85 GHz (C-band) three days after the detection of its radio revival. Following this successful detection, we then started monitoring XTE J1810–197 with the S110 and S45 receivers tuned to frequencies of 2.55 GHz (S-band) and 6 GHz receivers (C/X-band), respectively, to provide complementary frequency coverage to the JBO monitoring campaign. Since the S-band observations were less frequent and had a pulse profile with a lower signal-to-noise ratio (S/N) than the C/X-band observations, and were also recorded within 30 min of the C/X-band data, we decided to not include the S-band data in this study. Nonetheless we used the S-band data to verify the polarization calibration of the C/X-band data.

The circular-basis S110 and S60 receivers sent signals to the PSRIX backend⁶¹ and provided, respectively, 400 MHz and 500 MHz bandwidths split into 128 frequency channels. The PSRIX backend was configured to provide averaged pulse profiles folded over 10 s subintegrations across 1,024 phase bins using the rotational period determined from the most recently published timing ephemeris¹¹,³⁵. The data taken with the S45 linear-feed receiver were recorded with the PSRIX2 backend⁶² in PSRFITS search mode with 131 µs time resolution and 4 GHz of bandwidth centred at 6 GHz and split into 4,096 frequency channels. The data were folded offline across 2,048 phase bins using the latest ephemeris from ref. 35 to form a time-averaged pulse profile. The data were also folded across 16,384 phase bins to produce single-pulse archives. Twenty-five observations were recorded with the S45 receiver and the PSRIX2 backend, between 12 December 2018 and 18 June 2020.

All Effelsberg observations are usually between 10 to 20 min long. A 2 min scan of a polarized pulse noise diode recorded before each observation is used to calibrate the data in polarization with PSRCHIVE⁶³. We routinely observed the planetary nebula NGC 7027 to flux-calibrate the Effelsberg observations.

All JBO and Effelsberg data were then corrected for Faraday rotation using the rotation measure of 74.44 rad m⁻² value from ref. 15 and provide calibrated Stokes parameters I, Q, U, V⁶⁴.
<table><caption>**Table 1 | Parameter values for the two preferred models of free precession**</caption><thead><tr><th>Parameter</th><th>Value of model A</th><th>Value of model B</th></tr></thead><tbody><tr><td>Viewing angle, ζ (deg)</td><td>169.07±0.22</td><td>168.50<sup>+0.26</sup><sub>-0.21</sub></td></tr><tr><td>Initial wobble angle, θ<sub>0</sub> (deg)</td><td>19.52±0.38</td><td>30.52<sup>+0.66</sup><sub>-0.74</sub></td></tr><tr><td>Angle between the magnetic and symmetry axis, χ (deg)</td><td>173.34±0.13</td><td>171.37<sup>+0.20</sup><sub>-0.16</sub></td></tr><tr><td>Initial phase of the precession, Φ<sub>0</sub> (deg)</td><td>45<sup>+5</sup><sub>-4</sub></td><td>108±3</td></tr><tr><td>Constant ellipticity of the NS, ε<sub>0</sub></td><td>(1.24±0.03)×10<sup>-7</sup></td><td>(9.17±0.14)×10<sup>-8</sup></td></tr><tr><td>Initial ellipticity of the NS, ε<sub>1</sub></td><td>(2.37±0.05)×10<sup>-6</sup></td><td>(1.58±0.03)×10<sup>-6</sup></td></tr><tr><td>Ellipticity relaxation timescale, τ<sub>ε</sub> (days)</td><td>19.55±0.35</td><td>36.43±0.46</td></tr><tr><td>Wobble angle decay timescale, τ<sub>δ</sub> (days)</td><td>74.30±0.26</td><td>—</td></tr><tr><td>Frictional coupling timescale, τ<sub>c</sub> (s)</td><td>—</td><td>2.49±0.01</td></tr><tr><td>Ratio between the moment of inertia of the crust and the core, κ</td><td>—</td><td>&lt;0.01</td></tr><tr><td>Start time of the precession, T<sub>0</sub> (MJD)</td><td>58,444.5±0.4</td><td>58,444.9±0.5</td></tr><tr><td>Initial twist parameter, n<sub>0</sub></td><td>0.040<sup>+0.008</sup><sub>-0.040</sub></td><td>0.062<sup>+0.015</sup><sub>-0.062</sub></td></tr><tr><td>Twisted magnetic field relaxation timescale, τ<sub>t</sub> (days)</td><td>&gt;1,500</td><td>&gt;1,500</td></tr></tbody></table>

Model A denotes the phenomenological model of decaying wobble angle with relaxing ellipticity and model B denotes the model of frictional crust-core coupling with relaxing ellipticity. Both models include the effects of an decaying, eastward-twisted magnetic field.

magnetic field), we can compare the radio and X-ray derived geometry. We find our measurement of $\alpha$ is consistent with the results from the different analyses of X-ray data (if we take into consideration we use the historic convention of the PA increasing clockwise on the sky$^{53}$ or if we refer to the opposite magnetic pole, that is $\alpha = 180^\circ - \alpha$). The X-ray analyses, however, provide different and inconsistent estimates for $\zeta$ originating from the different assumptions behind their complex modelling. Our results rule out the modelling of the 2003 outburst by ref. 50 and strongly favour the analysis of ref. 51 based on, for example, thermal quiescent emission, a dipolar model for the magnetar surface temperature (instead of the uniform temperature assumed in ref. 50) and inclusion of light deflection from General Relativity. Combined with the Very Long Baseline Array estimate of the magnetar's distance of $2.5^{+0.4}_{-0.3}$ kpc (ref. 54), our accurate measurement of $\zeta$ could eventually constrain the compactness (mass $M$ over radius $R$) of the magnetar. Based on the modelling by ref. 51 assuming $M = 1.4 M_\odot$ our constraint on $\zeta$ indicates $R \gtrsim 15$ km.

Neutron stars with asymmetric deformation along an axis different from the rotation axis are expected to emit gravitational waves (GW) at both once and twice their rotation frequency$^{55}$. Numerous attempts have been made at detecting the GW emitted by normal and fast-spinning pulsars using ground-based interferometers such as the Laser Interferometer Gravitational-wave Observatory and Virgo. For slowly spinning NSs like magnetars and XTE J1810–197, the frequencies of the GW correspond to 0.18 and 0.36 Hz, falling in the frequency range ($10^{-4}$ Hz to 1 Hz) of the Laser Interferometer Space Antenna gravitational wave observatory. However, the characteristic strain amplitude of the GW at the Earth is a function of $P_s^{-2}$ (see, for example, ref. 56) and several orders of magnitude below the expected sensitivity of the Laser Interferometer Space Antenna gravitational wave observatory.

Freely precessing magnetars were recently proposed$^{6,7}$ to explain the periodicity observed in the activity of some repeating FRBs$^{57,58}$. We have shown here that free precession of XTE J1810–197, probably excited at the time of its 2018 X-ray outburst, has been damped within a timescale of a few months, with a corresponding decay in both the wobble angle and the ellipticity. These results argue against free precession

of magnetars as the mechanism behind repeating FRBs with activity periods of the order of months, at least for magnetars with $B_{surf}$ similar to XTE J1810–197.

Using our flux-calibrated single-pulse observations of XTE J1810–197 from Effelsberg, we investigated if the brightest single pulses (SPs) detected with peak flux density >200 Jy (still several orders of magnitude below the typical FRB luminosity$^5$) were occurring at some specific impact parameter $\beta$. However, we find no link between $\beta$ and the brightest SPs (Extended Data Fig. 6).

In summary, high-cadence radio and X-ray observations of magnetars, especially shortly after the detection of an outburst, are key to understanding the physics of free precession and could provide important information for testing the internal and magnetic structure of NSs.

## Methods

### Observations and data reduction

Following the detection of the radio revival of XTE J1810–197 on 8 December 2018, the Lovell Telescope at the Jodrell Bank Observatory (JBO) was used for a regular monitoring of the magnetar, at almost daily cadence during the first weeks of observations$^{11}$. The JBO data were acquired with the ROACH pulsar backend$^{59}$ tuned to a central frequency of 1,532 MHz (L-band). A 384 MHz bandwidth split over 768 channels was calibrated using a matrix template matching technique$^{60}$ and observations of the pulsars B0540+23, B0611+22 and B1737–30. A total of 36 observations, folded with the ephemeris from ref. 11, were recorded between 8 December 2018 and 23 January 2019, before the telescope went down for maintenance. One additional observation with the Lovell Telescope was made on 9 March 2019, during an extended period of maintenance. Most JBO observations are between 30 to 60 min long.

At the Effelsberg observatory, we observed XTE J1810–197 with the S60 receiver of the 100 m telescope tuned to a frequency of 4.85 GHz (C-band) three days after the detection of its radio revival. Following this successful detection, we then started monitoring XTE J1810–197 with the S110 and S45 receivers tuned to frequencies of 2.55 GHz (S-band) and 6 GHz receivers (C/X-band), respectively, to provide complementary frequency coverage to the JBO monitoring campaign. Since the S-band observations were less frequent and had a pulse profile with a lower signal-to-noise ratio (S/N) than the C/X-band observations, and were also recorded within 30 min of the C/X-band data, we decided to not include the S-band data in this study. Nonetheless we used the S-band data to verify the polarization calibration of the C/X-band data.

The circular-basis S110 and S60 receivers sent signals to the PSRIX backend$^{61}$ and provided, respectively, 400 MHz and 500 MHz bandwidths split into 128 frequency channels. The PSRIX backend was configured to provide averaged pulse profiles folded over 10 s subintegrations across 1,024 phase bins using the rotational period determined from the most recently published timing ephemeris$^{11,35}$. The data taken with the S45 linear-feed receiver were recorded with the PSRIX2 backend$^{62}$ in PSRFITS search mode with 131 µs time resolution and 4 GHz of bandwidth centred at 6 GHz and split into 4,096 frequency channels. The data were folded offline across 2,048 phase bins using the latest ephemeris from ref. 35 to form a time-averaged pulse profile. The data were also folded across 16,384 phase bins to produce single-pulse archives. Twenty-five observations were recorded with the S45 receiver and the PSRIX2 backend, between 12 December 2018 and 18 June 2020.

All Effelsberg observations are usually between 10 to 20 min long. A 2 min scan of a polarized pulse noise diode recorded before each observation is used to calibrate the data in polarization with PSRCHIVE$^{63}$. We routinely observed the planetary nebula NGC 7027 to flux-calibrate the Effelsberg observations.

All JBO and Effelsberg data were then corrected for Faraday rotation using the rotation measure of $74.44 \text{ rad m}^{-2}$ value from ref. 15 and provide calibrated Stokes parameters $I, Q, U, V^{64}$.
## Polarization of the pulse profiles

It has been previously reported¹⁵ that the PA of the linear polarization $L = \sqrt{Q^2 + U^2}$ of XTE J1810-197 pulse profiles recorded since the 2018 outburst shows complex structure and deviates from the simple 'S'-shaped curve predicted by the RVM¹⁶. The RVM is the standard geometrical model used by pulsar astronomers to interpret the PA and infer the viewing geometry of a pulsar, assuming a dipolar magnetic field.

Although the emission of the known radio magnetars has been shown to be usually highly linearly polarized, see, for example, refs. 13,45,46, we noticed some dips in the otherwise nearly 100% linearly polarized pulse profiles from multiple epochs at various phase ranges. This depolarization even occurs at different phase ranges for quasi-simultaneous observations at different frequencies (see, for example, the pulse profiles for MJD 58,474 in Supplementary Fig. 1). It is well known⁶⁵ that the superposition of different polarization mode can lead to depolarization and bias in the PA of the average pulse profile. To mitigate this, we excluded the phase ranges where the fractional linear polarization of the pulse drops below 80%. A two-dimensional histogram of the observed PA values of the linear polarization of the single-pulse profiles (Supplementary Fig. 8) shows that the depolarization in $L$ originates from the superposition of two different polarization modes and a PA branching, leading to some bias in the PA values after summation of the pulses. PA branching has been previously observed in the magnetar Swift J1818.0-1607 (ref. 22) and in a normal rotation-powered pulsar⁶⁶. This frequency-dependent effect, attributed to birefringence in the star's magnetosphere, was observed in the central part of the main component of the pulse starting around MJD 58,466 until it significantly decreased near MJD 58,500 (ref. 67). For the subsequent modelling of the PA data, these phase ranges were therefore excluded from the analysis.

## Comparing free and forced precession

The RVM, assuming the International Astronomical Union definition of the PA increasing counter-clockwise on the sky, can be written as:

$$ \tan(\psi_0 - \psi) = \frac{\sin \alpha \sin(\phi - \phi_0)}{\sin \zeta \cos \alpha - \cos \zeta \sin \alpha \cos(\phi - \phi_0)}, \quad (1) $$

with the viewing angle $\zeta = \alpha + \beta$ (ref. 17). Here $\alpha$ is the magnetic inclination angle and $\beta$ is the angle between the magnetic axis and our line of sight at the closest approach. The measured PA at any given pulse longitude $\phi$ is given by $\psi = \frac{1}{2} \arctan\left(\frac{U}{Q}\right)$. The parameter $\phi_0$ is the pulse longitude where the PA sweep is the steepest, that is, when $\psi = \psi_0$ (the PA of the rotation axis projected on the plane of the sky).

In contrast to the data recorded in 2006-2008 where an inter-pulse preceded the main pulse by about 100° (refs. 13,14), RVM fits to the PA of the individual epochs recorded after 2018 can only poorly constrain $\alpha$ and $\beta$. Therefore, we performed simultaneous RVM fits to all epochs to model the temporal PA variations as either free or forced precession. For the forced precession scenario, we set a common $\alpha$ parameter across all epochs whereas for the free-precession scenario, $\zeta$ is the common parameter across all epochs. For each epoch, the three remaining RVM parameters $\beta$, $\phi_0$ and $\psi_0$ are modelled independently.

We developed a set of codes based on model RVM⁶⁸ to simultaneously apply the RVM to the Stokes Q and U data of our dataset (instead of the traditional PA values) and explore the parameter space of the free and forced precession models with the nested sampling software PolyChord v.1.20.1 (ref. 33).

The analysis includes a scaling factor parameter per frequency band referred to as EFAC that acts as a multiplier factor to the Stokes Q and U off-pulse standard deviations. The dimensionality of the problem $N_{\text{dim}}$ is therefore $1 + 3 \times N_{\text{epochs}} + N_{\text{EFAC}}$, with $N_{\text{epochs}}$ and $N_{\text{EFAC}}$ being the number of epochs and the number of EFAC parameters included in the

analysis, respectively. All RVM parameters are sampled using $10 \times N_{\text{dim}}$ live points (with a minimum of 1,000 live points) from Gaussian priors with standard deviation of 5° (3° for $\beta$) and mean values derived from previous runs. We set the PolyChord $n_{\text{repeats}}$ parameter (that is the slice sampling chain length used to produce a new live point) to $5 \times N_{\text{dim}}$ to ensure we get reliable log-evidence ($\log z$) estimates from the nested sampling analysis. Each EFAC parameter is sampled from log-uniform prior in the $\log_{10}$ range [-0.4, 0.6]. The log-Bayes factor comes directly from the subtraction of $\log z$ between any two models. A log-Bayes factor >3 brings conclusive support for the model with the highest $\log z$⁶⁹.

In the following analyses, we selected the PA with phase bins that satisfy $L > 3\sigma_N$, where $\sigma_N$ is the off-pulse standard deviation, outside of the excluded phase ranges as described in the previous section.

We started by analysing the datasets from the two different radio bands independently (L and C/X bands) and then combined all the PA data from the two bands. The results are reported in Supplementary Table 1.

The log-Bayes factor are 54 for the JBO L-band and 9 for the Effelsberg C/X-band data in support of the free-precession interpretation of the PA. They thus provide very strong support the free-precession interpretation of the data. Unsurprisingly, the combined dataset also shows unarguable support for the free-precession scenario with a log-Bayes factor of 75.

As the free-precession model is the preferred model to explain the observed PA variations, we show in Supplementary Fig. 9 the posterior distributions for the parameter $\zeta$ from each analysis. The $\zeta$ posteriors from the two frequency bands and the combined dataset are consistent with each other. The posterior from the combined dataset gives a narrower distribution, with $\zeta = 161.0^\circ \pm 0.7^\circ$, arguing for XTE J1810-197 as being a near-aligned rotator. We show in Supplementary Fig. 10 the posterior distributions of the two EFAC parameters included in the analysis of the combined dataset. The mean values of the EFAC parameters are 1.12 and 1.44 for the L and C/X bands, respectively, indicating a good modelling of the PA.

## Twisted magnetic fields

It is assumed²⁵ that the magnetic field of a magnetar can become twisted following some abrupt crustal motion of the star. The following slow and gradual untwisting of the magnetic field lines is presumed to power the observed outburst. The perturbation of the twisted field lines on the PA $\psi$ can be written as²⁰:

$$ \Delta\psi_{\text{twist}} = -\frac{8}{9}\lambda\sin^2\theta_{\text{obs}} \quad (2) $$

where $\theta_{\text{obs}}$ is the colatitude of the line of sight defined by

$$ \cos \theta_{\text{obs}} = \cos \alpha \cos \zeta + \sin \alpha \sin \zeta \cos \phi \quad (3) $$

and $\lambda = \pm\sqrt{\frac{35}{16}(1-n)}$. Here $n \in [0, 1]$ defines the twist of the magnetic field lines. A positive $\lambda$ indicates an eastward twist of the magnetic field lines in the star's Southern hemisphere. Conversely, a negative $\lambda$ corresponds to a westward twist. When $n=1, \lambda=0$ and the twist perturbation becomes null. As an example, we show in Supplementary Fig. 11 the predictions from the twisted RVM for a geometry representative of the results derived in the previous section.

Following the toy model by ref. 70, we describe the evolution of $n$ during the outburst as an exponential decay of the form $n(t) = 1 - (1 - n_0)e^{-t/\tau_t}$, where $n_0$ is the initial twist value at $t = T_0$, the beginning of the 2018 outburst also assumed to be the start of the precession. $T_0$ is constrained by the X-ray detection of the outburst¹². The parameter $\tau_t$ represents the exponential decay timescale for the relaxation of the twisted magnetic field lines.
## Polarization of the pulse profiles

It has been previously reported¹⁵ that the PA of the linear polarization $L = \sqrt{Q^2 + U^2}$ of XTE J1810-197 pulse profiles recorded since the 2018 outburst shows complex structure and deviates from the simple 'S'-shaped curve predicted by the RVM¹⁶. The RVM is the standard geometrical model used by pulsar astronomers to interpret the PA and infer the viewing geometry of a pulsar, assuming a dipolar magnetic field.

Although the emission of the known radio magnetars has been shown to be usually highly linearly polarized, see, for example, refs. 13,45,46, we noticed some dips in the otherwise nearly 100% linearly polarized pulse profiles from multiple epochs at various phase ranges. This depolarization even occurs at different phase ranges for quasi-simultaneous observations at different frequencies (see, for example, the pulse profiles for MJD 58,474 in Supplementary Fig. 1). It is well known⁶⁵ that the superposition of different polarization mode can lead to depolarization and bias in the PA of the average pulse profile. To mitigate this, we excluded the phase ranges where the fractional linear polarization of the pulse drops below 80%. A two-dimensional histogram of the observed PA values of the linear polarization of the single-pulse profiles (Supplementary Fig. 8) shows that the depolarization in $L$ originates from the superposition of two different polarization modes and a PA branching, leading to some bias in the PA values after summation of the pulses. PA branching has been previously observed in the magnetar Swift J1818.0-1607 (ref. 22) and in a normal rotation-powered pulsar⁶⁶. This frequency-dependent effect, attributed to birefringence in the star's magnetosphere, was observed in the central part of the main component of the pulse starting around MJD 58,466 until it significantly decreased near MJD 58,500 (ref. 67). For the subsequent modelling of the PA data, these phase ranges were therefore excluded from the analysis.

## Comparing free and forced precession

The RVM, assuming the International Astronomical Union definition of the PA increasing counter-clockwise on the sky, can be written as:

$$ \tan(\psi_0 - \psi) = \frac{\sin \alpha \sin(\phi - \phi_0)}{\sin \zeta \cos \alpha - \cos \zeta \sin \alpha \cos(\phi - \phi_0)}, \quad (1) $$

with the viewing angle $\zeta = \alpha + \beta$ (ref. 17). Here $\alpha$ is the magnetic inclination angle and $\beta$ is the angle between the magnetic axis and our line of sight at the closest approach. The measured PA at any given pulse longitude $\phi$ is given by $\psi = \frac{1}{2} \arctan\left(\frac{U}{Q}\right)$. The parameter $\phi_0$ is the pulse longitude where the PA sweep is the steepest, that is, when $\psi = \psi_0$ (the PA of the rotation axis projected on the plane of the sky).

In contrast to the data recorded in 2006-2008 where an inter-pulse preceded the main pulse by about 100° (refs. 13,14), RVM fits to the PA of the individual epochs recorded after 2018 can only poorly constrain $\alpha$ and $\beta$. Therefore, we performed simultaneous RVM fits to all epochs to model the temporal PA variations as either free or forced precession. For the forced precession scenario, we set a common $\alpha$ parameter across all epochs whereas for the free-precession scenario, $\zeta$ is the common parameter across all epochs. For each epoch, the three remaining RVM parameters $\beta$, $\phi_0$ and $\psi_0$ are modelled independently.

We developed a set of codes based on modelRVM⁶⁸ to simultaneously apply the RVM to the Stokes Q and U data of our dataset (instead of the traditional PA values) and explore the parameter space of the free and forced precession models with the nested sampling software PolyChord v.1.20.1 (ref. 33).

The analysis includes a scaling factor parameter per frequency band referred to as EFAC that acts as a multiplier factor to the Stokes Q and U off-pulse standard deviations. The dimensionality of the problem $N_{\text{dim}}$ is therefore $1 + 3 \times N_{\text{epochs}} + N_{\text{EFAC}}$, with $N_{\text{epochs}}$ and $N_{\text{EFAC}}$ being the number of epochs and the number of EFAC parameters included in the

analysis, respectively. All RVM parameters are sampled using $10 \times N_{\text{dim}}$ live points (with a minimum of 1,000 live points) from Gaussian priors with standard deviation of 5° (3° for $\beta$) and mean values derived from previous runs. We set the PolyChord $n_{\text{repeats}}$ parameter (that is the slice sampling chain length used to produce a new live point) to $5 \times N_{\text{dim}}$ to ensure we get reliable log-evidence ($\log z$) estimates from the nested sampling analysis. Each EFAC parameter is sampled from log-uniform prior in the $\log_{10}$ range [-0.4, 0.6]. The log-Bayes factor comes directly from the subtraction of $\log z$ between any two models. A log-Bayes factor >3 brings conclusive support for the model with the highest $\log z$⁶⁹.

In the following analyses, we selected the PA with phase bins that satisfy $L > 3\sigma_N$, where $\sigma_N$ is the off-pulse standard deviation, outside of the excluded phase ranges as described in the previous section.

We started by analysing the datasets from the two different radio bands independently (L and C/X bands) and then combined all the PA data from the two bands. The results are reported in Supplementary Table 1.

The log-Bayes factor are 54 for the JBO L-band and 9 for the Effelsberg C/X-band data in support of the free-precession interpretation of the PA. They thus provide very strong support the free-precession interpretation of the data. Unsurprisingly, the combined dataset also shows unarguable support for the free-precession scenario with a log-Bayes factor of 75.

As the free-precession model is the preferred model to explain the observed PA variations, we show in Supplementary Fig. 9 the posterior distributions for the parameter $\zeta$ from each analysis. The $\zeta$ posteriors from the two frequency bands and the combined dataset are consistent with each other. The posterior from the combined dataset gives a narrower distribution, with $\zeta = 161.0^{\circ} \pm 0.7^{\circ}$, arguing for XTE J1810-197 as being a near-aligned rotator. We show in Supplementary Fig. 10 the posterior distributions of the two EFAC parameters included in the analysis of the combined dataset. The mean values of the EFAC parameters are 1.12 and 1.44 for the L and C/X bands, respectively, indicating a good modelling of the PA.

## Twisted magnetic fields

It is assumed²⁵ that the magnetic field of a magnetar can become twisted following some abrupt crustal motion of the star. The following slow and gradual untwisting of the magnetic field lines is presumed to power the observed outburst. The perturbation of the twisted field lines on the PA $\psi$ can be written as²⁰:

$$ \Delta\psi_{\text{twist}} = -\frac{8}{9}\lambda\sin^2\theta_{\text{obs}} \quad (2) $$

where $\theta_{\text{obs}}$ is the colatitude of the line of sight defined by

$$ \cos \theta_{\text{obs}} = \cos \alpha \cos \zeta + \sin \alpha \sin \zeta \cos \phi \quad (3) $$

and $\lambda = \pm\sqrt{\frac{35}{16}(1-n)}$. Here $n \in [0, 1]$ defines the twist of the magnetic field lines. A positive $\lambda$ indicates an eastward twist of the magnetic field lines in the star's Southern hemisphere. Conversely, a negative $\lambda$ corresponds to a westward twist. When $n=1, \lambda=0$ and the twist perturbation becomes null. As an example, we show in Supplementary Fig. 11 the predictions from the twisted RVM for a geometry representative of the results derived in the previous section.

Following the toy model by ref. 70, we describe the evolution of $n$ during the outburst as an exponential decay of the form $n(t) = 1 - (1 - n_0)e^{-t/\tau_t}$, where $n_0$ is the initial twist value at $t = T_0$, the beginning of the 2018 outburst also assumed to be the start of the precession. $T_0$ is constrained by the X-ray detection of the outburst¹². The parameter $\tau_t$ represents the exponential decay timescale for the relaxation of the twisted magnetic field lines.
The results of the model comparison between free and forced precession including the twist perturbations (three additional parameters) are reported in Supplementary Table 2. They support with a log-Bayes factor of 322 the relaxation of an eastward twist of the magnetic field lines in the free-precession model as the preferred model. Extended Data Fig. 1 shows the waveform for the twist parameter *n*. These results suggest a very large twist of the magnetic field lines at the time of the outburst with *n* = 0.05<sup>+0.16</sup><sub>-0.05</sub> at a 99% confidence interval with τ<sub>t</sub> > 730 days. The parameters *n*, *ζ* and *β* are covariant and with the twisted RVM, the value of the viewing angle has shifted to *ζ* = 167.3° ± 0.3°. The mean values of the EFAC parameters have decreased to 1.08 and 1.43 for the L and C/X bands, respectively.

### Emission height

The radio emission height $h_{em}$ can be inferred from the offset $\Delta\phi_{A/R}$ between the phase of the RVM inflection point (for example, $\phi_0$) and the centre of the pulse profile due to aberration and retardation effects<sup>71</sup>. It can be written as

$$h_{em} = \frac{\Delta\phi_{A/R}}{4} R_{LC} \quad (4)$$

where $R_{LC} = cP_s/2\pi$ is the light cylinder radius and *c* is the speed of light. The phases of the RVM inflection points are determined from the results of the free-precession modelling and the centre of the pulse is determined at the 4% and 0.5% total intensity levels for the JBO and Effelsberg data, respectively. We were able to derive emission heights for the Effelsberg data only from the first six epochs due to a drop of S/N for the subsequent epochs. Supplementary Fig. 12 shows the emission heights derived from the JBO and Effelsberg data. They show an increasing emission height in the range of 10,000 to 45,000 km with larger $|\beta|$ values. Interestingly, these results are also consistent with the concept of a radius-to-frequency mapping<sup>71,72</sup>, with higher radio frequency emission emitted at lower altitudes. Most of the normal pulsar population have emission heights less than 1,000 km (ref. 73) but magnetars tend to have higher emission altitudes (Fig. 14 of ref. 74).

The radio emission altitude is therefore well above the neutron star surface. The ratio $h_{em}/R$ where *R* is the neutron star radius, typically of size 12 km, is larger than 1,000 in all observations. At these heights, the strength of the magnetic field significantly decreases from its value at the surface, at least by a factor $(h_{em}/R)^3$ for the dipole and a factor $(h_{em}/R)^{\ell+2}$ for a multipole of order $\ell \ge 2$, the dipole being a $\ell = 1$ multipole. The surface dipole magnetic field strength is estimated to be around $B_{surf} \approx 2.7 \times 10^{10}$ T (ref. 44). At the radio emission site it is a billion times lower, amounting to a few Tesla or less. Even if quadrupolar surface magnetic fields are dominant at the surface, remaining dominant at an altitude of thousands of stellar radii would require a strength at the surface of at least 1,000 times larger than the dipole. Therefore, a physical process should maintain such a strong multipolar field for a time at least equal to the age of the magnetar. This is a stringent constraint on the magnetic field evolution and decay, difficult to explain with our current understanding of neutron star magnetic fields. Ref. 75 extensively discusses on the impact of multipolar fields on the neutron star spin-down luminosity. It has recently been shown<sup>76</sup> that adding a multipolar magnetic field component to the dipole does not alter the RVM PA expectations. However their results rely on an axisymmetric configuration where the multipole is oriented with respect to the dipole in such a way that the planes containing the field lines remain unchanged. If the orientation of this multipole would be random compared to the dipole, which is a more realistic configuration, and if this multipole would remain dominant at the radio emission heights derived above (a highly unlikely case) the RVM PA would be significantly altered. As a consequence, the impact of a multipole component can be confidently ignored when studying the PA evolution in radio.

### Correlation between spin-down and $\alpha$

The spin-down rate of neutron stars (hereafter referred as $\dot{v}$ similarly to the pulsar timing literature) is thought to vary approximately with $\sin^2\alpha$ (for example, ref. 47). We investigate here this correlation. We decided on using the Spearman correlation coefficient due to the non-Gaussian distribution of the data. To apply it, we first estimated $\dot{v}$ at the epochs of our dataset by linearly interpolating the $\dot{v}$ measurements from ref. 35 before computing the Spearman correlation coefficient $\rho$. The $\alpha$ measurements outside of the $\dot{v}$ measurements window were discarded. We find $\rho = -0.03$ with a *P* value of 0.86, where the null hypothesis is that our two datasets are linearly uncorrelated. This result, shown in Supplementary Fig. 13 supports the interpretation that we see no apparent correlation between the spin-down and the variation of geometry or that this variation is not contributing significantly to the spin-down compared to the intrinsic rotational instability commonly observed in magnetars.

### Reprocessing of the Effelsberg 2006 archival data

XTE J1810–197 was observed in the summer of 2006 with Effelsberg as part of a multifrequency observing campaign with the Lovell and Westerbork telescopes to study the magnetar’s polarization<sup>13</sup>. In this work, the PA data from the two main pulse components observed in 2006 was fitted independently to the RVM, but assuming a common $\zeta$. For the observation recorded on MJD 53,934 at 8.4 GHz (session 3), a RVM fit gave $\zeta = 83^\circ$ using the historical RVM convention<sup>53</sup>. This result would translate into $\zeta = 180 - \zeta_o = 97^\circ$ when referring to the same IAU definition of the PA as used in this work. The 2006 result is therefore in stark contrast to our current value of $\zeta$. We therefore decided to revisit the 2006 RVM results. We chose to focus on the Effelsberg 8.4 GHz observations as they exhibit the profiles with the highest S/N<sup>13</sup>. We selected the four epochs recorded at 8.4 GHz with the highest S/N where both components were observed, that is MJDs 53,926, 53,934, 53,938 and 53,944.

We performed separate RVM fits to these four epochs using MODELRVM, including in each case the four free parameters of the RVM plus an EFAC parameter. We used all PA data points with $L > 4\sigma_N$, except for MJD 53,944 where high level of linear depolarization is observed in the second component (see also Fig. 9 of ref. 13) and the corresponding phase range was excluded from the fit. Extended Data Fig. 2 shows the four selected polarimetric pulse profiles used to fit the PA data and Extended Data Fig. 3 shows the RVM posterior distributions for each epoch. Supplementary Table 3 summarizes the results of the modelling. We find that all modelled epochs from 2006 provide results on $\zeta$ consistent with our new dataset. We also find the angle $\beta$ to be constant between all four epochs within our error bars. We interpret the discrepancy between these new results and the results published in ref. 13 as being due to the exclusion of the PA data when depolarization is observed in the second pulse component.

We also applied the twisted RVM to each of the four epochs to try to measure the twist parameter for these observations taken two years before the radio disappearance of the magnetar. Unfortunately, the model could not constrain the twist parameter for any of the epochs and we attribute this to the narrowness of the pulse profile components in 2006.

Data from the Parkes radio telescope recorded in 2006 were also used to derive the geometry of XTE J1810–197 (ref. 14). Correcting for the likely use of the previous older convention for the RVM<sup>53</sup>, their results (no uncertainties were provided) would translate to $\alpha \approx 176^\circ$ and $\beta \approx -4^\circ$, in broad agreement with our reprocessing of the Effelsberg 2006 data.

### Theoretical models

**Free precession of neutron stars.** To set the body into free precession, a rotating NS must have some deformations misaligned with the
The results of the model comparison between free and forced precession including the twist perturbations (three additional parameters) are reported in Supplementary Table 2. They support with a log-Bayes factor of 322 the relaxation of an eastward twist of the magnetic field lines in the free-precession model as the preferred model. Extended Data Fig. 1 shows the waveform for the twist parameter *n*. These results suggest a very large twist of the magnetic field lines at the time of the outburst with *n* = 0.05<sup>+0.16</sup><sub>-0.05</sub> at a 99% confidence interval with τ<sub>t</sub> > 730 days. The parameters *n*, *ζ* and *β* are covariant and with the twisted RVM, the value of the viewing angle has shifted to *ζ* = 167.3° ± 0.3°. The mean values of the EFAC parameters have decreased to 1.08 and 1.43 for the L and C/X bands, respectively.

## Emission height

The radio emission height $h_{em}$ can be inferred from the offset $\Delta\phi_{A/R}$ between the phase of the RVM inflection point (for example, $\phi_0$) and the centre of the pulse profile due to aberration and retardation effects<sup>71</sup>. It can be written as

$$h_{em} = \frac{\Delta\phi_{A/R}}{4} R_{LC} \quad (4)$$

where $R_{LC} = cP_s/2\pi$ is the light cylinder radius and *c* is the speed of light. The phases of the RVM inflection points are determined from the results of the free-precession modelling and the centre of the pulse is determined at the 4% and 0.5% total intensity levels for the JBO and Effelsberg data, respectively. We were able to derive emission heights for the Effelsberg data only from the first six epochs due to a drop of S/N for the subsequent epochs. Supplementary Fig. 12 shows the emission heights derived from the JBO and Effelsberg data. They show an increasing emission height in the range of 10,000 to 45,000 km with larger $|\beta|$ values. Interestingly, these results are also consistent with the concept of a radius-to-frequency mapping<sup>71,72</sup>, with higher radio frequency emission emitted at lower altitudes. Most of the normal pulsar population have emission heights less than 1,000 km (ref. 73) but magnetars tend to have higher emission altitudes (Fig. 14 of ref. 74).

The radio emission altitude is therefore well above the neutron star surface. The ratio $h_{em}/R$ where *R* is the neutron star radius, typically of size 12 km, is larger than 1,000 in all observations. At these heights, the strength of the magnetic field significantly decreases from its value at the surface, at least by a factor $(h_{em}/R)^3$ for the dipole and a factor $(h_{em}/R)^{\ell+2}$ for a multipole of order $\ell \ge 2$, the dipole being a $\ell = 1$ multipole. The surface dipole magnetic field strength is estimated to be around $B_{surf} \approx 2.7 \times 10^{10}$ T (ref. 44). At the radio emission site it is a billion times lower, amounting to a few Tesla or less. Even if quadrupolar surface magnetic fields are dominant at the surface, remaining dominant at an altitude of thousands of stellar radii would require a strength at the surface of at least 1,000 times larger than the dipole. Therefore, a physical process should maintain such a strong multipolar field for a time at least equal to the age of the magnetar. This is a stringent constraint on the magnetic field evolution and decay, difficult to explain with our current understanding of neutron star magnetic fields. Ref. 75 extensively discusses on the impact of multipolar fields on the neutron star spin-down luminosity. It has recently been shown<sup>76</sup> that adding a multipolar magnetic field component to the dipole does not alter the RVM PA expectations. However their results rely on an axisymmetric configuration where the multipole is oriented with respect to the dipole in such a way that the planes containing the field lines remain unchanged. If the orientation of this multipole would be random compared to the dipole, which is a more realistic configuration, and if this multipole would remain dominant at the radio emission heights derived above (a highly unlikely case) the RVM PA would be significantly altered. As a consequence, the impact of a multipole component can be confidently ignored when studying the PA evolution in radio.

## Correlation between spin-down and $\alpha$

The spin-down rate of neutron stars (hereafter referred as $\dot{v}$ similarly to the pulsar timing literature) is thought to vary approximately with $\sin^2\alpha$ (for example, ref. 47). We investigate here this correlation. We decided on using the Spearman correlation coefficient due to the non-Gaussian distribution of the data. To apply it, we first estimated $\dot{v}$ at the epochs of our dataset by linearly interpolating the $\dot{v}$ measurements from ref. 35 before computing the Spearman correlation coefficient $\rho$. The $\alpha$ measurements outside of the $\dot{v}$ measurements window were discarded. We find $\rho = -0.03$ with a *P* value of 0.86, where the null hypothesis is that our two datasets are linearly uncorrelated. This result, shown in Supplementary Fig. 13 supports the interpretation that we see no apparent correlation between the spin-down and the variation of geometry or that this variation is not contributing significantly to the spin-down compared to the intrinsic rotational instability commonly observed in magnetars.

## Reprocessing of the Effelsberg 2006 archival data

XTE J1810–197 was observed in the summer of 2006 with Effelsberg as part of a multifrequency observing campaign with the Lovell and Westerbork telescopes to study the magnetar’s polarization<sup>13</sup>. In this work, the PA data from the two main pulse components observed in 2006 was fitted independently to the RVM, but assuming a common $\zeta$. For the observation recorded on MJD 53,934 at 8.4 GHz (session 3), a RVM fit gave $\zeta = 83^\circ$ using the historical RVM convention<sup>53</sup>. This result would translate into $\zeta = 180 - \zeta_o = 97^\circ$ when referring to the same IAU definition of the PA as used in this work. The 2006 result is therefore in stark contrast to our current value of $\zeta$. We therefore decided to revisit the 2006 RVM results. We chose to focus on the Effelsberg 8.4 GHz observations as they exhibit the profiles with the highest S/N<sup>13</sup>. We selected the four epochs recorded at 8.4 GHz with the highest S/N where both components were observed, that is MJDs 53,926, 53,934, 53,938 and 53,944.

We performed separate RVM fits to these four epochs using MODELRVM, including in each case the four free parameters of the RVM plus an EFAC parameter. We used all PA data points with $L > 4\sigma_N$, except for MJD 53,944 where high level of linear depolarization is observed in the second component (see also Fig. 9 of ref. 13) and the corresponding phase range was excluded from the fit. Extended Data Fig. 2 shows the four selected polarimetric pulse profiles used to fit the PA data and Extended Data Fig. 3 shows the RVM posterior distributions for each epoch. Supplementary Table 3 summarizes the results of the modelling. We find that all modelled epochs from 2006 provide results on $\zeta$ consistent with our new dataset. We also find the angle $\beta$ to be constant between all four epochs within our error bars. We interpret the discrepancy between these new results and the results published in ref. 13 as being due to the exclusion of the PA data when depolarization is observed in the second pulse component.

We also applied the twisted RVM to each of the four epochs to try to measure the twist parameter for these observations taken two years before the radio disappearance of the magnetar. Unfortunately, the model could not constrain the twist parameter for any of the epochs and we attribute this to the narrowness of the pulse profile components in 2006.

Data from the Parkes radio telescope recorded in 2006 were also used to derive the geometry of XTE J1810–197 (ref. 14). Correcting for the likely use of the previous older convention for the RVM<sup>53</sup>, their results (no uncertainties were provided) would translate to $\alpha \approx 176^\circ$ and $\beta \approx -4^\circ$, in broad agreement with our reprocessing of the Effelsberg 2006 data.

## Theoretical models

**Free precession of neutron stars.** To set the body into free precession, a rotating NS must have some deformations misaligned with the
centrifugal bulge. We can write the moment of inertia tensor of the NS as a sum of a spherical and two quadrupolar parts$^{77,78}$

$$I = I_0 \delta + \Delta I_{\omega} (\mathbf{n}_{\omega} \mathbf{n}_{\omega} - \delta/3) + \Delta I_{d} (\mathbf{n}_{d} \mathbf{n}_{d} - \delta/3). \quad (5)$$

The first term on the right-hand side is the spherical part of the non-rotating undeformed star with $\delta$ being the unit tensor. The second term is the change due to centrifugal forces, and has the unit vector $\mathbf{n}_{\omega}$ as the symmetric axis along the star's angular velocity $\omega$. The third term is the change due to crustal shear stresses or magnetic fields, which is assumed to be axisymmetric along a unit vector $\mathbf{n}_{d}$ for simplicity.

The angular momentum is then related to the angular velocity via

$$\mathbf{L} = (I_0 + 2/3 \Delta I_{\omega} - 1/3 \Delta I_{d}) \omega + \Delta I_{d} \omega_3 \mathbf{n}_{d}. \quad (6)$$

This shows that $\mathbf{L}$, $\omega$, and $\mathbf{n}_{d}$ are coplanar. As the angular momentum is conserved in free precession, this plane must rotate around $\mathbf{L}$. Taking the components of $\mathbf{L}$, we obtain

$$L_1 = (I_0 + 2/3 \Delta I_{\omega} - 1/3 \Delta I_{d}) \omega_1 \equiv I_1 \omega_1, \quad (7)$$

$$L_2 = (I_0 + 2/3 \Delta I_{\omega} - 1/3 \Delta I_{d}) \omega_2 \equiv I_1 \omega_2, \quad (8)$$

$$L_3 = (I_0 + 2/3 \Delta I_{\omega} + 2/3 \Delta I_{d}) \omega_3 \equiv I_3 \omega_3. \quad (9)$$

Despite the triaxiality of $I$, the components of the angular momentum are formally equivalent to a rigid symmetric rotator with $I_3 - I_1 = \Delta I_d$. We define the ellipticity of the deformed NS as

$$\epsilon \equiv \frac{I_3 - I_1}{I_1} = \frac{\Delta I_d}{I_1}. \quad (10)$$

In the body frame, the equation of motion of a freely precessing body can be described by the Euler equation

$$\dot{\mathbf{L}} + \mathbf{\omega} \times \mathbf{L} = 0, \quad (11)$$

where the over dot denotes the derivative with respect to the time $t$. The equation of motion only involves $\mathbf{L}$ and $\mathbf{\omega}$. Thus, equations (7)–(9) indicate that the precession of the triaxial body is formally equivalent to that of a rigid symmetric top$^{78}$. According to equation (11), the components of the angular velocities satisfy

$$\dot{\omega}_1 = -\epsilon \omega_2 \omega_3, \quad (12)$$

$$\dot{\omega}_2 = \epsilon \omega_1 \omega_3, \quad (13)$$

$$\dot{\omega}_3 = 0. \quad (14)$$

The solutions of equations (12)–(14) are

$$\omega_1 = a \cos(\omega_p t + \beta_0), \quad (15)$$

$$\omega_2 = a \sin(\omega_p t + \beta_0), \quad (16)$$

$$\omega_3 = b, \quad (17)$$

where $a$ and $b$ are constants satisfying $a^2 + b^2 = \omega^2$, $\omega_p = \epsilon\omega_3$ is the angular frequency and $\beta_0$ is an initial phase.

We introduce standard Euler angles to describe the orientation of the NS in the inertial frame, with the polar axis along the angular momentum $\mathbf{L}$. Let $\eta$ and $\theta$ denote the azimuthal and polar angles of the

deformation axis $\mathbf{n}_{d}$, and $\lambda$ represent a rotation about $\mathbf{n}_{d}$. As $\omega_3$ and $\mathbf{L}$ are constant, the wobble angle $\theta$ is also constant during the precession. We label the angle between $\mathbf{\omega}$ and $\mathbf{L}$ as $\theta_1$, which satisfies

$$\tan(\theta + \theta_1) = \frac{I_3}{I_1} \tan \theta = (1 + \epsilon) \tan \theta. \quad (18)$$

Expanding the above equation to the first order of $\epsilon$, we obtain

$$\theta_1 = \epsilon \sin \theta \cos \theta. \quad (19)$$

This angle is much smaller than $\theta$ because $\epsilon$ is extremely small. Therefore, we can neglect $\theta_1$ when evaluating the geometry of the star. To get the motion in the inertial frame, we decompose the angular velocity into

$$\mathbf{\omega} = \dot{\eta} \mathbf{n}_{\mathbf{L}} + \dot{\lambda} \mathbf{n}_{\mathbf{d}}, \quad (20)$$

where $\mathbf{n}_{\mathbf{L}}$ is the unit vector along $\mathbf{L}$ and $\dot{\eta}$ and $\dot{\lambda}$ are the time derivatives of $\eta$ and $\lambda$, respectively. Substituting this equation into equation (6), we get

$$\mathbf{L} = I_1 \dot{\eta}, \quad (21)$$

$$\dot{\lambda} = -\epsilon \omega_3 = -\omega_p. \quad (22)$$

Thus, as viewed from the inertia frame, the deformation axis rotates at a rate of $\dot{\eta}$ about $\mathbf{L}$ in a cone of half-angle $\theta$. We refer this angular frequency as the inertial precession frequency. Superimposed upon $\dot{\eta}\mathbf{n}_{\mathbf{L}}$ is a rotation about the deformation axis at a rate of $\dot{\lambda} = -\omega_p$. We refer to $\omega_p$ as the free-precession frequency and define the free-precession period as

$$P_f = \frac{2\pi}{|\omega_p|} \approx \frac{P_s}{|\epsilon \cos \theta|}, \quad (23)$$

where $P_s = 2\pi/\omega$ is the spin period.

The magnetic inclination $\alpha$ between $\mathbf{L}$ and the magnetic dipole $\mathbf{m}$ can be represented as

$$\cos \alpha = \sin \theta \sin \lambda \sin \chi + \cos \theta \cos \chi, \quad (24)$$

where $\chi$ is the angle between $\mathbf{m}$ and $\mathbf{n}_{d}$, and the precession phase $\lambda$ is

$$\lambda = \arctan \frac{\omega_1}{\omega_2} = \frac{\pi}{2} - \omega_p t - \beta_0. \quad (25)$$

The precession phase $\lambda$ changes with time as long as $\chi \neq 0$, which leads to the periodic variations of $\alpha$ during free precession. However, for XTE J1810–197, the variations of $\alpha$ damp away in less than a precession period. The simple free precession cannot explain this phenomenon. In the following, we investigate damped-precession models based on the free-precession solution.

**Free precession with relaxing ellipticity model.** In this model, we expand on the geometry of a freely precessing NS to consider a decreasing ellipticity $\epsilon$ of the NS such that $P_f$ becomes increasingly longer. The star gradually relaxes to a spherical shape. As a result, the precession phase is frozen at a certain epoch and $\alpha$ turns into constant. We assume

$$I_3 = I_1 [1 + \epsilon(t)], \quad (26)$$

with $\epsilon < 0$. In this case, the Euler equations are

$$\dot{\omega}_1 = -\epsilon \omega_2 \omega_3, \quad (27)$$
centrifugal bulge. We can write the moment of inertia tensor of the NS as a sum of a spherical and two quadrupolar parts$^{77,78}$

$$I = I_0 \delta + \Delta I_{\omega} (\mathbf{n}_{\omega} \mathbf{n}_{\omega} - \delta/3) + \Delta I_{d} (\mathbf{n}_{d} \mathbf{n}_{d} - \delta/3). \quad (5)$$

The first term on the right-hand side is the spherical part of the non-rotating undeformed star with $\delta$ being the unit tensor. The second term is the change due to centrifugal forces, and has the unit vector $\mathbf{n}_{\omega}$ as the symmetric axis along the star's angular velocity $\omega$. The third term is the change due to crustal shear stresses or magnetic fields, which is assumed to be axisymmetric along a unit vector $\mathbf{n}_{d}$ for simplicity.

The angular momentum is then related to the angular velocity via

$$\mathbf{L} = (I_0 + 2/3 \Delta I_{\omega} - 1/3 \Delta I_{d}) \omega + \Delta I_{d} \omega_3 \mathbf{n}_{d}. \quad (6)$$

This shows that $\mathbf{L}$, $\omega$, and $\mathbf{n}_{d}$ are coplanar. As the angular momentum is conserved in free precession, this plane must rotate around $\mathbf{L}$. Taking the components of $\mathbf{L}$, we obtain

$$L_1 = (I_0 + 2/3 \Delta I_{\omega} - 1/3 \Delta I_{d}) \omega_1 \equiv I_1 \omega_1, \quad (7)$$

$$L_2 = (I_0 + 2/3 \Delta I_{\omega} - 1/3 \Delta I_{d}) \omega_2 \equiv I_1 \omega_2, \quad (8)$$

$$L_3 = (I_0 + 2/3 \Delta I_{\omega} + 2/3 \Delta I_{d}) \omega_3 \equiv I_3 \omega_3. \quad (9)$$

Despite the triaxiality of $I$, the components of the angular momentum are formally equivalent to a rigid symmetric rotator with $I_3 - I_1 = \Delta I_d$. We define the ellipticity of the deformed NS as

$$\epsilon \equiv \frac{I_3 - I_1}{I_1} = \frac{\Delta I_d}{I_1}. \quad (10)$$

In the body frame, the equation of motion of a freely precessing body can be described by the Euler equation

$$\dot{\mathbf{L}} + \mathbf{\omega} \times \mathbf{L} = 0, \quad (11)$$

where the over dot denotes the derivative with respect to the time $t$. The equation of motion only involves $\mathbf{L}$ and $\mathbf{\omega}$. Thus, equations (7)–(9) indicate that the precession of the triaxial body is formally equivalent to that of a rigid symmetric top$^{78}$. According to equation (11), the components of the angular velocities satisfy

$$\dot{\omega}_1 = -\epsilon \omega_2 \omega_3, \quad (12)$$

$$\dot{\omega}_2 = \epsilon \omega_1 \omega_3, \quad (13)$$

$$\dot{\omega}_3 = 0. \quad (14)$$

The solutions of equations (12)–(14) are

$$\omega_1 = a \cos(\omega_p t + \beta_0), \quad (15)$$

$$\omega_2 = a \sin(\omega_p t + \beta_0), \quad (16)$$

$$\omega_3 = b, \quad (17)$$

where $a$ and $b$ are constants satisfying $a^2 + b^2 = \omega^2$, $\omega_p = \epsilon\omega_3$ is the angular frequency and $\beta_0$ is an initial phase.

We introduce standard Euler angles to describe the orientation of the NS in the inertial frame, with the polar axis along the angular momentum $\mathbf{L}$. Let $\eta$ and $\theta$ denote the azimuthal and polar angles of the

deformation axis $\mathbf{n}_{d}$, and $\lambda$ represent a rotation about $\mathbf{n}_{d}$. As $\omega_3$ and $\mathbf{L}$ are constant, the wobble angle $\theta$ is also constant during the precession. We label the angle between $\mathbf{\omega}$ and $\mathbf{L}$ as $\theta_1$, which satisfies

$$\tan(\theta + \theta_1) = \frac{I_3}{I_1} \tan \theta = (1 + \epsilon) \tan \theta. \quad (18)$$

Expanding the above equation to the first order of $\epsilon$, we obtain

$$\theta_1 = \epsilon \sin \theta \cos \theta. \quad (19)$$

This angle is much smaller than $\theta$ because $\epsilon$ is extremely small. Therefore, we can neglect $\theta_1$ when evaluating the geometry of the star. To get the motion in the inertial frame, we decompose the angular velocity into

$$\mathbf{\omega} = \dot{\eta} \mathbf{n}_{\mathbf{L}} + \dot{\lambda} \mathbf{n}_{\mathbf{d}}, \quad (20)$$

where $\mathbf{n}_{\mathbf{L}}$ is the unit vector along $\mathbf{L}$ and $\dot{\eta}$ and $\dot{\lambda}$ are the time derivatives of $\eta$ and $\lambda$, respectively. Substituting this equation into equation (6), we get

$$\mathbf{L} = I_1 \dot{\eta}, \quad (21)$$

$$\dot{\lambda} = -\epsilon \omega_3 = -\omega_p. \quad (22)$$

Thus, as viewed from the inertia frame, the deformation axis rotates at a rate of $\dot{\eta}$ about $\mathbf{L}$ in a cone of half-angle $\theta$. We refer this angular frequency as the inertial precession frequency. Superimposed upon $\dot{\eta}\mathbf{n}_{\mathbf{L}}$ is a rotation about the deformation axis at a rate of $\dot{\lambda} = -\omega_p$. We refer to $\omega_p$ as the free-precession frequency and define the free-precession period as

$$P_f = \frac{2\pi}{|\omega_p|} \approx \frac{P_s}{|\epsilon \cos \theta|}, \quad (23)$$

where $P_s = 2\pi/\omega$ is the spin period.

The magnetic inclination $\alpha$ between $\mathbf{L}$ and the magnetic dipole $\mathbf{m}$ can be represented as

$$\cos \alpha = \sin \theta \sin \lambda \sin \chi + \cos \theta \cos \chi, \quad (24)$$

where $\chi$ is the angle between $\mathbf{m}$ and $\mathbf{n}_{d}$, and the precession phase $\lambda$ is

$$\lambda = \arctan \frac{\omega_1}{\omega_2} = \frac{\pi}{2} - \omega_p t - \beta_0. \quad (25)$$

The precession phase $\lambda$ changes with time as long as $\chi \neq 0$, which leads to the periodic variations of $\alpha$ during free precession. However, for XTE J1810–197, the variations of $\alpha$ damp away in less than a precession period. The simple free precession cannot explain this phenomenon. In the following, we investigate damped-precession models based on the free-precession solution.

**Free precession with relaxing ellipticity model.** In this model, we expand on the geometry of a freely precessing NS to consider a decreasing ellipticity $\epsilon$ of the NS such that $P_f$ becomes increasingly longer. The star gradually relaxes to a spherical shape. As a result, the precession phase is frozen at a certain epoch and $\alpha$ turns into constant. We assume

$$I_3 = I_1 [1 + \epsilon(t)], \quad (26)$$

with $\epsilon < 0$. In this case, the Euler equations are

$$\dot{\omega}_1 = -\epsilon \omega_2 \omega_3, \quad (27)$$
$$ \dot{\omega}_2 = \epsilon\omega_1\omega_3, \qquad (28) $$

$$ \dot{\omega}_3 = -\frac{\dot{\epsilon}\omega_3}{1+\epsilon}. \qquad (29) $$

We can neglect the change of $\omega_3$ and $\theta$ when studying the geometry because the ellipticity is very small. In this approximation, the components of $\omega$ can be expressed as

$$ \omega_1(t) = a \cos[g(t) + \beta_0], \qquad (30) $$

$$ \omega_2(t) = a \sin[g(t) + \beta_0], \qquad (31) $$

$$ \omega_3(t) = b, \qquad (32) $$

where

$$ g(t) = b \int_{0}^{t} \epsilon(t) dt. \qquad (33) $$

The angle $\lambda$ is

$$ \lambda = \arctan \frac{\omega_1}{\omega_2} = \frac{\pi}{2} - g(t) - \beta_0. \qquad (34) $$

For simplicity, we assume that the ellipticity changes in an exponential form

$$ \epsilon(t) = \epsilon_0 e^{-t/\tau_c}, \qquad (35) $$

where $\tau_c$ is the relaxation timescale of the ellipticity. In this parametrization, the function $g(t)$ is

$$ g(t) = b\epsilon_0\tau - be^{-t/\tau_c}\epsilon_0\tau. \qquad (36) $$

**Phenomenological model of decaying wobble angle.** The decay of the wobble angle can also make the precession damp away, which is closely related to the internal couplings between the core and the crust. For simplicity, we first take a phenomenological model in which the wobble angle decays in an exponential form

$$ \theta(t) = \theta_0 e^{-t/\tau_c}, \qquad (37) $$

and the precession phase evolves in the free-precession form

$$ \dot{\lambda}(t) = -\epsilon\omega \cos\theta(t), \qquad (38) $$

where $\tau_0$ is the decaying timescale of the wobble angle.

**Phenomenological model of decaying wobble angle and relaxing ellipticity.** The phenomenological model of decaying wobble angle is unable to fit the data well. Therefore, we combined the decay of the wobble angle with the relaxing ellipticity. We considered two cases for the ellipticity decrease to represent a possible constant contribution to $\epsilon$,

$$ \epsilon(t) = \epsilon_1 e^{-t/\tau_c}, \qquad (39) $$

and

$$ \epsilon(t) = \epsilon_0 + \epsilon_1 e^{-t/\tau_c}, \qquad (40) $$

where $\epsilon_0$ and $\epsilon_1$ are constant, and $\tau_c$ is the relaxing timescale of the ellipticity. In this model, the precession phase is

$$ \lambda(t) = \frac{\pi}{2} - \int_{0}^{t} \epsilon(t) \omega \cos \theta(t) dt - \beta_0. \qquad (41) $$

**Model of crust-core frictional coupling combined with relaxing ellipticity.** In this model, we replace the phenomenological description of the decay of the wobble angle with a frictional coupling between the core and the crust of the NS. Following refs. 79–81, we take a two-component model for the NS. The component which is coupled directly with the external torque $\mathbf{N}_{\text{ext}}$ consists of the crust and the charged fluids in the core. The other component, which contains most of the moment of inertia, is the fluid in the core of the star. We do not consider any superfluid-pinning to the crust, because pinning will inevitably make the precession very fast$^{25,80}$, which is contradictory to our observations. The two components are coupled through an internal torque $\mathbf{N}_{\text{int}}$. We label the two components with ‘c’ for the crust and ‘f’ for the fluid in the core.

Standing in the corotating frame of the crust, the Euler equations describing the motion of the angular frequency vector $\mathbf{\omega}_c$ and $\mathbf{\Omega}_f$ can be written as

$$ \dot{\mathbf{L}}_c + \mathbf{\omega}_c \times \mathbf{L}_c = \mathbf{N}_{\text{ext}} + \mathbf{N}_{\text{int}}, \qquad (42) $$

$$ \dot{\mathbf{L}}_f + \mathbf{\Omega}_f \times \mathbf{L}_f = -\mathbf{N}_{\text{int}}, \qquad (43) $$

where $\mathbf{L}_c$ and $\mathbf{L}_f$ are the angular momentum of the two components. We ignore $\mathbf{N}_{\text{ext}}$ since it only has large effects on the orientation of the star in spin-down timescale (it does play an important role on the spin evolution).

In general, the precession of the fluid core is driven by an internal torque depending on the rotational velocity lags between the solid and the fluid. With this velocity dependent torque, the angular momentum of the system is conserved but the kinetic energy is dissipated. The precession of the crust will be damped with decaying wobble angle. We take the internal torque as

$$ \mathbf{N}_{\text{int}} = K(\mathbf{\Omega}_f - \mathbf{\omega}_c), \qquad (44) $$

where K is a constant depending on the microscopic physics of the frictional coupling. Then the Euler equations for the components of $\mathbf{\omega}_c$ are

$$ \dot{\omega}_1 = -\epsilon\omega_2\omega_3 - \frac{1}{\tau_c}(\omega_1 - \Omega_1), \qquad (45) $$

$$ \dot{\omega}_2 = \epsilon\omega_1\omega_3 - \frac{1}{\tau_c}(\omega_2 - \Omega_2), \qquad (46) $$

$$ \dot{\omega}_3 = -\frac{1}{(1+\epsilon)\tau_c}(\omega_3 - \Omega_3) - \frac{u_3\dot{\epsilon}}{1+\epsilon}, \qquad (47) $$

where $\omega_1, \omega_2$ and $\omega_3$ are the components of $\mathbf{\omega}_c, \mathbf{\Omega}_1, \mathbf{\Omega}_2$ and $\mathbf{\Omega}_3$ are the components of $\mathbf{\Omega}_f$, and $\tau_c = I_c/K$ is the timescale for the frictional coupling. The Euler equations for the components of $\mathbf{\Omega}_f$ are

$$ \dot{\Omega}_1 = -(\omega_2\Omega_3 - \omega_3\Omega_2) + \frac{K}{\tau_c}(\omega_1 - \Omega_1), \qquad (48) $$

$$ \dot{\Omega}_2 = -(\omega_3\Omega_1 - \omega_1\Omega_3) + \frac{K}{\tau_c}(\omega_2 - \Omega_2), \qquad (49) $$

$$ \dot{\Omega}_3 = -(\omega_1\Omega_2 - \omega_2\Omega_1) + \frac{K}{\tau_c}(\omega_3 - \Omega_3), \qquad (50) $$
$$ \dot{\omega}_2 = \epsilon \omega_1 \omega_3, \qquad (28) $$

$$ \dot{\omega}_3 = -\frac{\dot{\epsilon}}{1+\epsilon} \omega_3. \qquad (29) $$

We can neglect the change of $\omega_3$ and $\theta$ when studying the geometry because the ellipticity is very small. In this approximation, the components of $\omega$ can be expressed as

$$ \omega_1(t) = a \cos[g(t) + \beta_0], \qquad (30) $$

$$ \omega_2(t) = a \sin[g(t) + \beta_0], \qquad (31) $$

$$ \omega_3(t) = b, \qquad (32) $$

where

$$ g(t) = b \int_{0}^{t} \epsilon(t) dt. \qquad (33) $$

The angle $\lambda$ is

$$ \lambda = \arctan \frac{\omega_1}{\omega_2} = \frac{\pi}{2} - g(t) - \beta_0. \qquad (34) $$

For simplicity, we assume that the ellipticity changes in an exponential form

$$ \epsilon(t) = \epsilon_0 e^{-t/\tau_c}, \qquad (35) $$

where $\tau_c$ is the relaxation timescale of the ellipticity. In this parametrization, the function $g(t)$ is

$$ g(t) = b \epsilon_0 \tau - b e^{-t/\tau_c} \epsilon_0 \tau. \qquad (36) $$

**Phenomenological model of decaying wobble angle.** The decay of the wobble angle can also make the precession damp away, which is closely related to the internal couplings between the core and the crust. For simplicity, we first take a phenomenological model in which the wobble angle decays in an exponential form

$$ \theta(t) = \theta_0 e^{-t/\tau_c}, \qquad (37) $$

and the precession phase evolves in the free-precession form

$$ \dot{\lambda}(t) = -\epsilon\omega\cos\theta(t), \qquad (38) $$

where $\tau_0$ is the decaying timescale of the wobble angle.

**Phenomenological model of decaying wobble angle and relaxing ellipticity.** The phenomenological model of decaying wobble angle is unable to fit the data well. Therefore, we combined the decay of the wobble angle with the relaxing ellipticity. We considered two cases for the ellipticity decrease to represent a possible constant contribution to $\epsilon$,

$$ \epsilon(t) = \epsilon_1 e^{-t/\tau_c}, \qquad (39) $$

and

$$ \epsilon(t) = \epsilon_0 + \epsilon_1 e^{-t/\tau_c}, \qquad (40) $$

where $\epsilon_0$ and $\epsilon_1$ are constant, and $\tau_c$ is the relaxing timescale of the ellipticity. In this model, the precession phase is

$$ \lambda(t) = \frac{\pi}{2} - \int_{0}^{t} \epsilon(t) \omega \cos \theta(t) dt - \beta_0. \qquad (41) $$

**Model of crust-core frictional coupling combined with relaxing ellipticity.** In this model, we replace the phenomenological description of the decay of the wobble angle with a frictional coupling between the core and the crust of the NS. Following refs. 79–81, we take a two-component model for the NS. The component which is coupled directly with the external torque $\mathbf{N}_{\text{ext}}$ consists of the crust and the charged fluids in the core. The other component, which contains most of the moment of inertia, is the fluid in the core of the star. We do not consider any superfluid-pinning to the crust, because pinning will inevitably make the precession very fast$^{25,80}$, which is contradictory to our observations. The two components are coupled through an internal torque $\mathbf{N}_{\text{int}}$. We label the two components with ‘c’ for the crust and ‘f’ for the fluid in the core.

Standing in the corotating frame of the crust, the Euler equations describing the motion of the angular frequency vector $\mathbf{\omega}_c$ and $\mathbf{\Omega}_f$ can be written as

$$ \dot{\mathbf{L}}_c + \mathbf{\omega}_c \times \mathbf{L}_c = \mathbf{N}_{\text{ext}} + \mathbf{N}_{\text{int}}, \qquad (42) $$

$$ \dot{\mathbf{L}}_f + \mathbf{\Omega}_f \times \mathbf{L}_f = -\mathbf{N}_{\text{int}}, \qquad (43) $$

where $\mathbf{L}_c$ and $\mathbf{L}_f$ are the angular momentum of the two components. We ignore $\mathbf{N}_{\text{ext}}$ since it only has large effects on the orientation of the star in spin-down timescale (it does play an important role on the spin evolution).

In general, the precession of the fluid core is driven by an internal torque depending on the rotational velocity lags between the solid and the fluid. With this velocity dependent torque, the angular momentum of the system is conserved but the kinetic energy is dissipated. The precession of the crust will be damped with decaying wobble angle. We take the internal torque as

$$ \mathbf{N}_{\text{int}} = K(\mathbf{\Omega}_f - \mathbf{\omega}_c), \qquad (44) $$

where K is a constant depending on the microscopic physics of the frictional coupling. Then the Euler equations for the components of $\mathbf{\omega}_c$ are

$$ \dot{\omega}_1 = -\epsilon\omega_2\omega_3 - \frac{1}{\tau_c}(\omega_1 - \Omega_1), \qquad (45) $$

$$ \dot{\omega}_2 = \epsilon\omega_1\omega_3 - \frac{1}{\tau_c}(\omega_2 - \Omega_2), \qquad (46) $$

$$ \dot{\omega}_3 = -\frac{1}{(1+\epsilon)\tau_c}(\omega_3 - \Omega_3) - \frac{u_3 \dot{\epsilon}}{1+\epsilon}, \qquad (47) $$

where $\omega_1$, $\omega_2$ and $\omega_3$ are the components of $\mathbf{\omega}_c$, $\mathbf{\Omega}_1$, $\mathbf{\Omega}_2$ and $\mathbf{\Omega}_3$ are the components of $\mathbf{\Omega}_f$, and $\tau_c = I_c/K$ is the timescale for the frictional coupling. The Euler equations for the components of $\mathbf{\Omega}_f$ are

$$ \dot{\Omega}_1 = -(\omega_2\Omega_3 - \omega_3\Omega_2) + \frac{K}{\tau_c}(\omega_1 - \Omega_1), \qquad (48) $$

$$ \dot{\Omega}_2 = -(\omega_3\Omega_1 - \omega_1\Omega_3) + \frac{K}{\tau_c}(\omega_2 - \Omega_2), \qquad (49) $$

$$ \dot{\Omega}_3 = -(\omega_1\Omega_2 - \omega_2\Omega_1) + \frac{K}{\tau_c}(\omega_3 - \Omega_3), \qquad (50) $$
where $\kappa = I_c/I_f$ is the ratio between the moment of inertia of the crust and the core.

In principle, one can integrate equations (45)–(50) with appropriate initial conditions. However, it is very slow to integrate those six equations directly. To do parameter estimation more efficiently, we first make the following transformations

$$ \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix} = \begin{pmatrix} \cos \gamma - \sin \gamma & 0 \\ \sin \gamma & \cos \gamma & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{pmatrix}, \quad (51) $$

and

$$ \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix} = \begin{pmatrix} \cos \gamma - \sin \gamma & 0 \\ \sin \gamma & \cos \gamma & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \Omega_1 \\ \Omega_2 \\ \Omega_3 \end{pmatrix}. \quad (52) $$

Here the angle $\gamma$ is defined as

$$ \gamma(t) = -\epsilon\omega_3 dt. \quad (53) $$

Then the Euler equations can be written as

$$ \dot{u}_1 = -\frac{1}{\tau_c}(u_1 - v_1), \quad (54) $$

$$ \dot{u}_2 = -\frac{1}{\tau_c}(u_2 - v_2), \quad (55) $$

$$ \dot{u}_3 = -\frac{1}{\tau_c(1+\epsilon)}(u_3 - v_3) - \frac{u_3\dot{\epsilon}}{1+\epsilon}, \quad (56) $$

$$ \dot{v}_1 = \frac{\kappa}{\tau_c}(u_1 - v_1) - (u_2v_3 - u_3v_2) + \epsilon u_3 v_2, \quad (57) $$

$$ \dot{v}_2 = \frac{\kappa}{\tau_c}(u_2 - v_2) - (u_3v_1 - u_1v_3) - \epsilon u_3 v_1, \quad (58) $$

$$ \dot{v}_3 = \frac{\kappa}{\tau_c}(u_3 - v_3) - (u_1v_2 - u_2v_1). \quad (59) $$

Further, we take

$$ \epsilon\delta = u - v. \quad (60) $$

For XTE J1810–197, the damping timescale of the precession $\tau_d$ is on the order of $P_f$. According to the Bondi–Gold relation$^{79}$,

$$ \frac{P_s}{P_f} \approx \frac{\tau_c}{\tau_d}, \quad (61) $$

we quickly notice that $\tau_c \approx P_s$. The crust and the fluid core are so strongly coupled that $u$ and $v$ become nearly aligned in the timescale of $\tau_c$, with a difference on the order of $\epsilon$. Besides, the relaxing timescale of the ellipticity $\tau_c \gg \tau_c$. Thus, we can neglect the higher order contributions and set $\dot{\delta} = 0$. In this approximation, the Euler equations are

$$ \dot{u}_1 = -\frac{\epsilon}{\tau_c}\delta_1, \quad (62) $$

$$ \dot{u}_2 = -\frac{\epsilon}{\tau_c}\delta_2, \quad (63) $$

$$ \dot{u}_3 = -\frac{\epsilon}{\tau_c}\delta_3, \quad (64) $$

$$ 0 = -\frac{1+\kappa}{\tau_c}\delta_1 + u_3\delta_2 - u_2\delta_3 - u_2u_3, \quad (65) $$

$$ 0 = -u_3\delta_1 - \frac{1+\kappa}{\tau_c}\delta_2 + u_1\delta_3 + u_3u_1, \quad (66) $$

$$ 0 = u_2\delta_1 - u_1\delta_2 - \frac{1+\kappa}{\tau_c}\delta_3. \quad (67) $$

We have three differential equations and three algebraic equations, which can be solved very fast. Once we integrate out $u$ and $\delta$, the angular frequencies $\omega_c$ and $\Omega_f$ can be obtained by a rotation transformation. Then the precession phase can be obtained from the relation $\lambda = \arctan \frac{\omega_1}{\omega_2}$. We considered again the same description for the decay of $\epsilon$ as in model 3, with equations (39) and (40).

## Running the free-precession models

To perform the comparison between the different free-precession models described above, we again apply the RVM to the same Stokes Q and U data as used in the comparison between free and forced precession. Across all epochs, $\zeta$ remains the common RVM parameter. All priors for the parameters in the free-precession models are described in Supplementary Table 4. In all models, we assume $P_s$ to be constant after checking that taking into account $P_s$ does not affect the results. For each model, we also compared the Bayesian log-evidence log $Z$ with and without including the toy model for a twisted magnetosphere described previously. We report log $Z$ from all models in Extended Data Table 1. The preferred model with a Bayes factor of 45 is the phenomenological model of wobble angle decay with relaxing ellipticity to non-zero value. The posteriors and covariances between the parameters are shown in Supplementary Fig. 14. The next preferred model is the model based on frictional crust–core coupling combined with relaxing ellipticity to a non-zero value. The posteriors and covariances between the parameters of this model are shown in Supplementary Fig. 15.

## Data availability

The polarimetric profiles can be downloaded via the Max Planck Digital Library: https://keeper.mpdl.mpg.de/d/1ab5be1fcf974027bf07/.

## Code availability

The modelling codes can be downloaded from https://github.com/gdesvignes/XTE_J1810-197.

## References

1. Duncan, R. C. & Thompson, C. Formation of very strongly magnetized neutron stars: implications for gamma-ray bursts. Astrophys. J. 392, L9 (1992).

2. Beloborodov, A. M. Untwisting magnetospheres of neutron stars. Astrophys. J. 703, 1044–1060 (2009).

3. Olausen, S. A. & Kaspi, V. M. The McGill Magnetar Catalog. Astrophys. J. Suppl. Ser. 212, 6 (2014).

4. CHIME/FRB Collaboration. A bright millisecond-duration radio burst from a Galactic magnetar. Nature 587, 54–58 (2020).

5. Cordes, J. M. & Chatterjee, S. Fast radio bursts: an extragalactic enigma. Annu. Rev. Astron. Astrophys. 57, 417–465 (2019).

6. Zanazzi, J. J. & Lai, D. Periodic fast radio bursts with neutron star free precession. Astrophys. J. 892, L15 (2020).

7. Levin, Y., Beloborodov, A. M. & Bransgrove, A. Precessing flaring magnetar as a source of repeating FRB 180916.J0158+65. Astrophys. J. 895, L30 (2020).
where $\kappa = I_c/I_f$ is the ratio between the moment of inertia of the crust and the core.

In principle, one can integrate equations (45)–(50) with appropriate initial conditions. However, it is very slow to integrate those six equations directly. To do parameter estimation more efficiently, we first make the following transformations

$$ \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix} = \begin{pmatrix} \cos \gamma - \sin \gamma & 0 \\ \sin \gamma & \cos \gamma & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{pmatrix}, \quad (51) $$

and

$$ \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix} = \begin{pmatrix} \cos \gamma - \sin \gamma & 0 \\ \sin \gamma & \cos \gamma & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \Omega_1 \\ \Omega_2 \\ \Omega_3 \end{pmatrix}. \quad (52) $$

Here the angle $\gamma$ is defined as

$$ \gamma(t) = -\epsilon\omega_3 dt. \quad (53) $$

Then the Euler equations can be written as

$$ \dot{u}_1 = -\frac{1}{\tau_c} (u_1 - v_1), \quad (54) $$

$$ \dot{u}_2 = -\frac{1}{\tau_c} (u_2 - v_2), \quad (55) $$

$$ \dot{u}_3 = -\frac{1}{\tau_c(1+\epsilon)} (u_3 - v_3) - \frac{u_3 \dot{\epsilon}}{1+\epsilon}, \quad (56) $$

$$ \dot{v}_1 = \frac{\kappa}{\tau_c} (u_1 - v_1) - (u_2v_3 - u_3v_2) + \epsilon u_3 v_2, \quad (57) $$

$$ \dot{v}_2 = \frac{\kappa}{\tau_c} (u_2 - v_2) - (u_3v_1 - u_1v_3) - \epsilon u_3 v_1, \quad (58) $$

$$ \dot{v}_3 = \frac{\kappa}{\tau_c} (u_3 - v_3) - (u_1v_2 - u_2v_1). \quad (59) $$

Further, we take

$$ \epsilon\delta = u - v. \quad (60) $$

For XTE J1810–197, the damping timescale of the precession $\tau_d$ is on the order of $P_f$. According to the Bondi–Gold relation$^{79}$,

$$ \frac{P_s}{P_f} \approx \frac{\tau_c}{\tau_d}, \quad (61) $$

we quickly notice that $\tau_c \approx P_s$. The crust and the fluid core are so strongly coupled that $u$ and $v$ become nearly aligned in the timescale of $\tau_c$, with a difference on the order of $\epsilon$. Besides, the relaxing timescale of the ellipticity $\tau_c \gg \tau_c$. Thus, we can neglect the higher order contributions and set $\dot{\delta} = 0$. In this approximation, the Euler equations are

$$ \dot{u}_1 = -\frac{\epsilon}{\tau_c} \delta_1, \quad (62) $$

$$ \dot{u}_2 = -\frac{\epsilon}{\tau_c} \delta_2, \quad (63) $$

$$ \dot{u}_3 = -\frac{\epsilon}{\tau_c} \delta_3, \quad (64) $$

$$ 0 = -\frac{1+\kappa}{\tau_c} \delta_1 + u_3 \delta_2 - u_2 \delta_3 - u_2 u_3, \quad (65) $$

$$ 0 = -u_3 \delta_1 - \frac{1+\kappa}{\tau_c} \delta_2 + u_1 \delta_3 + u_3 u_1, \quad (66) $$

$$ 0 = u_2 \delta_1 - u_1 \delta_2 - \frac{1+\kappa}{\tau_c} \delta_3. \quad (67) $$

We have three differential equations and three algebraic equations, which can be solved very fast. Once we integrate out $u$ and $\delta$, the angular frequencies $\omega_c$ and $\Omega_f$ can be obtained by a rotation transformation. Then the precession phase can be obtained from the relation $\lambda = \arctan \frac{\omega_1}{\omega_2}$. We considered again the same description for the decay of $\epsilon$ as in model 3, with equations (39) and (40).

## Running the free-precession models

To perform the comparison between the different free-precession models described above, we again apply the RVM to the same Stokes Q and U data as used in the comparison between free and forced precession. Across all epochs, $\zeta$ remains the common RVM parameter. All priors for the parameters in the free-precession models are described in Supplementary Table 4. In all models, we assume $P_s$ to be constant after checking that taking into account $P_s$ does not affect the results. For each model, we also compared the Bayesian log-evidence log $Z$ with and without including the toy model for a twisted magnetosphere described previously. We report log $Z$ from all models in Extended Data Table 1. The preferred model with a Bayes factor of 45 is the phenomenological model of wobble angle decay with relaxing ellipticity to non-zero value. The posteriors and covariances between the parameters are shown in Supplementary Fig. 14. The next preferred model is the model based on frictional crust–core coupling combined with relaxing ellipticity to a non-zero value. The posteriors and covariances between the parameters of this model are shown in Supplementary Fig. 15.

## Data availability

The polarimetric profiles can be downloaded via the Max Planck Digital Library: https://keeper.mpdl.mpg.de/d/1ab5be1fcf974027bf07/.

## Code availability

The modelling codes can be downloaded from https://github.com/gdesvignes/XTE_J1810-197.

## References

1. Duncan, R. C. & Thompson, C. Formation of very strongly magnetized neutron stars: implications for gamma-ray bursts. Astrophys. J. **392**, L9 (1992).

2. Beloborodov, A. M. Untwisting magnetospheres of neutron stars. Astrophys. J. **703**, 1044–1060 (2009).

3. Olausen, S. A. & Kaspi, V. M. The McGill Magnetar Catalog. Astrophys. J. Suppl. Ser. **212**, 6 (2014).

4. CHIME/FRB Collaboration. A bright millisecond-duration radio burst from a Galactic magnetar. Nature **587**, 54–58 (2020).

5. Cordes, J. M. & Chatterjee, S. Fast radio bursts: an extragalactic enigma. Annu. Rev. Astron. Astrophys. **57**, 417–465 (2019).

6. Zanazzi, J. J. & Lai, D. Periodic fast radio bursts with neutron star free precession. Astrophys. J. **892**, L15 (2020).

7. Levin, Y., Beloborodov, A. M. & Bransgrove, A. Precessing flaring magnetar as a source of repeating FRB 180916.J0158+65. Astrophys. J. **895**, L30 (2020).
8. Ibrahim, A. I. et al. Discovery of a transient magnetar: XTE J1810-197. Astrophys. J. 609, L21-L24 (2004).

9. Camilo, F. et al. Transient pulsed radio emission from a magnetar. Nature 442, 892-895 (2006).

10. Camilo, F. et al. Radio disappearance of the magnetar XTE J1810-197 and continued X-ray timing. Astrophys. J. 820, 110 (2016).

11. Levin, L. et al. Spin frequency evolution and pulse profile variations of the recently re-activated radio magnetar XTE J1810-197. Mon. Not. R. Astron. Soc. 488, 5251-5258 (2019).

12. Gotthelf, E. V. et al. The 2018 X-ray and radio outburst of magnetar XTE J1810-197. Astrophys. J. 874, L25 (2019).

13. Kramer, M., Stappers, B. W., Jessner, A., Lyne, A. G. & Jordan, C. A. Polarized radio emission from a magnetar. Mon. Not. R. Astron. Soc. 377, 107-119 (2007).

14. Camilo, F. et al. Polarized radio emission from the magnetar XTE J1810-197. Astrophys. J. 659, L37-L40 (2007).

15. Dai, S. et al. Wideband polarized radio emission from the newly revived magnetar XTE J1810-197. Astrophys. J. 874, L14 (2019).

16. Radhakrishnan, V. & Cooke, D. J. Magnetic poles and the polarization structure of pulsar radiation. Astrophys. J. 3, 225-229 (1969).

17. Komesaroff, M. M. Possible mechanism for the pulsar radio emission. Nature 225, 612-614 (1970).

18. Desvignes, G. et al. Radio emission from a pulsar's magnetic pole revealed by general relativity. Science 365, 1013-1017 (2019).

19. Thompson, C., Lyutikov, M. & Kulkarni, S. R. Electrodynamics of magnetars: implications for the persistent X-ray emission and spin-down of the soft gamma repeaters and anomalous X-ray pulsars. Astrophys. J. 574, 332-355 (2002).

20. Tong, H., Wang, P. F., Wang, H. G. & Yan, Z. Rotating vector model for magnetars. Mon. Not. R. Astron. Soc. 502, 1549-1556 (2021).

21. Pétri, J. Polarized emission from an off-centred dipole. Mon. Not. R. Astron. Soc. 466, L73-L77 (2017).

22. Lower, M. E., Johnston, S., Shannon, R. M., Bailes, M. & Camilo, F. The dynamic magnetosphere of Swift J1818.0-1607. Mon. Not. R. Astron. Soc. 502, 127-139 (2021).

23. Dyks, J. Artefacts of circumpolar cartography in radio pulsar polarization. Mon. Not. R. Astron. Soc. 495, L118-L122 (2020).

24. Damour, T. & Ruffini, R. Certain new verifications of general relativity made possible by the discovery of a pulsar belonging to a binary system. C. R. Acad. Sci. Paris, Série A 279, 971-973 (1974).

25. Jones, D. I. & Andersson, N. Freely precessing neutron stars: model and observations. Mon. Not. R. Astron. Soc. 324, 811-824 (2001).

26. Shabanova, T. V., Lyne, A. G. & Urama, J. O. Evidence for free precession in the pulsar B1642-03. Astrophys. J. 552, 321-325 (2001).

27. Stairs, I. H., Lyne, A. G. & Shemar, S. L. Evidence for free precession in a pulsar. Nature 406, 484-486 (2000).

28. Lyne, A., Hobbs, G., Kramer, M., Stairs, I. & Stappers, B. Switched magnetospheric regulation of pulsar spin-down. Science 329, 408 (2010).

29. Stairs, I. H. et al. Mode switching and oscillations in PSR B1828-11. Mon. Not. R. Astron. Soc. 485, 3230-3240 (2019).

30. Makishima, K. et al. Possible evidence for free precession of a strongly magnetized neutron star in the magnetar 4U 0142+61. Phys. Rev. Lett. 112, 171102 (2014).

31. Makishima, K. et al. Evidence for a 36 ks phase modulation in the hard X-ray pulses from the magnetar 1E 1547.0-5408. Publ. Astron. Soc. Jpn 68, S12 (2016).

32. Makishima, K. et al. Discovery of 40.5 ks hard X-ray pulse-phase modulations from SGR 1900+14. Astrophys. J. 923, 63 (2021).

33. Handley, W. J., Hobson, M. P. & Lasenby, A. N. POLYCHORD: next-generation nested sampling. Mon. Not. R. Astron. Soc. 453, 4384-4398 (2015).

34. Maan, Y., Surnis, M. P., Chandra Joshi, B. & Bagchi, M. Magnetar XTE J1810-197: spectro-temporal evolution of average radio emission. Astrophys. J. 931, 67 (2022).

35. Caleb, M. et al. Radio and X-ray observations of giant pulses from XTE J1810-197. Mon. Not. R. Astron. Soc. 510, 1996-2010 (2022).

36. Ushomirsky, G., Cutler, C. & Bildsten, L. Deformations of accreting neutron star crusts and gravitational wave emission. Mon. Not. R. Astron. Soc. 319, 902 (2000).

37. Haskell, B., Jones, D. I. & Andersson, N. Mountains on neutron stars: accreted vs. non-accreted crusts. Mon. Not. R. Astron. Soc. 373, 1423-1439 (2006).

38. Johnson-McDaniel, N. K. & Owen, B. J. Maximum elastic deformations of relativistic stars. Phys. Rev. D 88, 044004 (2013).

39. Horowitz, C. J. & Kadau, K. The breaking strain of neutron star crust and gravitational waves. Phys. Rev. Lett. 102, 191102 (2009).

40. Haskell, B., Samuelsson, L., Glampedakis, K. & Andersson, N. Modelling magnetically deformed neutron stars. Mon. Not. R. Astron. Soc. 385, 531-542 (2008).

41. Akgun, T. & Wasserman, I. Toroidal magnetic fields in Type II superconducting neutron stars. Mon. Not. R. Astron. Soc. 383, 1551-1580 (2008).

42. Lander, S. K. & Jones, D. I. Magnetic fields in axisymmetric neutron stars. Mon. Not. R. Astron. Soc. 395, 2162 (2009).

43. Mastrano, A., Suvorov, A. G. & Melatos, A. Neutron star deformation due to poloidal-toroidal magnetic fields of arbitrary multipole order: a new analytic approach. Mon. Not. R. Astron. Soc. 447, 3475 (2015).

44. Güver, T., Özel, F., Gögüs, E. & Kouveliotou, C. The magnetar nature and the outburst mechanism of a transient anomalous X-ray pulsar. Astrophys. J. 667, L73-L76 (2007).

45. Levin, L. et al. A radio-loud magnetar in X-ray quiescence. Astrophys. J. 721, L33-L37 (2010).

46. Camilo, F. et al. Revival of the magnetar PSR J1622-4950: observations with MeerKAT, Parkes, XMM-Newton, Swift, Chandra, and NuSTAR. Astrophys. J. 856, 180 (2018).

47. Li, J., Spitkovsky, A. & Tchekhovskoy, A. Resistive solutions for pulsar magnetospheres. Astrophys. J. 746, 60 (2012).

48. Gangadhara, R. T. Circular polarization in pulsars due to curvature radiation. Astrophys. J. 710, 29-44 (2010).

49. Thompson, C. et al. Physical mechanisms for the variable spin-down and light curve of SGR 1900+14. Astrophys. J. 543, 340-350 (2000).

50. Albano, A. et al. A unified timing and spectral model for the anomalous X-ray pulsars XTE J1810-197 and CXOU J164710.2-455216. Astrophys. J. 722, 788-802 (2010).

51. Bernardini, F. et al. Emission geometry, radiation pattern and magnetic topology of the magnetar XTE J1810-197 in its quiescent state. Mon. Not. R. Astron. Soc. 418, 638-647 (2011).

52. Borghese, A. et al. The X-ray evolution and geometry of the 2018 outburst of XTE J1810-197. Mon. Not. R. Astron. Soc. 504, 5244-5257 (2021).

53. Everett, J. E. & Weisberg, J. M. Emission beam geometry of selected pulsars derived from average pulse polarization data. Astrophys. J. 553, 341-357 (2001).

54. Ding, H. et al. A magnetar parallax. Mon. Not. R. Astron. Soc. 498, 3736-3743 (2020).

55. Zimmermann, M. & Szedenits, E. Gravitational waves from rotating and precessing rigid bodies: simple models and applications to pulsars. Phys. Rev. D 20, 351-355 (1979).

56. Abbott, B. et al. Upper limits on gravitational wave emission from 78 radio pulsars. Phys. Rev. D 76, 042001 (2007).

57. Chime/Frb Collaboration. Periodic activity from a fast radio burst source. Nature 582, 351-355 (2020).

58. Rajwade, K. M. et al. Possible periodic activity in the repeating FRB 121102. Mon. Not. R. Astron. Soc. 495, 3551-3558 (2020).
8. Ibrahim, A. I. et al. Discovery of a transient magnetar: XTE J1810-197. Astrophys. J. 609, L21-L24 (2004).

9. Camilo, F. et al. Transient pulsed radio emission from a magnetar. Nature 442, 892-895 (2006).

10. Camilo, F. et al. Radio disappearance of the magnetar XTE J1810-197 and continued X-ray timing. Astrophys. J. 820, 110 (2016).

11. Levin, L. et al. Spin frequency evolution and pulse profile variations of the recently re-activated radio magnetar XTE J1810-197. Mon. Not. R. Astron. Soc. 488, 5251-5258 (2019).

12. Gotthelf, E. V. et al. The 2018 X-ray and radio outburst of magnetar XTE J1810-197. Astrophys. J. 874, L25 (2019).

13. Kramer, M., Stappers, B. W., Jessner, A., Lyne, A. G. & Jordan, C. A. Polarized radio emission from a magnetar. Mon. Not. R. Astron. Soc. 377, 107-119 (2007).

14. Camilo, F. et al. Polarized radio emission from the magnetar XTE J1810-197. Astrophys. J. 659, L37-L40 (2007).

15. Dai, S. et al. Wideband polarized radio emission from the newly revived magnetar XTE J1810-197. Astrophys. J. 874, L14 (2019).

16. Radhakrishnan, V. & Cooke, D. J. Magnetic poles and the polarization structure of pulsar radiation. Astrophys. J. 3, 225-229 (1969).

17. Komesaroff, M. M. Possible mechanism for the pulsar radio emission. Nature 225, 612-614 (1970).

18. Desvignes, G. et al. Radio emission from a pulsar's magnetic pole revealed by general relativity. Science 365, 1013-1017 (2019).

19. Thompson, C., Lyutikov, M. & Kulkarni, S. R. Electrodynamics of magnetars: implications for the persistent X-ray emission and spin-down of the soft gamma repeaters and anomalous X-ray pulsars. Astrophys. J. 574, 332-355 (2002).

20. Tong, H., Wang, P. F., Wang, H. G. & Yan, Z. Rotating vector model for magnetars. Mon. Not. R. Astron. Soc. 502, 1549-1556 (2021).

21. Pétri, J. Polarized emission from an off-centred dipole. Mon. Not. R. Astron. Soc. 466, L73-L77 (2017).

22. Lower, M. E., Johnston, S., Shannon, R. M., Bailes, M. & Camilo, F. The dynamic magnetosphere of Swift J1818.0-1607. Mon. Not. R. Astron. Soc. 502, 127-139 (2021).

23. Dyks, J. Artefacts of circumpolar cartography in radio pulsar polarization. Mon. Not. R. Astron. Soc. 495, L118-L122 (2020).

24. Damour, T. & Ruffini, R. Certain new verifications of general relativity made possible by the discovery of a pulsar belonging to a binary system. C. R. Acad. Sci. Paris, Série A 279, 971-973 (1974).

25. Jones, D. I. & Andersson, N. Freely precessing neutron stars: model and observations. Mon. Not. R. Astron. Soc. 324, 811-824 (2001).

26. Shabanova, T. V., Lyne, A. G. & Urama, J. O. Evidence for free precession in the pulsar B1642-03. Astrophys. J. 552, 321-325 (2001).

27. Stairs, I. H., Lyne, A. G. & Shemar, S. L. Evidence for free precession in a pulsar. Nature 406, 484-486 (2000).

28. Lyne, A., Hobbs, G., Kramer, M., Stairs, I. & Stappers, B. Switched magnetospheric regulation of pulsar spin-down. Science 329, 408 (2010).

29. Stairs, I. H. et al. Mode switching and oscillations in PSR B1828-11. Mon. Not. R. Astron. Soc. 485, 3230-3240 (2019).

30. Makishima, K. et al. Possible evidence for free precession of a strongly magnetized neutron star in the magnetar 4U 0142+61. Phys. Rev. Lett. 112, 171102 (2014).

31. Makishima, K. et al. Evidence for a 36 ks phase modulation in the hard X-ray pulses from the magnetar 1E 1547.0-5408. Publ. Astron. Soc. Jpn 68, S12 (2016).

32. Makishima, K. et al. Discovery of 40.5 ks hard X-ray pulse-phase modulations from SGR 1900+14. Astrophys. J. 923, 63 (2021).

33. Handley, W. J., Hobson, M. P. & Lasenby, A. N. POLYCHORD: next-generation nested sampling. Mon. Not. R. Astron. Soc. 453, 4384-4398 (2015).

34. Maan, Y., Surnis, M. P., Chandra Joshi, B. & Bagchi, M. Magnetar XTE J1810-197: spectro-temporal evolution of average radio emission. Astrophys. J. 931, 67 (2022).

35. Caleb, M. et al. Radio and X-ray observations of giant pulses from XTE J1810-197. Mon. Not. R. Astron. Soc. 510, 1996-2010 (2022).

36. Ushomirsky, G., Cutler, C. & Bildsten, L. Deformations of accreting neutron star crusts and gravitational wave emission. Mon. Not. R. Astron. Soc. 319, 902 (2000).

37. Haskell, B., Jones, D. I. & Andersson, N. Mountains on neutron stars: accreted vs. non-accreted crusts. Mon. Not. R. Astron. Soc. 373, 1423-1439 (2006).

38. Johnson-McDaniel, N. K. & Owen, B. J. Maximum elastic deformations of relativistic stars. Phys. Rev. D 88, 044004 (2013).

39. Horowitz, C. J. & Kadau, K. The breaking strain of neutron star crust and gravitational waves. Phys. Rev. Lett. 102, 191102 (2009).

40. Haskell, B., Samuelsson, L., Glampedakis, K. & Andersson, N. Modelling magnetically deformed neutron stars. Mon. Not. R. Astron. Soc. 385, 531-542 (2008).

41. Akgun, T. & Wasserman, I. Toroidal magnetic fields in Type II superconducting neutron stars. Mon. Not. R. Astron. Soc. 383, 1551-1580 (2008).

42. Lander, S. K. & Jones, D. I. Magnetic fields in axisymmetric neutron stars. Mon. Not. R. Astron. Soc. 395, 2162 (2009).

43. Mastrano, A., Suvorov, A. G. & Melatos, A. Neutron star deformation due to poloidal-toroidal magnetic fields of arbitrary multipole order: a new analytic approach. Mon. Not. R. Astron. Soc. 447, 3475 (2015).

44. Güver, T., Özel, F., Gögüs, E. & Kouveliotou, C. The magnetar nature and the outburst mechanism of a transient anomalous X-ray pulsar. Astrophys. J. 667, L73-L76 (2007).

45. Levin, L. et al. A radio-loud magnetar in X-ray quiescence. Astrophys. J. 721, L33-L37 (2010).

46. Camilo, F. et al. Revival of the magnetar PSR J1622-4950: observations with MeerKAT, Parkes, XMM-Newton, Swift, Chandra, and NuSTAR. Astrophys. J. 856, 180 (2018).

47. Li, J., Spitkovsky, A. & Tchekhovskoy, A. Resistive solutions for pulsar magnetospheres. Astrophys. J. 746, 60 (2012).

48. Gangadhara, R. T. Circular polarization in pulsars due to curvature radiation. Astrophys. J. 710, 29-44 (2010).

49. Thompson, C. et al. Physical mechanisms for the variable spin-down and light curve of SGR 1900+14. Astrophys. J. 543, 340-350 (2000).

50. Albano, A. et al. A unified timing and spectral model for the anomalous X-ray pulsars XTE J1810-197 and CXOU J164710.2-455216. Astrophys. J. 722, 788-802 (2010).

51. Bernardini, F. et al. Emission geometry, radiation pattern and magnetic topology of the magnetar XTE J1810-197 in its quiescent state. Mon. Not. R. Astron. Soc. 418, 638-647 (2011).

52. Borghese, A. et al. The X-ray evolution and geometry of the 2018 outburst of XTE J1810-197. Mon. Not. R. Astron. Soc. 504, 5244-5257 (2021).

53. Everett, J. E. & Weisberg, J. M. Emission beam geometry of selected pulsars derived from average pulse polarization data. Astrophys. J. 553, 341-357 (2001).

54. Ding, H. et al. A magnetar parallax. Mon. Not. R. Astron. Soc. 498, 3736-3743 (2020).

55. Zimmermann, M. & Szedenits, E. Gravitational waves from rotating and precessing rigid bodies: simple models and applications to pulsars. Phys. Rev. D 20, 351-355 (1979).

56. Abbott, B. et al. Upper limits on gravitational wave emission from 78 radio pulsars. Phys. Rev. D 76, 042001 (2007).

57. Chime/Frb Collaboration. Periodic activity from a fast radio burst source. Nature 582, 351-355 (2020).

58. Rajwade, K. M. et al. Possible periodic activity in the repeating FRB 121102. Mon. Not. R. Astron. Soc. 495, 3551-3558 (2020).
59. Bassa, C. G. et al. LEAP: the Large European Array for Pulsars. *Mon. Not. R. Astron. Soc.* **456**, 2196–2209 (2016).

60. van Straten, W. High-fidelity radio astronomical polarimetry using a millisecond pulsar as a polarized reference source. Astrophys. J. Suppl. Ser. **204**, 13 (2013).

61. Lazarus, P. et al. Prospects for high-precision pulsar timing with the new Effelsberg PSRIX backend. Mon. Not. R. Astron. Soc. **458**, 868–880 (2016).

62. Desvignes, G. et al. Large magneto-ionic variations toward the Galactic Center magnetar, PSR J1745-2900. Astrophys. J. **852**, L12 (2018).

63. van Straten, W., Demorest, P. & Olsowski, S. Pulsar data analysis with PSRCHIVE. Astron. Res. Technol. **9**, 237–256 (2012).

64. Lorimer, D. R. & Kramer, M. *Handbook of Pulsar Astronomy* (Cambridge Univ. Press, 2005).

65. Manchester, R. N., Taylor, J. H. & Huguenin, G. R. Observations of pulsar radio emission. II. Polarization of individual pulses. Astrophys. J. **196**, 83–102 (1975).

66. Ilie, C. D., Weltevrede, P., Johnston, S. & Chen, T. The drifting sub-pulses of PSR B0031-07 and its synchronously modulated radio polarization. Mon. Not. R. Astron. Soc. **491**, 3385–3394 (2020).

67. Lower, M. E. et al. Linear to circular conversion in the polarized radio emission of a magnetar. Preprint at https://doi.org/10.48550/arXiv.2311.04195 (2023).

68. Desvignes, G. gdesvignes/modelRVM: modelRVM: a nested sampling tool to model the RVM. Zenodo https://doi.org/10.5281/ zenodo.3265755 (2019).

69. Jeffreys, H. *Theory of Probability* 3rd edn (Oxford, 1961).

70. Tong, H. & Huang, L. Magnetar outburst and spin-down glitch. Mon. Not. R. Astron. Soc. **497**, 2680–2687 (2020).

71. Blaskiewicz, M., Cordes, J. M. & Wasserman, I. A relativistic model of pulsar polarization. Astrophys. J. **370**, 643–669 (1991).

72. Cordes, J. M. Observational limits on the location of pulsar emission regions. Astrophys. J. **222**, 1006–1011 (1978).

73. Johnston, S. et al. The Thousand-Pulsar-Array programme on MeerKAT - XI. Application of the rotating vector model. Mon. Not. R. Astron. Soc. **520**, 4801–4814 (2023).

74. Philippov, A. & Kramer, M. Pulsar magnetospheres and their radiation. Annu. Rev. Astron. Astrophys. **60**, 495–558 (2022).

75. Pétri, J. The illusion of neutron star magnetic field estimates. Mon. Not. R. Astron. Soc. **485**, 4573–4587 (2019).

76. Qiu, J. L., Tong, H. & Wang, H. G. Rotating vector model and radius-to-frequency mapping in the presence of a multipole magnetic field. Astrophys. J. **958**, 78 (2023).

77. Alpar, M. & Pines, D. Gravitational radiation from a solid crust neutron star. Nature **314**, 334–336 (1985).

78. Cutler, C. & Jones, D. I. Gravitational wave damping of neutron star wobble. Phys. Rev. D **63**, 024002 (2001).

79. Bondi, H. & Gold, T. On the damping of free nutation of the Earth. Mon. Not. R. Astron. Soc. **115**, 41 (1955).

80. Shaham, J. Free precession of neutron stars: role of possible vortex pinning. Astrophys. J. **214**, 251–260 (1977).

81. Sedrakian, A., Wasserman, I. & Cordes, J. M. Precession of isolated neutron stars I: effects of imperfect pinning. Astrophys. J. **524**, 341 (1998).

## Acknowledgements

We thank M. Mickaliger for help with the Lovell observations and M. Lower and H. Tao for useful discussions. Based on observations with the 100m telescope of the Max-Planck-Institut für Radioastronomie at Effelsberg. Pulsar research at the Jodrell Bank Centre for Astrophysics and Jodrell Bank Observatory is supported by a consolidated grant from the UK Science and Technology Facilities Council. We acknowledge the use of the Max-Planck-Institut für

Radioastronomie Hercules cluster and the HPC system Raven hosted at the Max Planck Computing and Data Facility. G.D., M.K. and K.L. gratefully acknowledge support from European Research Council Synergy Grant “BlackHoleCam” Grant Agreement Number 610058. Y.G. and L.S. are supported by the National SKA Program of China (202 OSKA0120300) and the Max Planck Partner Group Program funded by the Max Planck Society. D.I.J. acknowledges support from the Science and Technology Facilities Council via grant number ST/R00045X/1. M.C. acknowledges support of an Australian Research Council Discovery Early Career Research Award (project number DE220100819) funded by the Australian Government and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions, through project number CE170100013. J.P. acknowledges support from the ANR grant number ANR-20-CE31-0010.

## Author contributions

G.D. drafted the main paper with contributions from co-authors and led the Effelsberg observations enabled by M.K. and R.K. with help from K.L. P.W. led the JBO data processing with help from M.C., L.L., A.G.L. and B.S. Y.G developed the models of free precession and wrote the Methods section on ‘Theoretical models’ with input from D.I.J. and L.S. J.P. contributed to the theoretical interpretation of the polarization data.

## Funding

Open access funding provided by Max Planck Society.

## Competing interests

The authors declare no competing interests.

## Additional information

**Extended data** is available for this paper at
https://doi.org/10.1038/s41550-024-02226-7.

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41550-024-02226-7.

**Correspondence and requests for materials** should be addressed to Gregory Desvignes.

**Peer review information** *Nature Astronomy* thanks the anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permissions information** is available at
www.nature.com/reprints.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024
59. Bassa, C. G. et al. LEAP: the Large European Array for Pulsars. *Mon. Not. R. Astron. Soc.* **456**, 2196–2209 (2016).

60. van Straten, W. High-fidelity radio astronomical polarimetry using a millisecond pulsar as a polarized reference source. Astrophys. J. Suppl. Ser. **204**, 13 (2013).

61. Lazarus, P. et al. Prospects for high-precision pulsar timing with the new Effelsberg PSRIX backend. Mon. Not. R. Astron. Soc. **458**, 868–880 (2016).

62. Desvignes, G. et al. Large magneto-ionic variations toward the Galactic Center magnetar, PSR J1745-2900. Astrophys. J. **852**, L12 (2018).

63. van Straten, W., Demorest, P. & Olsowski, S. Pulsar data analysis with PSRCHIVE. Astron. Res. Technol. **9**, 237–256 (2012).

64. Lorimer, D. R. & Kramer, M. *Handbook of Pulsar Astronomy* (Cambridge Univ. Press, 2005).

65. Manchester, R. N., Taylor, J. H. & Huguenin, G. R. Observations of pulsar radio emission. II. Polarization of individual pulses. Astrophys. J. **196**, 83–102 (1975).

66. Ilie, C. D., Weltevrede, P., Johnston, S. & Chen, T. The drifting sub-pulses of PSR B0031-07 and its synchronously modulated radio polarization. Mon. Not. R. Astron. Soc. **491**, 3385–3394 (2020).

67. Lower, M. E. et al. Linear to circular conversion in the polarized radio emission of a magnetar. Preprint at https://doi.org/10.48550/arXiv.2311.04195 (2023).

68. Desvignes, G. gdesvignes/modelRVM: modelRVM: a nested sampling tool to model the RVM. Zenodo https://doi.org/10.5281/ zenodo.3265755 (2019).

69. Jeffreys, H. *Theory of Probability* 3rd edn (Oxford, 1961).

70. Tong, H. & Huang, L. Magnetar outburst and spin-down glitch. Mon. Not. R. Astron. Soc. **497**, 2680–2687 (2020).

71. Blaskiewicz, M., Cordes, J. M. & Wasserman, I. A relativistic model of pulsar polarization. Astrophys. J. **370**, 643–669 (1991).

72. Cordes, J. M. Observational limits on the location of pulsar emission regions. Astrophys. J. **222**, 1006–1011 (1978).

73. Johnston, S. et al. The Thousand-Pulsar-Array programme on MeerKAT - XI. Application of the rotating vector model. Mon. Not. R. Astron. Soc. **520**, 4801–4814 (2023).

74. Philippov, A. & Kramer, M. Pulsar magnetospheres and their radiation. Annu. Rev. Astron. Astrophys. **60**, 495–558 (2022).

75. Pétri, J. The illusion of neutron star magnetic field estimates. Mon. Not. R. Astron. Soc. **485**, 4573–4587 (2019).

76. Qiu, J. L., Tong, H. & Wang, H. G. Rotating vector model and radius-to-frequency mapping in the presence of a multipole magnetic field. Astrophys. J. **958**, 78 (2023).

77. Alpar, M. & Pines, D. Gravitational radiation from a solid crust neutron star. Nature **314**, 334–336 (1985).

78. Cutler, C. & Jones, D. I. Gravitational wave damping of neutron star wobble. Phys. Rev. D **63**, 024002 (2001).

79. Bondi, H. & Gold, T. On the damping of free nutation of the Earth. Mon. Not. R. Astron. Soc. **115**, 41 (1955).

80. Shaham, J. Free precession of neutron stars: role of possible vortex pinning. Astrophys. J. **214**, 251–260 (1977).

81. Sedrakian, A., Wasserman, I. & Cordes, J. M. Precession of isolated neutron stars I: effects of imperfect pinning. Astrophys. J. **524**, 341 (1998).

## Acknowledgements

We thank M. Mickaliger for help with the Lovell observations and M. Lower and H. Tao for useful discussions. Based on observations with the 100m telescope of the Max-Planck-Institut für Radioastronomie at Effelsberg. Pulsar research at the Jodrell Bank Centre for Astrophysics and Jodrell Bank Observatory is supported by a consolidated grant from the UK Science and Technology Facilities Council. We acknowledge the use of the Max-Planck-Institut für

Radioastronomie Hercules cluster and the HPC system Raven hosted at the Max Planck Computing and Data Facility. G.D., M.K. and K.L. gratefully acknowledge support from European Research Council Synergy Grant "BlackHoleCam" Grant Agreement Number 610058. Y.G. and L.S. are supported by the National SKA Program of China (202 OSKA0120300) and the Max Planck Partner Group Program funded by the Max Planck Society. D.I.J. acknowledges support from the Science and Technology Facilities Council via grant number ST/R00045X/1. M.C. acknowledges support of an Australian Research Council Discovery Early Career Research Award (project number DE220100819) funded by the Australian Government and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions, through project number CE170100013. J.P. acknowledges support from the ANR grant number ANR-20-CE31-0010.

## Author contributions

G.D. drafted the main paper with contributions from co-authors and led the Effelsberg observations enabled by M.K. and R.K. with help from K.L. P.W. led the JBO data processing with help from M.C., L.L., A.G.L. and B.S. Y.G developed the models of free precession and wrote the Methods section on 'Theoretical models' with input from D.I.J. and L.S. J.P. contributed to the theoretical interpretation of the polarization data.

## Funding

Open access funding provided by Max Planck Society.

## Competing interests

The authors declare no competing interests.

## Additional information

**Extended data** is available for this paper at
https://doi.org/10.1038/s41550-024-02226-7.

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41550-024-02226-7.

**Correspondence and requests for materials** should be addressed to Gregory Desvignes.

**Peer review information** *Nature Astronomy* thanks the anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permissions information** is available at
www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024
<figure><img src="image_6.png" /></figure>

**Extended Data Fig. 1** | Predicted waveform for the temporal evolution of the twist parameter *n* included in the free precession modelling. The grey contours, from darker to lighter shade, show the 68%, 95% and 99.7% confidence levels for the value of *n*. The two vertical red dashed lines delimit the time window of our observations.
<figure><img src="image_6.png" /></figure>

**Extended Data Fig. 1** | Predicted waveform for the temporal evolution of the twist parameter *n* included in the free precession modelling. The grey contours, from darker to lighter shade, show the 68%, 95% and 99.7% confidence levels for the value of *n*. The two vertical red dashed lines delimit the time window of our observations.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 2** | Polarimetric pulse profiles recorded at a frequency of 8.4 GHz with the Effelsberg radio telescope on four occasions during the summer of 2006. In the bottom part of each panel, the black, red and blue lines represent the total intensity I, linear L and circular V polarisation, respectively. In the top part, the black data points indicate the PA of the linear polarisation

included in the modelling. The uncertainty in the PA is calculated given L and the off-pulse standard deviation of I$^{53}$. The dotted lines delimit the 95% confidence levels on Ψ derived from the posterior distributions (seen in Extended Data Fig. 3) assuming the RVM.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 2** | Polarimetric pulse profiles recorded at a frequency of 8.4 GHz with the Effelsberg radio telescope on four occasions during the summer of 2006. In the bottom part of each panel, the black, red and blue lines represent the total intensity I, linear L and circular V polarisation, respectively. In the top part, the black data points indicate the PA of the linear polarisation

included in the modelling. The uncertainty in the PA is calculated given L and the off-pulse standard deviation of I$^{53}$. The dotted lines delimit the 95% confidence levels on Ψ derived from the posterior distributions (seen in Extended Data Fig. 3) assuming the RVM.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 3** | Posterior distributions for the $\zeta$ and $\beta$ parameters from the RVM analysis of the Effelsberg archival data. The observations were recorded at 8.4 GHz on MJDs 53926, 53934, 53938 and 53944. The vertical dashed lines delimit the 95% confidence levels on $\zeta$ from our preferred free precession modelling of the 2018-2020 PA data with an eastward twisted magnetosphere.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 3** | Posterior distributions for the $\zeta$ and $\beta$ parameters from the RVM analysis of the Effelsberg archival data. The observations were recorded at 8.4 GHz on MJDs 53926, 53934, 53938 and 53944. The vertical dashed lines delimit the 95% confidence levels on $\zeta$ from our preferred free precession modelling of the 2018-2020 PA data with an eastward twisted magnetosphere.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 4** | Temporal evolution of the ellipticity and wobble angle as predicted by the two best models of free precession. The top and bottom panels show the prediction for the ellipticity $\epsilon$ and wobble angle $\theta$, respectively. Predictions (mean values) from the preferred model of relaxing ellipticity with frictional coupling are shown with dashed lines while predictions from the model

of wobble angle decay with relaxing ellipticity to a non-zero value are shown with dotted lines. The light blue rectangle delimit the period of time when the second X-ray outburst occurred. The red vertical dashed line indicates the time of our first radio observation.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 4** | Temporal evolution of the ellipticity and wobble angle as predicted by the two best models of free precession. The top and bottom panels show the prediction for the ellipticity $\epsilon$ and wobble angle $\theta$, respectively. Predictions (mean values) from the preferred model of relaxing ellipticity with frictional coupling are shown with dashed lines while predictions from the model

of wobble angle decay with relaxing ellipticity to a non-zero value are shown with dotted lines. The light blue rectangle delimit the period of time when the second X-ray outburst occurred. The red vertical dashed line indicates the time of our first radio observation.
<figure><img src="image_13.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 5** | Temporal evolution of the polarisation of the pulse profile as our line of sight crosses the magnetic pole. The panels (a) and (b) show the total intensity as a function of time and rotational phase while panels (c) and (d) show the circular polarisation. Panels (a) and (c) show the pulse profile

near the time of the first PA reversal, when $\beta$ changes sign as shown by the orange line. Panels (b) and (d) show the pulse profile at the time of the second crossing of the magnetic pole. The $\beta$-scale on the y-axis is obtained from a linear fit to the $\beta$ values shown in Fig. 2 in the limited MJD range of the plots.
<figure><img src="image_13.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 5** | Temporal evolution of the polarisation of the pulse profile as our line of sight crosses the magnetic pole. The panels (a) and (b) show the total intensity as a function of time and rotational phase while panels (c) and (d) show the circular polarisation. Panels (a) and (c) show the pulse profile

near the time of the first PA reversal, when $\beta$ changes sign as shown by the orange line. Panels (b) and (d) show the pulse profile at the time of the second crossing of the magnetic pole. The $\beta$-scale on the y-axis is obtained from a linear fit to the $\beta$ values shown in Fig. 2 in the limited MJD range of the plots.
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 6** | Temporal evolution of the maximum single pulse peak flux densities. In panel (a) we report the same *β* measurements as shown in Fig. 2. In panel (b), we show the maximum SP peak flux densities recorded at 6 GHz with black crosses.
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 6** | Temporal evolution of the maximum single pulse peak flux densities. In panel (a) we report the same *β* measurements as shown in Fig. 2. In panel (b), we show the maximum SP peak flux densities recorded at 6 GHz with black crosses.
<table><caption>Extended Data Table 1 | Comparison between the different models of free precession, including the effects of twisted magnetic fields</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Without twist</th><th colspan="3">With twist</th></tr><tr><th>N<sub>dim</sub></th><th>log Z</th><th>N<sub>dim</sub></th><th>λ &lt; 0<br>log Z</th><th>λ &gt; 0<br>log Z</th></tr></thead><tbody><tr><td>RE</td><td>133</td><td>115521.5(6)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>WAD</td><td>133</td><td>109073.3(6)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>WAD + RE</td><td>134</td><td>117348.9(6)</td><td>136</td><td>117340.2(6)</td><td>117681.2(6)</td></tr><tr><td><b>WAD + RE (to non-zero)</b></td><td>135</td><td>117558.4(6)</td><td><b>137</b></td><td>117556.8(6)</td><td><b>117800.4(6)</b></td></tr><tr><td>FCCC + RE</td><td>135</td><td>116746.9(6)</td><td>137</td><td>116746.6(6)</td><td>117258.6(6)</td></tr><tr><td>FCCC + RE (to non-zero)</td><td>136</td><td>117518.9(6)</td><td>138</td><td>117532.8(6)</td><td>117755.2(6)</td></tr></tbody></table>

The acronyms RE, WAD and FCCC refer to the models of relaxing ellipticity, phenomenological wobble angle decay and frictional crust-core coupling. The two least likely models (RE and WAD) were not run with the twisted RVM hence the missing values.
<table><caption>**Extended Data Table 1 | Comparison between the different models of free precession, including the effects of twisted magnetic fields**</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Without twist</th><th colspan="3">With twist</th></tr><tr><th>N<sub>dim</sub></th><th>log Z</th><th>N<sub>dim</sub></th><th>λ &lt; 0<br>log Z</th><th>λ &gt; 0<br>log Z</th></tr></thead><tbody><tr><td>RE</td><td>133</td><td>115521.5(6)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>WAD</td><td>133</td><td>109073.3(6)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>WAD + RE</td><td>134</td><td>117348.9(6)</td><td>136</td><td>117340.2(6)</td><td>117681.2(6)</td></tr><tr><td><b>WAD + RE (to non-zero)</b></td><td>135</td><td>117558.4(6)</td><td><b>137</b></td><td>117556.8(6)</td><td><b>117800.4(6)</b></td></tr><tr><td>FCCC + RE</td><td>135</td><td>116746.9(6)</td><td>137</td><td>116746.6(6)</td><td>117258.6(6)</td></tr><tr><td>FCCC + RE (to non-zero)</td><td>136</td><td>117518.9(6)</td><td>138</td><td>117532.8(6)</td><td>117755.2(6)</td></tr></tbody></table>

The acronyms RE, WAD and FCCC refer to the models of relaxing ellipticity, phenomenological wobble angle decay and frictional crust-core coupling. The two least likely models (RE and WAD) were not run with the twisted RVM hence the missing values.